diff --git a/collects/math/private/bigfloat/bigfloat-log-arithmetic.rkt b/collects/math/private/bigfloat/bigfloat-log-arithmetic.rkt index 6ed2a619cc..5479f6b671 100644 --- a/collects/math/private/bigfloat/bigfloat-log-arithmetic.rkt +++ b/collects/math/private/bigfloat/bigfloat-log-arithmetic.rkt @@ -2,7 +2,7 @@ (require "bigfloat-struct.rkt") -(provide bflog* bflog/ bflog+ bflog- bflog1-) +(provide bflog* bflog/ bflog+ bflog- bflog1- bflogb) (: bflog* (Bigfloat Bigfloat -> Bigfloat)) (define (bflog* log-x log-y) (bf+ log-x log-y)) @@ -26,3 +26,23 @@ (define (bflog1- log-x) (cond [(log-x . bf> . (bflog (bf 0.5))) (bflog (bf- (bfexpm1 log-x)))] [else (bflog1p (bf- (bfexp log-x)))])) + +(: bflogb (Bigfloat Bigfloat -> Bigfloat)) +(define (bflogb b x) + (cond [(bf= x 1.bf) 0.bf] + [(bf= b 1.bf) +nan.bf] + [(not (and (bf<= 0.bf b) (bf<= b +inf.bf) (bf<= 0.bf x) (bf<= x +inf.bf))) +nan.bf] + [(bf= b 0.bf) + (cond [(bf= x 0.bf) +inf.bf] + [(bf= x +inf.bf) -inf.bf] + [(bf<= x 1.bf) 0.bf] + [else -0.bf])] + [(bf= b +inf.bf) + (cond [(bf= x 0.bf) -inf.bf] + [(bf= x +inf.bf) +inf.bf] + [(bf<= 1.bf x) 0.bf] + [else -0.bf])] + [(bf= x 0.bf) (if (bf< b 1.bf) +inf.bf -inf.bf)] + [(bf= x +inf.bf) (if (bf< b 1.bf) -inf.bf +inf.bf)] + [else + (bf/ (bflog x) (bflog b))])) diff --git a/collects/math/private/flonum/flonum-exp.rkt b/collects/math/private/flonum/flonum-exp.rkt index 8d9cb03a05..48976bd3bb 100644 --- a/collects/math/private/flonum/flonum-exp.rkt +++ b/collects/math/private/flonum/flonum-exp.rkt @@ -3,9 +3,10 @@ (require racket/performance-hint "flonum-functions.rkt" "flonum-constants.rkt" - "flonum-polyfun.rkt") + "flonum-polyfun.rkt" + "../unsafe.rkt") -(provide flexpm1 flexpsqr flgauss flexp1p) +(provide flexpm1 flexpsqr flgauss flexp1p flexp2) (define expm1-poly-numer (make-flpolyfun @@ -33,6 +34,10 @@ [else (fl+ (fl* x 0.10281276702880859e1) (fl* x (fl/ (expm1-poly-numer x) (expm1-poly-denom x))))])) +;; Integer arguments for flexp2 +(: flexp2s (Vectorof Nonnegative-Flonum)) +(define flexp2s (build-vector (- 1024 -1074) (λ: ([n : Index]) (fl (expt 2 (- n 1074)))))) + (begin-encourage-inline (: flexpm1 (Float -> Float)) @@ -83,4 +88,11 @@ (cond [(fl= lg2x lg2x+1) (flexp (fl+ 1.0 x))] [else (fl* (flexp x) (flexp 1.0))])])) + (: flexp2 (Flonum -> Nonnegative-Flonum)) + (define (flexp2 x) + (cond [(fl<= x -1075.0) 0.0] + [(fl>= x 1024.0) +inf.0] + [(fl= x (flround x)) (unsafe-vector-ref flexp2s (fl->exact-integer (fl+ x 1074.0)))] + [else (flexpt 2.0 x)])) + ) ; begin-encourage-inline diff --git a/collects/math/private/flonum/flonum-log.rkt b/collects/math/private/flonum/flonum-log.rkt index a3290b804d..aec8ec8053 100644 --- a/collects/math/private/flonum/flonum-log.rkt +++ b/collects/math/private/flonum/flonum-log.rkt @@ -10,7 +10,9 @@ (provide fllog1p fllog+ lg1+ lg+ lg1- lg- lgsum - fllog-quotient) + fllog-quotient + fllog2 + fllogb) (begin-encourage-inline @@ -105,3 +107,89 @@ [else (fl+ (fllog x) (- (fllog y)))])] [(s . fl= . 0.0) -inf.0] [else +nan.0]))) + +(define log-max.0 (fllog +max.0)) +(define log2.0 (fllog 2.0)) + +(: fllog2* (Flonum -> Flonum)) +;; Computes log2(x) with a least 8 extra bits precision, which reduces the probability of rounding +;; error significantly. Assumes 0.0 < x < +inf.0 and x != 1.0. +(define (fllog2* x) + (let* ([log-x (fllog x)] + ;; Solve for x^(2^k) = +max.0 (k is basically the number of extra bits precision) + [k (fl/ (fllog (fl/ log-max.0 (flabs log-x))) log2.0)] + ;; We'll be operating on x^adj, which is huge + [adj (flexp2 (flceiling (- k 1.0)))] + [adj (if (fl>= x 1.0) adj (- adj))] + ;; Compute floor(log2(x^adj)) + [y2 (fltruncate (fl/ (fl* adj log-x) log2.0))] + ;; Compute "remainder" log2(x^adj/2^y2) (note: dividing by 2^y2 is exact) + [y1 (fl/ (fllog (fl/ (flexpt x adj) (flexp2 y2))) log2.0)]) + (fl+ (fl/ y2 adj) (fl/ y1 adj)))) + +(: fllog2 (Flonum -> Flonum)) +;; Largest observed error is 0.5006 ulps +(define (fllog2 x) + (cond [(fl<= x 0.0) (if (fl< x 0.0) +nan.0 -inf.0)] + [(fl< x +inf.0) (if (fl= x 1.0) 0.0 (fllog2* x))] + [(fl= x +inf.0) +inf.0] + [else +nan.0])) + +(: fllogb (Flonum Flonum -> Flonum)) +;; Largest observed error is 2.1 ulps, but is usually < 0.7 ulps +(define (fllogb b x) + (cond [(fl= x 1.0) 0.0] + [(fl= b 1.0) + ;; For x != 1, first limit wrt x: +inf.0 or -inf.0 + +nan.0] + [(fl= b 2.0) + ;; Using the more accurate `fllog2' ensures that exact cases have zero error + (fllog2 x)] + [(not (and (fl<= 0.0 b) (fl<= b +inf.0) (fl<= 0.0 x) (fl<= x +inf.0))) + ;; One or both is out of bounds or is +nan.0 + +nan.0] + [(fl= b 0.0) + (cond [(fl= x 0.0) + ;; First limit wrt x: +inf.0 + ;; First limit wrt b: 0.0 + ;; +inf.0 corrects left-inverse case (fllogb 0.0 (flexpt 0.0 +inf.0)) + ;; +inf.0 corrects right-inverse case (flexpt 0.0 (fllogb 0.0 0.0)) + +inf.0] + [(fl= x +inf.0) + ;; First limit wrt x: -inf.0 + ;; First limit wrt b: 0.0 + ;; -inf.0 corrects left-inverse case (fllogb 0.0 (flexpt 0.0 -inf.0)) + ;; -inf.0 corrects right-inverse case (flexpt 0.0 (fllogb 0.0 +inf.0)) + -inf.0] + [(fl<= x 1.0) 0.0] + [else -0.0])] + [(fl= b +inf.0) + (cond [(fl= x 0.0) + ;; First limit wrt x: -inf.0 + ;; First limit wrt b: -0.0 + ;; -inf.0 corrects left-inverse case (fllogb +inf.0 (flexpt +inf.0 -inf.0)) + ;; -inf.0 corrects right-inverse case (flexpt +inf.0 (fllogb +inf.0 0.0)) + -inf.0] + [(fl= x +inf.0) + ;; First limit wrt x: +inf.0 + ;; First limit wrt b: 0.0 + ;; +inf.0 corrects left-inverse case (fllogb +inf.0 (flexpt +inf.0 +inf.0)) + ;; +inf.0 corrects right-inverse case (flexpt +inf.0 (fllogb +inf.0 +inf.0)) + +inf.0] + [(fl<= 1.0 x) 0.0] + [else -0.0])] + [(fl= x 0.0) (if (fl< b 1.0) +inf.0 -inf.0)] + [(fl= x +inf.0) (if (fl< b 1.0) -inf.0 +inf.0)] + [else + (define log-b (fllog b)) + (define y (fl/ (fllog x) log-b)) + ;; One Newton iteration reduces error to <= 1 ulp (instead of <= 2 ulps) + (define numer (fl- x (flexpt b y))) + (define denom (fl* x log-b)) + (cond [(and (numer . fl> . -inf.0) (numer . fl< . +inf.0) + (denom . fl> . 0.0) (denom . fl< . +inf.0)) + (fl+ y (fl/ numer denom))] + [else + ;; Oh noes! We had overflows or underflows! + ;; Not a lot we can do without introducing more error, so just return y + y])])) diff --git a/collects/math/private/utils/flonum-tests.rkt b/collects/math/private/utils/flonum-tests.rkt index 40de732c83..127aebc014 100644 --- a/collects/math/private/utils/flonum-tests.rkt +++ b/collects/math/private/utils/flonum-tests.rkt @@ -1,9 +1,7 @@ #lang typed/racket/base (require racket/math - racket/flonum racket/list - typed/rackunit "../base/base-random.rkt" "../flonum/expansion/expansion-base.rkt" "../flonum/expansion/expansion-exp.rkt" @@ -12,21 +10,64 @@ "../flonum/flonum-constants.rkt" "../flonum/flonum-bits.rkt" "../flonum/flonum-error.rkt" + "../flonum/flonum-log.rkt" + "../bigfloat/bigfloat-log-arithmetic.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) +(provide + print-fp-test-progress? + ;; Unary flonum tests + test-flabs + test-flsqrt + test-fllog + test-flexp + test-flsin + test-flcos + test-fltan + test-flasin + test-flacos + test-flatan + test-fllog2 + ;; Binary flonum tests + test-fl+ + test-fl* + test-fl- + test-fl/ + test-flexpt + test-fllogb + ;; Unary flop/error tests + test-flsqr/error + test-flsqrt/error + test-flexp/error + ;; Binary flop/error tests + test-fl+/error + test-fl-/error + test-fl*/error + test-fl//error + ;; fl2 conversion test + test-fl2 + ;; Unary fl2 tests + test-fl2abs + test-fl2sqr + test-fl2sqrt + test-fl2exp + test-fl2expm1 + test-fl2log + test-fl2log1p + ;; Binary fl2 tests + test-fl2+ + test-fl2- + test-fl2* + test-fl2/ + ;; Comprehensive test + test-floating-point) ;; Allowable error for different kinds of functions, in ulps (define flonum-fun-ulps 0.5) +(define fllog2-ulps 1.0) +(define fllogb-ulps 2.5) (define flonum/error-fun-ulps 0.5) (define flexp/error-fun-ulps 3.0) (define fl2-conversion-ulps 0.5) @@ -70,15 +111,15 @@ (or (not (flonum? e)) (e . fl> . ulps))) xes)) -(: print-test-progress? (Parameterof Boolean)) -(define print-test-progress? (make-parameter #t)) +(: print-fp-test-progress? (Parameterof Boolean)) +(define print-fp-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)) + (when (and (print-fp-test-progress?) (i . > . 0) (i . <= . m)) (let* ([flush? (cond [(= i 1) (printf "~a: " name)] [else #f])] [flush? (cond [(= 0 (modulo i progress-chunk-size)) @@ -170,48 +211,52 @@ ;; =================================================================================================== ;; Flonum functions -(: flonum-error (Flonum Bigfloat -> Any)) +(define-type Flonum-Error (U Flonum (List Symbol Flonum Real))) +(define-type Unary-Flonum-Failure (List (List Symbol Flonum) Flonum-Error)) +(define-type Binary-Flonum-Failure (List (List Symbol Flonum Flonum) Flonum-Error)) + +(: flonum-error (Flonum Bigfloat -> Flonum-Error)) (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)) +(: unary-flonum-fun-error ((Flonum -> Flonum) (Bigfloat -> Bigfloat) Flonum -> Flonum-Error)) (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)))) + -> (Listof Unary-Flonum-Failure))) (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)]) + (for/list: : (Listof Unary-Flonum-Failure) ([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))) (current-max-ulp-error))) (: binary-flonum-fun-error - ((Flonum Flonum -> Flonum) (Bigfloat Bigfloat -> Bigfloat) Flonum Flonum -> Any)) + ((Flonum Flonum -> Flonum) (Bigfloat Bigfloat -> Bigfloat) Flonum Flonum -> Flonum-Error)) (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)))) + -> (Listof Binary-Flonum-Failure))) (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)]) + (for/list: : (Listof Binary-Flonum-Failure) ([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))) (current-max-ulp-error))) @@ -219,24 +264,27 @@ ;; =================================================================================================== ;; fl2 conversion -(: fl2-error (Flonum Flonum Real -> Any)) +(define-type Fl2-Error (U Flonum (List Symbol Flonum Real))) +(define-type Fl2-Failure (List (List 'fl2 Real) Fl2-Error)) + +(: fl2-error (Flonum Flonum Real -> Fl2-Error)) (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)) +(: fl2-conversion-error (Real -> Fl2-Error)) (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)))) +(: test-fl2-conversion (Integer -> (Listof Fl2-Failure))) (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)]) + (for/list: : (Listof Fl2-Failure) ([x (in-list xs)] + [i (in-naturals 1)]) (maybe-print-progress 'fl2 i m) (list (list 'fl2 x) (fl2-conversion-error x))) (current-max-ulp-error))) @@ -244,8 +292,11 @@ ;; =================================================================================================== ;; Flonum arithmetic with error +(define-type Unary-Fl/Error-Failure (List (List Symbol Flonum) Fl2-Error)) +(define-type Binary-Fl/Error-Failure (List (List Symbol Flonum Flonum) Fl2-Error)) + (: unary-flonum/error-fun-error ((Flonum -> (Values Flonum Flonum)) (Bigfloat -> Bigfloat) Flonum - -> Any)) + -> Fl2-Error)) (define (unary-flonum/error-fun-error f g x) (define-values (z2 z1) (f x)) (fl2-error z2 z1 (parameterize ([bf-precision 256]) @@ -254,7 +305,7 @@ (: binary-flonum/error-fun-error ((Flonum Flonum -> (Values Flonum Flonum)) (Bigfloat Bigfloat -> Bigfloat) Flonum Flonum - -> Any)) + -> Fl2-Error)) (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]) @@ -262,29 +313,29 @@ (: test-unary-flonum/error-fun (Symbol (Flonum -> (Values Flonum Flonum)) (Bigfloat -> Bigfloat) Integer - -> (Listof (List (List Symbol Flonum) Any)))) + -> (Listof Unary-Fl/Error-Failure))) (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)]) + (for/list: : (Listof Unary-Fl/Error-Failure) ([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))) (current-max-ulp-error))) (: test-binary-flonum/error-fun (Symbol (Flonum Flonum -> (Values Flonum Flonum)) (Bigfloat Bigfloat -> Bigfloat) Integer - -> (Listof (List (List Symbol Flonum Flonum) Any)))) + -> (Listof Binary-Fl/Error-Failure))) (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)]) + (for/list: : (Listof Binary-Fl/Error-Failure) ([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))) (current-max-ulp-error))) @@ -292,8 +343,11 @@ ;; =================================================================================================== ;; Flonum expansions +(define-type Unary-Fl2-Failure (List (List Symbol Flonum Flonum) Fl2-Error)) +(define-type Binary-Fl2-Failure (List (List Symbol Flonum Flonum Flonum Flonum) Fl2-Error)) + (: unary-fl2-fun-error ((Flonum Flonum -> (Values Flonum Flonum)) (Bigfloat -> Bigfloat) - Flonum Flonum -> Any)) + Flonum Flonum -> Fl2-Error)) (define (unary-fl2-fun-error f g x2 x1) (define-values (z2 z1) (f x2 x1)) (fl2-error z2 z1 (parameterize ([bf-precision 256]) @@ -301,13 +355,13 @@ (: test-unary-fl2-fun (Symbol (Flonum Flonum -> (Values Flonum Flonum)) (Bigfloat -> Bigfloat) Integer - -> (Listof (List (List Symbol Flonum Flonum) Any)))) + -> (Listof Unary-Fl2-Failure))) (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)]) + (for/list: : (Listof Unary-Fl2-Failure) ([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))) @@ -316,7 +370,7 @@ (: binary-fl2-fun-error ((Flonum Flonum Flonum Flonum -> (Values Flonum Flonum)) (Bigfloat Bigfloat -> Bigfloat) Flonum Flonum Flonum Flonum - -> Any)) + -> Fl2-Error)) (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]) @@ -324,17 +378,16 @@ (: 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)))) + Integer -> (Listof Binary-Fl2-Failure))) (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)]) + (for/list: : (Listof Binary-Fl2-Failure) ([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)) @@ -342,151 +395,141 @@ (current-max-ulp-error))) ;; =================================================================================================== +;; Flonum tests -(: test-fpu-arith (Natural -> Any)) -(define (test-fpu-arith n) - (parameterize ([current-max-ulp-error flonum-fun-ulps]) - (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) - '()))) +(define-syntax-rule (define-unary-flonum-test test-name flop bfop mn mx ulps) + (begin (: test-name (Natural -> (Listof Unary-Flonum-Failure))) + (define (test-name n) + (parameterize ([current-max-ulp-error ulps]) + (test-unary-flonum-fun 'flop flop bfop n mn mx))))) -(: test-fpu-trig (Natural -> Any)) -(define (test-fpu-trig n) - (parameterize ([current-max-ulp-error flonum-fun-ulps]) - (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) - '()))) +(define-syntax-rule (define-binary-flonum-test test-name flop bfop ulps) + (begin (: test-name (Natural -> (Listof Binary-Flonum-Failure))) + (define (test-name n) + (parameterize ([current-max-ulp-error ulps]) + (test-binary-flonum-fun 'flop flop bfop n))))) -(: test-fpu-non-trig (Natural -> Any)) -(define (test-fpu-non-trig n) - (parameterize ([current-max-ulp-error flonum-fun-ulps]) - (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) - '()))) +(define-unary-flonum-test test-flabs flabs bfabs -inf.0 +inf.0 flonum-fun-ulps) +(define-unary-flonum-test test-flsqrt flsqrt bfsqrt 0.0 +inf.0 flonum-fun-ulps) +(define-unary-flonum-test test-fllog fllog bflog 0.0 +inf.0 flonum-fun-ulps) +(define-unary-flonum-test test-flexp flexp bfexp -746.0 710.0 flonum-fun-ulps) +(define-unary-flonum-test test-flsin flsin bfsin -inf.0 +inf.0 flonum-fun-ulps) +(define-unary-flonum-test test-flcos flcos bfcos -inf.0 +inf.0 flonum-fun-ulps) +(define-unary-flonum-test test-fltan fltan bftan -inf.0 +inf.0 flonum-fun-ulps) +(define-unary-flonum-test test-flasin flasin bfasin -1.0 1.0 flonum-fun-ulps) +(define-unary-flonum-test test-flacos flacos bfacos -1.0 1.0 flonum-fun-ulps) +(define-unary-flonum-test test-flatan flatan bfatan -inf.0 +inf.0 flonum-fun-ulps) +(define-unary-flonum-test test-fllog2 fllog2 bflog2 0.0 +inf.0 fllog2-ulps) +(define-binary-flonum-test test-fl+ fl+ bf+ flonum-fun-ulps) +(define-binary-flonum-test test-fl- fl- bf- flonum-fun-ulps) +(define-binary-flonum-test test-fl* fl* bf* flonum-fun-ulps) +(define-binary-flonum-test test-fl/ fl/ bf/ flonum-fun-ulps) +(define-binary-flonum-test test-flexpt flexpt bfexpt flonum-fun-ulps) +(define-binary-flonum-test test-fllogb fllogb bflogb fllogb-ulps) -(: test-fpu-arith/error (Natural -> Any)) -(define (test-fpu-arith/error n) - (parameterize ([current-max-ulp-error flonum/error-fun-ulps]) - (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) - '()))) +;; =================================================================================================== +;; Flonum/error tests -(: test-fpu-non-trig/error (Natural -> Any)) -(define (test-fpu-non-trig/error n) - (parameterize ([current-max-ulp-error flexp/error-fun-ulps]) - (check-equal? (test-unary-flonum/error-fun 'flexp/error flexp/error bfexp n) - '()))) +(define-syntax-rule (define-unary-flop/error-test test-name flop bfop ulps) + (begin (: test-name (Natural -> (Listof Unary-Fl/Error-Failure))) + (define (test-name n) + (parameterize ([current-max-ulp-error ulps]) + (test-unary-flonum/error-fun 'flop flop bfop n))))) -(: test-fpu-arith/fl2 (Natural -> Any)) -(define (test-fpu-arith/fl2 n) +(define-syntax-rule (define-binary-flop/error-test test-name flop bfop ulps) + (begin (: test-name (Natural -> (Listof Binary-Fl/Error-Failure))) + (define (test-name n) + (parameterize ([current-max-ulp-error ulps]) + (test-binary-flonum/error-fun 'flop flop bfop n))))) + +(define-unary-flop/error-test test-flsqr/error flsqr/error bfsqr flonum/error-fun-ulps) +(define-unary-flop/error-test test-flsqrt/error flsqrt/error bfsqrt flonum/error-fun-ulps) +(define-unary-flop/error-test test-flexp/error flexp/error bfexp flexp/error-fun-ulps) +(define-binary-flop/error-test test-fl+/error fl+/error bf+ flonum/error-fun-ulps) +(define-binary-flop/error-test test-fl-/error fl-/error bf- flonum/error-fun-ulps) +(define-binary-flop/error-test test-fl*/error fl*/error bf* flonum/error-fun-ulps) +(define-binary-flop/error-test test-fl//error fl//error bf/ flonum/error-fun-ulps) + +;; =================================================================================================== +;; fl2 tests + +(: test-fl2 (Natural -> (Listof Fl2-Failure))) +(define (test-fl2 n) (parameterize ([current-max-ulp-error fl2-conversion-ulps]) - (check-equal? (test-fl2-conversion n) - '())) - (parameterize ([current-max-ulp-error unary-fl2-fun-ulps]) - (check-equal? (test-unary-fl2-fun 'fl2abs fl2abs bfabs n) - '()) - (check-equal? (test-unary-fl2-fun 'fl2sqr fl2sqr bfsqr n) - '())) - (parameterize ([current-max-ulp-error binary-fl2-fun-ulps]) - (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-binary-fl2-fun 'fl2/ fl2/ bf/ n) - '()))) + (test-fl2-conversion n))) -(: test-fpu-non-trig/fl2 (Natural -> Any)) -(define (test-fpu-non-trig/fl2 n) - (parameterize ([current-max-ulp-error unary-fl2-fun-ulps]) - (check-equal? (test-unary-fl2-fun 'fl2sqrt fl2sqrt bfsqrt n) - '())) - (parameterize ([current-max-ulp-error fl2exp-fun-ulps]) - (check-equal? (test-unary-fl2-fun 'fl2exp fl2exp bfexp n) - '())) - (parameterize ([current-max-ulp-error fl2exp-fun-ulps]) - (check-equal? (test-unary-fl2-fun 'fl2expm1 fl2expm1 bfexpm1 n) - '())) - (parameterize ([current-max-ulp-error fl2log-fun-ulps]) - (check-equal? (test-unary-fl2-fun 'fl2log fl2log bflog n) - '())) - (parameterize ([current-max-ulp-error fl2log-fun-ulps]) - (check-equal? (test-unary-fl2-fun 'fl2log1p fl2log1p bflog1p n) - '())) - ) +(define-syntax-rule (define-unary-fl2-test test-name fl2op bfop ulps) + (begin (: test-name (Natural -> (Listof Unary-Fl2-Failure))) + (define (test-name n) + (parameterize ([current-max-ulp-error ulps]) + (test-unary-fl2-fun 'fl2op fl2op bfop 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-non-trig/error n) - (test-fpu-arith/fl2 n) - (test-fpu-non-trig/fl2 n)) +(define-syntax-rule (define-binary-fl2-test test-name fl2op bfop ulps) + (begin (: test-name (Natural -> (Listof Binary-Fl2-Failure))) + (define (test-name n) + (parameterize ([current-max-ulp-error ulps]) + (test-binary-fl2-fun 'fl2op fl2op bfop n))))) -(for*: ([x2 (list -inf.0 -max.0 -1.0 -min.0 -0.0 0.0 +min.0 +1.0 +max.0 +inf.0 +nan.0)] - [x1 (list -inf.0 -max.0 -1.0 -min.0 -0.0 0.0 +min.0 +1.0 +max.0 +inf.0 +nan.0)]) - (define n - (count (λ: ([b : Boolean]) b) - (map (λ: ([f : (Flonum Flonum -> Boolean)]) - (f x2 x1)) - (list fl2rational? fl2infinite? fl2nan?)))) - (unless (= n 1) (printf "x2 = ~v x1 = ~v~n" x2 x1))) +(define-unary-fl2-test test-fl2abs fl2abs bfabs unary-fl2-fun-ulps) +(define-unary-fl2-test test-fl2sqr fl2sqr bfsqr unary-fl2-fun-ulps) +(define-unary-fl2-test test-fl2sqrt fl2sqrt bfsqrt unary-fl2-fun-ulps) +(define-unary-fl2-test test-fl2exp fl2exp bfexp fl2exp-fun-ulps) +(define-unary-fl2-test test-fl2expm1 fl2expm1 bfexpm1 fl2exp-fun-ulps) +(define-unary-fl2-test test-fl2log fl2log bflog fl2log-fun-ulps) +(define-unary-fl2-test test-fl2log1p fl2log1p bflog1p fl2log-fun-ulps) +(define-binary-fl2-test test-fl2+ fl2+ bf+ binary-fl2-fun-ulps) +(define-binary-fl2-test test-fl2- fl2- bf- binary-fl2-fun-ulps) +(define-binary-fl2-test test-fl2* fl2* bf* binary-fl2-fun-ulps) +(define-binary-fl2-test test-fl2/ fl2/ bf/ binary-fl2-fun-ulps) -#| -Tests to add +;; =================================================================================================== +;; Comprehensive test -(for*: ([x2 (list -inf.0 -max.0 -1.0 -min.0 -0.0 0.0 +min.0 +1.0 +max.0 +inf.0 +nan.0)] - [x1 (list -inf.0 -max.0 -1.0 -min.0 -0.0 0.0 +min.0 +1.0 +max.0 +inf.0 +nan.0)]) - (define n - (count (λ: ([b : Boolean]) b) - (map (λ: ([f : (Flonum Flonum -> Boolean)]) - (f x2 x1)) - (list fl2rational? fl2infinite? fl2nan?)))) - (unless (= n 1) (printf "x2 = ~v x1 = ~v~n" x2 x1))) - -fl2= -fl2> -fl2< -fl2>= -fl2<= - -(fl2step x2 x1 n/2) twice = (fl2step x2 x1 n) - -|# - -(check-true (let-values ([(y2 y1) (fl+/error +max.hi +max.lo)]) - (fl2= y2 y1 +max.hi +max.lo))) - -(check-true (let*-values ([(y2 y1) (fl2next +max.hi +max.lo)]) - (fl2infinite? y2 y1))) +(: test-floating-point (Natural -> (Listof (U Unary-Flonum-Failure + Binary-Flonum-Failure + Unary-Fl/Error-Failure + Binary-Fl/Error-Failure + Fl2-Failure + Unary-Fl2-Failure + Binary-Fl2-Failure)))) +(define (test-floating-point n) + (append + ;; Hardware implementation tests + (test-flabs n) + (test-fl+ n) + (test-fl* n) + (test-fl- n) + (test-fl/ n) + (test-flsqrt n) + (test-fllog n) + (test-flexp n) + (test-flexpt n) + (test-flsin n) + (test-flcos n) + (test-fltan n) + (test-flasin n) + (test-flacos n) + (test-flatan n) + (test-fllog2 n) + (test-fllogb n) + ;; Derived tests + (test-fl+/error n) + (test-fl-/error n) + (test-fl*/error n) + (test-flsqr/error n) + (test-fl//error n) + (test-fl2 n) + (test-fl2abs n) + (test-fl2+ n) + (test-fl2- n) + (test-fl2* n) + (test-fl2sqr n) + (test-fl2/ n) + (test-flsqrt/error n) + (test-fl2sqrt n) + (test-flexp/error n) + (test-fl2exp n) + (test-fl2expm1 n) + (test-fl2log n) + (test-fl2log1p n) + )) diff --git a/collects/math/scribblings/math-flonum.scrbl b/collects/math/scribblings/math-flonum.scrbl index a7ef6c80fc..aa79fd5b31 100644 --- a/collects/math/scribblings/math-flonum.scrbl +++ b/collects/math/scribblings/math-flonum.scrbl @@ -1,12 +1,15 @@ #lang scribble/manual @(require scribble/eval + scribble/core + scribble/html-properties racket/sandbox (for-label racket/base racket/vector racket/list math plot (only-in typed/racket/base -> - Flonum Integer Index Real Boolean Any Listof Vectorof FlVector)) + Flonum Integer Index Real Boolean Any Listof Vectorof FlVector + Nonnegative-Flonum)) "utils.rkt") @(define untyped-eval (make-untyped-math-eval)) @@ -164,6 +167,82 @@ But see @racket[flexpt1p], which is more accurate still. Like @racket[(flexpt (+ 1.0 x) y)], but accurate for any @racket[x] and @racket[y]. } +@defproc[(flexp2 [x Flonum]) Nonnegative-Flonum]{ +Equivalent to @racket[(flexpt 2.0 x)], but faster when @racket[x] is an integer. +} + +@defproc[(fllog2 [x Flonum]) Flonum]{ +Computes the base-2 log of @racket[x] more accurately than @racket[(/ (fllog x) (fllog 2.0))]. +In particular, @racket[(fllog2 x)] is correct for any power of two @racket[x]. +@examples[#:eval untyped-eval + (fllog2 4.5) + (/ (fllog (flexp2 -1066.0)) (fllog 2.0)) + (fllog2 (flexp2 -1066.0))] +Maximum observed error is 0.5006 @tech{ulps}, but is almost always no more than 0.5 (i.e. it is +almost always @italic{correct}). +} + +@defproc[(fllogb [b Flonum] [x Flonum]) Flonum]{ +Computes the base-@racket[b] log of @racket[x] more accurately than @racket[(/ (fllog x) (fllog b))], +and handles limit values correctly. +@examples[#:eval untyped-eval + (plot3d (contour-intervals3d (λ (b x) (fllogb (fl b) (fl x))) 0 4 0 4) + #:x-label "b" #:y-label "x")] + +Maximum observed error is 2.1 @tech{ulps}, but is usually less than 0.7 (i.e. near rounding error). + +Except possibly at limit values (such as @racket[0.0] and @racket[+inf.0], and @racket[b = 1.0]) +and except when the inner expression underflows or overflows, @racket[fllogb] approximately meets +these identities for @racket[b > 0.0]: +@itemlist[@item{Left inverse: @racket[(fllogb b (flexpt b y)) = y]} + @item{Right inverse: @racket[(flexpt b (fllogb b x)) = x] when @racket[x > 0.0]}] + +Unlike with @racket[flexpt], there is no standard for @racket[fllogb]'s behavior at limit values. +Fortunately, deriving the following rules (applied in order) is not prohibitively difficult. +@centered{ +@tabular[#:style + (style 'plain + (list (table-columns (list (style 'plain (list 'left)) + (style 'plain (list 'center)) + (style 'plain (list 'right)))) + (attributes '((width . "90%"))))) + (list (list @bold{Case} @bold{Condition} @bold{Value}) + (list @racket[(fllogb b 1.0)] "" @racket[0.0]) + (list @racket[(fllogb 1.0 x)] "" @racket[+nan.0]) + (list @racket[(fllogb b x)] + @nested{@racket[b < 0.0] or @racket[x < 0.0]} + @racket[+nan.0]) + (list @italic{Double limits} 'cont 'cont) + (list @racket[(fllogb 0.0 0.0)] "" @racket[+inf.0]) + (list @racket[(fllogb 0.0 +inf.0)] "" @racket[-inf.0]) + (list @racket[(fllogb +inf.0 0.0)] "" @racket[-inf.0]) + (list @racket[(fllogb +inf.0 +inf.0)] "" @racket[+inf.0]) + (list @italic{Limits with respect to @racket[b]} 'cont 'cont) + (list @racket[(fllogb 0.0 x)] @nested{@racket[x < 1.0]} @racket[0.0]) + (list @racket[(fllogb 0.0 x)] @nested{@racket[x > 1.0]} @racket[-0.0]) + (list @racket[(fllogb +inf.0 x)] @nested{@racket[x > 1.0]} @racket[0.0]) + (list @racket[(fllogb +inf.0 x)] @nested{@racket[x < 1.0]} @racket[-0.0]) + (list @italic{Limits with respect to @racket[x]} 'cont 'cont) + (list @racket[(fllogb b 0.0)] @nested{@racket[b < 1.0]} @racket[+inf.0]) + (list @racket[(fllogb b 0.0)] @nested{@racket[b > 1.0]} @racket[-inf.0]) + (list @racket[(fllogb b +inf.0)] @nested{@racket[b > 1.0]} @racket[+inf.0]) + (list @racket[(fllogb b +inf.0)] @nested{@racket[b < 1.0]} @racket[-inf.0]))]} +Most of these rules are derived by taking limits of the mathematical base-@racket[b] log function. +Except for @racket[(fllogb 1.0 x)], when doing so gives rise to ambiguities, they are resolved using +@racket[flexpt]'s behavior, which follows the IEEE 754 and C99 standards for @tt{pow}. + +For example, consider @racket[(fllogb 0.0 0.0)]. +Taking an interated limit, we get ∞ if the outer limit is with respect to @racket[x], or 0 if the +outer limit is with respect to @racket[b]. +This would normally mean @racket[(fllogb 0.0 0.0) = +nan.0]. + +However, choosing @racket[+inf.0] ensures that these additional left-inverse and right-inverse +identities hold: +@racketblock[(fllogb 0.0 (flexpt 0.0 +inf.0)) = +inf.0 + (flexpt 0.0 (fllogb 0.0 0.0)) = 0.0] +Further, choosing @racket[0.0] does not ensure that any additional identities hold. +} + @defproc[(make-flexpt [x Real]) (Flonum -> Flonum)]{ Equivalent to @racket[(λ (y) (flexpt x y))] when @racket[x] is a flonum, but much more accurate for large @racket[y] when @racket[x] cannot be exactly represented diff --git a/collects/math/scribblings/math-utils.scrbl b/collects/math/scribblings/math-utils.scrbl index 50569399c4..f911cbe43d 100644 --- a/collects/math/scribblings/math-utils.scrbl +++ b/collects/math/scribblings/math-utils.scrbl @@ -16,11 +16,15 @@ @defmodule[math/utils] +@section[#:tag "utils:parallel"]{Parallelization} + @defparam[max-math-threads num Positive-Integer]{ The maximum number of threads a parallelized @racketmodname[math] function will use. The default value is @racket[(max 1 (processor-count))]. } +@section[#:tag "utils:dft"]{Discrete Fourier Transform Conventions} + @defparam[dft-convention lst (List Real Real)]{ A parameter controlling the convention used for scaling discrete Fourier transforms, such as those performed by @racket[array-fft]. The default value is @racket['(1 -1)], which represents the convention @@ -42,4 +46,54 @@ on @tt{Fourier}}, from which this excellent idea was stolen. Returns the convention used for inverse Fourier transforms, given the current convention. } +@section[#:tag "utils:fpu-test"]{Floating-Point Compliance Testing} + +@defproc[(test-floating-point [n Natural]) (Listof (List Any Any))]{ +Runs a comprehensive test of the system's IEEE 754 (floating-point) compliance, and reports +unexpected inaccuracies and errors. + +In each test, a function is applied to some carefully chosen values, as well as @racket[n] additional +random values. +Its corresponding @tech{bigfloat} function is applied to the same values, and the answers are +compared. +Each test returns a list of failures, which are appended and returned. + +Each failure in a failure list is formatted +@racketblock[(list (list name args ...) reason)] +where @racket[name] is the name of a function, such as @racket['fl+], @racket[args ...] are the +arguments it was applied to, and @racket[reason] is the reason for the failure. + +If @racket[reason] is a flonum, the failure was due to inaccuracy. For example, +@racketblock[(list (list 'fl+ 4.5 2.3) 0.76)] +means the result of @racket[(fl+ 4.5 2.3)] was off by @racket[0.76] @tech{ulps}. + +The threshold for reporting unexpected inaccuracy depends on the function tested. +All the arithmetic and irrational functions exported by @racketmodname[racket/flonum], for example, +must have no more than @racket[0.5] ulps error in order to be compliant. + +Two other possible failure reasons are +@racketblock[(list 'different-zero 0.0 -0.0) + (list 'different-zero -0.0 0.0)] +The first zero is the answer returned by the function, and the second zero is the expected answer. + +Other possible failure reasons have the form +@racketblock[(list 'not-fl2? x y)] +meaning that the result @racket[(values x y)] is not a valid flonum expansion. +Such reasons are only given for failures of functions whose names begin with @tt{fl2} or contain +@tt{/error}. +These functions are currently undocumented, but are used to implement many +@racketmodname[math/flonum], @racketmodname[math/special-functions], and +@racketmodname[math/distributions] functions. + +Tests of functions that operate on and return flonum expansions are the strictest tests, requiring +hardware arithmetic to be perfectly IEEE 754 compliant. +They reliably fail on seemingly innocuous noncompliant behavior, such as computing intermediate +results with 80-bit precision. +} + +@defparam[print-fp-test-progress? print? Boolean]{ +When @racket[(print-fp-test-progress?)] is @racket[#t], floating-point tests print and flush a +representation of their progress as they run. The default value is @racket[#t]. +} + @(close-eval untyped-eval) diff --git a/collects/math/tests/flonum-tests.rkt b/collects/math/tests/flonum-tests.rkt index d3ed5bf3a6..d9b9745b91 100644 --- a/collects/math/tests/flonum-tests.rkt +++ b/collects/math/tests/flonum-tests.rkt @@ -144,8 +144,56 @@ (check-equal? (flulp-error x (flnext x)) 1.0 (format "x = ~a" x))) +;; =================================================================================================== +;; fllog2 + +;; Make sure it's exact for exact powers of two +(for ([x (in-range -1074.0 1024.0)]) + (define y (flexp2 x)) + (define x0 (fllog2 y)) + (check-equal? x0 x)) + +;; =================================================================================================== +;; fl2 tests + +(for* ([x2 (list -inf.0 -max.0 -1.0 -min.0 -0.0 0.0 +min.0 +1.0 +max.0 +inf.0 +nan.0)] + [x1 (list -inf.0 -max.0 -1.0 -min.0 -0.0 0.0 +min.0 +1.0 +max.0 +inf.0 +nan.0)]) + (define n + (count (λ (b) b) + (map (λ (f) (f x2 x1)) + (list fl2rational? fl2infinite? fl2nan?)))) + (unless (= n 1) (printf "x2 = ~v x1 = ~v~n" x2 x1))) + +#| +Tests to add + +(for*: ([x2 (list -inf.0 -max.0 -1.0 -min.0 -0.0 0.0 +min.0 +1.0 +max.0 +inf.0 +nan.0)] + [x1 (list -inf.0 -max.0 -1.0 -min.0 -0.0 0.0 +min.0 +1.0 +max.0 +inf.0 +nan.0)]) + (define n + (count (λ: ([b : Boolean]) b) + (map (λ: ([f : (Flonum Flonum -> Boolean)]) + (f x2 x1)) + (list fl2rational? fl2infinite? fl2nan?)))) + (unless (= n 1) (printf "x2 = ~v x1 = ~v~n" x2 x1))) + +fl2= +fl2> +fl2< +fl2>= +fl2<= + +(fl2step x2 x1 n/2) twice = (fl2step x2 x1 n) +|# + +(check-true (let-values ([(y2 y1) (fl+/error +max.hi +max.lo)]) + (fl2= y2 y1 +max.hi +max.lo))) + +(check-true (let*-values ([(y2 y1) (fl2next +max.hi +max.lo)]) + (fl2infinite? y2 y1))) + ;; =================================================================================================== ;; FPU testing -(parameterize ([print-test-progress? #f]) - (test-fpu 1000)) +(check-equal? (parameterize ([print-fp-test-progress? #f]) + (test-floating-point 1000)) + empty)