| 1 | /** |
| 2 | * This file has no copyright assigned and is placed in the Public Domain. |
| 3 | * This file is part of the mingw-w64 runtime package. |
| 4 | * No warranty is given; refer to the file DISCLAIMER.PD within this package. |
| 5 | */ |
| 6 | #include <math.h> |
| 7 | #include <errno.h> |
| 8 | #include "fastmath.h" |
| 9 | |
| 10 | /* atanh (x) = 0.5 * log ((1.0 + x)/(1.0 - x)) */ |
| 11 | long double atanhl (long double x) |
| 12 | { |
| 13 | long double z; |
| 14 | if (isnan (x)) |
| 15 | return x; |
| 16 | z = fabsl (x); |
| 17 | if (z == 1.0L) |
| 18 | { |
| 19 | errno = ERANGE; |
| 20 | return (x > 0 ? INFINITY : -INFINITY); |
| 21 | } |
| 22 | if ( z > 1.0L) |
| 23 | { |
| 24 | errno = EDOM; |
| 25 | return nanl(""); |
| 26 | } |
| 27 | /* Rearrange formula to avoid precision loss for small x. |
| 28 | atanh(x) = 0.5 * log ((1.0 + x)/(1.0 - x)) |
| 29 | 	 = 0.5 * log1p ((1.0 + x)/(1.0 - x) - 1.0) |
| 30 | = 0.5 * log1p ((1.0 + x - 1.0 + x) /(1.0 - x)) |
| 31 | = 0.5 * log1p ((2.0 * x ) / (1.0 - x)) */ |
| 32 | z = 0.5L * __fast_log1pl ((z + z) / (1.0L - z)); |
| 33 | return copysignl(z, x); //ensure 0.0 -> 0.0 and -0.0 -> -0.0. |
| 34 | } |