Thread: Square root program

You should never use the == or != operators on floating point values. The reason is that they will not be exactly equal to each other, thus your program loops forever. It is better to check if the numbers are within a certain amount of each other.

Code:

if(abs(guess - oldguess) < 0.01)

Code:

if(abs(guess - oldguess) < 0.01)

Will do a partially decent job for your typical cases, but it's far from ideal. For example it will cause the function to fail miserably at calculating the sqrt of 1.0E-40 for example. Worse still, for very large numbers that epsilon wont be sufficient to prevent an infinite loop due to rounding.

One way to prevent infinite looping due to rounding whilst still making sure that the function obtains maximum accuracy is to perform two iterations per loop:

Code:

    while (guess != oldguess) {
        oldguess = guess;
        guess = (num/guess + guess) * 0.5;
        oldguess = guess;
        guess = (num/guess + guess) * 0.5;
    }

Accuracy problem gone!

I think a more obvious problem is that it will go into an infinite loop for NANs. In that case, the while loop evaluates NAN != NAN which is always true. For that, a special test for NANs should be performed beforehand.

However that's still not enough because if num is infinity then guess and oldGuess turn into a NAN during the loop and it still ends up looping forever, so we need a test for infinity as well.

The square root of a negative number is also supposed to be a NAN, not zero.
Thus I suggest this:

Code:

double sqrt(double num)
{
    if (num != num) // Test for NANs
        return num;

    if (num < 0)
        return std::numeric_limits<double>::quiet_NaN();

    if (num == std::numeric_limits<double>::infinity())
        return num;

    double guess = num/2, oldguess = 0;

    while(guess != oldguess) {
        oldguess = guess;
        guess = (num/guess + guess) * 0.5;
        oldguess = guess;
        guess = (num/guess + guess) * 0.5;
    }

    return guess;
}

There are ways of getting a more accurate initial guess too, but I think I've covered enough already.


Having said that, I was also lucky enough for it to not go into an infinite loop when running to code myself over millions of numbers, and I'm increasing by the minimum amount possible using the equivalent of nextafter. Of course I was looking at small numbers, and the problem may be more likely with larger numbers.