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