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 | float atanhf(float x) |
| 5 | { |
Szabolcs Nagy | 482ccd2 | 2012-12-11 23:06:20 +0100 | [diff] [blame] | 6 | union {float f; uint32_t i;} u = {.f = x}; |
| 7 | unsigned s = u.i >> 31; |
Szabolcs Nagy | f4d9bfb | 2013-09-05 16:57:46 +0000 | [diff] [blame] | 8 | float_t y; |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 9 | |
Szabolcs Nagy | 482ccd2 | 2012-12-11 23:06:20 +0100 | [diff] [blame] | 10 | /* |x| */ |
| 11 | u.i &= 0x7fffffff; |
Szabolcs Nagy | f4d9bfb | 2013-09-05 16:57:46 +0000 | [diff] [blame] | 12 | y = u.f; |
Szabolcs Nagy | 482ccd2 | 2012-12-11 23:06:20 +0100 | [diff] [blame] | 13 | |
| 14 | if (u.i < 0x3f800000 - (1<<23)) { |
Szabolcs Nagy | f4d9bfb | 2013-09-05 16:57:46 +0000 | [diff] [blame] | 15 | if (u.i < 0x3f800000 - (32<<23)) { |
| 16 | /* handle underflow */ |
| 17 | if (u.i < (1<<23)) |
| 18 | FORCE_EVAL((float)(y*y)); |
| 19 | } else { |
| 20 | /* |x| < 0.5, up to 1.7ulp error */ |
| 21 | y = 0.5f*log1pf(2*y + 2*y*y/(1-y)); |
| 22 | } |
Szabolcs Nagy | 482ccd2 | 2012-12-11 23:06:20 +0100 | [diff] [blame] | 23 | } else { |
Szabolcs Nagy | f4d9bfb | 2013-09-05 16:57:46 +0000 | [diff] [blame] | 24 | /* avoid overflow */ |
| 25 | y = 0.5f*log1pf(2*(y/(1-y))); |
Szabolcs Nagy | 482ccd2 | 2012-12-11 23:06:20 +0100 | [diff] [blame] | 26 | } |
Szabolcs Nagy | f4d9bfb | 2013-09-05 16:57:46 +0000 | [diff] [blame] | 27 | return s ? -y : y; |
Rich Felker | b69f695 | 2012-03-13 01:17:53 -0400 | [diff] [blame] | 28 | } |