From: kovarik@mcmail.cis.McMaster.CA (Zdislav V. Kovarik)
Newsgroups: comp.lang.fortran,sci.math.num-analysis
Subject: Re: Optimal evaluation of 1/SQRT(x)
Date: 11 Sep 1995 21:37:27 -0400

In article <431uap$bt2@noc.usfca.edu>,
Prof. Loren Meissner <meissner@usfca.edu> wrote:
:>jenhel@stud.unit.no (Jens Helmers) wrote:
[ deleted ]
:>> I want to evaluate 1/SQRT(x) for a very large number of arguments.
[ deleted ]
:>> Is it possible to speed up these evaluations by using a special 
:>>algorithm?
[ deleted ]
:>Special algorithm for doing the reciprocal SqRt all at once? Probably
:>not. 

From the ancient history of numerical algorithms:

Apply Newton's Method to  F(r) = r^(-2) - A  ,  the root being  A^(-1/2):

after simplification,  

                r_(n+1) = (1/2) * r_n * (3 - a * (r_n)^2);

... will converge as long as  0 < r_0 < (3/A)^(1/2)
    i.e. r_0 > 0 and A * (r_0)^2 < 3  (plenty of room),

... the constants to be stored are  (1/2) and 3;

... only addition and multiplication is required (which was an advantage 
    when division was time-consuming compared with multiplication).

Quadratic convergence follows from general theory.

Even if it doesn't help (the co-processor may beat this method by not
having to use the main chip's time), it's of historical interest. 

Remark: you can have a division-free improvement method for finding the 
reciprocal, too, if you apply Newton's method to  F(R) = (1/R) - A and 
simplify. Applied to matrices, we have Hotelling's method for the 
improvement of inverses:  R_(n+1) = R_n * (2 - A * R_n) (this time, order 
of factors is essential since matrices do not commute in general).

Have fun, ZVK (Slavek).