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