[gambit-list] Re: Implementation of read-ieee-float64

Bradley Lucier lucier at math.purdue.edu
Fri Dec 9 14:51:24 EST 2005


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


-------------- next part --------------
A non-text attachment was scrubbed...
Name: expt.scm
Type: application/octet-stream
Size: 2836 bytes
Desc: not available
URL: <http://mailman.iro.umontreal.ca/pipermail/gambit-list/attachments/20051209/235a85cb/attachment.obj>


More information about the Gambit-list mailing list