Recursive Newton Square Root Function Only Terminates for Perfect Squares

1

I wrote a program that looks very similar to other recursive newton square root functions I have seen on the web. For some reason this one only works with perfect squares and I can't seem to find a reason why. I have tried passing 3 variables(epsilon), setting x=a, setting a = to the equation passed in line 9 then passing abs(a*a-x). I tried to describe it the best I could it's a slightly new topic for me and I am just not sure if this is only capable of finding perfect roots or if my code/equation is incorrect.

 1  #include <cmath>
 2  #include <iostream>
 3  using namespace std;
 4
 5  double newtroot(double x, double a) {
 6      if (fabs(a*a - x) <= DBL_EPSILON)
 7          return a;
 8      else 
 9          return newtroot(x, fabs(a*a + x)/(2*a));
10  }
11
12  int main() {
13      cout << newtroot(9, 1) << endl;
14      system("pause");
15      return 0;
16  }

EDIT: The function does not only work for perfect squares but it only returns correctly for perfect squares. If it is not a perfect square a ends up being the correct value (checked in the debugger) but the recursion never stops. I figured it must be something with the comparison in line 6 so I tried replacing DBL_EPSILON with a and the incorrect value is returned.
This error also shows on line 6 when a non-perfect square is entered:

Unhandled exception at 0x00007FFE8E9C06F0 (ucrtbased.dll) in RecursionProgrammingExcercisesMurphyT.exe: 0xC00000FD: Stack overflow (parameters: 0x0000000000000001, 0x00000013B2603FE8). occurred

recursion
square-root
perfect-square
asked on Stack Overflow Nov 8, 2017 by Trevor Murphy • edited Feb 18, 2018 by greybeard

1 Answer

0

Not detecting the base case is due to an overly optimistic interpretation:
DBL_EPSILON is not "a value suitable to terminate improving any approximation," but
Difference between 1 and the least value greater than 1 that is representable:
it needs to be scaled to limit relative error.
return (fabs(approx*approx - x) <= 2*x*DBL_EPSILON) ? approx
        : newtroot(x, fabs(approx*approx + x)/(2*approx));
When an approximation is fast enough to tolerate numerical errors, a promising approach to terminate it is to check that the new approximate value is neither the current nor the previous one.

/*! approximate the value of x**.5 */
double newtroot(double x, double approx, double previous) {
    double next = fabs(approx*approx + x)/(2*approx);
    return next == approx || next == previous
        ? approx : newtroot(x, next, approx);
}
answered on Stack Overflow Feb 18, 2018 by greybeard

User contributions licensed under CC BY-SA 3.0