OLD | NEW |
1 #include "libm.h" | 1 #include "libm.h" |
2 | 2 |
3 /* tanh(x) = (exp(x) - exp(-x))/(exp(x) + exp(-x)) | 3 /* tanh(x) = (exp(x) - exp(-x))/(exp(x) + exp(-x)) |
4 * = (exp(2*x) - 1)/(exp(2*x) - 1 + 2) | 4 * = (exp(2*x) - 1)/(exp(2*x) - 1 + 2) |
5 * = (1 - exp(-2*x))/(exp(-2*x) - 1 + 2) | 5 * = (1 - exp(-2*x))/(exp(-2*x) - 1 + 2) |
6 */ | 6 */ |
7 double tanh(double x) | 7 double tanh(double x) { |
8 { | 8 union { |
9 » union {double f; uint64_t i;} u = {.f = x}; | 9 double f; |
10 » uint32_t w; | 10 uint64_t i; |
11 » int sign; | 11 } u = {.f = x}; |
12 » double_t t; | 12 uint32_t w; |
| 13 int sign; |
| 14 double_t t; |
13 | 15 |
14 » /* x = |x| */ | 16 /* x = |x| */ |
15 » sign = u.i >> 63; | 17 sign = u.i >> 63; |
16 » u.i &= (uint64_t)-1/2; | 18 u.i &= (uint64_t)-1 / 2; |
17 » x = u.f; | 19 x = u.f; |
18 » w = u.i >> 32; | 20 w = u.i >> 32; |
19 | 21 |
20 » if (w > 0x3fe193ea) { | 22 if (w > 0x3fe193ea) { |
21 » » /* |x| > log(3)/2 ~= 0.5493 or nan */ | 23 /* |x| > log(3)/2 ~= 0.5493 or nan */ |
22 » » if (w > 0x40340000) { | 24 if (w > 0x40340000) { |
23 » » » /* |x| > 20 or nan */ | 25 /* |x| > 20 or nan */ |
24 » » » /* note: this branch avoids raising overflow */ | 26 /* note: this branch avoids raising overflow */ |
25 » » » t = 1 - 0/x; | 27 t = 1 - 0 / x; |
26 » » } else { | 28 } else { |
27 » » » t = expm1(2*x); | 29 t = expm1(2 * x); |
28 » » » t = 1 - 2/(t+2); | 30 t = 1 - 2 / (t + 2); |
29 » » } | 31 } |
30 » } else if (w > 0x3fd058ae) { | 32 } else if (w > 0x3fd058ae) { |
31 » » /* |x| > log(5/3)/2 ~= 0.2554 */ | 33 /* |x| > log(5/3)/2 ~= 0.2554 */ |
32 » » t = expm1(2*x); | 34 t = expm1(2 * x); |
33 » » t = t/(t+2); | 35 t = t / (t + 2); |
34 » } else if (w >= 0x00100000) { | 36 } else if (w >= 0x00100000) { |
35 » » /* |x| >= 0x1p-1022, up to 2ulp error in [0.1,0.2554] */ | 37 /* |x| >= 0x1p-1022, up to 2ulp error in [0.1,0.2554] */ |
36 » » t = expm1(-2*x); | 38 t = expm1(-2 * x); |
37 » » t = -t/(t+2); | 39 t = -t / (t + 2); |
38 » } else { | 40 } else { |
39 » » /* |x| is subnormal */ | 41 /* |x| is subnormal */ |
40 » » /* note: the branch above would not raise underflow in [0x1p-102
3,0x1p-1022) */ | 42 /* note: the branch above would not raise underflow in [0x1p-1023,0x1p-1022) |
41 » » FORCE_EVAL((float)x); | 43 */ |
42 » » t = x; | 44 FORCE_EVAL((float)x); |
43 » } | 45 t = x; |
44 » return sign ? -t : t; | 46 } |
| 47 return sign ? -t : t; |
45 } | 48 } |
OLD | NEW |