[fpc-pascal] Re: fast sqrt routine

[Thomas Schatzl]

 >Andrew Bennett
 >Sat, 11 Mar 2006 17:24:51 -0800
 >
 >On Sat, 11 Mar 2006 13:44:52 +0100, Pianoman wrote:
 >
 >>Hi, I need to perform very fast squareroot (sqrt) operations on double
 >> type nukmbers.
 >>I tried this
 >>function mysqrt(a:double):double;
 >> [...]
 >>yn:=(y*y+a)/(2*y);
 >>until abs(yn-y) < 10e-16;
 >>mysqrt:=yn;
 >>end;
 >>but it is quite slow (FPC uses Heron iterations) which is 10 times 
 >>faster
 >>than this code but I need even faster sqrt routin.
 >>Any ideas for optimizing the function written above would be 
 >>appreciated.
 >>>Pianoman
 >>Once upon a time, I did this in single precision:
 >
 >1) I made a crude 1st guess by manufacturing, in
 >Assembler, a linear fit mantissa - 2 different fits
 >depending if exponent was odd or even - and the
 >appropriate exponent. The worst case error of this
 >is easily calculated. Process takes 2 floating ops and
 >some integer stuff taking out and putting back the
 >exponent.
 >
 >2) I then explicitly coded the iterations required to
 >converge in this worst case: just 3 floating ops per
 >(or is it 2 plus integer exponent manipulation?)
 >and it didn't take very many. No test. Saving the test
 >saved more time than was wasted doing unnecessary
 >iterations, at least on the machine I did it on. This
 >may not be so in double precision ...

For some reason I recently have been interested in fast sqrt computation 
as well, and found the following (C) pseudo-code, which looks quite fast.

It approximates 1/sqrt(x), which can be done without a divide, just 
multiplies. Maybe it is of help, don't know the Heron iteration thing 
you mention...

(See especially section "B.  sqrt(x) by Reciproot Iteration", 
http://www.netlib.org/fdlibm/e_sqrt.c)

Regards,
   Thomas