fast sqrt routine

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

fast sqrt routine

pianoman-3
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
----- Original Message -----
From: <[hidden email]>
To: <[hidden email]>
Sent: Saturday, March 11, 2006 12:00 PM
Subject: fpc-pascal Digest, Vol 19, Issue 17


> Send fpc-pascal mailing list submissions to
> [hidden email]
>
> To subscribe or unsubscribe via the World Wide Web, visit
> http://lists.freepascal.org/mailman/listinfo/fpc-pascal
> or, via email, send a message with subject or body 'help' to
> [hidden email]
>
> You can reach the person managing the list at
> [hidden email]
>
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of fpc-pascal digest..."
>
>
> Today's Topics:
>
>    1.  problem setting up serial port (Marc Santhoff)
>    2. Re:  problem setting up serial port (Sebastian G?nther)
>    3. Re:  problem setting up serial port (Michael Van Canneyt)
>    4. Re:  problem setting up serial port (Marc Santhoff)
>    5. Re:  problem setting up serial port (Marc Santhoff)
>    6.  makeskel --update (Vincent Snijders)
>    7. Re:  makeskel --update (Michael Van Canneyt)
>    8. Re:  problem setting up serial port (Marc Santhoff)
>
>
> ----------------------------------------------------------------------
>
> Message: 1
> Date: Fri, 10 Mar 2006 15:51:46 +0100
> From: Marc Santhoff <[hidden email]>
> Subject: [fpc-pascal] problem setting up serial port
> To: FPC-Pascal users discussions <[hidden email]>
> Message-ID: <[hidden email]>
> Content-Type: text/plain
>
> Hi,
>
> I'm trying to use a serial port at 1200baud, 7N2 with this code
> sequence:
>
> <snip>
> fCom := fpOpen('/dev/cuaa0', O_RDWR OR O_NOCTTY );//OR O_NONBLOCK);
> if (fCom < 0) then begin
> writeln(stderr, 'Couldn''t open serial port.');
> halt(1);
> end;
>
> r := 0;
> fillchar(tios, sizeof(tios), #0);
>
> cfsetispeed(tios, B1200);
> cfsetospeed(tios, B1200);
>
> tios.c_cflag := CREAD or CLOCAL or CS7 or CSTOPB;
> tios.c_oflag := 0;
> tios.c_iflag := IGNBRK OR IGNPAR;
> tios.c_lflag := 0;
> r := tcsetattr(fCom, TCSANOW, tios);
> if (r = -1) then begin
> writeln(stderr, 'tcsetattr failed!');
> fpClose(fCom);
> halt(1);
> end;
> </snip>
>
> and all I get when running it is: 'tcsetattr failed!'. The rights are
> okay and running as root does make no difference. The hardware is
> working, connections are checked okay ... another program doing mostly
> the same is working as it should.
>
> Maybe I have looked at it for too long, where is the problem?
>
> TIA,
> Marc
>
>
>
> ------------------------------
>
> Message: 2
> Date: Fri, 10 Mar 2006 16:46:06 +0100
> From: Sebastian G?nther <[hidden email]>
> Subject: Re: [fpc-pascal] problem setting up serial port
> To: FPC-Pascal users discussions <[hidden email]>
> Message-ID: <[hidden email]>
> Content-Type: text/plain; charset=ISO-8859-15; format=flowed
>
> Marc Santhoff schrieb:
> >
> > Maybe I have looked at it for too long, where is the problem?
>
> just by the way, there is a unit called 'serial' in Free Pascal. Perhaps
> you might want to use this one, or at least have a look at the source
> code...
>
>
> Regards,
> Sebastian
>
> ------------------------------
>
> Message: 3
> Date: Fri, 10 Mar 2006 17:06:37 +0100 (Romance Standard Time)
> From: Michael Van Canneyt <[hidden email]>
> Subject: Re: [fpc-pascal] problem setting up serial port
> To: FPC-Pascal users discussions <[hidden email]>
> Message-ID: <[hidden email]>
> Content-Type: text/plain; charset="iso-8859-1"
>
>
>
> On Fri, 10 Mar 2006, Sebastian Günther wrote:
>
> > Marc Santhoff schrieb:
> >>
> >> Maybe I have looked at it for too long, where is the problem?
> >
> > just by the way, there is a unit called 'serial' in Free Pascal. Perhaps
you
> > might want to use this one, or at least have a look at the source
code...

>
> Better yet, the german Magazine Toolbox distributes a unit by Rainer
> Reush, with a  'TSerial' component.
> AFAIK: Cross-platform and cross-delphi-lazarus...
>
> Michael.
>
> ------------------------------
>
> Message: 4
> Date: Fri, 10 Mar 2006 17:11:32 +0100
> From: Marc Santhoff <[hidden email]>
> Subject: Re: [fpc-pascal] problem setting up serial port
> To: [hidden email]
> Message-ID: <[hidden email]>
> Content-Type: text/plain; charset=ISO-8859-15
>
> Am Freitag, den 10.03.2006, 16:46 +0100 schrieb Sebastian Günther:
> > Marc Santhoff schrieb:
> > >
> > > Maybe I have looked at it for too long, where is the problem?
> >
> > just by the way, there is a unit called 'serial' in Free Pascal. Perhaps
> > you might want to use this one, or at least have a look at the source
> > code...
>
> Yes, I know and used it at first, but since it does not give result
> codes back I had to rip it up to create a test program for detecting why
> it doesn't work.
>
> Thank you anyways,
> Marc
>
>
>
> ------------------------------
>
> Message: 5
> Date: Fri, 10 Mar 2006 21:01:46 +0100
> From: Marc Santhoff <[hidden email]>
> Subject: Re: [fpc-pascal] problem setting up serial port
> To: [hidden email]
> Message-ID: <[hidden email]>
> Content-Type: text/plain; charset=ISO-8859-15
>
> Am Freitag, den 10.03.2006, 17:06 +0100 schrieb Michael Van Canneyt:
> >
> > On Fri, 10 Mar 2006, Sebastian Günther wrote:
> >
> > > Marc Santhoff schrieb:
> > >>
> > >> Maybe I have looked at it for too long, where is the problem?
> > >
> > > just by the way, there is a unit called 'serial' in Free Pascal.
Perhaps you
> > > might want to use this one, or at least have a look at the source
code...

> >
> > Better yet, the german Magazine Toolbox distributes a unit by Rainer
> > Reush, with a  'TSerial' component.
> > AFAIK: Cross-platform and cross-delphi-lazarus...
>
> Yes, but it is not open source, I'd need something under BSD or LGPL
> style licence.
>
> And my code should work ... ;)
>
> Marc
>
>
>
> ------------------------------
>
> Message: 6
> Date: Fri, 10 Mar 2006 22:18:24 +0100
> From: Vincent Snijders <[hidden email]>
> Subject: [fpc-pascal] makeskel --update
> To: FPC-Pascal users discussions <[hidden email]>
> Message-ID: <[hidden email]>
> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>
> Hi,
>
> If I use the makeskel --update parameter, what is the best way to merge
> the resulting file to the existing file?
>
> Vincent
>
> ------------------------------
>
> Message: 7
> Date: Fri, 10 Mar 2006 22:50:20 +0100 (Romance Standard Time)
> From: Michael Van Canneyt <[hidden email]>
> Subject: Re: [fpc-pascal] makeskel --update
> To: FPC-Pascal users discussions <[hidden email]>
> Message-ID: <Pine.WNT.4.58.0603102250040.1376@home-office>
> Content-Type: TEXT/PLAIN; charset=US-ASCII
>
>
>
> On Fri, 10 Mar 2006, Vincent Snijders wrote:
>
> > Hi,
> >
> > If I use the makeskel --update parameter, what is the best way to merge
> > the resulting file to the existing file?
>
> Copy and paste.
>
> Michael.
>
> ------------------------------
>
> Message: 8
> Date: Sat, 11 Mar 2006 00:51:52 +0100
> From: Marc Santhoff <[hidden email]>
> Subject: Re: [fpc-pascal] problem setting up serial port
> To: [hidden email]
> Message-ID: <[hidden email]>
> Content-Type: text/plain
>
> Hi again,
>
> from trying and decoding error messages I know the call to tcsetattr
> fails with an invalid input. The problem is caused by trying to set the
> field termios.c_ospeed, if left out the call is successfull. So I went
> the way trough the sources and have found something.
>
> On FreeBSD 4.11 termios is defined like this:
>
> /usr/include/termios.h -> sys/termios.h
>
> struct termios {
> tcflag_t c_iflag; /* input flags */
> tcflag_t c_oflag; /* output flags */
> tcflag_t c_cflag; /* control flags */
> tcflag_t c_lflag; /* local flags */
> cc_t c_cc[NCCS]; /* control chars */
> speed_t c_ispeed; /* input speed */
> speed_t c_ospeed; /* output speed */
> };
>
>
> The fpc source file "" defines it as:
>
> <fpc-2.0.2>/share/src/fpc-2.0.2/rtl/freebsd/termios.inc
>
> type
>   Termios = packed record
>     c_iflag,
>     c_oflag,
>     c_cflag,
>     c_lflag  : longint;
>     c_line   : char;
>     c_cc     : array[0..NCCS-1] of byte;
>    {$IFDEF BSD}
>     c_ispeed,
>     c_ospeed : longint;
>    {$endif}
>   end;
>   TTermios=Termios;
>
> Note the line in pascal saying "c_line: char" which is missing in the
> systems declaration. This is an error, insn't it?
>
> Marc
>
>
>
> ------------------------------
>
> _______________________________________________
> fpc-pascal maillist  -  [hidden email]
> http://lists.freepascal.org/mailman/listinfo/fpc-pascal
>
> End of fpc-pascal Digest, Vol 19, Issue 17
> ******************************************
>
>
> __________ Informacia od NOD32 1.1438 (20060310) __________
>
> Tato sprava bola preverena antivirusovym systemom NOD32.
> http://www.eset.sk
>
>

_______________________________________________
fpc-pascal maillist  -  [hidden email]
http://lists.freepascal.org/mailman/listinfo/fpc-pascal
Reply | Threaded
Open this post in threaded view
|

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