OLD | NEW |
1 #include "libm.h" | 1 #include "libm.h" |
2 | 2 |
3 /* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */ | 3 /* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */ |
4 double atanh(double x) | 4 double atanh(double x) { |
5 { | 5 union { |
6 » union {double f; uint64_t i;} u = {.f = x}; | 6 double f; |
7 » unsigned e = u.i >> 52 & 0x7ff; | 7 uint64_t i; |
8 » unsigned s = u.i >> 63; | 8 } u = {.f = x}; |
9 » double_t y; | 9 unsigned e = u.i >> 52 & 0x7ff; |
| 10 unsigned s = u.i >> 63; |
| 11 double_t y; |
10 | 12 |
11 » /* |x| */ | 13 /* |x| */ |
12 » u.i &= (uint64_t)-1/2; | 14 u.i &= (uint64_t)-1 / 2; |
13 » y = u.f; | 15 y = u.f; |
14 | 16 |
15 » if (e < 0x3ff - 1) { | 17 if (e < 0x3ff - 1) { |
16 » » if (e < 0x3ff - 32) { | 18 if (e < 0x3ff - 32) { |
17 » » » /* handle underflow */ | 19 /* handle underflow */ |
18 » » » if (e == 0) | 20 if (e == 0) |
19 » » » » FORCE_EVAL((float)y); | 21 FORCE_EVAL((float)y); |
20 » » } else { | 22 } else { |
21 » » » /* |x| < 0.5, up to 1.7ulp error */ | 23 /* |x| < 0.5, up to 1.7ulp error */ |
22 » » » y = 0.5*log1p(2*y + 2*y*y/(1-y)); | 24 y = 0.5 * log1p(2 * y + 2 * y * y / (1 - y)); |
23 » » } | 25 } |
24 » } else { | 26 } else { |
25 » » /* avoid overflow */ | 27 /* avoid overflow */ |
26 » » y = 0.5*log1p(2*(y/(1-y))); | 28 y = 0.5 * log1p(2 * (y / (1 - y))); |
27 » } | 29 } |
28 » return s ? -y : y; | 30 return s ? -y : y; |
29 } | 31 } |
OLD | NEW |