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