In the implementation of the withdrawn SRFI-56, read-ieee-float64 has the following code
(define (mantissa expn b2 b3 b4 b5 b6 b7 b8) (case expn ; recognize special literal exponents ((255) #f) ; won't handle NaN and +/- Inf ((0) ; denormalized (exact->inexact (* (expt 2.0 (- 1 (+ 1023 52))) (combine b2 b3 b4 b5 b6 b7 b8)))) (else (exact->inexact (* (expt 2.0 (- expn (+ 1023 52))) (combine (+ b2 16) b3 b4 b5 b6 b7 b8))))))
Gambit 4.0b15 fails some srfi-56 tests because the first instance of 2.0 is not 2; unfortunately, Gambit-C calculates
(expt 2.0 -1074)
as
(/ (expt 2.0 1074))
Since the exponent range of IEEE arithmetic is not symmetric, (expt 2.0 1074) is +inf.0, yet (expt 2.0 -1074) is not 0.
If both instances of 2.0 are changed to 2, then Gambit passes all the tests. I think that this is what you may have meant in the first place, because you take the trouble to call exact->inexact on what is already an inexact number.
Brad
Afficher les réponses par date
On 12/7/05, Bradley Lucier lucier@math.purdue.edu wrote:
Gambit 4.0b15 fails some srfi-56 tests because the first instance of 2.0 is not 2; unfortunately, Gambit-C calculates
(expt 2.0 -1074)
as
(/ (expt 2.0 1074))
Hi, thanks for looking into this and getting it working with Gambit. The (expt 2.0 ...) form was there for the sake of Bigloo, for which
(expt N -M) => 0
for all exact integers N, M > 1. This is clearly a bug on Bigloo's part, but it seemed safe to compensate for it. In particular, every other tested implementation correctly evaluates (expt 2.0 -1024).
[Or so I thought. Scheme48 also gives 0.0 for this.]
Choosing between the lesser of two evil^H^H^H^Hbugs, I'd consider Bigloo's behavior the more egregious here, so let's call (expt 2 ...) the correct reference form.
Congratulations, Gambit now passes all tests!
[And possibly the latest version of Scheme48 does as well.]
-- Alex
On Dec 8, 2005, at 9:12 PM, Alex Shinn wrote:
On 12/7/05, Bradley Lucier lucier@math.purdue.edu wrote:
Gambit 4.0b15 fails some srfi-56 tests because the first instance of 2.0 is not 2; unfortunately, Gambit-C calculates
(expt 2.0 -1074)
as
(/ (expt 2.0 1074))
Hi, thanks for looking into this and getting it working with Gambit. The (expt 2.0 ...) form was there for the sake of Bigloo, for which
(expt N -M) => 0
for all exact integers N, M > 1.
This may not be unreasonable for an implementation without rationals, like bigloo, but perhaps it would be better to do an inexact expt in this case.
Choosing between the lesser of two evil^H^H^H^Hbugs,
I'm working on a new expt that trusts libm's pow function more ...
Brad
On 9-Dec-05, at 12:21 AM, Bradley Lucier wrote:
On Dec 8, 2005, at 9:12 PM, Alex Shinn wrote:
On 12/7/05, Bradley Lucier lucier@math.purdue.edu wrote:
Gambit 4.0b15 fails some srfi-56 tests because the first instance of 2.0 is not 2; unfortunately, Gambit-C calculates
(expt 2.0 -1074)
as
(/ (expt 2.0 1074))
Hi, thanks for looking into this and getting it working with Gambit. The (expt 2.0 ...) form was there for the sake of Bigloo, for which
(expt N -M) => 0
for all exact integers N, M > 1.
This may not be unreasonable for an implementation without rationals, like bigloo, but perhaps it would be better to do an inexact expt in this case.
Choosing between the lesser of two evil^H^H^H^Hbugs,
I'm working on a new expt that trusts libm's pow function more ...
Brad
Brad, would it be sound to compute (expt X -Y) as (expt (/ X) Y) ? This would solve the denormalized result case but I'm wondering how this compares precision-wise to (/ (expt X Y)).
Marc
On Dec 9, 2005, at 8:27 AM, Marc Feeley wrote:
Brad, would it be sound to compute (expt X -Y) as (expt (/ X) Y) ? This would solve the denormalized result case but I'm wondering how this compares precision-wise to (/ (expt X Y)).
Neither is great if X or Y is a floating-point number and the other is real. We really need to use the built-in pow, which is more accurate than what we can do with the repeated squaring algorithm.
I include *UNTESTED* code below; since I have a cold today I won't test it soon.
For example, with this code you get (expt is the current expt)
(compile-file "expt.scm")
#t
(load "expt")
"/Users/lucier/programs/gambc40b15/expt/expt.o13"
(exp (* (+ (expt 2 54) 7) (log (- 1 (expt 2. -53)))))
.13533528323661256
(expt (- 1 (expt 2. -53)) (+ (expt 2 54) 7))
.13533528122258084
(##my-expt (- 1 (expt 2. -53)) (+ (expt 2 54) 7))
.1353352832366126
(expt 2. -1074)
0.
(##my-expt 2. -1074)
5e-324
etc. So the way you do expt when y is an integer can lead to really big relative errors, even away from the underflow threshold.
Brad
Marc:
I bit the bullet and did the optimal expt, you'll be amazed. It still is not tested extensively, but at least the code is straightforward. So now we have things like:
(##my-expt 8/27 1/3)
2/3
In the process, I discovered that ##exact-int.nth-root was buggy as hell and fixed it (and tested it). I trust it now.
The new routines begin with ##my-...
If other people want to test ##my-expt, that would help. Put it in a sibling directory to the lib directory in the gambc40b15 directory to compile it.
I removed srfi-56 and Alex Shinn from the CC list ;-).
Brad
On 9-Dec-05, at 9:39 PM, Bradley Lucier wrote:
Marc:
I bit the bullet and did the optimal expt, you'll be amazed.
Hey... I'm always amazed when I read your code ;-)
It still is not tested extensively, but at least the code is straightforward. So now we have things like:
(##my-expt 8/27 1/3)
2/3
In the process, I discovered that ##exact-int.nth-root was buggy as hell and fixed it (and tested it). I trust it now.
The new routines begin with ##my-...
If other people want to test ##my-expt, that would help. Put it in a sibling directory to the lib directory in the gambc40b15 directory to compile it.
I removed srfi-56 and Alex Shinn from the CC list ;-).
Brad
<expt.scm>
Great! In the meantime I will add "pow" to gambit.h so that ##flonum.pow can be inlined. Let me know when your testing is done.
Marc
On Dec 9, 2005, at 9:16 PM, Marc Feeley wrote:
Let me know when your testing is done.
Well, I don't know if it's done, but I'm done with it.
Just do a code review, clean up the two routines to match your coding standards, change ##my-whatever to ##whatever, and install it. If there are errors, they're too subtle for me to find right now. I did test that the new ratnum exponent code works and that for a lot of random arguments the new expt has similar accuracy to the old expt.
Brad
On Dec 9, 2005, at 9:16 PM, Marc Feeley wrote:
Let me know when your testing is done.
Well, I don't know if it's done, but I'm done with it.
Just do a code review, clean up the two routines to match your coding standards, change ##my-whatever to ##whatever, and install it. If there are errors, they're too subtle for me to find right now. I did test that the new ratnum exponent code works and that for a lot of random arguments the new expt has similar accuracy to the old expt.
Brad