Redo some of the floating-point parameters. Test them.
authorTaylor R Campbell <campbell@mumble.net>
Tue, 20 Nov 2018 06:55:35 +0000 (06:55 +0000)
committerTaylor R Campbell <campbell@mumble.net>
Tue, 20 Nov 2018 06:57:09 +0000 (06:57 +0000)
My initial attempt at these went wrong in various ways, oops.

src/relnotes/flonum
src/runtime/arith.scm
src/runtime/runtime.pkg
tests/microcode/test-flonum-except.scm
tests/runtime/test-arith.scm

index bedc759ac9a684ac7291ad17f1067189e769aa7e..e852e3d9c627c9625d505d7fd811bb506c5a102f 100644 (file)
@@ -1,8 +1,10 @@
 New flonum-related definitions:
 
-. flo:radix - Floating-point radix.  This is essentially always 2 but
-  useful for assertions if you want to future-proof code that assumes
-  it with a noisy failure in case we ever violate the assumption.
+. flo:radix - Floating-point radix as an integer.  This is essentially
+  always 2 but useful for assertions if you want to future-proof code
+  that assumes it with a noisy failure in case we ever violate the
+  assumption.
+- flo:radix. - Floating-point radix as a flonum.
 . flo:error-bound - Greatest possible relative error in rounding to nearest.
 . flo:log-error-bound = (log flo:error-bound)
 . flo:ulp-of-one - Distance from 1 to next greater floating-point number.
@@ -10,21 +12,41 @@ New flonum-related definitions:
 . flo:log-ulp-of-one = (log flo:ulp-of-one)
 . (flo:ulp x) - Distance from x to the next floating-point number
   larger in magnitude with the same sign.
-. flo:normal-exponent-max-base-2 - Exponent of greatest power of two
-  that is a floating-point number.
-. flo:normal-exponent-min-base-2 - Exponent of least positive power of
-  two that is a _normal_ floating-point number.
-. flo:normal-exponent-max-base-e, flo:normal-exponent-min-base-e,
-  flo:normal-exponent-max-base-10, flo:normal-exponent-min-base-10 -
-  Exponents of the greatest and least positive powers of e and 10 that
-  are normal floating-point numbers.
-. flo:subnormal-exponent-min-base-2 - Exponent of least positive power
-  of two that is a nonzero, subnormal floating-point number.  This
-  number is the least positive floating-point number.
-. flo:subnormal-exponent-min-base-e, flo:subnormal-exponent-min-base-10 -
-  Exponents of least positive powers of e and 10.
+. flo:normal-exponent-max - Exponent of greatest integer power of
+  flo:radix that is a finite floating-point number, as an exact
+  integer.  Note that there are floating-point numbers greater than
+  (expt flo:radix. flo:normal-exponent-max) -- the greatest one is just
+  below flo:radix times that.  But there are none that are a greater
+  _integer_ power of flo:radix.
+. flo:normal-exponent-min - Exponent of least positive integer power of
+  flo:radix that is a normal floating-point number, as an exact
+  integer.  (expt flo:radix. flo:normal-exponent-min) is also named
+  flo:smallest-positive-normal.
+. flo:subnormal-exponent-min - Exponent of least positive integer power
+  of flo:radix that is a nonzero, subnormal floating-point number, as
+  an exact integer.  (expt flo:radix. flo:subnormal-exponent-min) is
+  the least positive floating-point number, also named
+  flo:smallest-positive-subnormal.
 . flo:smallest-positive-subnormal - Smallest positive subnormal.
 . flo:smallest-positive-normal - Smallest positive normal.
 . flo:largest-positive-normal - Largest positive normal.
+. flo:least-subnormal-exponent-base-2 - Least flonum input x for which
+  (expt 2. x) gives a nonzero result.
+. flo:least-subnormal-exponent-base-e - Least flonum input x for which
+  (exp x) gives a nonzero result.
+. flo:least-subnormal-exponent-base-10 - Least flonum input x for which
+  (expt 10. x) gives a nonzero result.
+. flo:least-normal-exponent-base-2 - Least flonum input x for which
+  (expt 2. x) gives a normal result.
+. flo:least-normal-exponent-base-e - Least flonum input x for which
+  (exp x) gives a normal result.
+. flo:least-normal-exponent-base-10 - Least flonum input x for which
+  (expt 10. x) gives a normal result.
+. flo:greatest-normal-exponent-base-2 - Greatest flonum input x for
+  which (expt 2. x) gives a finite result.
+. flo:greatest-normal-exponent-base-e - Greatest flonum input x for
+  which (exp x) gives a finite result.
+. flo:greatest-normal-exponent-base-10 - Greatest flonum input x for
+  which (expt 10. x) gives a finite result.
 . (flo:ldexp x e) = x * 2^e
 . (flo:scalbn x e) = x * b^e, where b is flo:radix
index a0fd4485393f8be78b024309b1ddc770cd8b46c2..d1883d3014bdfcc2fe4f64cf15afce654c55cd24 100644 (file)
@@ -93,22 +93,26 @@ USA.
 (define rec:pi/2 (flo:* 2. (flo:atan2 1. 1.)))
 (define rec:pi (flo:* 2. rec:pi/2))
 
-(define flo:radix 2)
+(define-integrable flo:radix 2)
+(define-integrable flo:radix. 2.)
 (define flo:ulp-of-one)
 (define flo:error-bound)
 (define flo:log-error-bound)
-(define flo:normal-exponent-max-base-2)
-(define flo:normal-exponent-min-base-2)
-(define flo:normal-exponent-max-base-e)
-(define flo:normal-exponent-min-base-e)
-(define flo:normal-exponent-max-base-10)
-(define flo:normal-exponent-min-base-10)
-(define flo:subnormal-exponent-min-base-2)
-(define flo:subnormal-exponent-min-base-e)
-(define flo:subnormal-exponent-min-base-10)
+(define flo:subnormal-exponent-min)
+(define flo:normal-exponent-min)
+(define flo:normal-exponent-max)
 (define flo:smallest-positive-subnormal)
 (define flo:smallest-positive-normal)
 (define flo:largest-positive-normal)
+(define flo:least-subnormal-exponent-base-2)
+(define flo:least-subnormal-exponent-base-e)
+(define flo:least-subnormal-exponent-base-10)
+(define flo:least-normal-exponent-base-2)
+(define flo:least-normal-exponent-base-e)
+(define flo:least-normal-exponent-base-10)
+(define flo:greatest-normal-exponent-base-2)
+(define flo:greatest-normal-exponent-base-e)
+(define flo:greatest-normal-exponent-base-10)
 (define flo:significand-digits-base-2)
 (define flo:significand-digits-base-10)
 (define int:flonum-integer-limit)
@@ -130,33 +134,100 @@ USA.
   (set! flo:ulp-of-one microcode-id/floating-epsilon)
   (set! flo:error-bound (flo:/ flo:ulp-of-one 2.))
   (set! flo:log-error-bound (flo:log flo:error-bound))
-  (set! flo:normal-exponent-max-base-2 microcode-id/floating-exponent-max)
-  (set! flo:normal-exponent-min-base-2 microcode-id/floating-exponent-min)
-  (set! flo:subnormal-exponent-min-base-2
-       (int:- flo:normal-exponent-min-base-2
+  (set! flo:normal-exponent-max microcode-id/floating-exponent-max)
+  (set! flo:normal-exponent-min microcode-id/floating-exponent-min)
+  (set! flo:subnormal-exponent-min
+       (int:- flo:normal-exponent-min
               (int:- flo:significand-digits-base-2 1)))
-  (set! flo:normal-exponent-max-base-e
-       (flo:log (flo:expt 2. (int:->flonum flo:normal-exponent-max-base-2))))
-  (set! flo:normal-exponent-min-base-e
-       (flo:log (flo:expt 2. (int:->flonum flo:normal-exponent-min-base-2))))
-  (set! flo:subnormal-exponent-min-base-e
-       (flo:+
-        flo:normal-exponent-min-base-e
-        (flo:log
-         (flo:expt 2.
-                   (flo:- 0. (int:->flonum flo:significand-digits-base-2))))))
-  (set! flo:normal-exponent-max-base-10
-       (flo:/ flo:normal-exponent-max-base-e (flo:log 10.)))
-  (set! flo:normal-exponent-min-base-10
-       (flo:/ flo:normal-exponent-min-base-e (flo:log 10.)))
-  (set! flo:subnormal-exponent-min-base-10
-       (flo:/ flo:subnormal-exponent-min-base-e (flo:log 10.)))
   (set! flo:smallest-positive-subnormal
-       (flo:ldexp 1. flo:subnormal-exponent-min-base-2))
+       (flo:scalbn 1. flo:subnormal-exponent-min))
   (set! flo:smallest-positive-normal
-       (flo:ldexp 1. flo:normal-exponent-min-base-2))
+       (flo:scalbn 1. flo:normal-exponent-min))
   (set! flo:largest-positive-normal
-       (flo:ldexp (flo:nextafter 2. 0.) flo:normal-exponent-max-base-2))
+       (flo:scalbn (flo:nextafter flo:radix. 0.) flo:normal-exponent-max))
+  (define (rejigger x worst good?)
+    ;; We could use bisection but the initial approximations should be
+    ;; good enough that this search should take at most a couple steps.
+    (define (good x i)
+      (if (fix:>= i 10)
+         (error "Floating-point parameters are hosed, can't find 'em!"))
+      (let ((x* (flo:nextafter x worst)))
+       (if (good? x*) (good x* (fix:+ i 1)) x)))
+    (define (bad x i)
+      (if (fix:>= i 10)
+         (error "Floating-point parameters are hosed, can't find 'em!"))
+      (let ((x* (flo:nextafter x (flo:negate worst))))
+       (if (good? x*) x* (bad x* (fix:+ i 1)))))
+    (if (good? x) (good x 0) (bad x 0)))
+  (set! flo:least-subnormal-exponent-base-2
+       (rejigger (flo:- (int:->flonum flo:subnormal-exponent-min) 1.)
+                 (flo:negate flo:largest-positive-normal)
+                 (lambda (x) (not (flo:zero? (flo:expt 2. x))))))
+  (set! flo:least-subnormal-exponent-base-e
+       (rejigger (flo:- (flo:log flo:smallest-positive-subnormal) flo:log2)
+                 (flo:negate flo:largest-positive-normal)
+                 (lambda (x) (not (flo:zero? (flo:exp x))))))
+  (set! flo:least-subnormal-exponent-base-10
+       (rejigger (flo:/ flo:least-subnormal-exponent-base-e (flo:log 10.))
+                 (flo:negate flo:largest-positive-normal)
+                 (lambda (x) (not (flo:zero? (flo:expt 10. x))))))
+  (set! flo:least-normal-exponent-base-2
+       (int:->flonum flo:normal-exponent-min))
+  (set! flo:least-normal-exponent-base-e
+       (rejigger (flo:log flo:smallest-positive-normal)
+                 (flo:negate flo:largest-positive-normal)
+                 (lambda (x) (flo:normal? (exp x)))))
+  (set! flo:least-normal-exponent-base-10
+       (rejigger (flo:/ flo:least-normal-exponent-base-e (flo:log 10.))
+                 (flo:negate flo:largest-positive-normal)
+                 (lambda (x) (flo:normal? (flo:expt 10. x)))))
+  (let ((b flo:radix.)
+       (emax (int:->flonum flo:normal-exponent-max)))
+    ;; Let eps = ulp(1)/2 = 1/(2 b^{p - 1}).  The largest positive
+    ;; floating-point number m is
+    ;;
+    ;; b^emax (b - 1 + (b^{p - 1} - 1)/b^{p - 1})
+    ;; b^emax (b - 1 + (2 b^{p - 1} - 2)/(2 b^{p - 1}))
+    ;; = b^emax (b - 1 + 1 - 2 eps)
+    ;; = b^emax (b - 2 eps).
+    ;;
+    ;; What _would_ be the next larger floating-point number is m' =
+    ;; b^{emax + 1}; anything up to (and excluding) halfway from m to
+    ;; m' gets rounded to m, while anything halfway to m' or above is
+    ;; rounded to +infinity.
+    ;;
+    ;; We seek the greatest floating-point x with fl(e^x) < (m + m')/2.
+    ;; We can start by computing an approximation to the log of
+    ;;
+    ;; (b*b^emax + b^emax (b - 2 eps))/2
+    ;; = b^emax (b + b - 2 eps)/2
+    ;; = b^emax (2b - 2 eps)/2
+    ;; = b^emax 2b (1 - 2 eps/b)/2
+    ;; = b^{emax + 1} (1 - 2 eps/b)
+    ;;
+    ;; by
+    ;;
+    ;; log(b^{emax + 1} (1 - 2 eps/b))
+    ;; = log(b^{emax + 1}) + log(1 - 2 eps/b)
+    ;; = (emax + 1) log(b) + log1p(-2 eps/b)
+    ;; = (emax + 1) log(b) + log1p(-ulp(1)/b)
+    ;;
+    ;; It will be off by at most a handful of rounding errors, so we
+    ;; rejigger this estimate with a short search for the correct one.
+    ;;
+    (set! flo:greatest-normal-exponent-base-e
+         (rejigger (flo:+ (flo:* (flo:+ emax 1.) (flo:log b))
+                          (flo:log1p (flo:/ (flo:negate flo:ulp-of-one) b)))
+                   flo:largest-positive-normal
+                   (lambda (x) (flo:finite? (flo:exp x)))))
+    (set! flo:greatest-normal-exponent-base-2
+         (rejigger (flo:/ flo:greatest-normal-exponent-base-e (flo:log 2.))
+                   flo:largest-positive-normal
+                   (lambda (x) (flo:finite? (flo:expt 2. x)))))
+    (set! flo:greatest-normal-exponent-base-10
+         (rejigger (flo:/ flo:greatest-normal-exponent-base-e (flo:log 10.))
+                   flo:largest-positive-normal
+                   (lambda (x) (flo:finite? (flo:expt 10. x))))))
   unspecific)
 
 (define (initialize-package!)
@@ -2083,7 +2154,7 @@ USA.
 
 (define (log1pexp x)
   (guarantee-real x 'log1pexp)
-  (cond ((<= x flo:subnormal-exponent-min-base-e) 0.)
+  (cond ((< x flo:least-subnormal-exponent-base-e) 0.)
        ((<= x flo:log-error-bound) (exp x))
        ((<= x 18) (log1p (exp x)))
        ((<= x 33.3) (+ x (exp (- x))))
@@ -2212,7 +2283,7 @@ USA.
 
 (define (logistic x)
   (guarantee-real x 'logistic)
-  (cond ((<= x flo:subnormal-exponent-min-base-e)
+  (cond ((< x flo:least-subnormal-exponent-base-2)
         ;; e^x/(1 + e^x) < e^x < smallest positive float.  (XXX Should
         ;; raise inexact and underflow here.)
         0.)
index 7ee5cadbd979a5df89a9f7013ff3d22efbdca662..16d53848e0d3b197040ad25df0ca05c6761dfdb0 100644 (file)
@@ -3359,22 +3359,26 @@ USA.
          exact-nonnegative-integer?
          exact-positive-integer?
          flo:error-bound
+         flo:greatest-normal-exponent-base-10
+         flo:greatest-normal-exponent-base-2
+         flo:greatest-normal-exponent-base-e
          flo:largest-positive-normal
+         flo:least-normal-exponent-base-10
+         flo:least-normal-exponent-base-2
+         flo:least-normal-exponent-base-e
+         flo:least-subnormal-exponent-base-10
+         flo:least-subnormal-exponent-base-2
+         flo:least-subnormal-exponent-base-e
          flo:log-error-bound
-         flo:normal-exponent-max-base-10
-         flo:normal-exponent-max-base-2
-         flo:normal-exponent-max-base-e
-         flo:normal-exponent-min-base-10
-         flo:normal-exponent-min-base-2
-         flo:normal-exponent-min-base-e
+         flo:normal-exponent-max
+         flo:normal-exponent-min
          flo:radix
+         flo:radix.
          flo:significand-digits-base-10
          flo:significand-digits-base-2
          flo:smallest-positive-normal
          flo:smallest-positive-subnormal
-         flo:subnormal-exponent-min-base-10
-         flo:subnormal-exponent-min-base-2
-         flo:subnormal-exponent-min-base-e
+         flo:subnormal-exponent-min
          flo:ulp-of-one
          flonum-printer:engineering-output
          flonum-printer:normal-output
index 3cbaa4271b8f9361c5357b422a331723ddb9fe3a..8bf8502feeb455969760c9b958ad7be4e8c90e6a 100644 (file)
@@ -32,7 +32,7 @@ USA.
   x)
 
 (define (flo:subnormal? x)
-  (flo:< (flo:abs x) (flo:ldexp 1. flo:normal-exponent-min-base-2)))
+  (flo:< (flo:abs x) flo:smallest-positive-normal))
 
 (define assert-flonum
   (predicate-assertion flo:flonum? "object"))
@@ -208,21 +208,21 @@ USA.
 
 ;; XXX Check rounding modes.
 (define-overflow-flag-test 'flonum-multiply
-  (applicator flo:* 2. (flo:ldexp 1. flo:normal-exponent-max-base-2)))
+  (applicator flo:* flo:radix. (flo:scalbn 1. flo:normal-exponent-max)))
 (define-overflow-trap-test 'flonum-multiply
-  (applicator flo:* 2. (flo:ldexp 1. flo:normal-exponent-max-base-2)))
+  (applicator flo:* flo:radix. (flo:scalbn 1. flo:normal-exponent-max)))
 
 ;;; IEEE 754-2008, Sec. 7.5
 
 (define-underflow-flag-test 'flonum-multiply
-  (applicator flo:* .50000001 (flo:ldexp 1. flo:normal-exponent-min-base-2)))
+  (applicator flo:* .50000001 (flo:scalbn 1. flo:normal-exponent-min)))
 (define-underflow-trap-test 'flonum-multiply
-  (applicator flo:* .50000001 (flo:ldexp 1. flo:normal-exponent-min-base-2)))
+  (applicator flo:* .50000001 (flo:scalbn 1. flo:normal-exponent-min)))
 
 ;;; IEEE 754-2008, Sec. 7.6
 
 (define-inexact-flag-test 'flonum-multiply
-  (applicator flo:* .50000001 (flo:ldexp 1. flo:normal-exponent-min-base-2)))
+  (applicator flo:* .50000001 (flo:scalbn 1. flo:normal-exponent-min)))
 
 ;;; Miscellaneous inexact results
 
index 31f87ee4c4c0371df24148c52c209611af21a7ba..24332b01c664c5aa92aa3b97575a26c8d5817ce2 100644 (file)
@@ -35,6 +35,15 @@ USA.
   (assert-true (flo:flonum? object))
   (assert-true (flo:nan? object)))
 
+(define (assert-zero object)
+  (assert-= object 0))
+
+(define (assert-zero+ object)
+  (assert-eqv object 0.))
+
+(define (assert-zero- object)
+  (assert-eqv object -0.))
+
 (define (assert-inf- object)
   (assert-eqv object (flo:-inf.0)))
 
@@ -47,12 +56,29 @@ USA.
 (define assert-integer
   (predicate-assertion integer? "integer"))
 
+(define assert-exact-integer
+  (predicate-assertion exact-integer? "integer"))
+
 (define assert-not-integer
   (predicate-assertion not-integer? "not integer"))
 
+(define assert-flonum
+  (predicate-assertion flo:flonum? "flonum"))
+
 (define assert-real
   (predicate-assertion real? "real number"))
 
+(define assert-normal
+  (predicate-assertion flo:normal? "normal floating-point number"))
+
+(define (flo:subnormal? x)
+  (and (flo:finite? x)
+       (not (flo:zero? x))
+       (not (flo:normal? x))))
+
+(define assert-subnormal
+  (predicate-assertion flo:subnormal? "subnormal floating-point number"))
+
 (define (with-expected-failure xfail? body)
   (case xfail?
     ((xfail) (expect-failure body))
@@ -621,4 +647,96 @@ USA.
           (xfail? (if (<= 3 (vector-length v)) (vector-ref v 2) #f)))
       (with-expected-failure xfail?
         (lambda ()
-          (assert-<= (relerr t (atan x)) 1e-15))))))
\ No newline at end of file
+          (assert-<= (relerr t (atan x)) 1e-15))))))
+
+(define-test 'radix
+  (lambda ()
+    (assert-exact-integer flo:radix)
+    (assert-flonum flo:radix.)
+    (assert-= flo:radix flo:radix.)))
+
+(define-test 'error-ulp
+  (lambda ()
+    (assert-flonum flo:ulp-of-one)
+    (assert-flonum flo:error-bound)
+    (assert-> flo:ulp-of-one 0)
+    (assert-< flo:ulp-of-one 1)
+    (assert-= (/ flo:ulp-of-one 2) flo:error-bound)))
+
+(define-test 'exponents
+  (lambda ()
+    (assert-exact-integer flo:normal-exponent-max)
+    (assert-exact-integer flo:normal-exponent-min)
+    (assert-exact-integer flo:subnormal-exponent-min)))
+
+(define-test 'extremes
+  (lambda ()
+    (assert-flonum flo:smallest-positive-subnormal)
+    (assert-flonum flo:smallest-positive-normal)
+    (assert-flonum flo:largest-positive-normal)
+    (assert-eqv flo:smallest-positive-subnormal
+                (flo:scalbn 1. flo:subnormal-exponent-min))
+    (assert-eqv flo:smallest-positive-normal
+                (flo:scalbn 1. flo:normal-exponent-min))
+    (assert-eqv flo:largest-positive-normal
+                (flo:scalbn (- flo:radix. flo:ulp-of-one)
+                            flo:normal-exponent-max))
+    (assert-subnormal flo:smallest-positive-subnormal)
+    (assert-zero+ (flo:nextafter flo:smallest-positive-subnormal (flo:-inf.0)))
+    (assert-normal flo:smallest-positive-normal)
+    (assert-subnormal
+     (flo:nextafter flo:smallest-positive-normal (flo:-inf.0)))
+    (assert-normal flo:largest-positive-normal)
+    (assert-inf+ (flo:nextafter flo:largest-positive-normal (flo:+inf.0)))))
+
+(define-test 'least-subnormal-exponents
+  (lambda ()
+    (assert-flonum flo:least-subnormal-exponent-base-2)
+    (assert-flonum flo:least-subnormal-exponent-base-e)
+    (assert-flonum flo:least-subnormal-exponent-base-10)
+    (assert-subnormal (expt 2. flo:least-subnormal-exponent-base-2))
+    (assert-subnormal (exp flo:least-subnormal-exponent-base-e))
+    (assert-subnormal (expt 10. flo:least-subnormal-exponent-base-10))
+    (assert-zero+
+     (expt 2.
+           (flo:nextafter flo:least-subnormal-exponent-base-2 (flo:-inf.0))))
+    (assert-zero+
+     (exp (flo:nextafter flo:least-subnormal-exponent-base-e (flo:-inf.0))))
+    (assert-zero+
+     (expt 10.
+           (flo:nextafter flo:least-subnormal-exponent-base-10
+                          (flo:-inf.0))))))
+
+(define-test 'least-normal-exponents
+  (lambda ()
+    (assert-flonum flo:least-normal-exponent-base-2)
+    (assert-flonum flo:least-normal-exponent-base-e)
+    (assert-flonum flo:least-normal-exponent-base-10)
+    (assert-normal (expt 2. flo:least-normal-exponent-base-2))
+    (assert-normal (exp flo:least-normal-exponent-base-e))
+    (assert-normal (expt 10. flo:least-normal-exponent-base-10))
+    (assert-subnormal
+     (expt 2.
+           (flo:nextafter flo:least-normal-exponent-base-2 (flo:-inf.0))))
+    (assert-subnormal
+     (exp (flo:nextafter flo:least-normal-exponent-base-e (flo:-inf.0))))
+    (assert-subnormal
+     (expt 10.
+           (flo:nextafter flo:least-normal-exponent-base-10 (flo:-inf.0))))))
+
+(define-test 'greatest-normal-exponents
+  (lambda ()
+    (assert-flonum flo:greatest-normal-exponent-base-2)
+    (assert-flonum flo:greatest-normal-exponent-base-e)
+    (assert-flonum flo:greatest-normal-exponent-base-10)
+    (assert-normal (expt 2. flo:greatest-normal-exponent-base-2))
+    (assert-normal (exp flo:greatest-normal-exponent-base-e))
+    (assert-normal (expt 10. flo:greatest-normal-exponent-base-10))
+    (assert-inf+
+     (expt 2.
+           (flo:nextafter flo:greatest-normal-exponent-base-2 (flo:+inf.0))))
+    (assert-inf+
+     (exp (flo:nextafter flo:greatest-normal-exponent-base-e (flo:+inf.0))))
+    (assert-inf+
+     (expt 10.
+           (flo:nextafter flo:greatest-normal-exponent-base-2 (flo:+inf.0))))))
\ No newline at end of file