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
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);
}
User contributions licensed under CC BY-SA 3.0