Hi Leonard, On Wednesday, 2008-08-06 00:47:00 +0300, Leonard Mada wrote:
>> It seems to be one of the most difficult functions. To get an idea of >> that function, have a look at the file toms708.c of the R-project. > > There is some hidden beauty in that file. ;) > > Without being of much help, I just want to draw attention to some > features I spotted during a peak in the code: > - often, (1-x) is decomposed into more stable terms, > - from the simple 0.5 - x -0.5 > [and x - 0.5 - 0.5 (for x-1)], > - to more complex terms like: > return w * (0.5 - 1.0 / w + 0.5); > [this is actually w - 1] Indeed, nice tricks. Though the second one also already has some cancellation error for w=0.99999 and the result (in C/C++) bitwise is not identical to w-1, but that's after the 20th decimal digit in this case, neglectable. > The code was well written by numerical analysts, and optimised for > accuracy over a large domain of values. Is it feasible to use this > library Unfortunately not, R is licensed under GPL. > or a similar library (if one is found)? The license would have to be compatible with LGPLv3 and pass Sun Legal department. > Rewriting this code isn't going to be easy. I'd say for a non-totally-expert on that field a true rewrite is nearly impossible. Eike -- PGP/OpenPGP/GnuPG encrypted mail preferred in all private communication. Key ID: 0x293C05FD - 997A 4C60 CE41 0149 0DB3 9E96 2F1A D073 293C 05FD
pgp0XBukXGCia.pgp
Description: PGP signature
