Precise flonum tests (error usually must be <= 0.5 ulp), and prerequisite
additions/changes * More accurate `flulp-error' * Added `flonum->fields', `fields->flonum', `flonum->sig+exp', `sig+exp->flonum' (currently undocumented) * Correctly rounded, robust `bigfloat->fl2' and `fl2' * Correctly rounded, robust `fl+/error', `fl-/error', `fl*/error', `flsqr/error', `fl//error' * Much faster but slightly less accurate fl2 ops (shamelessly stolen from crlibm, which is LGPL) * Added `fl2ulp', `fl2ulp-error', `fl2?' (which detects overlap), `+max-fl2-subnormal.0' (which was tricky), `fl2abs' * Added deterministic and randomized flonum op tests (against MPFR) * Added deterministic and randomized flonum/error op tests (against MPFR) * Added deterministic and randomized fl2 op tests (against MPFR) * Exposed FPU tests in `math/utils' (currently undocumented)
This commit is contained in:
Neil Toronto
parent
fc02d40a66
commit
e55a31480e
|
|
@ -11,7 +11,7 @@
|
|||
"private/flonum/flonum-factorial.rkt"
|
||||
"private/flonum/flonum-log1pmx.rkt"
|
||||
"private/flonum/flonum-polyfun.rkt"
|
||||
"private/flonum/flonum-syntax.rkt"
|
||||
"private/flonum/flonum-error.rkt"
|
||||
"private/flonum/expansion/expansion-base.rkt"
|
||||
"private/flonum/expansion/expansion-exp.rkt"
|
||||
"private/flonum/expansion/expansion-log.rkt"
|
||||
|
|
@ -28,7 +28,7 @@
|
|||
"private/flonum/flonum-factorial.rkt"
|
||||
"private/flonum/flonum-log1pmx.rkt"
|
||||
"private/flonum/flonum-polyfun.rkt"
|
||||
"private/flonum/flonum-syntax.rkt"
|
||||
"private/flonum/flonum-error.rkt"
|
||||
"private/flonum/expansion/expansion-base.rkt"
|
||||
"private/flonum/expansion/expansion-exp.rkt"
|
||||
"private/flonum/expansion/expansion-log.rkt"
|
||||
|
|
|
|||
|
|
@ -1,6 +1,7 @@
|
|||
#lang typed/racket/base
|
||||
|
||||
(require (only-in "mpfr.rkt" 1ary-funs 1ary-preds 1ary2-funs 2ary-funs)
|
||||
(require racket/flonum
|
||||
(only-in "mpfr.rkt" 1ary-funs 1ary-preds 1ary2-funs 2ary-funs)
|
||||
"../base/base-random.rkt"
|
||||
"utils.rkt")
|
||||
|
||||
|
|
@ -139,11 +140,18 @@
|
|||
(: bigfloat->fl2 (Bigfloat -> (Values Flonum Flonum)))
|
||||
(define (bigfloat->fl2 x)
|
||||
(define x2 (bigfloat->flonum x))
|
||||
(values x2 (bigfloat->flonum (bf- x (flonum->bigfloat x2)))))
|
||||
(cond [(rational? x2)
|
||||
(let ([x2 (+ x2 (bigfloat->flonum (bf- x (flonum->bigfloat x2))))])
|
||||
(cond [(rational? x2)
|
||||
(values x2 (bigfloat->flonum (bf- x (flonum->bigfloat x2))))]
|
||||
[else
|
||||
(values x2 0.0)]))]
|
||||
[else (values x2 0.0)]))
|
||||
|
||||
(: fl2->bigfloat (Flonum Flonum -> Bigfloat))
|
||||
(define (fl2->bigfloat x2 x1)
|
||||
(bf+ (flonum->bigfloat x1) (flonum->bigfloat x2)))
|
||||
(cond [(fl= x1 0.0) (bf x2)]
|
||||
[else (bf+ (flonum->bigfloat x1) (flonum->bigfloat x2))]))
|
||||
|
||||
(provide
|
||||
;; Parameters
|
||||
|
|
|
|||
|
|
@ -6,81 +6,124 @@ Arithmetic based on:
|
|||
Jonathan Richard Shewchuk
|
||||
Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates
|
||||
Discrete & Computational Geometry 18(3):305–363, October 1997
|
||||
|
||||
Other parts shamelessly stolen from crlibm (which is LGPL)
|
||||
|#
|
||||
|
||||
(require "../flonum-functions.rkt"
|
||||
"../flonum-syntax.rkt")
|
||||
"../flonum-bits.rkt"
|
||||
"../flonum-error.rkt"
|
||||
"../flonum-constants.rkt"
|
||||
"../utils.rkt")
|
||||
|
||||
(provide fl2 fl2->real
|
||||
fl2+ fl2- fl2*split-fl fl2* fl2/
|
||||
flsqrt/error fl2sqrt)
|
||||
(provide fl2? fl2 fl2->real
|
||||
fl2ulp fl2ulp-error
|
||||
+max-fl2-subnormal.0 -max-fl2-subnormal.0
|
||||
fl2abs fl2+ fl2- fl2*split-fl fl2* fl2sqr fl2/
|
||||
fl2sqrt flsqrt/error)
|
||||
|
||||
(: floverlapping? (Flonum Flonum -> Boolean))
|
||||
(define (floverlapping? x2 x1)
|
||||
(define-values (s2 e2) (flonum->sig+exp (flabs x2)))
|
||||
(define-values (s1 e1) (flonum->sig+exp (flabs x1)))
|
||||
(define-values (n1 n2)
|
||||
(if (e2 . > . e1)
|
||||
(values s1 (arithmetic-shift s2 (- e2 e1)))
|
||||
(values (arithmetic-shift s1 (- e1 e2)) s2)))
|
||||
(not (= (bitwise-ior n1 n2)
|
||||
(bitwise-xor n1 n2))))
|
||||
|
||||
(: fl2? (Flonum Flonum -> Boolean))
|
||||
(define (fl2? x2 x1)
|
||||
(cond [(flrational? x2)
|
||||
(cond [(flrational? x1)
|
||||
(cond [((flabs x2) . < . (flabs x1)) #f]
|
||||
[else (not (floverlapping? x2 x1))])]
|
||||
[else #f])]
|
||||
[else (fl= x1 0.0)]))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; Conversion
|
||||
|
||||
(: fl2 (case-> (Real -> (Values Flonum Flonum))
|
||||
(Flonum Flonum -> (Values Flonum Flonum))))
|
||||
(define fl2
|
||||
(case-lambda
|
||||
[(x)
|
||||
(cond [(flonum? x) (values x 0.0)]
|
||||
[(single-flonum? x) (values (fl x) 0.0)]
|
||||
[else (let ([x2 (fl x)])
|
||||
(values x2 (fl (- x (inexact->exact x2)))))])]
|
||||
[(x y)
|
||||
(fast-fl+/error x y)]))
|
||||
(: fl2 (Real -> (Values Flonum Flonum)))
|
||||
(define (fl2 x)
|
||||
(cond [(flonum? x) (values x 0.0)]
|
||||
[(single-flonum? x) (values (fl x) 0.0)]
|
||||
[else (let ([x2 (fl x)])
|
||||
(cond [(flrational? x2)
|
||||
(let ([x2 (fl+ x2 (fl (- x (inexact->exact x2))))])
|
||||
(cond [(flrational? x2)
|
||||
(values x2 (fl (- x (inexact->exact x2))))]
|
||||
[else (values x2 0.0)]))]
|
||||
[else (values x2 0.0)]))]))
|
||||
|
||||
(: fl2->real (Flonum Flonum -> Real))
|
||||
(define (fl2->real x2 x1)
|
||||
(cond [(and (x1 . fl> . -inf.0) (x1 . fl< . +inf.0)
|
||||
(x2 . fl> . -inf.0) (x2 . fl< . +inf.0))
|
||||
(cond [(and (flrational? x2) (flrational? x1))
|
||||
(+ (inexact->exact x2) (inexact->exact x1))]
|
||||
[else (fl+ x1 x2)]))
|
||||
|
||||
(: fl3->fl2 (Flonum Flonum Flonum -> (Values Flonum Flonum)))
|
||||
(define (fl3->fl2 e3 e2 e1)
|
||||
(values e3 (fl+ e2 e1)))
|
||||
|
||||
(: fl4->fl2 (Flonum Flonum Flonum Flonum -> (Values Flonum Flonum)))
|
||||
(define (fl4->fl2 e4 e3 e2 e1)
|
||||
(values e4 (fl+ e3 (fl+ e2 e1))))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; Error
|
||||
|
||||
(: fl2ulp (Flonum Flonum -> Flonum))
|
||||
(define (fl2ulp x2 x1)
|
||||
(cond [(fl= x2 0.0) 0.0]
|
||||
[else (flmax +min.0 (fl* (flulp x2) (flexpt 2.0 -52.0)))]))
|
||||
|
||||
(: fl2ulp-error (Flonum Flonum Real -> Flonum))
|
||||
(define (fl2ulp-error x2 x1 r)
|
||||
(define x (fl2->real x2 x1))
|
||||
(define-values (r2 r1) (fl2 r))
|
||||
(cond [(eqv? x r) 0.0]
|
||||
[(and (fl= x2 0.0) (fl= r2 0.0)) 0.0]
|
||||
[(and (fl= x2 +inf.0) (fl= r2 +inf.0)) 0.0]
|
||||
[(and (fl= x2 -inf.0) (fl= r2 -inf.0)) 0.0]
|
||||
[(zero? r) +inf.0]
|
||||
[(and (rational? x) (flrational? r2))
|
||||
(flabs (fl (/ (- (inexact->exact x) (inexact->exact r))
|
||||
(inexact->exact (flmax +min.0 (fl2ulp r2 r1))))))]
|
||||
[else +inf.0]))
|
||||
|
||||
(define +max-fl2-subnormal.0 1.0020841800044863e-292) ; (flprev (flexpt 2.0 -970.0))
|
||||
(define -max-fl2-subnormal.0 (- +max-fl2-subnormal.0))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; Absolute value
|
||||
|
||||
(: fl2abs (case-> (Flonum -> (Values Flonum Flonum))
|
||||
(Flonum Flonum -> (Values Flonum Flonum))))
|
||||
(define (fl2abs x2 [x1 0.0])
|
||||
(define x (fl+ x2 x1))
|
||||
(cond [(not (flrational? x)) (values (flabs x2) 0.0)]
|
||||
[(fl= x 0.0) (values 0.0 0.0)]
|
||||
[(fl> x 0.0) (values x2 x1)]
|
||||
[else (values (- x2) (- x1))]))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; Addition and subtraction
|
||||
|
||||
(: raw-fl2+fl (Flonum Flonum Flonum -> (Values Flonum Flonum Flonum)))
|
||||
(define (raw-fl2+fl e2 e1 b)
|
||||
(let*-values ([(Q h1) (fast-fl+/error b e1)]
|
||||
[(h3 h2) (fast-fl+/error Q e2)])
|
||||
(values h3 h2 h1)))
|
||||
|
||||
(: raw-fl2+ (Flonum Flonum Flonum Flonum -> (Values Flonum Flonum Flonum Flonum)))
|
||||
(define (raw-fl2+ e2 e1 f2 f1)
|
||||
(let*-values ([(h3 h2 h1) (raw-fl2+fl e2 e1 f1)]
|
||||
[(h4 h3 h2) (raw-fl2+fl h3 h2 f2)])
|
||||
(values h4 h3 h2 h1)))
|
||||
|
||||
(: fl2+ (case-> (Flonum Flonum Flonum -> (Values Flonum Flonum))
|
||||
(Flonum Flonum Flonum Flonum -> (Values Flonum Flonum))))
|
||||
(define fl2+
|
||||
(case-lambda
|
||||
[(e2 e1 b)
|
||||
(let-values ([(h3 h2 h1) (raw-fl2+fl e2 e1 b)])
|
||||
(fl3->fl2 h3 h2 h1))]
|
||||
[(x2 x1 y2 y1)
|
||||
(let*-values ([(e4 e3 e2 e1) (raw-fl2+ x2 x1 y2 y1)])
|
||||
(fl4->fl2 e4 e3 e2 e1))]))
|
||||
(define (fl2+ x2 x1 y2 [y1 0.0])
|
||||
(define r (fl+ x2 y2))
|
||||
(cond [(not (flrational? r)) (values r 0.0)]
|
||||
[(and (fl= x2 0.0) (fl= y2 0.0)) (values r 0.0)]
|
||||
[else
|
||||
(define s (if ((flabs x2) . fl> . (flabs y2))
|
||||
(fl+ (fl+ (fl+ (fl- x2 r) y2) y1) x1)
|
||||
(fl+ (fl+ (fl+ (fl- y2 r) x2) x1) y1)))
|
||||
(define z2 (fl+ r s))
|
||||
(values z2 (fl+ (fl- r z2) s))]))
|
||||
|
||||
(: fl2- (case-> (Flonum Flonum Flonum -> (Values Flonum Flonum))
|
||||
(Flonum Flonum Flonum Flonum -> (Values Flonum Flonum))))
|
||||
(define fl2-
|
||||
(case-lambda
|
||||
[(e2 e1 b)
|
||||
(let-values ([(h3 h2 h1) (raw-fl2+fl e2 e1 (- b))])
|
||||
(fl3->fl2 h3 h2 h1))]
|
||||
[(x2 x1 y2 y1)
|
||||
(let*-values ([(e4 e3 e2 e1) (raw-fl2+ x2 x1 (- y2) (- y1))])
|
||||
(fl4->fl2 e4 e3 e2 e1))]))
|
||||
(define (fl2- x2 x1 y2 [y1 0.0])
|
||||
(fl2+ x2 x1 (- y2) (- y1)))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; Multiplication and division
|
||||
|
|
@ -119,56 +162,104 @@ Discrete & Computational Geometry 18(3):305–363, October 1997
|
|||
|
||||
(: fl2* (case-> (Flonum Flonum Flonum -> (Values Flonum Flonum))
|
||||
(Flonum Flonum Flonum Flonum -> (Values Flonum Flonum))))
|
||||
(define fl2*
|
||||
(define (fl2* x2 x1 y2 [y1 0.0])
|
||||
(define z (fl* x2 y2))
|
||||
(cond [(and (flrational? z) (not (fl= z 0.0)) (not (flsubnormal? z)))
|
||||
(define dx (near-pow2 x2))
|
||||
(define dy (near-pow2 y2))
|
||||
(define d (fl* dx dy))
|
||||
(let* ([x2 (fl/ x2 dx)]
|
||||
[x1 (fl/ x1 dx)]
|
||||
[y2 (fl/ y2 dy)]
|
||||
[y1 (fl/ y1 dy)]
|
||||
[up (fl* x2 (fl+ 1.0 (flexpt 2.0 27.0)))]
|
||||
[vp (fl* y2 (fl+ 1.0 (flexpt 2.0 27.0)))]
|
||||
[u1 (fl+ (fl- x2 up) up)]
|
||||
[v1 (fl+ (fl- y2 vp) vp)]
|
||||
[u2 (fl- x2 u1)]
|
||||
[v2 (fl- y2 v1)]
|
||||
[m2 (fl* x2 y2)]
|
||||
[m1 (fl+ (fl+ (fl+ (fl+ (fl+ (fl- (fl* u1 v1) m2)
|
||||
(fl* u1 v2))
|
||||
(fl* u2 v1))
|
||||
(fl* u2 v2))
|
||||
(fl* x2 y1))
|
||||
(fl* x1 y2))]
|
||||
[z2 (fl+ m2 m1)]
|
||||
[z1 (fl+ (fl- m2 z2) m1)]
|
||||
[z2 (fl* d z2)])
|
||||
(values z2 (if (flrational? z2) (fl* d z1) 0.0)))]
|
||||
[else
|
||||
(values z 0.0)]))
|
||||
|
||||
(: fl2sqr (case-> (Flonum -> (Values Flonum Flonum))
|
||||
(Flonum Flonum -> (Values Flonum Flonum))))
|
||||
;; Derived from fl2*
|
||||
(define fl2sqr
|
||||
(case-lambda
|
||||
[(e2 e1 b)
|
||||
(let-values ([(b-hi b-lo) (flsplit b)])
|
||||
(fl2*split-fl e2 e1 b-hi b-lo))]
|
||||
[(x2 x1 y2 y1)
|
||||
(let*-values ([(x2-hi x2-lo) (flsplit x2)]
|
||||
[(x1-hi x1-lo) (flsplit x1)]
|
||||
[(y2-hi y2-lo) (flsplit y2)]
|
||||
[(y1-hi y1-lo) (flsplit y1)]
|
||||
[(a2 a1) (split-fl2*split-fl x2-hi x2-lo x1-hi x1-lo y1-hi y1-lo)]
|
||||
[(b2 b1) (split-fl2*split-fl x2-hi x2-lo x1-hi x1-lo y2-hi y2-lo)])
|
||||
(fl2+ a2 a1 b2 b1))]))
|
||||
[(x) (flsqr/error x)]
|
||||
[(x2 x1)
|
||||
(define z (fl* x2 x2))
|
||||
(cond [(and (flrational? z) (not (fl= z 0.0)) (not (flsubnormal? z)))
|
||||
(define d (near-pow2 x2))
|
||||
(define d^2 (fl* d d))
|
||||
(let* ([x2 (fl/ x2 d)]
|
||||
[x1 (fl/ x1 d)]
|
||||
[up (fl* x2 (fl+ 1.0 (flexpt 2.0 27.0)))]
|
||||
[u1 (fl+ (fl- x2 up) up)]
|
||||
[u2 (fl- x2 u1)]
|
||||
[m2 (fl* x2 x2)]
|
||||
[m1 (fl+ (fl+ (fl+ (fl- (fl* u1 u1) m2)
|
||||
(fl* 2.0 (fl* u1 u2)))
|
||||
(fl* u2 u2))
|
||||
(fl* 2.0 (fl* x2 x1)))]
|
||||
[z2 (fl+ m2 m1)]
|
||||
[z1 (fl+ (fl- m2 z2) m1)]
|
||||
[z2 (fl* d^2 z2)])
|
||||
(values z2 (if (flrational? z2) (fl* d^2 z1) 0.0)))]
|
||||
[else
|
||||
(values z 0.0)])]))
|
||||
|
||||
(: fl2/ (case-> (Flonum Flonum Flonum -> (Values Flonum Flonum))
|
||||
(Flonum Flonum Flonum Flonum -> (Values Flonum Flonum))))
|
||||
(define fl2/
|
||||
(case-lambda
|
||||
[(x2 x1 y)
|
||||
(let*-values ([(a2 a1) (fast-fl//error x1 y)]
|
||||
[(b2 b1) (fast-fl//error x2 y)])
|
||||
(fl2+ a2 a1 b2 b1))]
|
||||
[(x2 x1 y2 y1)
|
||||
;; Compute three "digits" (flonums) of two-flonum long division; the third ensures the result is
|
||||
;; correctly rounded
|
||||
(let*-values ([(z2) (fl/ x2 y2)]
|
||||
[(w2 w1) (fl2* y2 y1 z2)]
|
||||
[(x2 x1) (fl2- x2 x1 w2 w1)]
|
||||
[(z1) (fl/ x2 y2)]
|
||||
[(w2 w1) (fl2* y2 y1 z1)]
|
||||
[(x2 x1) (fl2- x2 x1 w2 w1)])
|
||||
(fl3->fl2 z2 z1 (/ x2 y2)))]))
|
||||
(define (fl2/ x2 x1 y2 [y1 0.0])
|
||||
(define z (fl/ x2 y2))
|
||||
(cond [(and (flrational? z) (not (fl= z 0.0)) (flrational? y2))
|
||||
(define d (near-pow2/div x2 y2))
|
||||
(let*-values ([(x2 x1) (values (fl/ x2 d) (fl/ x1 d))]
|
||||
[(y2 y1) (values (fl/ y2 d) (fl/ y1 d))]
|
||||
[(c2) (fl/ x2 y2)]
|
||||
[(u2 u1) (fl*/error c2 y2)]
|
||||
[(c1) (fl/ (fl- (fl+ (fl- (fl- x2 u2) u1) x1) (fl* c2 y1)) y2)]
|
||||
[(z2) (fl+ c2 c1)])
|
||||
(values z2 (if (flrational? z2) (fl+ (fl- c2 z2) c1) 0.0)))]
|
||||
[else
|
||||
(values z 0.0)]))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; Square roots
|
||||
|
||||
(: flsqrt/error (Flonum -> (Values Flonum Flonum)))
|
||||
(: fl2sqrt (case-> (Flonum -> (Values Flonum Flonum))
|
||||
(Flonum Flonum -> (Values Flonum Flonum))))
|
||||
;; One-flonum estimate followed by one Newton's method iteration
|
||||
;; This could be a little faster if `y' were split only once
|
||||
(define (flsqrt/error x)
|
||||
(let*-values ([(y) (flsqrt x)]
|
||||
[(z2 z1) (fast-flsqr/error y)]
|
||||
[(dy2 dy1) (fl2+ (- z2) (- z1) x)]
|
||||
[(dy2 dy1) (fl2/ dy2 dy1 y)])
|
||||
(fl2+ (* 0.5 dy2) (* 0.5 dy1) y)))
|
||||
(define (fl2sqrt x2 [x1 0.0])
|
||||
(cond [(and (flrational? x2) (not (fl= x2 0.0)))
|
||||
(define-values (d^2 d)
|
||||
(cond [(x2 . fl<= . +max-fl2-subnormal.0) (values (flexpt 2.0 -104.0)
|
||||
(flexpt 2.0 -52.0))]
|
||||
[(x2 . fl> . 1e300) (values (flexpt 2.0 104.0)
|
||||
(flexpt 2.0 52.0))]
|
||||
[else (values 1.0 1.0)]))
|
||||
(let*-values ([(x2 x1) (values (fl/ x2 d^2) (fl/ x1 d^2))]
|
||||
[(y) (flsqrt (fl+ x2 x1))]
|
||||
[(z2 z1) (fast-flsqr/error y)]
|
||||
[(dy2 dy1) (fl2- x2 x1 z2 z1)]
|
||||
[(dy2 dy1) (fl2/ dy2 dy1 y)]
|
||||
[(y2 y1) (fl2+ (fl* 0.5 dy2) (fl* 0.5 dy1) y)]
|
||||
[(y2) (fl* y2 d)])
|
||||
(values y2 (if (flrational? y2) (fl* y1 d) 0.0)))]
|
||||
[else
|
||||
(values (flsqrt x2) 0.0)]))
|
||||
|
||||
(: fl2sqrt (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
(define (fl2sqrt x2 x1)
|
||||
(let*-values ([(y) (flsqrt (fl+ x1 x2))]
|
||||
[(z2 z1) (fast-flsqr/error y)]
|
||||
[(dy2 dy1) (fl2- x2 x1 z2 z1)]
|
||||
[(dy2 dy1) (fl2/ dy2 dy1 y)])
|
||||
(fl2+ (* 0.5 dy2) (* 0.5 dy1) y)))
|
||||
(: flsqrt/error (Flonum -> (Values Flonum Flonum)))
|
||||
(define (flsqrt/error x) (fl2sqrt x 0.0))
|
||||
|
|
|
|||
|
|
@ -17,7 +17,7 @@ Note to self: exp(x+y)-1 = (exp(x)-1) * (exp(y)-1) + (exp(x)-1) + (exp(y)-1)
|
|||
"../../../base.rkt"
|
||||
"../../../bigfloat.rkt"
|
||||
"../flonum-functions.rkt"
|
||||
"../flonum-syntax.rkt"
|
||||
"../flonum-error.rkt"
|
||||
"../flonum-exp.rkt"
|
||||
"../flonum-constants.rkt"
|
||||
"expansion-base.rkt"
|
||||
|
|
|
|||
|
|
@ -8,7 +8,7 @@ perform one Newton step.
|
|||
|#
|
||||
|
||||
(require "../flonum-functions.rkt"
|
||||
"../flonum-syntax.rkt"
|
||||
"../flonum-error.rkt"
|
||||
"../flonum-log.rkt"
|
||||
"expansion-base.rkt"
|
||||
"expansion-exp.rkt")
|
||||
|
|
|
|||
|
|
@ -4,6 +4,8 @@
|
|||
racket/performance-hint)
|
||||
|
||||
(provide flonum->bit-field bit-field->flonum
|
||||
flonum->fields fields->flonum
|
||||
flonum->sig+exp sig+exp->flonum
|
||||
flonum->ordinal ordinal->flonum
|
||||
flstep flnext flprev flonums-between
|
||||
flulp)
|
||||
|
|
@ -20,6 +22,48 @@
|
|||
[else
|
||||
(raise-argument-error 'bit-field->flonum "Integer in [0 .. #xffffffffffffffff]" i)]))
|
||||
|
||||
(define implicit-leading-one (arithmetic-shift 1 52))
|
||||
(define max-significand (- implicit-leading-one 1))
|
||||
(define max-exponent 2047)
|
||||
(define max-signed-exponent 1023)
|
||||
(define min-signed-exponent -1074)
|
||||
|
||||
(: flonum->fields (Flonum -> (Values (U 0 1) Index Natural)))
|
||||
(define (flonum->fields x)
|
||||
(define n (flonum->bit-field x))
|
||||
(values (if (zero? (bitwise-bit-field n 63 64)) 0 1)
|
||||
(assert (bitwise-bit-field n 52 63) index?)
|
||||
(bitwise-bit-field n 0 52)))
|
||||
|
||||
(: fields->flonum (Integer Integer Integer -> Flonum))
|
||||
(define (fields->flonum s e m)
|
||||
(cond [(not (or (= s 0) (= s 1)))
|
||||
(raise-argument-error 'fields->flonum "(U 0 1)" 0 s e m)]
|
||||
[(or (e . < . 0) (e . > . max-exponent))
|
||||
(raise-argument-error 'fields->flonum (format "Natural <= ~e" max-exponent) 1 s e m)]
|
||||
[(or (m . < . 0) (m . > . max-significand))
|
||||
(raise-argument-error 'fields->flonum (format "Natural <= ~a" max-significand) 2 s e m)]
|
||||
[else
|
||||
(bit-field->flonum (bitwise-ior (arithmetic-shift s 63)
|
||||
(arithmetic-shift e 52)
|
||||
m))]))
|
||||
|
||||
(: flonum->sig+exp (Flonum -> (Values Integer Fixnum)))
|
||||
(define (flonum->sig+exp x)
|
||||
(define-values (s e m) (flonum->fields x))
|
||||
(let-values ([(sig exp) (if (= e 0)
|
||||
(values m -1074)
|
||||
(values (bitwise-ior m implicit-leading-one)
|
||||
(assert (- e 1075) fixnum?)))])
|
||||
(values (if (zero? s) sig (- sig)) exp)))
|
||||
|
||||
(: sig+exp->flonum (Integer Integer -> Flonum))
|
||||
(define (sig+exp->flonum sig exp)
|
||||
(cond [(= sig 0) 0.0]
|
||||
[(exp . > . max-signed-exponent) (if (sig . < . 0) -inf.0 +inf.0)]
|
||||
[(exp . < . min-signed-exponent) (if (sig . < . 0) -0.0 0.0)]
|
||||
[else (real->double-flonum (* sig (expt 2 exp)))]))
|
||||
|
||||
(: flonum->ordinal (Flonum -> Integer))
|
||||
(define (flonum->ordinal x)
|
||||
(cond [(x . fl< . 0.0) (- (flonum->bit-field (fl- 0.0 x)))]
|
||||
|
|
|
|||
184
collects/math/private/flonum/flonum-error.rkt
Normal file
184
collects/math/private/flonum/flonum-error.rkt
Normal file
|
|
@ -0,0 +1,184 @@
|
|||
#lang racket/base
|
||||
|
||||
(provide flsplit
|
||||
fast-mono-fl+/error
|
||||
fast-mono-fl-/error
|
||||
fast-fl+/error
|
||||
fast-fl-/error
|
||||
fast-fl*/error
|
||||
fast-flsqr/error
|
||||
fast-fl//error
|
||||
fl+/error
|
||||
fl-/error
|
||||
fl*/error
|
||||
flsqr/error
|
||||
fl//error)
|
||||
|
||||
(module untyped-defs racket/base
|
||||
(require (for-syntax racket/base)
|
||||
"flonum-functions.rkt")
|
||||
|
||||
(provide (all-defined-out))
|
||||
|
||||
;(: flsplit (Flonum -> (Values Flonum Flonum)))
|
||||
;; Splits a flonum into a two flonums `hi' and `lo' with 26 bits precision each, such that
|
||||
;; |hi| >= |lo| and hi + lo = a. (The extra sign bit accounts for the missing bit.)
|
||||
;; This function returns (values +nan.0 +nan.0) for |a| >= 1.3393857490036326e+300.
|
||||
(define-syntax-rule (flsplit a-expr)
|
||||
(let ([a a-expr])
|
||||
(let* ([c (fl* a (fl+ 1.0 (flexpt 2.0 27.0)))]
|
||||
[x2 (fl- c (fl- c a))])
|
||||
(values x2 (fl- a x2)))))
|
||||
|
||||
;; =================================================================================================
|
||||
;; Fast monotone addition and subtraction
|
||||
|
||||
;(: fast-mono-fl+/error (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
;; Returns a+b and its rounding error
|
||||
;; Assumes |a| >= |b|
|
||||
(define-syntax-rule (fast-mono-fl+/error a-expr b-expr)
|
||||
(let ([a a-expr] [b b-expr])
|
||||
(let ([x2 (+ a b)])
|
||||
(values x2 (- b (- x2 a))))))
|
||||
|
||||
;(: fast-mono-fl-/error (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
;; Returns a+b and its rounding error
|
||||
;; Assumes |a| >= |b|
|
||||
(define-syntax-rule (fast-mono-fl-/error a-expr b-expr)
|
||||
(let ([a a-expr] [b b-expr])
|
||||
(let ([x2 (- a b)])
|
||||
(values x2 (- (- a x2) b)))))
|
||||
|
||||
;; =================================================================================================
|
||||
;; Fast arithmetic that returns rounding error
|
||||
|
||||
;(: fast-fl+/error (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
;; Returns a+b and its rounding error
|
||||
(define-syntax-rule (fast-fl+/error a-expr b-expr)
|
||||
(let ([a a-expr] [b b-expr])
|
||||
(let* ([x2 (fl+ a b)]
|
||||
[v (fl- x2 a)])
|
||||
(values x2 (fl+ (fl- a (fl- x2 v)) (fl- b v))))))
|
||||
|
||||
;(: fast-fl-/error (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
;; Returns a-b and its rounding error
|
||||
(define-syntax-rule (fast-fl-/error a-expr b-expr)
|
||||
(let ([a a-expr] [b b-expr])
|
||||
(let* ([x2 (fl- a b)]
|
||||
[v (fl- x2 a)])
|
||||
(values x2 (fl- (fl- a (fl- x2 v)) (fl+ b v))))))
|
||||
|
||||
;(: fast-fl*/error (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
;; Returns a*b and its rounding error
|
||||
(define-syntax-rule (fast-fl*/error a-expr b-expr)
|
||||
(let ([a a-expr] [b b-expr])
|
||||
(let*-values ([(x2) (fl* a b)]
|
||||
[(a2 a1) (flsplit a)]
|
||||
[(b2 b1) (flsplit b)])
|
||||
(values x2 (- (fl- (fl- (fl- (fl- x2 (fl* a2 b2))
|
||||
(fl* a1 b2))
|
||||
(fl* a2 b1))
|
||||
(fl* a1 b1)))))))
|
||||
|
||||
;(: fast-flfma/error (Flonum Flonum Flonum -> (Values Flonum Flonum)))
|
||||
;; Returns a*b+c and its rounding error
|
||||
(define-syntax-rule (fast-flfma/error a-expr b-expr c-expr)
|
||||
(let*-values ([(y2 y1) (fast-fl*/error a-expr b-expr)]
|
||||
[(h0 h1) (fast-fl+/error c-expr y1)]
|
||||
[(h3 h2) (fast-fl+/error h0 y2)])
|
||||
(values h3 (fl+ h2 h1))))
|
||||
|
||||
#;; If we had hardware fused multiply-add:
|
||||
(define-syntax-rule (fast-fl*/error a-expr b-expr)
|
||||
(let ([a a-expr] [b b-expr])
|
||||
(let ([x2 (fl* a b)])
|
||||
(values x2 (flfma a b (- x2))))))
|
||||
|
||||
;(: fast-flsqr/error (Flonum -> (Values Flonum Flonum)))
|
||||
;; Returns a*a and its rounding error
|
||||
(define-syntax-rule (fast-flsqr/error a-expr)
|
||||
(let ([a a-expr])
|
||||
(let*-values ([(x2) (fl* a a)]
|
||||
[(a2 a1) (flsplit a)])
|
||||
(values x2 (- (fl- (fl- (fl- x2 (fl* a2 a2))
|
||||
(fl* 2.0 (fl* a2 a1)))
|
||||
(fl* a1 a1)))))))
|
||||
|
||||
;(: fast-fl//error (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
;; Returns a/b and its rounding error
|
||||
(define-syntax-rule (fast-fl//error a-expr b-expr)
|
||||
(let ([a a-expr] [b b-expr])
|
||||
(let*-values ([(x2) (fl/ a b)]
|
||||
[(w2 w1) (fast-fl*/error x2 b)])
|
||||
(fast-mono-fl+/error x2 (fl/ (fl- (fl- a w2) w1) b)))))
|
||||
|
||||
) ; module
|
||||
|
||||
(require (submod "." untyped-defs))
|
||||
|
||||
(module typed-defs typed/racket/base
|
||||
(require racket/performance-hint
|
||||
(submod ".." untyped-defs)
|
||||
"flonum-functions.rkt"
|
||||
"utils.rkt")
|
||||
|
||||
(provide (all-defined-out))
|
||||
|
||||
;; =================================================================================================
|
||||
;; Function versions of the above that are well-defined for the largest domain, and return 0.0 as
|
||||
;; the second argument whenever the first isn't rational
|
||||
|
||||
(begin-encourage-inline
|
||||
|
||||
(: fl+/error (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
(define (fl+/error a b)
|
||||
(let-values ([(x2 x1) (fast-fl+/error a b)])
|
||||
(values x2 (if (flrational? x2) x1 0.0))))
|
||||
|
||||
(: fl-/error (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
(define (fl-/error a b)
|
||||
(let-values ([(x2 x1) (fast-fl-/error a b)])
|
||||
(values x2 (if (flrational? x2) x1 0.0))))
|
||||
|
||||
(: fl*/error (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
(define (fl*/error a b)
|
||||
(let ([x2 (fl* a b)])
|
||||
(values x2 (if (and (flrational? x2) (not (flsubnormal? x2)))
|
||||
(let*-values ([(da db) (values (near-pow2 a) (near-pow2 b))]
|
||||
[(d) (fl* da db)]
|
||||
[(a2 a1) (flsplit (fl/ a da))]
|
||||
[(b2 b1) (flsplit (fl/ b db))])
|
||||
(fl* d (- (fl- (fl- (fl- (fl- (fl/ x2 d) (fl* a2 b2))
|
||||
(fl* a1 b2))
|
||||
(fl* a2 b1))
|
||||
(fl* a1 b1)))))
|
||||
0.0))))
|
||||
|
||||
(: flsqr/error (Flonum -> (Values Flonum Flonum)))
|
||||
(define (flsqr/error a)
|
||||
(let ([x2 (fl* a a)])
|
||||
(values x2 (if (and (flrational? x2) (not (flsubnormal? x2)))
|
||||
(let*-values ([(d) (near-pow2 a)]
|
||||
[(d^2) (fl* d d)]
|
||||
[(a2 a1) (flsplit (fl/ a d))])
|
||||
(fl* d^2 (- (fl- (fl- (fl- (fl/ x2 d^2) (fl* a2 a2))
|
||||
(fl* 2.0 (fl* a1 a2)))
|
||||
(fl* a1 a1)))))
|
||||
0.0))))
|
||||
|
||||
(: fl//error (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
(define (fl//error a b)
|
||||
(let ([x2 (fl/ a b)])
|
||||
(values x2 (if (and (flrational? x2) (flrational? b))
|
||||
(let* ([d (near-pow2/div a b)]
|
||||
[a (fl/ a d)]
|
||||
[b (fl/ b d)])
|
||||
(let-values ([(w2 w1) (fl*/error x2 b)])
|
||||
(fl/ (fl- (fl- a w2) w1) b)))
|
||||
0.0))))
|
||||
|
||||
) ; begin-encourage-inline
|
||||
|
||||
) ; module
|
||||
|
||||
(require (submod "." typed-defs))
|
||||
|
|
@ -8,7 +8,7 @@
|
|||
"flonum-functions.rkt"
|
||||
"flonum-log.rkt"
|
||||
"flonum-more-functions.rkt"
|
||||
"flonum-syntax.rkt")
|
||||
"flonum-error.rkt")
|
||||
|
||||
(provide flfactorial
|
||||
flbinomial
|
||||
|
|
|
|||
|
|
@ -83,13 +83,15 @@
|
|||
|
||||
(: flulp-error (Flonum Real -> Flonum))
|
||||
(define (flulp-error x r)
|
||||
(define r.0 (fl r))
|
||||
(cond [(eqv? x r) 0.0]
|
||||
[(and (fl= x 0.0) (zero? r)) 0.0]
|
||||
[(and (fl= x 0.0) (fl= r.0 0.0)) 0.0]
|
||||
[(and (fl= x +inf.0) (fl= r.0 +inf.0)) 0.0]
|
||||
[(and (fl= x -inf.0) (fl= r.0 -inf.0)) 0.0]
|
||||
[(zero? r) +inf.0]
|
||||
[(and (rational? x) (rational? r))
|
||||
(flabs (real->double-flonum
|
||||
(/ (- (inexact->exact x) (inexact->exact r))
|
||||
(inexact->exact (flulp x)))))]
|
||||
[(and (flrational? x) (flrational? r.0))
|
||||
(flabs (fl (/ (- (inexact->exact x) (inexact->exact r))
|
||||
(inexact->exact (flmax +min.0 (flulp r.0))))))]
|
||||
[else +inf.0]))
|
||||
|
||||
;; ===================================================================================================
|
||||
|
|
|
|||
|
|
@ -5,7 +5,7 @@
|
|||
"flonum-functions.rkt"
|
||||
"flonum-constants.rkt"
|
||||
"flonum-exp.rkt"
|
||||
"flonum-syntax.rkt"
|
||||
"flonum-error.rkt"
|
||||
"flvector.rkt")
|
||||
|
||||
(provide fllog1p fllog+
|
||||
|
|
|
|||
|
|
@ -5,7 +5,7 @@
|
|||
"flonum-constants.rkt"
|
||||
"flonum-exp.rkt"
|
||||
"flonum-log.rkt"
|
||||
"flonum-syntax.rkt"
|
||||
"flonum-error.rkt"
|
||||
"flvector.rkt")
|
||||
|
||||
(provide flsqrt1pm1
|
||||
|
|
|
|||
|
|
@ -1,105 +0,0 @@
|
|||
#lang racket/base
|
||||
|
||||
(require (for-syntax racket/base)
|
||||
racket/flonum)
|
||||
|
||||
(provide flsplit
|
||||
fast-mono-fl+/error fast-fl+/error fl+/error
|
||||
fast-mono-fl-/error fast-fl-/error fl-/error
|
||||
fast-fl*/error fast-flsqr/error fl*/error
|
||||
fast-fl//error)
|
||||
|
||||
;(: flsplit (Flonum -> (Values Flonum Flonum)))
|
||||
;; Splits a flonum into a two flonums `hi' and `lo' with 26 bits precision each, such that
|
||||
;; |hi| >= |lo| and hi + lo = a. (The extra sign bit accounts for the missing bit.)
|
||||
(define-syntax-rule (flsplit a-expr)
|
||||
(let ([a a-expr])
|
||||
(let* ([s (if ((flabs a) . fl< . 1e300) 1.0 (flexpt 2.0 52.0))]
|
||||
[a (fl/ a s)]
|
||||
[c (fl* a (fl+ 1.0 (flexpt 2.0 27.0)))]
|
||||
[a-hi (fl- c (fl- c a))]
|
||||
[a-lo (fl- a a-hi)])
|
||||
(values (fl* a-hi s)
|
||||
(fl* a-lo s)))))
|
||||
|
||||
;(: fast-mono-fl+/error (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
;; Returns a+b and its rounding error
|
||||
;; Assumes |a| >= |b|
|
||||
(define-syntax-rule (fast-mono-fl+/error a-expr b-expr)
|
||||
(let ([a a-expr] [b b-expr])
|
||||
(let ([x (+ a b)])
|
||||
(values x (- b (- x a))))))
|
||||
|
||||
;(: fast-mono-fl-/error (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
;; Returns a+b and its rounding error
|
||||
;; Assumes |a| >= |b|
|
||||
(define-syntax-rule (fast-mono-fl-/error a-expr b-expr)
|
||||
(let ([a a-expr] [b b-expr])
|
||||
(let ([x (- a b)])
|
||||
(values x (- (- a x) b)))))
|
||||
|
||||
;(: fast-fl+/error (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
;; Returns a+b and its rounding error
|
||||
(define-syntax-rule (fast-fl+/error a-expr b-expr)
|
||||
(let ([a a-expr] [b b-expr])
|
||||
(let* ([x (fl+ a b)]
|
||||
[v (fl- x a)])
|
||||
(values x (fl+ (fl- a (fl- x v)) (fl- b v))))))
|
||||
|
||||
;(: fast-fl-/error (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
;; Returns a-b and its rounding error
|
||||
(define-syntax-rule (fast-fl-/error a-expr b-expr)
|
||||
(let ([a a-expr] [b b-expr])
|
||||
(let* ([x (fl- a b)]
|
||||
[v (fl- x a)])
|
||||
(values x (fl- (fl- a (fl- x v)) (fl+ b v))))))
|
||||
|
||||
;(: fast-fl*/error (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
;; Returns a*b and its rounding error
|
||||
(define-syntax-rule (fast-fl*/error a-expr b-expr)
|
||||
(let ([a a-expr] [b b-expr])
|
||||
(let*-values ([(x) (fl* a b)]
|
||||
[(a-hi a-lo) (flsplit a)]
|
||||
[(b-hi b-lo) (flsplit b)])
|
||||
(values x (- (fl- (fl- (fl- (fl- x (fl* a-hi b-hi))
|
||||
(fl* a-lo b-hi))
|
||||
(fl* a-hi b-lo))
|
||||
(fl* a-lo b-lo)))))))
|
||||
|
||||
(define-syntax-rule (fast-flsqr/error a-expr)
|
||||
(let ([a a-expr])
|
||||
(let*-values ([(x) (fl* a a)]
|
||||
[(a-hi a-lo) (flsplit a)])
|
||||
(values x (- (fl- (fl- (fl- x (fl* a-hi a-hi))
|
||||
(fl* 2.0 (fl* a-hi a-lo)))
|
||||
(fl* a-lo a-lo)))))))
|
||||
|
||||
;(: fast-fl//error (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
;; Returns a/b and its rounding error
|
||||
(define-syntax-rule (fast-fl//error a-expr b-expr)
|
||||
(let ([a a-expr] [b b-expr])
|
||||
(let*-values ([(q-hi) (fl/ a b)]
|
||||
[(q0 q1) (fast-fl*/error q-hi b)])
|
||||
(values q-hi (fl/ (fl+ (- q1) (fl- a q0)) b)))))
|
||||
|
||||
#;; If we had a fused multiply-add:
|
||||
(define (fast-fl*/error a b)
|
||||
(let ([p (* a b)])
|
||||
(values p (flfma a b (- p)))))
|
||||
|
||||
(define-syntax-rule (define-slow-wrapper name fast-flop/error op)
|
||||
(define-syntax-rule (name a-expr b-expr)
|
||||
(let ([a a-expr] [b b-expr])
|
||||
(let-values ([(x e) (fast-flop/error a b)])
|
||||
(if (and (x . > . -inf.0) (x . < . +inf.0)
|
||||
(e . > . -inf.0) (e . < . +inf.0))
|
||||
(values x e)
|
||||
(let* ([v (op (inexact->exact a) (inexact->exact b))]
|
||||
[x (real->double-flonum v)])
|
||||
(if (and (x . > . -inf.0) (x . < . +inf.0))
|
||||
(values x (real->double-flonum (- v (inexact->exact x))))
|
||||
(values x 0.0))))))))
|
||||
|
||||
(define-slow-wrapper fl+/error fast-fl+/error +)
|
||||
(define-slow-wrapper fl-/error fast-fl-/error -)
|
||||
(define-slow-wrapper fl*/error fast-fl*/error *)
|
||||
22
collects/math/private/flonum/utils.rkt
Normal file
22
collects/math/private/flonum/utils.rkt
Normal file
|
|
@ -0,0 +1,22 @@
|
|||
#lang racket/base
|
||||
|
||||
(provide near-pow2
|
||||
near-pow2/div)
|
||||
|
||||
(module untyped-defs racket/base
|
||||
(require racket/flonum)
|
||||
|
||||
(provide (all-defined-out))
|
||||
|
||||
(define-syntax-rule (near-pow2 a)
|
||||
(flexpt 2.0 (flmax -511.0 (flmin 511.0 (flround (fl/ (fllog (flabs a)) (fllog 2.0)))))))
|
||||
|
||||
(define-syntax-rule (near-pow2/div a b)
|
||||
;; Clamping both values makes this work properly when a or b is infinite or zero
|
||||
(let ([ea (flmax -511.0 (flmin 511.0 (fl/ (fllog (flabs a)) (fllog 2.0))))]
|
||||
[eb (flmax -511.0 (flmin 511.0 (fl/ (fllog (flabs b)) (fllog 2.0))))])
|
||||
(flexpt 2.0 (flround (fl* 0.5 (fl+ ea eb))))))
|
||||
|
||||
) ; module
|
||||
|
||||
(require (submod "." untyped-defs))
|
||||
418
collects/math/private/utils/flonum-tests.rkt
Normal file
418
collects/math/private/utils/flonum-tests.rkt
Normal file
|
|
@ -0,0 +1,418 @@
|
|||
#lang typed/racket/base
|
||||
|
||||
(require racket/math
|
||||
racket/flonum
|
||||
racket/list
|
||||
typed/rackunit
|
||||
"../base/base-random.rkt"
|
||||
"../flonum/expansion/expansion-base.rkt"
|
||||
"../flonum/flonum-functions.rkt"
|
||||
"../flonum/flonum-constants.rkt"
|
||||
"../flonum/flonum-bits.rkt"
|
||||
"../flonum/flonum-error.rkt"
|
||||
"../distributions/dist-struct.rkt"
|
||||
"../distributions/geometric-dist.rkt"
|
||||
"../../bigfloat.rkt")
|
||||
|
||||
(provide print-test-progress?
|
||||
test-fpu-arith
|
||||
test-fpu-trig
|
||||
test-fpu-non-trig
|
||||
test-fpu-arith/error
|
||||
test-fpu-arith/fl2
|
||||
test-fpu-non-trig/fl2
|
||||
test-fpu)
|
||||
|
||||
;; Allowable error for different kinds of functions, in ulps
|
||||
(define flonum-fun-ulps 0.5)
|
||||
(define flonum/error-fun-ulps 0.5)
|
||||
(define unary-fl2-fun-ulps 1.0)
|
||||
(define binary-fl2-fun-ulps 8.0)
|
||||
|
||||
;; ===================================================================================================
|
||||
;; Helpers
|
||||
|
||||
(: different-zero? (Real Real -> Boolean))
|
||||
(define (different-zero? x y)
|
||||
(or (and (eqv? x -0.0) (eqv? y 0.0))
|
||||
(and (eqv? x 0.0) (eqv? y -0.0))))
|
||||
|
||||
(: fl2->real* (Flonum Flonum -> Real))
|
||||
;; Like `fl2->real', but returns signed flonum zeros
|
||||
(define (fl2->real* x2 x1)
|
||||
(define x.0 (fl+ x2 x1))
|
||||
(cond [(zero? x.0) x2]
|
||||
[else (fl2->real x2 x1)]))
|
||||
|
||||
(: bigfloat->real* (Bigfloat -> Real))
|
||||
;; Like `bigfloat->real*', but returns a signed infinity or a signed flonum zero if conversion would
|
||||
;; overflow or underflow a flonum
|
||||
(define (bigfloat->real* x)
|
||||
(define x.0 (bigfloat->flonum x))
|
||||
(cond [(fl= x.0 0.0) x.0]
|
||||
[(flrational? x.0) (bigfloat->real x)]
|
||||
[else x.0]))
|
||||
|
||||
(: filter/ulp-error (All (A B) ((Listof (List A (U B Flonum)))
|
||||
Flonum -> (Listof (List A (U B Flonum))))))
|
||||
(define (filter/ulp-error xes ulps)
|
||||
(filter (λ: ([xe : (List A (U B Flonum))])
|
||||
(define e (second xe))
|
||||
(or (not (flonum? e)) (e . fl> . ulps)))
|
||||
xes))
|
||||
|
||||
(: print-test-progress? (Parameterof Boolean))
|
||||
(define print-test-progress? (make-parameter #t))
|
||||
|
||||
(define progress-chunk-size 200)
|
||||
(define progress-superchunk-chunks 5)
|
||||
|
||||
(: maybe-print-progress (Symbol Integer Natural -> Void))
|
||||
(define (maybe-print-progress name i m)
|
||||
(when (and (print-test-progress?) (i . > . 0) (i . <= . m))
|
||||
(let* ([flush? (cond [(= i 1) (printf "~a: " name)]
|
||||
[else #f])]
|
||||
[flush? (cond [(= 0 (modulo i progress-chunk-size))
|
||||
(cond [(= 0 (modulo i (* progress-superchunk-chunks
|
||||
progress-chunk-size)))
|
||||
(printf "* ~a " i)]
|
||||
[else (printf "*")])]
|
||||
[else flush?])]
|
||||
[flush? (cond [(= i m) (printf "* ~a~n" m)]
|
||||
[else flush?])])
|
||||
(when flush? (flush-output)))))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; Test case generation
|
||||
|
||||
;; Deteriministic test cases
|
||||
|
||||
(define standard-xs
|
||||
(list
|
||||
;; Test the sign of the return value of `flexpt'
|
||||
-1001.0 -10.0 -0.1 +0.1 +10.0 +1001.0
|
||||
;; Test squaring
|
||||
(- (flsqrt +min.0)) (flsqrt +min.0)
|
||||
(- (flsqrt +max-subnormal.0)) (flsqrt +max-subnormal.0)
|
||||
(- (flsqrt +max-fl2-subnormal.0)) (flsqrt +max-fl2-subnormal.0)
|
||||
(- (flsqrt +max.0)) (flsqrt +max.0)
|
||||
;; Test exp limits
|
||||
(fllog +min.0) (fllog +max-subnormal.0) (fllog +max-fl2-subnormal.0) (fllog +max.0)
|
||||
;; Standard special values
|
||||
-inf.0 -max.0 -1.0 -max-fl2-subnormal.0 -max-subnormal.0 -min.0 -0.0
|
||||
+inf.0 +max.0 +1.0 +max-fl2-subnormal.0 +max-subnormal.0 +min.0 +0.0
|
||||
+nan.0))
|
||||
|
||||
(define standard-rs
|
||||
(append standard-xs
|
||||
(list +10 +1 +1/7 +1/10 +1/13
|
||||
-10 -1 -1/7 -1/10 -1/13
|
||||
0)))
|
||||
|
||||
(: product (All (A B) ((Listof A) (Listof B) -> (Values (Listof A) (Listof B)))))
|
||||
(define (product as bs)
|
||||
(define abs
|
||||
(append*
|
||||
(for/list: : (Listof (Listof (Pair A B))) ([a (in-list as)])
|
||||
(for/list: : (Listof (Pair A B)) ([b (in-list bs)])
|
||||
(cons a b)))))
|
||||
(values (map (inst car A B) abs) (map (inst cdr A B) abs)))
|
||||
|
||||
;; Random test cases
|
||||
|
||||
(define min-subnormal-ord (flonum->ordinal -max-subnormal.0))
|
||||
(define max-subnormal-ord (+ 1 (flonum->ordinal +max-subnormal.0)))
|
||||
|
||||
(define min-fl2-subnormal-ord (flonum->ordinal -max-fl2-subnormal.0))
|
||||
(define max-fl2-subnormal-ord (+ 1 (flonum->ordinal +max-fl2-subnormal.0)))
|
||||
|
||||
(: sample-flonum (case-> (Integer -> (Listof Flonum))
|
||||
(Integer Flonum Flonum -> (Listof Flonum))))
|
||||
(define (sample-flonum n [mn -inf.0] [mx +inf.0])
|
||||
(define min-ord (flonum->ordinal mn))
|
||||
(define max-ord (+ 1 (flonum->ordinal mx)))
|
||||
(let ([min-subnormal-ord (max min-ord min-subnormal-ord)]
|
||||
[max-subnormal-ord (min max-ord max-subnormal-ord)]
|
||||
[min-fl2-subnormal-ord (max min-ord min-fl2-subnormal-ord)]
|
||||
[max-fl2-subnormal-ord (min max-ord max-fl2-subnormal-ord)])
|
||||
(build-list
|
||||
n (λ (_)
|
||||
(define r (random))
|
||||
(ordinal->flonum
|
||||
(cond [(and (min-subnormal-ord . < . max-subnormal-ord) (r . < . 0.1))
|
||||
(random-integer min-subnormal-ord max-subnormal-ord)]
|
||||
[(and (min-fl2-subnormal-ord . < . max-fl2-subnormal-ord) (r . < . 0.2))
|
||||
(random-integer min-subnormal-ord max-subnormal-ord)]
|
||||
[else
|
||||
(random-integer min-ord max-ord)]))))))
|
||||
|
||||
(define denom-dist (geometric-dist 1e-32))
|
||||
|
||||
(: sample-rational (Integer -> (Listof Exact-Rational)))
|
||||
(define (sample-rational n)
|
||||
(map (λ: ([f1 : Flonum] [d : Integer])
|
||||
(+ (inexact->exact f1)
|
||||
(* (if ((random) . > . 0.5) -1 1)
|
||||
(/ (random-natural (+ 1 d)) d)
|
||||
(expt 2 (- (exact-round (/ (fllog (flabs f1)) (fllog 2.0))) 52)))))
|
||||
(sample-flonum n)
|
||||
(map (λ: ([x : Flonum]) (+ 1 (exact-floor x))) (sample denom-dist n))))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; Flonum functions
|
||||
|
||||
(: flonum-error (Flonum Bigfloat -> Any))
|
||||
(define (flonum-error z z0.bf)
|
||||
(define z0 (bigfloat->real* z0.bf))
|
||||
(cond [(different-zero? z z0) (list 'different-zero? z z0)]
|
||||
[else (flulp-error z z0)]))
|
||||
|
||||
(: unary-flonum-fun-error ((Flonum -> Flonum) (Bigfloat -> Bigfloat) Flonum -> Any))
|
||||
(define (unary-flonum-fun-error f g x)
|
||||
(flonum-error (f x) (parameterize ([bf-precision 53])
|
||||
(g (bf x)))))
|
||||
|
||||
(: test-unary-flonum-fun
|
||||
(Symbol (Flonum -> Flonum) (Bigfloat -> Bigfloat) Integer Flonum Flonum
|
||||
-> (Listof (List (List Symbol Flonum) Any))))
|
||||
(define (test-unary-flonum-fun name f g n mn mx)
|
||||
(define xs (append standard-xs (sample-flonum n mn mx)))
|
||||
(define m (length xs))
|
||||
(filter/ulp-error
|
||||
(for/list: : (Listof (List (List Symbol Flonum) Any)) ([x (in-list xs)]
|
||||
[i (in-naturals 1)])
|
||||
(maybe-print-progress name i m)
|
||||
(list (list name x) (unary-flonum-fun-error f g x)))
|
||||
flonum-fun-ulps))
|
||||
|
||||
(: binary-flonum-fun-error
|
||||
((Flonum Flonum -> Flonum) (Bigfloat Bigfloat -> Bigfloat) Flonum Flonum -> Any))
|
||||
(define (binary-flonum-fun-error f g x y)
|
||||
(flonum-error (f x y) (parameterize ([bf-precision 53])
|
||||
(g (bf x) (bf y)))))
|
||||
|
||||
(: test-binary-flonum-fun
|
||||
(Symbol (Flonum Flonum -> Flonum) (Bigfloat Bigfloat -> Bigfloat) Integer
|
||||
-> (Listof (List (List Symbol Flonum Flonum) Any))))
|
||||
(define (test-binary-flonum-fun name f g n)
|
||||
(define-values (pre-xs pre-ys) (product standard-xs standard-xs))
|
||||
(define xs (append pre-xs (sample-flonum n)))
|
||||
(define ys (append pre-ys (sample-flonum n)))
|
||||
(define m (length xs))
|
||||
(filter/ulp-error
|
||||
(for/list: : (Listof (List (List Symbol Flonum Flonum) Any)) ([x (in-list xs)]
|
||||
[y (in-list ys)]
|
||||
[i (in-naturals 1)])
|
||||
(maybe-print-progress name i m)
|
||||
(list (list name x y) (binary-flonum-fun-error f g x y)))
|
||||
flonum-fun-ulps))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; fl2 conversion
|
||||
|
||||
(: fl2-error (Flonum Flonum Real -> Any))
|
||||
(define (fl2-error x2 x1 x)
|
||||
(cond [(not (fl2? x2 x1)) (list 'not-fl2? x2 x1)]
|
||||
[(different-zero? x2 x) (list 'different-zero? x2 x)]
|
||||
[else (fl2ulp-error x2 x1 x)]))
|
||||
|
||||
(: fl2-conversion-error (Real -> Any))
|
||||
(define (fl2-conversion-error x)
|
||||
(define-values (x2 x1) (fl2 x))
|
||||
(fl2-error x2 x1 x))
|
||||
|
||||
(: test-fl2-conversion (Integer -> (Listof (List (List 'fl2 Real) Any))))
|
||||
(define (test-fl2-conversion n)
|
||||
(define xs (append standard-rs (sample-rational n)))
|
||||
(define m (length xs))
|
||||
(filter/ulp-error
|
||||
(for/list: : (Listof (List (List 'fl2 Real) Any)) ([x (in-list xs)]
|
||||
[i (in-naturals 1)])
|
||||
(maybe-print-progress 'fl2 i m)
|
||||
(list (list 'fl2 x) (fl2-conversion-error x)))
|
||||
flonum/error-fun-ulps))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; Flonum arithmetic with error
|
||||
|
||||
(: unary-flonum/error-fun-error ((Flonum -> (Values Flonum Flonum)) (Bigfloat -> Bigfloat) Flonum
|
||||
-> Any))
|
||||
(define (unary-flonum/error-fun-error f g x)
|
||||
(define-values (z2 z1) (f x))
|
||||
(fl2-error z2 z1 (parameterize ([bf-precision 256])
|
||||
(bigfloat->real* (g (bf x))))))
|
||||
|
||||
(: binary-flonum/error-fun-error ((Flonum Flonum -> (Values Flonum Flonum))
|
||||
(Bigfloat Bigfloat -> Bigfloat)
|
||||
Flonum Flonum
|
||||
-> Any))
|
||||
(define (binary-flonum/error-fun-error f g x y)
|
||||
(define-values (z2 z1) (f x y))
|
||||
(fl2-error z2 z1 (parameterize ([bf-precision 256])
|
||||
(bigfloat->real* (g (bf x) (bf y))))))
|
||||
|
||||
(: test-unary-flonum/error-fun
|
||||
(Symbol (Flonum -> (Values Flonum Flonum)) (Bigfloat -> Bigfloat) Integer
|
||||
-> (Listof (List (List Symbol Flonum) Any))))
|
||||
(define (test-unary-flonum/error-fun name f g n)
|
||||
(define xs (append standard-xs (sample-flonum n)))
|
||||
(define m (length xs))
|
||||
(filter/ulp-error
|
||||
(for/list: : (Listof (List (List Symbol Flonum) Any)) ([x (in-list xs)]
|
||||
[i (in-naturals 1)])
|
||||
(maybe-print-progress name i m)
|
||||
(list (list name x) (unary-flonum/error-fun-error f g x)))
|
||||
flonum/error-fun-ulps))
|
||||
|
||||
(: test-binary-flonum/error-fun
|
||||
(Symbol (Flonum Flonum -> (Values Flonum Flonum)) (Bigfloat Bigfloat -> Bigfloat) Integer
|
||||
-> (Listof (List (List Symbol Flonum Flonum) Any))))
|
||||
(define (test-binary-flonum/error-fun name f g n)
|
||||
(define-values (pre-xs pre-ys) (product standard-xs standard-xs))
|
||||
(define xs (append pre-xs (sample-flonum n)))
|
||||
(define ys (append pre-ys (sample-flonum n)))
|
||||
(define m (length xs))
|
||||
(filter/ulp-error
|
||||
(for/list: : (Listof (List (List Symbol Flonum Flonum) Any)) ([x (in-list xs)]
|
||||
[y (in-list ys)]
|
||||
[i (in-naturals 1)])
|
||||
(maybe-print-progress name i m)
|
||||
(list (list name x y) (binary-flonum/error-fun-error f g x y)))
|
||||
flonum/error-fun-ulps))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; Flonum expansions
|
||||
|
||||
(: unary-fl2-fun-error ((Flonum Flonum -> (Values Flonum Flonum)) (Bigfloat -> Bigfloat)
|
||||
Flonum Flonum -> Any))
|
||||
(define (unary-fl2-fun-error f g x2 x1)
|
||||
(define-values (z2 z1) (f x2 x1))
|
||||
(fl2-error z2 z1 (parameterize ([bf-precision 256])
|
||||
(bigfloat->real* (g (bf (fl2->real* x2 x1)))))))
|
||||
|
||||
(: test-unary-fl2-fun
|
||||
(Symbol (Flonum Flonum -> (Values Flonum Flonum)) (Bigfloat -> Bigfloat) Integer
|
||||
-> (Listof (List (List Symbol Flonum Flonum) Any))))
|
||||
(define (test-unary-fl2-fun name f g n)
|
||||
(define xs (append standard-rs (sample-rational n)))
|
||||
(define m (length xs))
|
||||
(filter/ulp-error
|
||||
(for/list: : (Listof (List (List Symbol Flonum Flonum) Any)) ([x (in-list xs)]
|
||||
[i (in-naturals 1)])
|
||||
(maybe-print-progress name i m)
|
||||
(define-values (x2 x1) (fl2 x))
|
||||
(list (list name x2 x1) (unary-fl2-fun-error f g x2 x1)))
|
||||
unary-fl2-fun-ulps))
|
||||
|
||||
(: binary-fl2-fun-error ((Flonum Flonum Flonum Flonum -> (Values Flonum Flonum))
|
||||
(Bigfloat Bigfloat -> Bigfloat)
|
||||
Flonum Flonum Flonum Flonum
|
||||
-> Any))
|
||||
(define (binary-fl2-fun-error f g x2 x1 y2 y1)
|
||||
(define-values (z2 z1) (f x2 x1 y2 y1))
|
||||
(fl2-error z2 z1 (parameterize ([bf-precision 256])
|
||||
(bigfloat->real* (g (bf (fl2->real* x2 x1)) (bf (fl2->real* y2 y1)))))))
|
||||
|
||||
(: test-binary-fl2-fun
|
||||
(Symbol (Flonum Flonum Flonum Flonum -> (Values Flonum Flonum)) (Bigfloat Bigfloat -> Bigfloat)
|
||||
Integer -> (Listof (List (List Symbol Flonum Flonum Flonum Flonum) Any))))
|
||||
(define (test-binary-fl2-fun name f g n)
|
||||
(define-values (pre-xs pre-ys) (product standard-rs standard-rs))
|
||||
(define xs (append pre-xs (sample-rational n)))
|
||||
(define ys (append pre-ys (sample-rational n)))
|
||||
(define m (length xs))
|
||||
(filter/ulp-error
|
||||
(for/list: : (Listof (List (List Symbol Flonum Flonum Flonum Flonum) Any)
|
||||
) ([x (in-list xs)]
|
||||
[y (in-list ys)]
|
||||
[i (in-naturals 1)])
|
||||
(maybe-print-progress name i m)
|
||||
(define-values (x2 x1) (fl2 x))
|
||||
(define-values (y2 y1) (fl2 y))
|
||||
(list (list name x2 x1 y2 y1) (binary-fl2-fun-error f g x2 x1 y2 y1)))
|
||||
binary-fl2-fun-ulps))
|
||||
|
||||
;; ===================================================================================================
|
||||
|
||||
(: test-fpu-arith (Natural -> Any))
|
||||
(define (test-fpu-arith n)
|
||||
(check-equal? (test-unary-flonum-fun 'flabs flabs bfabs n -inf.0 +inf.0)
|
||||
'())
|
||||
(check-equal? (test-binary-flonum-fun 'fl+ fl+ bf+ n)
|
||||
'())
|
||||
(check-equal? (test-binary-flonum-fun 'fl- fl- bf- n)
|
||||
'())
|
||||
(check-equal? (test-binary-flonum-fun 'fl* fl* bf* n)
|
||||
'())
|
||||
(check-equal? (test-binary-flonum-fun 'fl/ fl/ bf/ n)
|
||||
'()))
|
||||
|
||||
(: test-fpu-trig (Natural -> Any))
|
||||
(define (test-fpu-trig n)
|
||||
(check-equal? (test-unary-flonum-fun 'flsin flsin bfsin n -inf.0 +inf.0)
|
||||
'())
|
||||
(check-equal? (test-unary-flonum-fun 'flcos flcos bfcos n -inf.0 +inf.0)
|
||||
'())
|
||||
(check-equal? (test-unary-flonum-fun 'fltan fltan bftan n -inf.0 +inf.0)
|
||||
'())
|
||||
(check-equal? (test-unary-flonum-fun 'flasin flasin bfasin n -1.0 1.0)
|
||||
'())
|
||||
(check-equal? (test-unary-flonum-fun 'flacos flacos bfacos n -1.0 1.0)
|
||||
'())
|
||||
(check-equal? (test-unary-flonum-fun 'flatan flatan bfatan n -inf.0 +inf.0)
|
||||
'()))
|
||||
|
||||
(: test-fpu-non-trig (Natural -> Any))
|
||||
(define (test-fpu-non-trig n)
|
||||
(check-equal? (test-unary-flonum-fun 'flsqrt flsqrt bfsqrt n 0.0 +inf.0)
|
||||
'())
|
||||
(check-equal? (test-unary-flonum-fun 'fllog fllog bflog n 0.0 +inf.0)
|
||||
'())
|
||||
(check-equal? (test-unary-flonum-fun 'flexp flexp bfexp n -746.0 710.0)
|
||||
'())
|
||||
(check-equal? (test-binary-flonum-fun 'flexpt flexpt bfexpt n)
|
||||
'()))
|
||||
|
||||
(: test-fpu-arith/error (Natural -> Any))
|
||||
(define (test-fpu-arith/error n)
|
||||
(check-equal? (test-binary-flonum/error-fun 'fl+/error fl+/error bf+ n)
|
||||
'())
|
||||
(check-equal? (test-binary-flonum/error-fun 'fl-/error fl-/error bf- n)
|
||||
'())
|
||||
(check-equal? (test-binary-flonum/error-fun 'fl*/error fl*/error bf* n)
|
||||
'())
|
||||
(check-equal? (test-unary-flonum/error-fun 'flsqr/error flsqr/error bfsqr n)
|
||||
'())
|
||||
(check-equal? (test-binary-flonum/error-fun 'fl//error fl//error bf/ n)
|
||||
'()))
|
||||
|
||||
(: test-fpu-arith/fl2 (Natural -> Any))
|
||||
(define (test-fpu-arith/fl2 n)
|
||||
(check-equal? (test-fl2-conversion n)
|
||||
'())
|
||||
(check-equal? (test-unary-fl2-fun 'fl2abs fl2abs bfabs n)
|
||||
'())
|
||||
(check-equal? (test-binary-fl2-fun 'fl2+ fl2+ bf+ n)
|
||||
'())
|
||||
(check-equal? (test-binary-fl2-fun 'fl2- fl2- bf- n)
|
||||
'())
|
||||
(check-equal? (test-binary-fl2-fun 'fl2* fl2* bf* n)
|
||||
'())
|
||||
(check-equal? (test-unary-fl2-fun 'fl2sqr fl2sqr bfsqr n)
|
||||
'())
|
||||
(check-equal? (test-binary-fl2-fun 'fl2/ fl2/ bf/ n)
|
||||
'()))
|
||||
|
||||
(: test-fpu-non-trig/fl2 (Natural -> Any))
|
||||
(define (test-fpu-non-trig/fl2 n)
|
||||
(check-equal? (test-unary-fl2-fun 'fl2sqrt fl2sqrt bfsqrt n)
|
||||
'()))
|
||||
|
||||
(: test-fpu (Natural -> Any))
|
||||
(define (test-fpu n)
|
||||
(test-fpu-arith n)
|
||||
(test-fpu-trig n)
|
||||
(test-fpu-non-trig n)
|
||||
(test-fpu-arith/error n)
|
||||
(test-fpu-arith/fl2 n)
|
||||
(test-fpu-non-trig/fl2 n))
|
||||
|
|
@ -1,20 +1,151 @@
|
|||
#lang racket
|
||||
|
||||
(require math/flonum
|
||||
math/base
|
||||
math/utils
|
||||
rackunit)
|
||||
|
||||
(for* ([x '(+inf.0 +nan.0 -inf.0)]
|
||||
[y '(+inf.0 +nan.0 -inf.0)])
|
||||
(cond [(eqv? x y)
|
||||
(check-eqv? (flulp-error x y) 0.0 (format "(flulp-error ~v ~v)" x y))]
|
||||
[else
|
||||
(check-eqv? (flulp-error x y) +inf.0 (format "(flulp-error ~v ~v)" x y))]))
|
||||
;; ===================================================================================================
|
||||
;; Test `flulp-error' heavily; it MUST be correct, or all the FPU tests are suspect
|
||||
|
||||
(check-equal? (flulp-error 55123.135123 55123.135123)
|
||||
0.0)
|
||||
;; Both arguments eqv?
|
||||
(check-equal? (flulp-error -inf.0 -inf.0) 0.0)
|
||||
(check-equal? (flulp-error -max.0 -max.0) 0.0)
|
||||
(check-equal? (flulp-error -1.0 -1.0) 0.0)
|
||||
(check-equal? (flulp-error -min.0 -min.0) 0.0)
|
||||
(check-equal? (flulp-error -0.0 -0.0) 0.0)
|
||||
(check-equal? (flulp-error +0.0 +0.0) 0.0)
|
||||
(check-equal? (flulp-error +min.0 +min.0) 0.0)
|
||||
(check-equal? (flulp-error +1.0 +1.0) 0.0)
|
||||
(check-equal? (flulp-error +max.0 +max.0) 0.0)
|
||||
(check-equal? (flulp-error +inf.0 +inf.0) 0.0)
|
||||
|
||||
(check-equal? (flulp-error 1.0 (flnext 1.0))
|
||||
1.0)
|
||||
;; Both arguments zero
|
||||
(check-equal? (flulp-error -0.0 +0.0) 0.0)
|
||||
(check-equal? (flulp-error +0.0 -0.0) 0.0)
|
||||
|
||||
(check-equal? (flulp-error +max.0 (flprev +max.0))
|
||||
1.0)
|
||||
;; LHS argument +inf.0
|
||||
(check-equal? (flulp-error +inf.0 -inf.0) +inf.0)
|
||||
(check-equal? (flulp-error +inf.0 -max.0) +inf.0)
|
||||
(check-equal? (flulp-error +inf.0 -1.0) +inf.0)
|
||||
(check-equal? (flulp-error +inf.0 -min.0) +inf.0)
|
||||
(check-equal? (flulp-error +inf.0 -0.0) +inf.0)
|
||||
(check-equal? (flulp-error +inf.0 +0.0) +inf.0)
|
||||
(check-equal? (flulp-error +inf.0 +min.0) +inf.0)
|
||||
(check-equal? (flulp-error +inf.0 +1.0) +inf.0)
|
||||
(check-equal? (flulp-error +inf.0 +max.0) +inf.0)
|
||||
|
||||
;; LHS argument -inf.0
|
||||
(check-equal? (flulp-error -inf.0 -max.0) +inf.0)
|
||||
(check-equal? (flulp-error -inf.0 -1.0) +inf.0)
|
||||
(check-equal? (flulp-error -inf.0 -min.0) +inf.0)
|
||||
(check-equal? (flulp-error -inf.0 -0.0) +inf.0)
|
||||
(check-equal? (flulp-error -inf.0 +0.0) +inf.0)
|
||||
(check-equal? (flulp-error -inf.0 +min.0) +inf.0)
|
||||
(check-equal? (flulp-error -inf.0 +1.0) +inf.0)
|
||||
(check-equal? (flulp-error -inf.0 +max.0) +inf.0)
|
||||
(check-equal? (flulp-error -inf.0 +inf.0) +inf.0)
|
||||
|
||||
;; RHS argument +inf.0
|
||||
(check-equal? (flulp-error -max.0 +inf.0) +inf.0)
|
||||
(check-equal? (flulp-error -1.0 +inf.0) +inf.0)
|
||||
(check-equal? (flulp-error -min.0 +inf.0) +inf.0)
|
||||
(check-equal? (flulp-error -0.0 +inf.0) +inf.0)
|
||||
(check-equal? (flulp-error +0.0 +inf.0) +inf.0)
|
||||
(check-equal? (flulp-error +min.0 +inf.0) +inf.0)
|
||||
(check-equal? (flulp-error +1.0 +inf.0) +inf.0)
|
||||
(check-equal? (flulp-error +max.0 +inf.0) +inf.0)
|
||||
|
||||
;; RHS argument -inf.0
|
||||
(check-equal? (flulp-error -max.0 -inf.0) +inf.0)
|
||||
(check-equal? (flulp-error -1.0 -inf.0) +inf.0)
|
||||
(check-equal? (flulp-error -min.0 -inf.0) +inf.0)
|
||||
(check-equal? (flulp-error -0.0 -inf.0) +inf.0)
|
||||
(check-equal? (flulp-error +0.0 -inf.0) +inf.0)
|
||||
(check-equal? (flulp-error +min.0 -inf.0) +inf.0)
|
||||
(check-equal? (flulp-error +1.0 -inf.0) +inf.0)
|
||||
(check-equal? (flulp-error +max.0 -inf.0) +inf.0)
|
||||
|
||||
;; RHS argument 0.0
|
||||
(check-equal? (flulp-error -max.0 +0.0) +inf.0)
|
||||
(check-equal? (flulp-error -1.0 +0.0) +inf.0)
|
||||
(check-equal? (flulp-error -min.0 +0.0) +inf.0)
|
||||
(check-equal? (flulp-error +min.0 +0.0) +inf.0)
|
||||
(check-equal? (flulp-error +1.0 +0.0) +inf.0)
|
||||
(check-equal? (flulp-error +max.0 +0.0) +inf.0)
|
||||
|
||||
;; RHS argument -0.0
|
||||
(check-equal? (flulp-error -max.0 -0.0) +inf.0)
|
||||
(check-equal? (flulp-error -1.0 -0.0) +inf.0)
|
||||
(check-equal? (flulp-error -min.0 -0.0) +inf.0)
|
||||
(check-equal? (flulp-error +min.0 -0.0) +inf.0)
|
||||
(check-equal? (flulp-error +1.0 -0.0) +inf.0)
|
||||
(check-equal? (flulp-error +max.0 -0.0) +inf.0)
|
||||
|
||||
;; Small errors
|
||||
(check-equal? (flulp-error +0.0 -min.0) 1.0)
|
||||
(check-equal? (flulp-error +0.0 +min.0) 1.0)
|
||||
(check-equal? (flulp-error -0.0 -min.0) 1.0)
|
||||
(check-equal? (flulp-error -0.0 +min.0) 1.0)
|
||||
(check-equal? (flulp-error -min.0 +min.0) 2.0)
|
||||
(check-equal? (flulp-error +min.0 -min.0) 2.0)
|
||||
|
||||
(define large-flulp-error-xys
|
||||
(list (list -1.0 -max.0)
|
||||
(list -1.0 +max.0)
|
||||
(list -min.0 -max.0)
|
||||
(list -min.0 +max.0)
|
||||
(list -1.0 +1.0)
|
||||
(list -max.0 +max.0)
|
||||
(list -max.0 -1.0)
|
||||
(list -max.0 -min.0)
|
||||
(list -max.0 +min.0)
|
||||
(list -max.0 +1.0)
|
||||
(list -1.0 -min.0)
|
||||
(list -1.0 +min.0)
|
||||
(list -min.0 -1.0)
|
||||
(list -min.0 +1.0)
|
||||
(list -0.0 -max.0)
|
||||
(list -0.0 -1.0)
|
||||
(list -0.0 +1.0)
|
||||
(list -0.0 +max.0)
|
||||
(list +0.0 -max.0)
|
||||
(list +0.0 -1.0)
|
||||
(list +0.0 +1.0)
|
||||
(list +0.0 +max.0)
|
||||
(list +min.0 -max.0)
|
||||
(list +min.0 -1.0)
|
||||
(list +min.0 +1.0)
|
||||
(list +min.0 +max.0)
|
||||
(list +1.0 -max.0)
|
||||
(list +1.0 -1.0)
|
||||
(list +1.0 -min.0)
|
||||
(list +1.0 +min.0)
|
||||
(list +1.0 +max.0)
|
||||
(list +max.0 -max.0)
|
||||
(list +max.0 -1.0)
|
||||
(list +max.0 -min.0)
|
||||
(list +max.0 +min.0)
|
||||
(list +max.0 +1.0)))
|
||||
|
||||
;; Large errors
|
||||
(for ([xy (in-list large-flulp-error-xys)])
|
||||
(match-define (list x y) xy)
|
||||
(check-true ((flulp-error x y) . >= . (expt 2 52))
|
||||
(format "x = ~a y = ~a" x y)))
|
||||
|
||||
(check-equal? (flulp-error 1.0 (flnext 1.0)) 1.0)
|
||||
(check-equal? (flulp-error +max.0 (flprev +max.0)) 1.0)
|
||||
|
||||
(for ([_ (in-range 1000)])
|
||||
(define s (random))
|
||||
(define e (fl (random-integer -1074 1024)))
|
||||
(define x (* s (flexpt 2.0 e)))
|
||||
(check-equal? (flulp-error x (flnext x)) 1.0
|
||||
(format "x = ~a" x)))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; FPU testing
|
||||
|
||||
(parameterize ([print-test-progress? #f])
|
||||
(test-fpu 1000))
|
||||
|
|
|
|||
|
|
@ -1,6 +1,8 @@
|
|||
#lang racket
|
||||
|
||||
(require "private/parameters.rkt")
|
||||
(require "private/parameters.rkt"
|
||||
"private/utils/flonum-tests.rkt")
|
||||
|
||||
(provide (all-from-out
|
||||
"private/parameters.rkt"))
|
||||
"private/parameters.rkt"
|
||||
"private/utils/flonum-tests.rkt"))
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user