racket icon indicating copy to clipboard operation
racket copied to clipboard

Inaccurate (expt flonum exact-int)

Open gambiteer opened this issue 7 months ago • 14 comments

What version of Racket are you using? Welcome to Racket v8.16.0.2 [cs].

What program did you run? Please include a short example that triggers the bug

Welcome to Racket v8.16.0.2 [cs].
> (expt 1.000000000000001 (expt 2 56))
5.540619075645279e+34
> (expt -1.000000000000001 (expt 2 56))                                                     
5.540619075645279e+34
> (expt 1.000000000000001 (+ 1 (expt 2 56)))                                                
5.5406190756452855e+34
> (expt -1.000000000000001 (+ 1 (expt 2 56)))
-5.5406190756452855e+34

What should have happened?

Roughly (from Gambit):

> (expt 1.000000000000001 (expt 2 56))
5.540622384393264e34
> (expt -1.000000000000001 (expt 2 56))                                                     
5.540622384393264e34
> (expt 1.000000000000001 (+ 1 (expt 2 56)))                                                
5.540622384393264e34
> (expt -1.000000000000001 (+ 1 (expt 2 56)))
-5.540622384393264e34

If you got an error message, please include it here.

Please include any other relevant details e.g., the operating system used or how you are running the code.

The system C library power function should compute this very accurately, (log 1.000000000000001) is easily computed to machine accuracy, multiplication by (expt 2 56) is exact, and then calculating exp is also very accurate:

> (exp (* (log 1.000000000000001) (expt 2 56)))
5.540622384393274e34
> (flexpt 1.000000000000001 (inexact (expt 2 56)))
5.540622384393264e34

The iterative method used in Racket's library is very inaccurate for this example (getting only five correct digits in the result). One can compute the correctly rounded results (here with my computable reals library, but also with MPFR, etc.):

> (computable->inexact (computable-expt (->computable 1.000000000000001) (expt 2 56))) 
5.540622384393264e34
> (computable->inexact (computable-expt (->computable 1.000000000000001) (+ 1 (expt 2 56))))
5.54062238439327e34

For (expt flonum exact-int) I would recommend computing (flonum-expt (flabs flonum) (->inexact exact-int)) and then changing the sign of the result if flonum is negative and exact-int is odd.

gambiteer avatar Mar 16 '25 00:03 gambiteer