5

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;
}
csciXAV_12
  • 309
  • 1
  • 5
  • 13

5 Answers5

6

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.

ale64bit
  • 6,232
  • 3
  • 24
  • 44
  • 2
    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. – M Oehm Jan 15 '15 at 07:26
  • @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`) – leemes Jan 15 '15 at 07:35
  • @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) – ale64bit Jan 15 '15 at 07:36
4

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;
}
Tony Delroy
  • 102,968
  • 15
  • 177
  • 252
  • 1
    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. – M Oehm Jan 15 '15 at 07:35
  • 1
    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 – MSalters Jan 15 '15 at 08:39
  • @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. – Tony Delroy Jan 15 '15 at 10:29
  • Well, the fix would be easy anyway. `sqrt(1/x) = 1/sqrt(x)`, so just add `if (x<1) return 1.0/squareroot(x);`. – MSalters Jan 15 '15 at 13:20
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.

barak manos
  • 29,648
  • 10
  • 62
  • 114
0

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

Álvaro Gómez
  • 306
  • 1
  • 6
  • Thanks, I just updated it to use the fabs() function and I am still getting an infinite loop. – csciXAV_12 Jan 15 '15 at 07:03
  • 4
    Apparently there are a few C++ programmers out there that don't know about `std::abs` overloaded function. – n. m. could be an AI Jan 15 '15 at 07:03
  • @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 ``. – M Oehm Jan 15 '15 at 07:09
  • @csciXAV_12 There is a overload for `std::abs` (and possibly `abs` in header ``. See http://en.cppreference.com/w/cpp/numeric/math/fabs. It is important to post an MCVE, and if possible, use `std::abs` explicitly. – juanchopanza Jan 15 '15 at 07:09
0

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;
}
denis-bu
  • 3,426
  • 1
  • 17
  • 11
  • 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) – MSalters Jan 15 '15 at 08:48