Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 1 | #include "libm.h" |
| 2 | |
Szabolcs Nagy | 482ccd2 | 2012-12-11 23:06:20 +0100 | [diff] [blame] | 3 | /* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */ |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 4 | double atanh(double x) |
| 5 | { |
Szabolcs Nagy | 482ccd2 | 2012-12-11 23:06:20 +0100 | [diff] [blame] | 6 | union {double f; uint64_t i;} u = {.f = x}; |
| 7 | unsigned e = u.i >> 52 & 0x7ff; |
| 8 | unsigned s = u.i >> 63; |
Szabolcs Nagy | f4d9bfb | 2013-09-05 16:57:46 +0000 | [diff] [blame] | 9 | double_t y; |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 10 | |
Szabolcs Nagy | 482ccd2 | 2012-12-11 23:06:20 +0100 | [diff] [blame] | 11 | /* |x| */ |
| 12 | u.i &= (uint64_t)-1/2; |
Szabolcs Nagy | f4d9bfb | 2013-09-05 16:57:46 +0000 | [diff] [blame] | 13 | y = u.f; |
Szabolcs Nagy | 482ccd2 | 2012-12-11 23:06:20 +0100 | [diff] [blame] | 14 | |
| 15 | if (e < 0x3ff - 1) { |
Szabolcs Nagy | f4d9bfb | 2013-09-05 16:57:46 +0000 | [diff] [blame] | 16 | if (e < 0x3ff - 32) { |
| 17 | /* handle underflow */ |
| 18 | if (e == 0) |
| 19 | FORCE_EVAL((float)y); |
| 20 | } else { |
| 21 | /* |x| < 0.5, up to 1.7ulp error */ |
| 22 | y = 0.5*log1p(2*y + 2*y*y/(1-y)); |
| 23 | } |
Szabolcs Nagy | 482ccd2 | 2012-12-11 23:06:20 +0100 | [diff] [blame] | 24 | } else { |
Szabolcs Nagy | f4d9bfb | 2013-09-05 16:57:46 +0000 | [diff] [blame] | 25 | /* avoid overflow */ |
| 26 | y = 0.5*log1p(2*(y/(1-y))); |
Szabolcs Nagy | 482ccd2 | 2012-12-11 23:06:20 +0100 | [diff] [blame] | 27 | } |
Szabolcs Nagy | f4d9bfb | 2013-09-05 16:57:46 +0000 | [diff] [blame] | 28 | return s ? -y : y; |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 29 | } |