Chez Scheme: Inaccurate (expt flonum exact-int)

[gambiteer]

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]

I chose these examples from

racket/pkgs/racket-test-core/tests/racket/number.rktl

Line 642 in d111b80

(test 5.540619075645279e+34 expt 1.000000000000001 (expt 2 56))

[soegaard]

@gambiteer

I know @bdeket has worked on complex expt recently (March 12).
Is v8.16.0.2 before or after this date?

[gambiteer]

One gets the same results with git mainline head.

[bdeket]

Changes from 12 march where only about fllonum-flonum operations.
Which are now improved:

#lang racket/base

(require math/bigfloat
         math/flonum)
(bf-precision 1024)

(expt 1.000000000000001 (expt 2 56))
(expt 1.000000000000001 (fl (expt 2 56)))
(bigfloat->flonum (bfexpt (bf 1.000000000000001) (bf (expt 2 56))))
;5.540619075645279e+34
;5.540622384393264e+34
;5.540622384393264e+34

[gambiteer]

Here's a way that might compute (expt flonum exact-int) more accurately when (inexact exact-int) does not equal inexact-int.
Here's a test program, run under Gambit with my computable reals package and the C function nextafter:

(do ((i 0 (fx+ i 1)))
    ((fx> i 16))
  (let ((arg (+ (expt 2 56) i)))
    (let ((gambit
           (expt 1.000000000000001 arg))
          (correct
           (computable->inexact (computable-expt (->computable 1.000000000000001) arg)))
          (new
           (let* ((inexact-arg (let ((possible-result (inexact arg)))
                                 (if (<= possible-result arg)
                                     possible-result
                                     (nextafter possible-result 0.)))) ;; next lower flonum
                  (rest-arg (inexact (- arg (exact inexact-arg)))))    ;; nonnegative flonum
             (* (flexpt 1.000000000000001 inexact-arg)
                (flexpt 1.000000000000001 rest-arg)))))
      (for-each display
                (list `(expt 1.000000000000001 ,arg) #\space
                      gambit: #\space
                      gambit #\space
                      correct: #\space
                      correct #\tab
                      new: #\space
                      new
                      #\newline)))))

with output

> (load "expt-test2")
(expt 1.000000000000001 72057594037927936) gambit 5.540622384393264e34 correct 5.540622384393264e34   new 5.540622384393264e34
(expt 1.000000000000001 72057594037927937) gambit 5.540622384393264e34 correct 5.54062238439327e34      new 5.54062238439327e34
(expt 1.000000000000001 72057594037927938) gambit 5.540622384393264e34 correct 5.5406223843932765e34    new 5.540622384393276e34
(expt 1.000000000000001 72057594037927939) gambit 5.540622384393264e34 correct 5.540622384393282e34     new 5.540622384393282e34
(expt 1.000000000000001 72057594037927940) gambit 5.540622384393264e34 correct 5.5406223843932885e34    new 5.5406223843932885e34
(expt 1.000000000000001 72057594037927941) gambit 5.540622384393264e34 correct 5.540622384393295e34     new 5.540622384393294e34
(expt 1.000000000000001 72057594037927942) gambit 5.540622384393264e34 correct 5.5406223843933005e34    new 5.5406223843933005e34
(expt 1.000000000000001 72057594037927943) gambit 5.540622384393264e34 correct 5.540622384393307e34     new 5.540622384393307e34
(expt 1.000000000000001 72057594037927944) gambit 5.540622384393264e34 correct 5.540622384393313e34     new 5.540622384393312e34
(expt 1.000000000000001 72057594037927945) gambit 5.540622384393362e34 correct 5.540622384393319e34     new 5.540622384393319e34
(expt 1.000000000000001 72057594037927946) gambit 5.540622384393362e34 correct 5.540622384393325e34     new 5.540622384393325e34
(expt 1.000000000000001 72057594037927947) gambit 5.540622384393362e34 correct 5.540622384393332e34     new 5.540622384393331e34
(expt 1.000000000000001 72057594037927948) gambit 5.540622384393362e34 correct 5.540622384393337e34     new 5.540622384393337e34
(expt 1.000000000000001 72057594037927949) gambit 5.540622384393362e34 correct 5.540622384393344e34     new 5.540622384393344e34
(expt 1.000000000000001 72057594037927950) gambit 5.540622384393362e34 correct 5.54062238439335e34      new 5.540622384393349e34
(expt 1.000000000000001 72057594037927951) gambit 5.540622384393362e34 correct 5.540622384393357e34     new 5.540622384393356e34
(expt 1.000000000000001 72057594037927952) gambit 5.540622384393362e34 correct 5.540622384393362e34     new 5.540622384393362e34

This might be the best way to go.

[soegaard]

Results from Chez Scheme:

Chez Scheme Version 10.1.0-pre-release.3
Copyright 1984-2024 Cisco Systems, Inc.

> (expt 1.000000000000001 (expt 2 56))
5.540619075645279e34
> (expt -1.000000000000001 (expt 2 56))
5.540619075645279e34
> (expt 1.000000000000001 (+ 1 (expt 2 56))) 
5.5406190756452855e34
> (expt -1.000000000000001 (+ 1 (expt 2 56)))
-5.5406190756452855e34

This matches the initial example.
So I think, we need to look at expt in the Chez Scheme source.

[gambiteer]

The code that computes this case of expt is here:

racket/racket/src/ChezScheme/s/5_3.ss

Line 1717 in d111b80

[else

[soegaard]

For easier experimentation, I have made a version of expt that runs in plain Racket.

#lang racket
(require (rename-in racket/flonum
                    [flexpt      $flexpt])
         (rename-in racket/base
                    [inexact? $inexactnum?])
         (only-in rnrs/base-6 inexact)
         (only-in racket/fixnum fx= fx- fx<=)
         (only-in math/flonum flsgn))

(define (bignum? x)      (and (integer? x) (not (fixnum? x))))
(define (ratnum? x)      (and (number? x) (not (flonum? x)) (not (integer? x))))
(define ($flinteger? x)  (and (flonum? x) (integer? x)))
(define ($exactnum? x)   (and (number? x) (exact? x)))
(define ($nan? x)        (nan? x))
(define ($flonum-sign x) (case (flsgn x) [(0.0) 0] [(1.0) 1] [(-1.0) -1]))
(define (ash int count) (arithmetic-shift int count))

(define $inexactnum-real-part real-part)

(define (exact-integer-fits-float? x) #f)  ; TODO

(define-syntax type-case
  (syntax-rules (else)
    [(_ expr
        [(pred1 pred2 ...) e1 e2 ...] ...
        [else ee1 ee2 ...])
     (let ([t expr])
       (cond
         [(or (pred1 t) (pred2 t) ...) e1 e2 ...]
         ...
         [else ee1 ee2 ...]))]))

(define max-float-exponent  1023)
(define min-float-exponent -1023)

(define floatable?
  ; x is "floatable" if it can be made inexact without overflow or underflow
  (lambda (x)
    (type-case x
               [(fixnum?) #t]
               [(bignum?) (fx<= (integer-length x) (constant max-float-exponent))]
               [(ratnum?) (fx<= (constant min-float-exponent)
                                (fx- (integer-length (numerator x))
                                     (integer-length (denominator x)))
                                (constant max-float-exponent))]
               [($exactnum?) (and (floatable? (real-part x))
                                  (floatable? (imag-part x)))]
               [else #t])))


(define-syntax $impoops        (syntax-rules () [(_ expr ...)  (error expr ...)]))
(define-syntax nonnumber-error (syntax-rules () [(_ name expr) (error "not a number" expr)]))
(define-syntax constant        (syntax-rules () [(_ value) value]))
(define-syntax negated-flonum? (syntax-rules () [(_ x) (fx= ($flonum-sign x) 1)]))


(define expt
   (lambda (x y)
      (type-case y
         [(fixnum? bignum?)
          (cond
            [(and (eq? y 0) (number? x)) 1]
            [(eq? x 0)
             (if (< y 0)
                 ($impoops 'expt "undefined for values ~s and ~s" x y)
                 0)]
            [(eq? x 1) 1]
            [(eq? x 2) (if (< y 0) (/ (ash 1 (- y))) (ash 1 y))]
            [(and (flonum? x) (exact-integer-fits-float? y)) ($flexpt x (inexact y))]
            [(and ($inexactnum? x) (exact-integer-fits-float? y)) (exp (* y (log x)))]
            [(not (number? x)) (nonnumber-error 'expt x)]
            [(ratnum? x)
             (if (< y 0)
                 (let ([y (- y)])
                   (/ (expt (denominator x) y) (expt (numerator x) y)))
                 (/ (expt (numerator x) y) (expt (denominator x) y)))]
            [else
             (displayln "this case")
             (let ()
               (define (f x n)
                 (displayln (list (list 'x x) (list 'n n)))
                 (let loop ([i (integer-length n)] [a 1])
                   (displayln (list (list 'i i) (list 'a a)))
                   (let ([a (if (bitwise-bit-set? n i)
                                (* a x)
                                a)])
                     (if (fx= i 0)
                         a
                         (loop (fx- i 1)
                               (* a a))))))
               (if (< y 0)
                   (if (or (fixnum? x) (bignum? x))
                       (/ (f x (- y)))
                       (f (/ x) (- y)))
                   (f x y)))])]
         [(flonum?)
          (type-case x
             [(flonum?)
              (if (and (fl< x 0.0) (not ($flinteger? y)))
                  (exp (* y (log x)))
                  ($flexpt x y))]
             [($inexactnum? $exactnum?) (exp (* y (log x)))]
             [(fixnum? bignum? ratnum?)
              (cond
               [(eq? x 0)
                (cond
                 [(fl< y 0.0) ($impoops 'expt "undefined for values ~s and ~s" x y)]
                 [(fl= y 0.0) 1.0]
                 [($nan? y) +nan.0]
                 [else 0])]
               [(eq? x 1) 1]
               [else
                (if (floatable? x)
                    (expt (inexact x) y)
                    (exp (* y (log x))))])]
             [else (nonnumber-error 'expt x)])]
         [($inexactnum?)
          (cond
           [(eq? x 0)
            (let ([r ($inexactnum-real-part y)])
              (cond
               [(fl> r 0.0) 0]
               [else
                ($impoops 'expt "undefined for values ~s and ~s" x y)]))]
           [(eq? x 1) 1]
           [(and (flonum? x) (fl= x 0.0) (not (negated-flonum? ($inexactnum-real-part y))))
            0.0]
           [else
            (unless (number? x) (nonnumber-error 'expt x))
            (exp (* y (log x)))])]
         [(ratnum? $exactnum?)
          (unless (number? x) (nonnumber-error 'expt x))
          (cond
             [(eqv? y 1/2) (sqrt x)]
             [(eq? x 0)
              (if (> (real-part y) 0)
                  0
                  ($impoops 'expt "undefined for values ~s and ~s" x y))]
             [(eq? x 1) 1]
             [(and (floatable? y)
                   (let ([y (inexact y)])
                     ;; Don't use this case if `(inexact y)` loses
                     ;; precision and becomes an an integer, in which
                     ;; case the result would be real (but should be
                     ;; non-real complex)
                     (and (not (and (flonum? y)
                                    ($flinteger? y)))
                          y)))
              => (lambda (y) (expt x y))]
             [else (exp (* y (log x)))])]
         [else (nonnumber-error 'expt y)])))

;     new expt                                     old expt               gambit
(list (expt  1.000000000000001 (expt 2 56))        5.540619075645279e+34  5.540622384393264e34)
(list (expt -1.000000000000001 (expt 2 56))        5.540619075645279e+34  5.540622384393264e34)
(list (expt  1.000000000000001 (+ 1 (expt 2 56)))  5.5406190756452855e+34 5.540622384393264e34)
(list (expt -1.000000000000001 (+ 1 (expt 2 56))) -5.5406190756452855e+34 -5.540622384393264e34)

[bdeket]

Some extra boundaries on the problem:
For (< (abs y) (expt 2 53)) the flonum-flonum path is taken
For (< (log +min.0 (flprev 1.)) (abs y)) the result will be 0. ±1. or ±inf.0 **
If y falls within these bounds x must be within:
(< (expt +min.0 (/ 1 (expt 2 53))) (abs x) (expt +max.0 (/ 1 (expt 2 53)))) = (< 0.9999999999999174 (abs x) 1.0000000000000788)
otherwise the result would again be 0. or ±inf.0

Splitting the result in the product of 2 flonum operations is in this region a lot better than the loop. Generally flulp-error < 1, but it can go up:
(expt 0.999999999999923 9113206036900121) has a flulp-error of 690 (9 bits). But the current implementation has an error of 34 bits.

** in practice upper bound is better (inexact->exact (+ (log +min.0 (flprev 1.)) 1e16)) due to rounding near zero, but all the extra results will be +min.0 (when (= x (flprev 1.))).
the upperbound for (< 1. x) of (log +max.0 (flnext 1.)) does not have the same problem

[gambiteer]

Interesting analysis. I don't quite understand:

Splitting the result in the product of 2 flonum operations is in this region a lot better than the loop. Generally flulp-error < 1, but it can go up:
(expt 0.999999999999923 9113206036900121) has a flulp-error of 690 (9 bits). But the current implementation has an error of 34 bits.

I did some experiments around these example args with the following test program. It appears that writing the exponent as a sum of integers that can both be converted exactly to flonums leads to very accurate answers near these example arguments. Test program below.

(define (test-range x y min-range max-range)
  (do ((i min-range (fx+ i 1)))
      ((fx> i max-range))
    (let ((arg (+ y i)))
      (let ((gambit
             (expt x arg))
            (correct
             (computable->inexact (computable-expt (->computable x) arg)))
            (new
             (let* ((inexact-arg (let ((possible-result (inexact arg)))
                                   (if (<= possible-result arg)
                                       possible-result
                                       (nextafter possible-result 0.)))) ;; next lower flonum
                    (rest-arg (inexact (- arg (exact inexact-arg)))))    ;; nonnegative flonum
               (* (flexpt x inexact-arg)
                  (flexpt x rest-arg)))))
        (for-each display
                  (list `(expt ,x ,arg) #\space
                        gambit: #\space
                        gambit #\tab
                        correct: #\space
                        correct #\tab
                        new: #\space
                        new
                        #\newline))))))

(test-range 0.999999999999923 9113206036900121 -4 4)

Result:

> (load "expt-test2")
(expt .999999999999923 9113206036900117) gambit 1.1282707411035222e-305 correct 1.1282707411034352e-305 new 1.1282707411034353e-305
(expt .999999999999923 9113206036900118) gambit 1.1282707411033483e-305 correct 1.1282707411033483e-305 new 1.1282707411033483e-305
(expt .999999999999923 9113206036900119) gambit 1.1282707411031744e-305 correct 1.1282707411032614e-305 new 1.1282707411032614e-305
(expt .999999999999923 9113206036900120) gambit 1.1282707411031744e-305 correct 1.1282707411031744e-305 new 1.1282707411031744e-305
(expt .999999999999923 9113206036900121) gambit 1.1282707411031744e-305 correct 1.1282707411030875e-305 new 1.1282707411030875e-305
(expt .999999999999923 9113206036900122) gambit 1.1282707411030006e-305 correct 1.1282707411030006e-305 new 1.1282707411030006e-305
(expt .999999999999923 9113206036900123) gambit 1.1282707411028267e-305 correct 1.1282707411029136e-305 new 1.1282707411029137e-305
(expt .999999999999923 9113206036900124) gambit 1.1282707411028267e-305 correct 1.1282707411028267e-305 new 1.1282707411028267e-305
(expt .999999999999923 9113206036900125) gambit 1.1282707411028267e-305 correct 1.1282707411027398e-305 new 1.1282707411027398e-305
"/home/lucier/text/courses/computation/computational-reals/src/expt-test2.scm"

[bdeket]

Yes you are right,
I had an error in my calculation for the second exponent (your rest-arg)
With the correct version I don't find any error bigger than 2ulp which is more in line with what I expected :) Since a 1×product should introduce at worst 0.5 flulp and flexpt should be < 1 flulp
Corrected, the worst I've found is ( 1.1 flulp ):
(expt 0.9999999999999614 18337817992015593) ; = 2.0150635414821977e-308

[gambiteer]

Here's some lightly tested code that I'll probably integrate into Gambit (so it has some Gambitisms). Maybe it will be helpful.

(include "lib/header.scm")

(macro-case-target
 ((C)
  (define nextafter
    (c-lambda (double double) double "nextafter")))
 (else
  #t)
 )


(define-prim (expt-flonum-exact-int x y)
  ;; Operations on known flonums are flonum specific.
  ;; So there are no operations for the compiler to test
  ;; for flonum arguments and inline.
  (declare (mostly-fixnum))
  (cond
   ;; The one really special case
   ((eqv? y 0)
    1)
   ;; The fast path
   ((##exact-int->flonum-exact? y)
    (flexpt x (inexact y)))
   ;; The next two cases simplify arguments about
   ;; later comparisons.
   ((flnan? x)
    x)
   ((flzero? x)
    (if (and (odd? y) (flnegative? (##flcopysign 1. x)))
        -0.
        0.))
   ;; Unfortunately, we need to deal with the case
   ;; of large y separately, because (inexact y)
   ;; might overflow to an infinity.
   ((or (<= y (- (expt 2 63)))
        (<= (expt 2 63) y))
    (cond ((fl= (flabs x) 1.)
           (if (and (odd? y) (fl= x -1.0))
               -1.
               1.))
          ((fl< 1. (flabs x))
           (if (positive? y)
               (if (and (odd? y) (flnegative? x))
                   -inf.0
                   +inf.0)
               (if (and (odd? y) (flnegative? x))
                   -0.
                   +0.)))
          (else
           ;; (fl< (flabs x) 1.)
           (if (positive? y)
               (if (and (odd? y) (flnegative? x))
                   -0.
                   +0.)
               (if (and (odd? y) (flnegative? x))
                   -inf.0
                   +inf.0)))))
   (else
    (macro-case-target
     ((C)
      ;; we decompose y into the sum of two integers
      ;; with the same sign, each of which can be
      ;; converted exactly to a flonum.
      (let* ((big-y
              (let ((inexact-y (inexact y)))
                (if (<= (flabs inexact-y) (abs y))
                    inexact-y
                    (nextafter inexact-y 0.))))
             (rest-y
              (inexact (- y (exact big-y)))))
        ;; the following multiplication is not
        ;; an infinity times a zero.
        (fl* (flexpt x big-y)
             (flexpt x rest-y))))
     (else
      ;; punt
      (flexpt x (inexact y)))))))

[gambiteer]

I committed a change to Gambit that achieves roughly the same thing without using nextafter (and so works on all Gambit backends): gambit/gambit@9f8adf2

[gambiteer]

Here's the final code to fix Gambit. There's a test before calling this code to return 1 if the exponent is exact 0. Perhaps you can do something similar for Racket. (##exact-int->flonum-exact? could be replaced by a less precise test, e.g., something like (< (abs x) (expt 2 53)).)

(define (flonum-exact-int-expt x y)
    ;; x is a flonum, y is a nonzero exact-int
    (cond
     ;; The fast path
     ((##exact-int->flonum-exact? y)
      (flexpt x (inexact y)))
     ;; singular cases
     ((flnan? x)         x)
     ((flzero? x)        (if (odd? y) x 0.))
     ((fl= (flabs x) 1.) (if (odd? y) x 1.))
     ;; extremely large exponents
     ((or (<= y (- (expt 2 63)))
          (<= (expt 2 63) y))
      ;; Only extreme values (+/- 0, infty) are possible
      ;; in IEEE double precision.
      ;; for future reference
      ;; (nextafter 1. +inf.0) => 1.0000000000000002
      ;; (nextafter 1. +.0   ) =>  .9999999999999999
      (if (eq? (fl< 1. (flabs x)) (positive? y))
          ;; result is an infinity
          (if (and (odd? y) (flnegative? x)) -inf.0 +inf.0)
          ;; result is a zero
          (if (and (odd? y) (flnegative? x)) -0. +0.)))
     (else
      (let* ((abs-big-y
              ;; remove the lowest 12 bits of (abs y)
              (arithmetic-shift (arithmetic-shift (abs y) -12) 12))
             (big-part-of-y
              (if (negative? y)
                  (- abs-big-y)
                  abs-big-y))
             (rest-of-y
              (- y big-part-of-y)))
        ;; y = big-part-of-y + rest-of-y, both terms on right are nonzero
        ;; because (abs y) <= 2^63 and it can't be converted exactly to a flonum,
        ;; so there are more than 53 upper bits of y.
        ;; big-part-of-y and rest-of-y can be converted to flonum without roundoff error.
        ;; y, big-part-of-y, and rest-of-y have the same sign, so the following
        ;; product is not an infinity times a zero.
        (fl* (flexpt x (inexact big-part-of-y))
             (flexpt x (inexact rest-of-y)))))))