The sqrt function to end all sqrt functions [RESOLVED]

[sail3005]

ok, i have found this equation, thanks to a link by jim:

An initial, simplistic application of Newton's method on f(x) = x^2-y yields x <= (x + y/x)*0.5. However, the divide is too expensive, which makes this approach unattractive. So this lead me to believe that you could use a newton's method on the reciprocal of x, simultaneously with the newton approximation on x. This leads to (a,b) <= ((a+y*b)/2,b*(2-a*b)), where a is supposed to converge to the square root and b to the reciprocal of the square root. Here is a practical implementation:


#include <stdio.h>

double a,b,t,y;

main() {
int c;

y = 0.1; /* Number to square root */

a = 0.5*(y+1.0); /* Initial guess */
if( a<1.0 ) b = 2.0; /* Initial reciprocal of guess */
else b = 0.5;

for(c=0;c<14;c++) {
a = 0.5*(a+y*b);
t = b*(2.0-a*b);
b = (t>0) ? t : b * 0.5;
printf("(a,b)=(%19.16f,%19.16f)\n",a,b);
}

}

/*

Output:

(a,b)=( 0.3750000000000000, 2.5000000000000000)
(a,b)=( 0.3125000000000000, 3.0468750000000000)
(a,b)=( 0.3085937500000000, 3.2289361953735350)
(a,b)=( 0.3157436847686768, 3.1659195913714560)
(a,b)=( 0.3161678219529112, 3.1628742879915540)
(a,b)=( 0.3162276253760334, 3.1622789545412300)
(a,b)=( 0.3162277604150782, 3.1622777161854930)
(a,b)=( 0.3162277660168138, 3.1622776601686200)
(a,b)=( 0.3162277660168379, 3.1622776601683800)
(a,b)=( 0.3162277660168380, 3.1622776601683790)
(a,b)=( 0.3162277660168380, 3.1622776601683790)
(a,b)=( 0.3162277660168380, 3.1622776601683790)
(a,b)=( 0.3162277660168380, 3.1622776601683790)
(a,b)=( 0.3162277660168380, 3.1622776601683790)

*/



Note odd way that the initial conditions are chosen. There is some work required in the initial condition to make sure the algorithm converges. Ideally, in practical implementations, the inital condition can be worked out by taking the exponent (base 4) of the initial double and divide by two for the "guessed" exponent. On most systems integer to floating point conversions can be done in hardware, and the floating point format's exponent is easy to extract. Further, a table look up with the high bits of the mantissa and the parity of the resultant exponent can help the accuracy of initial guesses as well.

This algorithm roughly doubles in the number of digits of accuracy on every iteration.

i have tested this, and it works perfectly! now, would someone please be kind enough, and tell me how it works?

[sail3005]

ok, i understand it now