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-324etc. 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
Attachment:
expt.scm
Description: Binary data