Newton Method Square Root Iterations

Question

I pulled this code from http://blog.shay.co/newtons-method/:

//a - the number to square root
//times - the number of iterations
public double Sqrt(double a, int times)
{
    if (a < 0)
        throw new Exception("Can not sqrt a negative number");
    double x = 1;
    while (times-- > 0)
        x = x / 2 + a / (2 * x);
    return x;
}

What is a good rule of thumb for the number of iterations on a number, if one exists? (For example, "Use n/2 iterations".)

Answer 1

What is a good rule of thumb for the number of iterations on a number, if one exists?

Newton's method has quadratic convergence, ie. at every step of the algorithm, the number of significant digits in the answer doubles. Thus the algorithm computes square roots upto D digits of precision in O(log D) time. Thus the number of iterations in your loop will depend upon the accuracy expected. So to have a fine grained control over the accuracy of the result, you can add a check in your code that returns the answer when the estimate is not outside the error bounds.

public double Sqrt(double a){
    if (a < 0)
        throw new Exception("Can not sqrt a negative number");
    double error = 0.00001;
    double x = 1;
    while (true){
        double val = x*x;
        if(abs(val-a) <= error)
             return x;
        x = x / 2 + a / (2 * x);
    }
}

Answer 2

If you want a guaranteed relative accuracy with a fixed number of iterations, then prepare the iteration by dividing a by 4 until the result is between 1/2 and 2. Or multiply if starting with a value smaller than 1. Remember the number of divisions, since

sqrt(a)=2^k*sqrt(4^(-k)*a)

requires to multiply the result with a like number of factors 2. Then the relative error for the square root of the reduced a will fall like

3*(1/6)^(2^m)

after m iterations, or faster, according to

(x[m]-sqrt(a))/(x[m]+sqrt(a)) = ( (1-sqrt(a))/(1+sqrt(a)) )^(2^m)

You can extract and manipulate the required powers of 4 and 2 fastest by using the access functions to the floating point exponent, frexp and ldexp in the C math library, similar methods in the Double class in Java.