Patch to improve elementary functions for complex and exact arguments
Marc: Here is a diff containing the following changes to _num.scm: new macros: macro-exp-1/2, macro-exp--1/2, macro-log-2 new function: ##flonum.full-precision, which checks that (abs x) is finite and above the largest subnormal number. In ##log: Define exact-log, which gives a finite answer for extremely large and extremely small exact arguments, and better relative accuracy for exact rationals near 1. It is slower than the old behavior of log on exact rationals near 1, so if you just want a quick and dirty approximation, use (log (exact->inexact x)). Define complex-log-magnitude, which computes (log (magnitude x)) carefully in one step instead of in two steps. in ##atan2: Return exact 0 if x and y are exact, x is positive, and y is 0. Carefully scale arguments before handling off to ##flonum.atan so more answers are finite. Remove some checks, put them in calling functions. in atan: add checks before calling ##atan and ##atan2. in ##sqrt: Define complex-sqrt-magnitude, which carefully calculates (sqrt (magnitude x)) in one step instead of two. For complex arguments, return an exact result where possible, e.g., (##sqrt +2i) => 1+i. I've included a test program below, and the results of that test program in old-results and new-results. Many previous +nan.0 and +inf. 0 results are now finite; finite results that differ in more than the last bit or so have more accuracy. Brad
Afficher les réponses par date
Thanks Brad! The patch is not committed to the repository. Note that your tests give slightly different results on my machine (as far as I can see probably only the least significant bit is different). Marc On 30-Oct-08, at 3:53 PM, Bradley Lucier wrote:
Marc:
Here is a diff containing the following changes to _num.scm:
new macros: macro-exp-1/2, macro-exp--1/2, macro-log-2
new function: ##flonum.full-precision, which checks that (abs x) is finite and above the largest subnormal number.
In ##log:
Define exact-log, which gives a finite answer for extremely large and extremely small exact arguments, and better relative accuracy for exact rationals near 1. It is slower than the old behavior of log on exact rationals near 1, so if you just want a quick and dirty approximation, use (log (exact->inexact x)).
Define complex-log-magnitude, which computes (log (magnitude x)) carefully in one step instead of in two steps.
in ##atan2:
Return exact 0 if x and y are exact, x is positive, and y is 0. Carefully scale arguments before handling off to ##flonum.atan so more answers are finite. Remove some checks, put them in calling functions.
in atan:
add checks before calling ##atan and ##atan2.
in ##sqrt:
Define complex-sqrt-magnitude, which carefully calculates (sqrt (magnitude x)) in one step instead of two. For complex arguments, return an exact result where possible, e.g., (##sqrt +2i) => 1+i.
I've included a test program below, and the results of that test program in old-results and new-results. Many previous +nan.0 and +inf.0 results are now finite; finite results that differ in more than the last bit or so have more accuracy.
Brad
<_num.scm-diff><test-log.scm><old-results.gz><new-results.gz>
On Oct 30, 2008, at 5:06 PM, Marc Feeley wrote:
Thanks Brad! The patch is not committed to the repository.
Do you mean "not" or "now"? I tried a "make update" on my linux box for testing and got
frying-pan:~/programs/gambc-v4_3_0-devel> make update D .hgignore D .hgtags Already on branch "master"
make etc.
Note that your tests give slightly different results on my machine (as far as I can see probably only the least significant bit is different).
I reconfigured on my PPC G5 to not use the fused multiply-add instruction, i.e., [descartes:~/programs/gambc-v4_2_8] lucier% ./configure CC='gcc - mcpu=970 -m64 -mno-fused-madd' and got the following differences myself. "new-results" are what I sent you before, "my-results are without the fused multiply-add instruction:
[descartes:~/programs/gambc-v4_2_8] lucier% diff new-results my- results 419c419 < -1.5707963267948966-0.i) ---
-1.5707963267948966-2.2737367544323206e-13i)
421c421 < -1.5707963267948966-0.i) ---
-1.5707963267948966+2.2737367544323206e-13i)
457c457 < .7853981633974483+363.1531368262252i) ---
.7853981633974483+363.15313682622525i)
459c459 < .7853981633974483-363.1531368262252i) ---
.7853981633974483-363.15313682622525i)
495c495 < -.7853981633974483+363.1531368262252i) ---
-.7853981633974483+363.15313682622525i)
497c497 < -.7853981633974483-363.1531368262252i) ---
-.7853981633974483-363.15313682622525i)
1051c1051 < 3295.836866004329-1.5707963267948966i) ---
3295.8368660043293-1.5707963267948966i)
1053c1053 < -725.6131264718905+1.5707963267948966i) ---
-725.6131264718906+1.5707963267948966i)
1055c1055 < -725.6131264718905-1.5707963267948966i) ---
-725.6131264718906-1.5707963267948966i)
1438c1438 < 3295.836866004329+3.141592653589793i) ---
3295.8368660043293+3.141592653589793i)
1476c1476 < -725.6131264718905) ---
-725.6131264718906)
1514c1514 < -725.6131264718905+3.141592653589793i) ---
-725.6131264718906+3.141592653589793i)
So that accounts for a few of the differences. I hope we don't have to disable fused-madd on PPC, it speeds some things up a lot. I believe that all the others are because of the following somewhat disheartening difference. On my x86-64 linux box running Ubuntu 8.4:
frying-pan:~> gsi Gambit v4.3.0
(atan 1) .7853981633974483 (sin (atan 1)) .7071067811865475 (cos (atan 1)) .7071067811865476
On my Mac G5 running 10.5.5:
[descartes:~/programs/gambc-v4_2_8] lucier% gsi/gsi Gambit v4.2.9
(atan 1) .7853981633974483 (sin (atan 1)) .7071067811865476 (cos (atan 1)) .7071067811865475
Which one is correct, I don't know. Perhaps there's a double- rounding going on somewhere. (I'm glad I finally tracked this down.) Brad
On Oct 30, 2008, at 6:16 PM, Bradley Lucier wrote:
I tried a "make update" on my linux box for testing and got
frying-pan:~/programs/gambc-v4_3_0-devel> make update D .hgignore D .hgtags Already on branch "master"
make
etc.
I see from the git browser that you need to do hg pull hg pull --tags once by hand before make update will work in the future; is that right? Brad
On 30-Oct-08, at 9:40 PM, Bradley Lucier wrote:
I see from the git browser that you need to do
hg pull hg pull --tags
once by hand before make update will work in the future; is that right?
This is a problem with the makefiles which was hopefully fixed 2 days ago, but unfortunately after v4.3.0 was released. So if you have v4.3.0 the "make update" won't work by itself. In the future a "make update" (with an installed v4.3.1 or later) should do the right thing. Currently you have to "git pull" and "git fetch --tags" manually before the "make update". Could a git guru check the logic of the "make update" in ./makefile and report any problems? Marc
On Oct 30, 2008, at 6:16 PM, Bradley Lucier wrote:
I believe that all the others are because of the following somewhat disheartening difference.
On my x86-64 linux box running Ubuntu 8.4:
frying-pan:~> gsi Gambit v4.3.0
(atan 1) .7853981633974483 (sin (atan 1)) .7071067811865475 (cos (atan 1)) .7071067811865476
On my Mac G5 running 10.5.5:
[descartes:~/programs/gambc-v4_2_8] lucier% gsi/gsi Gambit v4.2.9
(atan 1) .7853981633974483 (sin (atan 1)) .7071067811865476 (cos (atan 1)) .7071067811865475
Which one is correct, I don't know. Perhaps there's a double- rounding going on somewhere. (I'm glad I finally tracked this down.)
Well, I ran things through my computable reals package (with a quick implementation of computable-sin and computable-cos, so there could be bugs here) and got
[descartes:computation/computational-reals/src] lucier% gsi Gambit v4.3.0
(load "exact-reals.scm") "/Volumes/lucier/text/courses/computation/computational-reals/src/ exact-reals.scm" (computable->inexact (computable-/-by-integer computable-pi 4)) .7853981633974483 (computable->inexact (computable-sin (computable-/-by-integer computable-pi 4))) .7071067811865476 (computable->inexact (computable-sin (exact->computable (inexact- exact (computable->inexact (computable-/-by-integer computable-pi 4)))))) .7071067811865475 (computable->inexact (computable-cos (computable-/-by-integer computable-pi 4))) .7071067811865476 (computable->inexact (computable-cos (exact->computable (inexact- exact (computable->inexact (computable-/-by-integer computable-pi 4)))))) .7071067811865476
So the inexact (double-precision floating-point) approximations to the sine and cosine of exact pi/4 are equal, as one would expect, but the inexact approximations to the sine and cosine of the inexact approximation to pi/4 are not equal (it's obvious that they're not equal to infinite precision, but indeed they round to different floating-point numbers) and it appears that glibc on linux is getting it right and the math library on Mac OS X 10.5.5 is getting it wrong. I'm sure everyone on this list was just dying to hear how this mystery turned out ;-). Brad
On Oct 31, 2008, at 11:44 AM, Bradley Lucier wrote:
Well, I ran things through my computable reals package (with a quick implementation of computable-sin and computable-cos, so there could be bugs here) and got
[descartes:computation/computational-reals/src] lucier% gsi Gambit v4.3.0
(load "exact-reals.scm") "/Volumes/lucier/text/courses/computation/computational-reals/src/ exact-reals.scm" (computable->inexact (computable-/-by-integer computable-pi 4)) .7853981633974483 (computable->inexact (computable-sin (computable-/-by-integer computable-pi 4))) .7071067811865476 (computable->inexact (computable-sin (exact->computable (inexact- exact (computable->inexact (computable-/-by-integer computable-pi 4)))))) .7071067811865475 (computable->inexact (computable-cos (computable-/-by-integer computable-pi 4))) .7071067811865476 (computable->inexact (computable-cos (exact->computable (inexact- exact (computable->inexact (computable-/-by-integer computable-pi 4)))))) .7071067811865476
So the inexact (double-precision floating-point) approximations to the sine and cosine of exact pi/4 are equal, as one would expect, but the inexact approximations to the sine and cosine of the inexact approximation to pi/4 are not equal (it's obvious that they're not equal to infinite precision, but indeed they round to different floating-point numbers) and it appears that glibc on linux is getting it right and the math library on Mac OS X 10.5.5 is getting it wrong.
There's the further test:
(computable-< (exact->computable (inexact->exact . 7853981633974483)) (computable-/-by-integer computable-pi 4)) #t
So the floating point approximation to pi/4 is smaller than pi/4, so its cosine is bigger than its sine. Brad
participants (2)
-
Bradley Lucier -
Marc Feeley