OLD | NEW |
1 #include "libm.h" | 1 #include "libm.h" |
2 | 2 |
3 /* sinh(x) = (exp(x) - 1/exp(x))/2 | 3 /* sinh(x) = (exp(x) - 1/exp(x))/2 |
4 * = (exp(x)-1 + (exp(x)-1)/exp(x))/2 | 4 * = (exp(x)-1 + (exp(x)-1)/exp(x))/2 |
5 * = x + x^3/6 + o(x^5) | 5 * = x + x^3/6 + o(x^5) |
6 */ | 6 */ |
7 double sinh(double x) | 7 double sinh(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 » double t, h, absx; | 11 } u = {.f = x}; |
| 12 uint32_t w; |
| 13 double t, h, absx; |
12 | 14 |
13 » h = 0.5; | 15 h = 0.5; |
14 » if (u.i >> 63) | 16 if (u.i >> 63) |
15 » » h = -h; | 17 h = -h; |
16 » /* |x| */ | 18 /* |x| */ |
17 » u.i &= (uint64_t)-1/2; | 19 u.i &= (uint64_t)-1 / 2; |
18 » absx = u.f; | 20 absx = u.f; |
19 » w = u.i >> 32; | 21 w = u.i >> 32; |
20 | 22 |
21 » /* |x| < log(DBL_MAX) */ | 23 /* |x| < log(DBL_MAX) */ |
22 » if (w < 0x40862e42) { | 24 if (w < 0x40862e42) { |
23 » » t = expm1(absx); | 25 t = expm1(absx); |
24 » » if (w < 0x3ff00000) { | 26 if (w < 0x3ff00000) { |
25 » » » if (w < 0x3ff00000 - (26<<20)) | 27 if (w < 0x3ff00000 - (26 << 20)) |
26 » » » » /* note: inexact and underflow are raised by exp
m1 */ | 28 /* note: inexact and underflow are raised by expm1 */ |
27 » » » » /* note: this branch avoids spurious underflow *
/ | 29 /* note: this branch avoids spurious underflow */ |
28 » » » » return x; | 30 return x; |
29 » » » return h*(2*t - t*t/(t+1)); | 31 return h * (2 * t - t * t / (t + 1)); |
30 » » } | 32 } |
31 » » /* note: |x|>log(0x1p26)+eps could be just h*exp(x) */ | 33 /* note: |x|>log(0x1p26)+eps could be just h*exp(x) */ |
32 » » return h*(t + t/(t+1)); | 34 return h * (t + t / (t + 1)); |
33 » } | 35 } |
34 | 36 |
35 » /* |x| > log(DBL_MAX) or nan */ | 37 /* |x| > log(DBL_MAX) or nan */ |
36 » /* note: the result is stored to handle overflow */ | 38 /* note: the result is stored to handle overflow */ |
37 » t = 2*h*__expo2(absx); | 39 t = 2 * h * __expo2(absx); |
38 » return t; | 40 return t; |
39 } | 41 } |
OLD | NEW |