Added flexp2', fllog2', `fllogb'; refactored and documented flonum testing
Note: With this refactoring, `math/utils' no longer depends on `rackunit'. * (flexp2 x) computes (flexpt 2.0 x) but in about 1/3 the time for integer `x' using a lookup table. Written for exact argument reduction in `fllog2' after discovering that (flexpt 2.0 x) was the main performance bottleneck. * (fllog2 x) computes (/ (fllog x) (fllog 2.0)) with near perfect accuracy. Invented an algorithm to compute it with at least 8 extra bits before final rounding; quite pleased with the result. Needed `fllog2' to ensure (fllogb 2.0 x) would be exact when `x' is a power of two. * (fllogb b x) computes (/ (fllog x) (fllog b)) with better accuracy, and also handles limit values in a way that's consistent with the mathematical limits. When those are ambiguous, it's consistent with `flexpt', which follows IEEE 754 and C99. Otherwise returns +nan.0. See docs for details. * `bflogb' is currently just for testing `fllogb'. * Refactored FPU testing and documented it. So far, the only documented way to do it is by calling `test-floating-point', which runs a comprehensive deterministic+randomized suite of tests and returns a list representing failed tests. I'll document individual tests after I document flonum expansions and result/error functions like `fl+/error'. * Added `fllog2' and `fllogb' to the flonum tests.
This commit is contained in:
Neil Toronto
parent
e0a7329d86
commit
aed3b39546
|
|
@ -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))]))
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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])]))
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
))
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user