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)).
(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)
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) 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