From 83bf9b9fa283b9fe68e5df0feb0f03381f9333af Mon Sep 17 00:00:00 2001 From: Taylor R Campbell Date: Tue, 20 Nov 2018 06:55:35 +0000 Subject: [PATCH] Redo some of the floating-point parameters. Test them. My initial attempt at these went wrong in various ways, oops. --- src/relnotes/flonum | 54 +++++++--- src/runtime/arith.scm | 141 +++++++++++++++++++------ src/runtime/runtime.pkg | 22 ++-- tests/microcode/test-flonum-except.scm | 12 +-- tests/runtime/test-arith.scm | 120 ++++++++++++++++++++- 5 files changed, 282 insertions(+), 67 deletions(-) diff --git a/src/relnotes/flonum b/src/relnotes/flonum index bedc759ac..e852e3d9c 100644 --- a/src/relnotes/flonum +++ b/src/relnotes/flonum @@ -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 diff --git a/src/runtime/arith.scm b/src/runtime/arith.scm index a0fd44853..d1883d301 100644 --- a/src/runtime/arith.scm +++ b/src/runtime/arith.scm @@ -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.) diff --git a/src/runtime/runtime.pkg b/src/runtime/runtime.pkg index 7ee5cadbd..16d53848e 100644 --- a/src/runtime/runtime.pkg +++ b/src/runtime/runtime.pkg @@ -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 diff --git a/tests/microcode/test-flonum-except.scm b/tests/microcode/test-flonum-except.scm index 3cbaa4271..8bf8502fe 100644 --- a/tests/microcode/test-flonum-except.scm +++ b/tests/microcode/test-flonum-except.scm @@ -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 diff --git a/tests/runtime/test-arith.scm b/tests/runtime/test-arith.scm index 31f87ee4c..24332b01c 100644 --- a/tests/runtime/test-arith.scm +++ b/tests/runtime/test-arith.scm @@ -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 -- 2.25.1