>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