[fpc-pascal] Re: fast sqrt routine

[Andrew Bennett]

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;
>var y,yn:double;
>begin
>yn:=a;
>repeat
>y:=yn;
>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 ...

Best of luck.

Andrew Bennett, Avondale Vineyard, Nova Scotia.