Newton Raphson iteration trapped in infinite loop

Question

I'm quite a beginner in this topic, and couldn't find out the reason: sometimes the program works, sometimes not (after asking the question, it simply doensn't want to take in my answers, than I can write in as much as I want, it doesn't respond, just list out the numbers, I tiped in)

  #include <stdio.h>

float abszolut (float szam)
{
    float abszoluterteke;
    if (szam >=0)
         abszoluterteke = szam;
    else 
        abszoluterteke = -szam;
    return abszoluterteke;
}

float negyzetgyok (float szam)
{
    float pontossag = 0.000001;
    float tipp = 1;
    if (szam <0)
    {
        printf ("Megszakítás elfogadva! \nKöszönjük, hogy programunkat választotta!\n");
        return -1;
    }
    else
        {while (abszolut (tipp*tipp-szam) >= pontossag)
            tipp = (szam/tipp + tipp)/2;
        return tipp;
    }
}

int main (void)
{
    float alap, eredmeny;
    for (;;)
    {
        printf ("Melyik számnak szeretnéd meghatározni a négyzetgyökét ilyen módszerrel?\n");
        scanf ("%f", &alap);
        eredmeny = negyzetgyok (alap);
        if (eredmeny == -1)
            return 1;
        else
         printf ("A(z) %f négyzetgyöke megfelelő közelítéssel: %f\n", alap, eredmeny);




    }
    return 0;
}

The main goal was to ask in a number, find out it's square root using the Newton-Raphson method if its value is positive and aborting the process if its negative. — Soma Turi, Feb 28 '15 at 12:07
Can you provide data for which program works and data for which it fails? — VolAnd, Feb 28 '15 at 12:13
It works for 1 to 41 and doesn't works for anything above 41 — Taimour, Feb 28 '15 at 12:14
Read http://floating-point-gui.de/ and test the result of `scanf` — Basile Starynkevitch, Feb 28 '15 at 12:29

francis · Answer 1 · 2015-02-28T14:20:42.053

Change for abszolut (tipp*tipp-szam) >= pontossag*szam

The while loop must stop once tipp*tipp is close to szam. But IEEE floating point computations have a limited precision : about 7 digits for float and 15 digits for double.

So the error on float szam is about 0.0000001*szam. It's the same for tipp. Consequently, the error on tipp*tipp-szam is higher than 0.0000001*szam. If szam is large, this error will hardly become lower than 0.000001. Even if double precision is used, it is likely that while (abszolut (tipp*tipp-szam) >= pontossag) triggers an infinite loop for very large numbers.

On the other side, what happens if szam is very small, say 1e-10 ? The while loop prematurely exits and the square root of 1e-10 is computed as something about 1e-3, instead of 1e-5... The relative error is about 10000%... And using double does not change anything !

To avoid this, you can change for abszolut (tipp*tipp-szam) >= pontossag*szam.

Notice that both sides have the same dimension. If szam were in square feet, tipp would be in feet and pontossag, the precision, is dimensionless. It is a good practice to compare things having the same dimension.

If you keep noticing infinite loops, switch to double precision or increase pontossag.

To avoid infinite loop, add a counter int i; and exit the while loop if the number of iteration is 100. 100 should be sufficient, since your Newton-Raphson iteration features a quadratic convergence.

score 1 · Answer 2 · answered Feb 28 '15 at 14:56

There are a number of problems with your code.

The exit condition in your loop is flawed.
The problem with your square root algorithm is the use of the error limit pontossag. Your algorithm will give erroneous results for very small numbers and it will loop forever for numbers larger than 20 or so. To fix this, change the loop test from abszolut (tipp*tipp-szam) >= pontossag to abszolut (tipp*tipp-szam) >= pontossag*szam.

You aren't checking for all the problem cases.
If your computer uses IEEE 754 floating point, your algorithm happens to work. That's just luck. Never rely on luck when doing numerical programming. It's easy to input an infinity. For example, 3.5e38 (350000000000000000000000000000000000000) does the trick with single precision numbers (float). Your function negyzetgyok should check for infinity:

if (isinf (szam))
{
    return szam;
}

You can do much better than an initial guess of 1.0 as the square root.
An initial guess of 1.0 against 3.4e38 means a lot of needless looping. A fast and easy way to form a good initial guess is to take advantage of the fact that floating point numbers are represented internally as (1+fractional_part)*2^exponent. A good first guess is 1*2^(exponent/2). With single precision numbers,

int expo;
float tipp;
frexpf (szam, &expo);
tipp = ldexpf (1.0f, n/2);

You are using %f rather than %g to parse floating point numbers.
The %g format can parse anything that can be parsed with the %f format, plus a whole lot more.

You aren't checking the status of fscanf. Enter x when prompted to enter a number. The scanner will read that character, which stops the scan. The scanner will put that character (x) back into the input buffer, and return 0, indicating that nothing was scanned. The next time around, the scanner will read the character x again, once again put that character back into the input buffer, and once again return 0. Infinite loop! Always check the status of any of the scanf family of functions to see if the scanner scanned the number of items expected.

You are using fscanf.
There are a number of existing questions and answers at this site that address the many problems with using fscanf to read from a file. This is particularly so when reading human-generated input. People make mistakes. Ignoring that people do make mistakes on entering data is a programming error. A better approach would be to read a line into a buffer using frets and parse that line with sscanf.

Newton Raphson iteration trapped in infinite loop

2 Answers2

Linked