Re: fast sqrt routine

classic Classic list List threaded Threaded
1 message Options
Reply | Threaded
Open this post in threaded view
|

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
_______________________________________________
fpc-pascal maillist  -  [hidden email]
http://lists.freepascal.org/mailman/listinfo/fpc-pascal