Fixed some limits and accuracy issues with flonum expansions
This commit is contained in:
Neil Toronto
parent
57f233cff1
commit
564f589601
|
|
@ -10,16 +10,20 @@ Discrete & Computational Geometry 18(3):305–363, October 1997
|
|||
Other parts shamelessly stolen from crlibm (which is LGPL)
|
||||
|#
|
||||
|
||||
(require "../flonum-functions.rkt"
|
||||
(require racket/math
|
||||
"../flonum-functions.rkt"
|
||||
"../flonum-bits.rkt"
|
||||
"../flonum-error.rkt"
|
||||
"../flonum-constants.rkt"
|
||||
"../utils.rkt")
|
||||
|
||||
(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/
|
||||
(provide fl2? fl2zero? fl2rational? fl2positive? fl2negative? fl2infinite? fl2nan?
|
||||
fl2 fl2->real
|
||||
fl2ulp fl2ulp-error fl2step fl2next fl2prev
|
||||
+max.hi +max.lo -max.hi -max.lo
|
||||
+max-subnormal.hi -max-subnormal.hi
|
||||
fl2abs fl2+ fl2- fl2= fl2> fl2< fl2>= fl2<=
|
||||
fl2*split-fl fl2* fl2sqr fl2/
|
||||
fl2sqrt flsqrt/error)
|
||||
|
||||
(: floverlapping? (Flonum Flonum -> Boolean))
|
||||
|
|
@ -38,30 +42,58 @@ Other parts shamelessly stolen from crlibm (which is LGPL)
|
|||
(cond [(flrational? x2)
|
||||
(cond [(flrational? x1)
|
||||
(cond [((flabs x2) . < . (flabs x1)) #f]
|
||||
[else (not (floverlapping? x2 x1))])]
|
||||
[else (not (floverlapping? x2 x1))])]
|
||||
[else #f])]
|
||||
[else (fl= x1 0.0)]))
|
||||
[else
|
||||
(fl= x1 0.0)]))
|
||||
|
||||
(define-syntax-rule (define-simple-fl2-predicate fl2pred? flpred?)
|
||||
(begin
|
||||
(: fl2pred? (Flonum Flonum -> Boolean))
|
||||
(define (fl2pred? x2 x1)
|
||||
(flpred? (fl+ x2 x1)))))
|
||||
|
||||
(define-simple-fl2-predicate fl2zero? (λ (x) (fl= x 0.0)))
|
||||
(define-simple-fl2-predicate fl2positive? (λ (x) (fl> x 0.0)))
|
||||
(define-simple-fl2-predicate fl2negative? (λ (x) (fl< x 0.0)))
|
||||
(define-simple-fl2-predicate fl2rational? flrational?)
|
||||
(define-simple-fl2-predicate fl2nan? flnan?)
|
||||
(define-simple-fl2-predicate fl2infinite? flinfinite?)
|
||||
|
||||
;; ===================================================================================================
|
||||
;; Conversion
|
||||
|
||||
(: 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 (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
|
||||
(define x2 (fl x))
|
||||
(if (flinfinite? x2)
|
||||
(values x2 0.0)
|
||||
(let* ([x (- x (inexact->exact x2))]
|
||||
[x1 (fl x)]
|
||||
[x (- x (inexact->exact x1))])
|
||||
(let-values ([(x2 x1) (fl+/error x2 x1)])
|
||||
(values x2 (fl+ x1 (fl x))))))])]
|
||||
[(x2 x1)
|
||||
(if (and (fl= x2 0.0) (fl= x1 0.0))
|
||||
(values x2 0.0)
|
||||
(fl+/error x2 x1))]))
|
||||
|
||||
(: fl2eqv? (case-> (Flonum Flonum Flonum -> Boolean)
|
||||
(Flonum Flonum Flonum Flonum -> Boolean)))
|
||||
(define (fl2eqv? x2 x1 y2 [y1 0.0])
|
||||
(and (eqv? x2 y2) (fl= x1 y1)))
|
||||
|
||||
(: fl2->real (Flonum Flonum -> Real))
|
||||
(define (fl2->real x2 x1)
|
||||
(cond [(and (flrational? x2) (flrational? x1))
|
||||
(+ (inexact->exact x2) (inexact->exact x1))]
|
||||
[else (fl+ x1 x2)]))
|
||||
(if (flrational? x2)
|
||||
(+ (inexact->exact x2) (inexact->exact x1))
|
||||
x2))
|
||||
|
||||
(: fl4->fl2 (Flonum Flonum Flonum Flonum -> (Values Flonum Flonum)))
|
||||
(define (fl4->fl2 e4 e3 e2 e1)
|
||||
|
|
@ -73,7 +105,7 @@ Other parts shamelessly stolen from crlibm (which is LGPL)
|
|||
(: 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)))]))
|
||||
[else (flmax +min.0 (fl* (flulp x2) epsilon.0))]))
|
||||
|
||||
(: fl2ulp-error (Flonum Flonum Real -> Flonum))
|
||||
(define (fl2ulp-error x2 x1 r)
|
||||
|
|
@ -89,20 +121,47 @@ Other parts shamelessly stolen from crlibm (which is LGPL)
|
|||
(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))
|
||||
(define-values (+max.hi +max.lo)
|
||||
(values +max.0 (flprev (* 0.5 (flulp +max.0)))))
|
||||
|
||||
(define-values (-max.hi -max.lo)
|
||||
(values (- +max.hi) (- +max.lo)))
|
||||
|
||||
(: fl2step (Flonum Flonum Integer -> (Values Flonum Flonum)))
|
||||
(define (fl2step x2 x1 n)
|
||||
(let-values ([(x2 x1) (fast-fl+/error x2 x1)])
|
||||
(cond [(flnan? x2) (values +nan.0 0.0)]
|
||||
[(fl= x2 +inf.0) (fl+/error +max.hi (flstep +max.lo (+ n 1)))]
|
||||
[(fl= x2 -inf.0) (fl+/error -max.hi (flstep -max.lo (- n 1)))]
|
||||
[else (fl+/error x2 (flstep x1 n))])))
|
||||
|
||||
(: fl2next (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
(define (fl2next x2 x1) (fl2step x2 x1 1))
|
||||
|
||||
(: fl2prev (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
(define (fl2prev x2 x1) (fl2step x2 x1 -1))
|
||||
|
||||
(define +min-normal.hi (fl/ (flnext +max-subnormal.0) epsilon.0))
|
||||
|
||||
(define-values (+max-subnormal.hi +max-subnormal.lo)
|
||||
(fl2prev +min-normal.hi 0.0))
|
||||
|
||||
(define-values (-max-subnormal.hi -max-subnormal.lo)
|
||||
(values (- +max-subnormal.hi) (- +max-subnormal.lo)))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; 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))]))
|
||||
(define fl2abs
|
||||
(case-lambda
|
||||
[(x) (values (flabs x) 0.0)]
|
||||
[(x2 x1)
|
||||
(cond [(flnan? x2) (values +nan.0 0.0)]
|
||||
[(fl= x2 0.0) (values 0.0 0.0)]
|
||||
[(fl> x2 0.0) (values x2 x1)]
|
||||
[else (values (- x2) (- x1))])]))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; Addition and subtraction
|
||||
|
|
@ -125,6 +184,22 @@ Other parts shamelessly stolen from crlibm (which is LGPL)
|
|||
(define (fl2- x2 x1 y2 [y1 0.0])
|
||||
(fl2+ x2 x1 (- y2) (- y1)))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; Comparison
|
||||
|
||||
(define-syntax-rule (define-fl2-comparison name flcomp)
|
||||
(begin
|
||||
(: name (Flonum Flonum Flonum Flonum -> Boolean))
|
||||
(define (name x2 x1 y2 y1)
|
||||
(let-values ([(z2 z1) (fl2- x2 x1 y2 y1)])
|
||||
((fl+ z2 z1) . flcomp . 0.0)))))
|
||||
|
||||
(define-fl2-comparison fl2= fl=)
|
||||
(define-fl2-comparison fl2> fl>)
|
||||
(define-fl2-comparison fl2< fl<)
|
||||
(define-fl2-comparison fl2>= fl>=)
|
||||
(define-fl2-comparison fl2<= fl<=)
|
||||
|
||||
;; ===================================================================================================
|
||||
;; Multiplication and division
|
||||
|
||||
|
|
@ -164,10 +239,13 @@ Other parts shamelessly stolen from crlibm (which is LGPL)
|
|||
(Flonum Flonum Flonum Flonum -> (Values Flonum Flonum))))
|
||||
(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)))
|
||||
(cond [(fl= z 0.0) (values z 0.0)]
|
||||
[(flsubnormal? z) (values z 0.0)]
|
||||
[(and (flrational? x2) (flrational? y2) (z . fl>= . -inf.0) (z . fl<= . +inf.0))
|
||||
(define dx (near-pow2 x2))
|
||||
(define dy (near-pow2 y2))
|
||||
(define d (fl* dx dy))
|
||||
(define d? (and (d . fl> . 0.0) (d . fl< . +inf.0)))
|
||||
(let* ([x2 (fl/ x2 dx)]
|
||||
[x1 (fl/ x1 dx)]
|
||||
[y2 (fl/ y2 dy)]
|
||||
|
|
@ -187,8 +265,8 @@ Other parts shamelessly stolen from crlibm (which is LGPL)
|
|||
(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)))]
|
||||
[z2 (if d? (fl* z2 d) (fl* (fl* z2 dx) dy))])
|
||||
(values z2 (if (flrational? z2) (if d? (fl* z1 d) (fl* (fl* z1 dx) dy)) 0.0)))]
|
||||
[else
|
||||
(values z 0.0)]))
|
||||
|
||||
|
|
@ -200,11 +278,14 @@ Other parts shamelessly stolen from crlibm (which is LGPL)
|
|||
[(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)]
|
||||
(cond [(fl= z 0.0) (values z 0.0)]
|
||||
[(flsubnormal? z) (values z 0.0)]
|
||||
[(and (flrational? x2) (z . fl>= . -inf.0) (z . fl<= . +inf.0))
|
||||
(define dx (near-pow2 x2))
|
||||
(define d (fl* dx dx))
|
||||
(define d? (and (d . fl> . 0.0) (d . fl< . +inf.0)))
|
||||
(let* ([x2 (fl/ x2 dx)]
|
||||
[x1 (fl/ x1 dx)]
|
||||
[up (fl* x2 (fl+ 1.0 (flexpt 2.0 27.0)))]
|
||||
[u1 (fl+ (fl- x2 up) up)]
|
||||
[u2 (fl- x2 u1)]
|
||||
|
|
@ -215,8 +296,8 @@ Other parts shamelessly stolen from crlibm (which is LGPL)
|
|||
(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)))]
|
||||
[z2 (if d? (fl* z2 d) (fl* (fl* z2 dx) dx))])
|
||||
(values z2 (if (flrational? z2) (if d? (fl* z1 d) (fl* (fl* z1 dx) dx)) 0.0)))]
|
||||
[else
|
||||
(values z 0.0)])]))
|
||||
|
||||
|
|
@ -245,8 +326,8 @@ Other parts shamelessly stolen from crlibm (which is LGPL)
|
|||
(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))]
|
||||
(cond [(x2 . fl<= . +max-subnormal.hi) (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)]))
|
||||
|
|
|
|||
|
|
@ -88,6 +88,7 @@ negative.
|
|||
(cons y2 y1)))))
|
||||
|#
|
||||
|
||||
(: positive-ds (Vectorof (Pair Flonum Flonum)))
|
||||
(define positive-ds
|
||||
#((0.00048511293333186416 . 7.287053259455084e-21)
|
||||
(0.0005002765039165264 . 5.174812552711603e-20)
|
||||
|
|
@ -553,6 +554,7 @@ negative.
|
|||
(1.0134858165681662e+305 . -1.1960732643227e+288)
|
||||
(+inf.0 . -inf.0)))
|
||||
|
||||
(: negative-ds (Vectorof (Pair Flonum Flonum)))
|
||||
(define negative-ds
|
||||
#((-0.0004848777128822606 . -2.106959040416988e-20)
|
||||
(-0.0005000263524810411 . 1.7625505823099235e-20)
|
||||
|
|
|
|||
|
|
@ -72,14 +72,13 @@ Note to self: exp(x+y)-1 = (exp(x)-1) * (exp(y)-1) + (exp(x)-1) + (exp(y)-1)
|
|||
|
||||
;; See "expansion-exp-reduction.rkt" for details on the argument reduction here
|
||||
|
||||
(: flexp/error (Flonum -> (Values Flonum Flonum)))
|
||||
(define (flexp/error x)
|
||||
(let: loop : (Values Flonum Flonum) ([x : Flonum x] [n : Nonnegative-Fixnum 0])
|
||||
(: flexp/positive/error (Flonum -> (Values Flonum Flonum)))
|
||||
(define (flexp/positive/error x)
|
||||
(let: loop ([x x] [n : Nonnegative-Fixnum 0])
|
||||
(cond
|
||||
[(or ((flabs x) . fl< . exp-min) (n . fx> . 15)) ; even n > 5 should never happen
|
||||
(let-values ([(y2 y1) (flexpm1-small/error x)])
|
||||
(fl2+ y2 y1 1.0))]
|
||||
[(x . fl< . (- (fllog +max.0))) (values (flexp x) 0.0)]
|
||||
[(x . fl> . (fllog +max.0)) (values +inf.0 0.0)]
|
||||
[(rational? x)
|
||||
(let*-values ([(x d2 d1) (flexpm1-reduction x)]
|
||||
|
|
@ -87,14 +86,29 @@ Note to self: exp(x+y)-1 = (exp(x)-1) * (exp(y)-1) + (exp(x)-1) + (exp(y)-1)
|
|||
[(y2 y1) (loop x (fx+ n 1))])
|
||||
(fl2* d2 d1 y2 y1))]
|
||||
[else
|
||||
(values +nan.0 +nan.0)])))
|
||||
(values +nan.0 0.0)])))
|
||||
|
||||
(: flexp/error (Flonum -> (Values Flonum Flonum)))
|
||||
(define (flexp/error x)
|
||||
(cond [(x . fl> . 0.0) (flexp/positive/error x)]
|
||||
[(x . fl= . 0.0) (values 1.0 0.0)]
|
||||
[(x . fl> . (- (fllog +max.0)))
|
||||
(let*-values ([(y2 y1) (flexp/positive/error (- x))])
|
||||
(fl2/ 1.0 0.0 y2 y1))]
|
||||
[(x . fl>= . -inf.0) (values (flexp x) 0.0)]
|
||||
[else (values +nan.0 0.0)]))
|
||||
|
||||
(: fl2exp (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
(define (fl2exp x2 x1)
|
||||
(let*-values ([(a2 a1) (flexp/error x2)]
|
||||
[(b2 b1) (flexpm1-tiny/error x1)]
|
||||
[(b2 b1) (fl2+ b2 b1 1.0)])
|
||||
(fl2* a2 a1 b2 b1)))
|
||||
(cond [(x2 . fl> . -746.0)
|
||||
(let*-values ([(a2 a1) (flexp/error x2)]
|
||||
[(b2 b1) (flexpm1-tiny/error x1)]
|
||||
[(b2 b1) (fl2+ b2 b1 1.0)])
|
||||
(fl2* a2 a1 b2 b1))]
|
||||
[(x2 . fl>= . -inf.0)
|
||||
(values 0.0 0.0)]
|
||||
[else
|
||||
(values +nan.0 0.0)]))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; expm1
|
||||
|
|
@ -118,12 +132,8 @@ The computation of R+1 definitely loses precision, but it doesn't seem to matter
|
|||
(: flexpm1/error (Flonum -> (Values Flonum Flonum)))
|
||||
(define (flexpm1/error x)
|
||||
(cond
|
||||
[(x . < . -1.0)
|
||||
;; exp(x) is near zero here, so this is more accurate
|
||||
(let-values ([(y2 y1) (flexp/error x)])
|
||||
(fl2+ y2 y1 -1.0))]
|
||||
[else
|
||||
(let: loop : (Values Flonum Flonum) ([x : Flonum x] [n : Nonnegative-Fixnum 0])
|
||||
[(x . fl> . -1.0)
|
||||
(let: loop ([x x] [n : Nonnegative-Fixnum 0])
|
||||
(cond [(or ((flabs x) . fl< . expm1-min) (n . fx> . 15)) ; even n > 5 should never happen
|
||||
(flexpm1-small/error x)]
|
||||
[(x . fl< . (- (fllog +max.0))) (values -1.0 0.0)]
|
||||
|
|
@ -135,14 +145,29 @@ The computation of R+1 definitely loses precision, but it doesn't seem to matter
|
|||
[(y2 y1) (fl2* d2 d1 w2 w1)])
|
||||
(fl2+ y2 y1 r2 r1))]
|
||||
[else
|
||||
(values +nan.0 +nan.0)]))]))
|
||||
(values +nan.0 0.0)]))]
|
||||
[(x . fl> . (fllog (* 0.25 epsilon.0)))
|
||||
;; exp(x) is near zero here, so this is more accurate
|
||||
(let-values ([(y2 y1) (flexp/error x)])
|
||||
(fl2+ y2 y1 -1.0))]
|
||||
[else
|
||||
;; The high-order flonum is -1 here, so return -1 + exp(x) directly
|
||||
(values -1.0 (flexp x))]))
|
||||
|
||||
(define-values (expm1-max-hi expm1-max-lo)
|
||||
(values 709.782712893384 2.3691528222554853e-14))
|
||||
|
||||
(: fl2expm1 (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
(define (fl2expm1 x2 x1)
|
||||
(define x (fl+ x2 x1))
|
||||
(cond
|
||||
[((fl+ x1 x2) . < . -1.0)
|
||||
[(x . fl< . -1.0)
|
||||
(let-values ([(y2 y1) (fl2exp x2 x1)])
|
||||
(fl2+ y2 y1 -1.0))]
|
||||
[(x2 x1 . fl2>= . expm1-max-hi expm1-max-lo)
|
||||
(values +inf.0 0.0)]
|
||||
[(x . fl= . 0.0)
|
||||
(values x2 0.0)]
|
||||
[else
|
||||
(let*-values ([(a2 a1) (flexpm1/error x2)]
|
||||
[(b2 b1) (flexpm1-tiny/error x1)]
|
||||
|
|
|
|||
|
|
@ -50,7 +50,8 @@ A value of k that reduces any x to [0.5, 1.0] is
|
|||
(define (fl2log x2 x1)
|
||||
(define x (fl+ x1 x2))
|
||||
(cond [(x . fl<= . 0.0) (cond [(fl= x 0.0) (values -inf.0 0.0)]
|
||||
[else (values +nan.0 +nan.0)])]
|
||||
[else (values +nan.0 0.0)])]
|
||||
[(x . fl= . +inf.0) (values +inf.0 0.0)]
|
||||
[(or (x . fl< . 0.5) (x . fl> . 2.5))
|
||||
;; Reduce arguments
|
||||
(let*-values ([(k x2 x1) (fl2log-reduction x2 x1)]
|
||||
|
|
@ -80,10 +81,8 @@ A `k' that reduces any argument `x' to (-1/2,1/2) is
|
|||
k = round(log1p(x)/log(2))
|
||||
|#
|
||||
|
||||
(: fl2log1p-reduction (Flonum Flonum -> (Values Flonum Flonum Flonum)))
|
||||
(define (fl2log1p-reduction x2 x1)
|
||||
(define-values (a2 a1) (fl2+ x2 x1 1.0))
|
||||
(define y (fllog+ a2 a1))
|
||||
(: fl2log1p-reduction (Flonum Flonum Flonum Flonum Flonum -> (Values Flonum Flonum Flonum)))
|
||||
(define (fl2log1p-reduction x2 x1 a2 a1 y)
|
||||
(define k (flround (fl/ y (fllog 2.0))))
|
||||
(define 2^k (flexpt 2.0 k))
|
||||
(define-values (j2 j1) (fast-fl-/error (/ 1.0 2^k) 1.0))
|
||||
|
|
@ -96,9 +95,10 @@ A `k' that reduces any argument `x' to (-1/2,1/2) is
|
|||
(define-values (a2 a1) (fl2+ x2 x1 1.0))
|
||||
(define y (fllog+ a2 a1))
|
||||
(cond
|
||||
[(y . fl< . -0.5) (fl2log a2 a1)]
|
||||
[(or (y . fl< . -0.5) (a2 . fl> . 2.0)) (fl2log a2 a1)]
|
||||
[(fl= (fl+ x2 x1) 0.0) (values x2 0.0)]
|
||||
[(y . fl> . 0.5)
|
||||
(let*-values ([(k x2 x1) (fl2log1p-reduction x2 x1)]
|
||||
(let*-values ([(k x2 x1) (fl2log1p-reduction x2 x1 a2 a1 y)]
|
||||
[(y2 y1) (fl2log1p x2 x1)]
|
||||
[(z2 z1) (fl2* log2-hi log2-lo k)])
|
||||
(fl2+ y2 y1 z2 z1))]
|
||||
|
|
|
|||
|
|
@ -146,12 +146,15 @@
|
|||
(values x2 (if (and (flrational? x2) (not (flsubnormal? x2)))
|
||||
(let*-values ([(da db) (values (near-pow2 a) (near-pow2 b))]
|
||||
[(d) (fl* da db)]
|
||||
[(d?) (and (d . fl> . 0.0) (d . fl< . +inf.0))]
|
||||
[(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)))))
|
||||
[(b2 b1) (flsplit (fl/ b db))]
|
||||
[(x2) (if d? (fl/ x2 d) (fl/ (fl/ x2 da) db))]
|
||||
[(x1) (- (fl- (fl- (fl- (fl- x2 (fl* a2 b2))
|
||||
(fl* a1 b2))
|
||||
(fl* a2 b1))
|
||||
(fl* a1 b1)))])
|
||||
(if d? (fl* x1 d) (fl* (fl* x1 da) db)))
|
||||
0.0))))
|
||||
|
||||
(: flsqr/error (Flonum -> (Values Flonum Flonum)))
|
||||
|
|
@ -160,10 +163,13 @@
|
|||
(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)))))
|
||||
[(d^2?) (and (d^2 . fl> . 0.0) (d^2 . fl< . +inf.0))]
|
||||
[(a2 a1) (flsplit (fl/ a d))]
|
||||
[(x2) (if d^2? (fl/ x2 d^2) (fl/ (fl/ x2 d) d))]
|
||||
[(x1) (- (fl- (fl- (fl- x2 (fl* a2 a2))
|
||||
(fl* 2.0 (fl* a1 a2)))
|
||||
(fl* a1 a1)))])
|
||||
(if d^2? (fl* x1 d^2) (fl* (fl* x1 d) d)))
|
||||
0.0))))
|
||||
|
||||
(: fl//error (Flonum Flonum -> (Values Flonum Flonum)))
|
||||
|
|
|
|||
|
|
@ -9,7 +9,7 @@
|
|||
|
||||
(provide (all-from-out racket/flonum)
|
||||
fl
|
||||
flsubnormal? flrational? flnan? flinteger?
|
||||
flsubnormal? flrational? flinfinite? flnan? flinteger?
|
||||
flnext* flprev*
|
||||
flulp-error
|
||||
fleven? flodd? flsgn flhypot fllog/base
|
||||
|
|
@ -40,6 +40,10 @@
|
|||
(λ: ([x : Flonum])
|
||||
(and (x . fl> . -inf.0) (x . fl< . +inf.0))))
|
||||
|
||||
(define flinfinite?
|
||||
(λ: ([x : Flonum])
|
||||
(or (x . fl= . -inf.0) (x . fl= . +inf.0))))
|
||||
|
||||
(define flnan?
|
||||
(λ: ([x : Flonum])
|
||||
(not (and (x . fl>= . -inf.0) (x . fl<= . +inf.0)))))
|
||||
|
|
|
|||
|
|
@ -9,7 +9,7 @@
|
|||
(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)))))))
|
||||
(flexpt 2.0 (flmax -1023.0 (flmin 1023.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
|
||||
|
|
|
|||
|
|
@ -6,6 +6,8 @@
|
|||
typed/rackunit
|
||||
"../base/base-random.rkt"
|
||||
"../flonum/expansion/expansion-base.rkt"
|
||||
"../flonum/expansion/expansion-exp.rkt"
|
||||
"../flonum/expansion/expansion-log.rkt"
|
||||
"../flonum/flonum-functions.rkt"
|
||||
"../flonum/flonum-constants.rkt"
|
||||
"../flonum/flonum-bits.rkt"
|
||||
|
|
@ -26,8 +28,15 @@
|
|||
;; Allowable error for different kinds of functions, in ulps
|
||||
(define flonum-fun-ulps 0.5)
|
||||
(define flonum/error-fun-ulps 0.5)
|
||||
(define flexp/error-fun-ulps 2.0)
|
||||
(define fl2-conversion-ulps 0.5)
|
||||
(define unary-fl2-fun-ulps 1.0)
|
||||
(define binary-fl2-fun-ulps 8.0)
|
||||
(define fl2exp-fun-ulps 2.5)
|
||||
(define fl2log-fun-ulps 1.5)
|
||||
|
||||
(: current-max-ulp-error (Parameterof Nonnegative-Flonum))
|
||||
(define current-max-ulp-error (make-parameter 0.0))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; Helpers
|
||||
|
|
@ -94,13 +103,13 @@
|
|||
;; 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-subnormal.hi)) (flsqrt +max-subnormal.hi)
|
||||
(- (flsqrt +max.0)) (flsqrt +max.0)
|
||||
;; Test exp limits
|
||||
(fllog +min.0) (fllog +max-subnormal.0) (fllog +max-fl2-subnormal.0) (fllog +max.0)
|
||||
(fllog +min.0) (fllog +max-subnormal.0) (fllog +max-subnormal.hi) (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
|
||||
-inf.0 -max.0 -1.0 -max-subnormal.hi -max-subnormal.0 -min.0 -0.0
|
||||
+inf.0 +max.0 +1.0 +max-subnormal.hi +max-subnormal.0 +min.0 +0.0
|
||||
+nan.0))
|
||||
|
||||
(define standard-rs
|
||||
|
|
@ -123,8 +132,8 @@
|
|||
(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)))
|
||||
(define min-fl2-subnormal-ord (flonum->ordinal -max-subnormal.hi))
|
||||
(define max-fl2-subnormal-ord (+ 1 (flonum->ordinal +max-subnormal.hi)))
|
||||
|
||||
(: sample-flonum (case-> (Integer -> (Listof Flonum))
|
||||
(Integer Flonum Flonum -> (Listof Flonum))))
|
||||
|
|
@ -183,7 +192,7 @@
|
|||
[i (in-naturals 1)])
|
||||
(maybe-print-progress name i m)
|
||||
(list (list name x) (unary-flonum-fun-error f g x)))
|
||||
flonum-fun-ulps))
|
||||
(current-max-ulp-error)))
|
||||
|
||||
(: binary-flonum-fun-error
|
||||
((Flonum Flonum -> Flonum) (Bigfloat Bigfloat -> Bigfloat) Flonum Flonum -> Any))
|
||||
|
|
@ -205,7 +214,7 @@
|
|||
[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))
|
||||
(current-max-ulp-error)))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; fl2 conversion
|
||||
|
|
@ -230,7 +239,7 @@
|
|||
[i (in-naturals 1)])
|
||||
(maybe-print-progress 'fl2 i m)
|
||||
(list (list 'fl2 x) (fl2-conversion-error x)))
|
||||
flonum/error-fun-ulps))
|
||||
(current-max-ulp-error)))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; Flonum arithmetic with error
|
||||
|
|
@ -262,7 +271,7 @@
|
|||
[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))
|
||||
(current-max-ulp-error)))
|
||||
|
||||
(: test-binary-flonum/error-fun
|
||||
(Symbol (Flonum Flonum -> (Values Flonum Flonum)) (Bigfloat Bigfloat -> Bigfloat) Integer
|
||||
|
|
@ -278,7 +287,7 @@
|
|||
[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))
|
||||
(current-max-ulp-error)))
|
||||
|
||||
;; ===================================================================================================
|
||||
;; Flonum expansions
|
||||
|
|
@ -302,7 +311,7 @@
|
|||
(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))
|
||||
(current-max-ulp-error)))
|
||||
|
||||
(: binary-fl2-fun-error ((Flonum Flonum Flonum Flonum -> (Values Flonum Flonum))
|
||||
(Bigfloat Bigfloat -> Bigfloat)
|
||||
|
|
@ -330,83 +339,110 @@
|
|||
(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))
|
||||
(current-max-ulp-error)))
|
||||
|
||||
;; ===================================================================================================
|
||||
|
||||
(: 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)
|
||||
'()))
|
||||
(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)
|
||||
'())))
|
||||
|
||||
(: 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)
|
||||
'()))
|
||||
(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)
|
||||
'())))
|
||||
|
||||
(: 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)
|
||||
'()))
|
||||
(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)
|
||||
'())))
|
||||
|
||||
(: 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)
|
||||
'()))
|
||||
(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)
|
||||
'())))
|
||||
|
||||
(: 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)
|
||||
'())))
|
||||
|
||||
(: 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)
|
||||
'()))
|
||||
(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-fpu-non-trig/fl2 (Natural -> Any))
|
||||
(define (test-fpu-non-trig/fl2 n)
|
||||
(check-equal? (test-unary-fl2-fun 'fl2sqrt fl2sqrt bfsqrt 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)
|
||||
'()))
|
||||
)
|
||||
|
||||
(: test-fpu (Natural -> Any))
|
||||
(define (test-fpu n)
|
||||
|
|
@ -414,5 +450,43 @@
|
|||
(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))
|
||||
|
||||
(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)))
|
||||
|
||||
#|
|
||||
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)))
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user