square root algorithm C++

Question

I can not figure out why this algorithm enters an infinite loop if the number entered is over 12 digits. Can anyone see why it will never end? Thanks. I just updated the algorithm to use the fabs() function and still get an infinite loop.

double squareroot(double x)
{ /* computes the square root of x */
/* make sure x is not negative .. no math crimes allowed! */
assert( x >= 0 );
if (x==0) return 0;
/* the sqrt must be between xhi and xlo */
double xhi = x;
double xlo = 0;
double guess = x/2;
/* We stop when guess*guess-x is very small */
while (abs(guess*guess-x) > 0.00001 )
{
 if (guess*guess > x){
 xhi = guess;
 }
 else {
 xlo = guess;
 }
 guess = (xhi + xlo)/2;
}
return guess;
}

• M Oehm: abs is a function for integers. Did you mean fabs?
• Ankur: 12 digits of guess*guess ? that would be big..
• n. m. could be an AI: floating-point-gui.de
• M Oehm: Your fixed epsilon for the termination criterion should probably also be expressed in relation to the input value. Otherwise your calculation for 16.0e-8, say, will not even enter the first iteration.
• TNA: @M Oehm thats true for C but in C++, there is also a floting point overload for abs in <cmath>.

Answer 1

I believe you should use relative error for the termination, and not the absolute error.

while (abs((guess*guess-x) / guess) > 0.00001)

Otherwise it will take very long time (it's not an infinite loop) to compute square root of very long values.

http://en.wikipedia.org/wiki/Approximation_error

Cheers!

EDIT: moreover, as pointed below in the comments, it is worthy to check if the guess was already guessed in order to avoid infinite loop with some specific corner cases.

• M Oehm: The loop may well become infinite when the numbers are large: When xhi and xlo are adjacent, the midpoint will always be one of them. When xhi - xlo is greater than the chosen epsilon, you get an infinite loop.
• leemes: @MOehm Good point. This can be checked additionally, i.e. the loop should be exited when the new guess was already guessed before (guess == xhi || guess == xlo)
• ale64bit: @MOehm Agree. The answer of TonyD fixes this issue as leemes pointed out in his comment. (although it would be cool to get an actual x in order to make it happen :D)

Answer 2

I suggest waiting until you've got a stable answer, rather than fiddling with epsilon values:

double squareroot(double x)
{
 if (x < 1) return 1.0 / squareroot(x); // MSalter's general solution
 double xhi = x;
 double xlo = 0;
 double guess = x/2;
 while (guess * guess != x)
 {
 if (guess * guess > x)
 xhi = guess;
 else
 xlo = guess;
 double new_guess = (xhi + xlo) / 2;
 if (new_guess == guess)
 break; // not getting closer
 guess = new_guess;
 }
 return guess;
}

• M Oehm: Very good! The remedy to an infinite loop that occurs because the same guess is used over and over is of course to enforce a different guess.
• MSalters: Does this work for x=0.25? Because the algorithm appears to assume sqrt(x) lies in the range [0,x] while sqrt(0.25) > 0.25. The first guess is 0.125, the second guess 0.1875 etc. It appears xlo converges to xhi and the loop terminates when xlo=xhi=guess=0.25, which obviously is not 0.5
• Tony Delroy: @MSalters: nope... and it's worthy of mention; I'm not attempting to fix the algo - just dealing with the convergence / precision issue in the question.
• MSalters: Well, the fix would be easy anyway. sqrt(1/x) = 1/sqrt(x), so just add if (x<1) return 1.0/squareroot(x);.

Answer 3

This is not a direct answer to your question, but an alternative solution.

You can use Newton's method for finding roots:

assert(x >= 0);
if (x == 0)
 return 0;
double guess = x;
for (int i=0; i<NUM_OF_ITERATIONS; i++)
 guess -= (guess*guess-x)/(2*guess);
return guess;

24 iterations should get you a good enough approximation, but you can also check the absolute difference.

Answer 4

http://floating-point-gui.de/errors/comparison/

Heh, looks like you are not supposed to use abs() like that. There are some cases where it should stop but i won't since it's limited precision.

Instead use fabs()

http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm

• csciXAV_12: Thanks, I just updated it to use the fabs() function and I am still getting an infinite loop.
• n. m. could be an AI: Apparently there are a few C++ programmers out there that don't know about std::abs overloaded function.
• M Oehm: @n.m.: Touché. I came here via the algorithm tag and I didn't see the C++ tag. Just looked at the code and saw it was C. No streams, no namespaces, no <math>.
• juanchopanza: @csciXAV_12 There is a overload for std::abs (and possibly abs in header <cmath>. See en.cppreference.com/w/cpp/numeric/math/fabs. It is important to post an MCVE, and if possible, use std::abs explicitly.

Answer 5

I'd say that when numbers are big enough you can't use absolute epsilon value because it doesn't fit into precision.

Try to use relative comparison instead. Consider following function to check if 2 doubles are equal:

bool are_eql(double n1, double n2, double abs_e, double rel_e)
{
 // also add check that n1 and n2 are not NaN
 double diff = abs(n1 - n2);
 if (diff < abs_e)
 return true;
 double largest = max(abs(n1), abs(n2));
 if (diff < largest * rel_e)
 return true;
 return false;
}

• MSalters: The first lijne obviously should be abs(n1-n2). I suppose largest should be abs(max(n1,n2)) - what if both are negative? In your code, largest can be negative, in which case diff>largest (as +0.01 > -100.0)