1#lang typed/racket/base 2 3(require racket/performance-hint 4 "flonum-functions.rkt" 5 "flonum-constants.rkt" 6 "flonum-exp.rkt" 7 "flonum-log.rkt" 8 "flonum-error.rkt" 9 "flvector.rkt" 10 "utils.rkt") 11 12(provide flsqrt1pm1 13 flsinh flcosh fltanh 14 flasinh flacosh flatanh 15 make-flexpt flexpt+ flexpt1p 16 flfma) 17 18;; --------------------------------------------------------------------------------------------------- 19;; sqrt(1+x)-1 20 21(: flsqrt1pm1 (Float -> Float)) 22(define (flsqrt1pm1 x) 23 (cond [((flabs x) . fl> . 0.75) 24 (fl- (flsqrt (fl+ 1.0 x)) 1.0)] 25 [else 26 (flexpm1 (fl* 0.5 (fllog1p x)))])) 27 28;; --------------------------------------------------------------------------------------------------- 29;; Hyperbolic sine 30 31(: flsinh (Float -> Float)) 32(define (flsinh x) 33 (cond [(x . fl< . 0.0) 34 ;; Odd function 35 (- (flsinh (- x)))] 36 [(x . fl< . (flexpt 2.0 -26.0)) 37 ;; sinh(x) ~ x 38 x] 39 [(x . fl< . 18.5) 40 ;; sinh(x) = (exp(2*x) - 1) / (2*exp(x)) 41 (define y (flexpm1 x)) 42 (fl* 0.5 (fl+ y (fl/ y (fl+ y 1.0))))] 43 [(x . fl< . (fllog +max.0)) 44 ;; sinh(x) ~ exp(x) / 2 45 (fl* 0.5 (flexp x))] 46 [else 47 ;; sinh(x) ~ exp(x) / 2 = (exp(x/2) / 2) * exp(x/2) 48 (define y (flexp (fl* 0.5 x))) 49 (fl* (fl* 0.5 y) y)])) 50 51;; --------------------------------------------------------------------------------------------------- 52;; Hyperbolic cosine 53 54(: flcosh (Float -> Float)) 55(define (flcosh x) 56 ;; cosh(x) = cosh(-x) 57 (let ([x (flabs x)]) 58 (cond [(x . fl< . (flexpt 2.0 -26.0)) 59 ;; cosh(x) ~ 1 60 1.0] 61 [(x . fl< . (fl* 0.5 (fllog 2.0))) 62 ;; cosh(x) = 1 + (exp(x) - 1)^2 / (2*exp(x)) 63 (define y (flexpm1 x)) 64 (fl+ 1.0 (fl/ (fl* y y) (fl* 2.0 (fl+ 1.0 y))))] 65 [(x . fl< . 18.5) 66 ;; cosh(x) = (exp(x) + 1/exp(x)) / 2 67 (define y (flexp x)) 68 (fl+ (fl* 0.5 y) (fl/ 0.5 y))] 69 [(x . fl< . (fllog +max.0)) 70 ;; cosh(x) ~ exp(x) / 2 71 (fl* 0.5 (flexp x))] 72 [else 73 ;; cosh(x) ~ exp(x) / 2 = (exp(x/2) / 2) * exp(x/2) 74 (define y (flexp (fl* 0.5 x))) 75 (fl* (fl* 0.5 y) y)]))) 76 77;; --------------------------------------------------------------------------------------------------- 78;; Hyperbolic tangent 79 80(: fltanh (Float -> Float)) 81(define (fltanh x) 82 (cond [(x . fl< . 0.0) 83 ;; tanh(x) = -tanh(-x) 84 (- (fltanh (- x)))] 85 [(x . fl< . 1e-16) 86 ;; tanh(x) ~ x + x^2 87 (fl* x (fl+ 1.0 x))] 88 [(x . fl< . 0.5) 89 ;; tanh(x) = (exp(2*x) - 1) / (exp(2*x) + 1) 90 (define y (flexpm1 (fl* -2.0 x))) 91 (- (fl/ y (fl+ 2.0 y)))] 92 [(x . fl< . 19.5) 93 ;; tanh(x) = (exp(2*x) - 1) / (exp(2*x) + 1) 94 (define y (flexp (fl* 2.0 x))) 95 (fl/ (fl- y 1.0) (fl+ y 1.0))] 96 [(x . fl<= . +inf.0) 97 ;; tanh(x) ~ 1 98 1.0] 99 [else +nan.0])) 100 101;; --------------------------------------------------------------------------------------------------- 102;; Inverse hyperbolic sine 103 104(: flasinh (Float -> Float)) 105(define (flasinh x) 106 (cond [(x . fl< . 0.0) (- (flasinh (- x)))] 107 [(x . fl< . 2e-8) x] 108 [(x . fl< . 0.00018) 109 ;; Taylor series order 3 110 (fl* x (fl+ 1.0 (fl* (fl* #i-1/6 x) x)))] 111 [(x . fl< . 1.0) 112 ;; Standard definition, rearranged to preserve digits 113 (fllog1p (fl+ x (flsqrt1pm1 (fl* x x))))] 114 [(x . fl< . 3e3) 115 ;; Standard definition 116 (fllog (fl+ x (flsqrt (fl+ (fl* x x) 1.0))))] 117 [(x . fl< . 1e307) 118 ;; Laurent series in 1/x at 0+ order from -1 to 1 119 (fl+ (fllog (fl* x 2.0)) (fl/ 1.0 (fl* (fl* 4.0 x) x)))] 120 [else 121 ;; Laurent series, rearranged to not overflow 122 (fl+ (fllog x) (fllog 2.0))])) 123 124;; --------------------------------------------------------------------------------------------------- 125;; Inverse hyperbolic cosine 126 127(: flacosh (Float -> Float)) 128(define (flacosh x) 129 (cond [(x . fl< . 1.0) +nan.0] 130 [(x . fl< . 1.5) 131 ;; Standard definition, rearranged to preserve digits when x is near 1.0 132 (define y (fl- x 1.0)) 133 (fllog1p (fl+ y (flsqrt (fl+ (fl* y y) (fl* 2.0 y)))))] 134 [(x . fl< . 1e8) 135 ;; Standard definition 136 (fllog (fl+ x (flsqrt (fl- (fl* x x) 1.0))))] 137 [(x . fl< . 1e307) 138 ;; Laurent series in 1/x at 0+ order from -1 to 0 139 (fllog (fl* x 2.0))] 140 [else 141 ;; Laurent series, rearranged to avoid overflow 142 (fl+ (fllog x) (fllog 2.0))])) 143 144;; --------------------------------------------------------------------------------------------------- 145;; Inverse hyperbolic tangent 146 147(: flatanh (Float -> Float)) 148(define (flatanh x) 149 (cond [(x . fl< . 0.0) (- (flatanh (- x)))] 150 [(x . fl< . 1e-8) x] 151 [(x . fl< . 0.00015) 152 ;; Taylor series order 2 153 (fl+ x (fl* (fl* (fl* #i1/3 x) x) x))] 154 [(x . fl< . 0.5) 155 ;; Standard definition, rearranged to preserve digits when x is near 0.0 156 (fl* 0.5 (fl- (fllog1p x) (fllog1p (- x))))] 157 [(x . fl< . 1.0) 158 ;; Standard definition 159 (fl* 0.5 (fllog (fl/ (fl+ 1.0 x) (fl- 1.0 x))))] 160 [(x . fl= . 1.0) +inf.0] 161 [else +nan.0])) 162 163;; --------------------------------------------------------------------------------------------------- 164;; Exponential with high-precision bases 165 166(begin-encourage-inline 167 168 (: make-flexpt (Positive-Exact-Rational -> (Flonum -> Flonum))) 169 (define (make-flexpt b) 170 (define b-hi (fl b)) 171 (define b-lo (fl (- (/ (inexact->exact b-hi) b) 1))) 172 (cond [(fl= b-lo 0.0) (λ: ([x : Flonum]) (flexpt b-hi x))] 173 [else 174 (λ: ([x : Flonum]) 175 (fl/ (flexpt b-hi x) 176 (flexp (fl* x (fllog1p b-lo)))))])) 177 178 (: flexpt+ (Flonum Flonum Flonum -> Flonum)) 179 (define (flexpt+ a b y) 180 (define-values (x-hi x-lo) (fast-fl+/error a b)) 181 (fl/ (flexpt x-hi y) 182 (flexp (fl* y (fllog1p (- (/ x-lo x-hi))))))) 183 184 (: flexpt1p (Flonum Flonum -> Flonum)) 185 (define (flexpt1p x y) 186 (cond [(and (x . > . -0.5) (x . < . +inf.0)) 187 (define-values (a-hi a-lo) (fast-fl+/error 1.0 x)) 188 (fl/ (flexpt a-hi y) 189 (flexp (fl* y (fllog1p (- (/ a-lo a-hi))))))] 190 [else (flexpt (+ 1.0 x) y)])) 191 192 ) ; begin-encourage-inline 193 194;; --------------------------------------------------------------------------------------------------- 195;; Fused multiply-add 196 197(: slow-flfma (-> Flonum Flonum Flonum Flonum)) 198(define (slow-flfma a b c) 199 (define n (near-pow2 (max (flsqrt (abs a)) (flsqrt (abs b))))) 200 (define 1/n (/ 1.0 n)) 201 (* n n (fast-flfma (* a 1/n) (* b 1/n) (* c 1/n 1/n)))) 202 203(begin-encourage-inline 204 205 (: fast-flfma (-> Flonum Flonum Flonum Flonum)) 206 (define (fast-flfma a b c) 207 (let-values ([(d-hi d-lo) (fast-flfma/error a b c)]) 208 (+ d-hi d-lo))) 209 210 (: flfma (-> Flonum Flonum Flonum Flonum)) 211 (define (flfma a b c) 212 (let ([d (fast-flfma a b c)]) 213 (if (flrational? d) d (slow-flfma a b c)))) 214 215 ) ; begin-encourage-inline 216