My initial attempt at these went wrong in various ways, oops.
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.
. 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
(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)
(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!)
(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))))
(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.)
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
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"))
;; 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
(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)))
(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))
(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