| OLD | NEW |
| 1 #include "libm.h" | 1 #include "libm.h" |
| 2 | 2 |
| 3 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 | 3 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 |
| 4 long double atanhl(long double x) | 4 long double atanhl(long double x) { |
| 5 { | 5 return atanh(x); |
| 6 » return atanh(x); | |
| 7 } | 6 } |
| 8 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 | 7 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 |
| 9 /* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */ | 8 /* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */ |
| 10 long double atanhl(long double x) | 9 long double atanhl(long double x) { |
| 11 { | 10 union ldshape u = {x}; |
| 12 » union ldshape u = {x}; | 11 unsigned e = u.i.se & 0x7fff; |
| 13 » unsigned e = u.i.se & 0x7fff; | 12 unsigned s = u.i.se >> 15; |
| 14 » unsigned s = u.i.se >> 15; | |
| 15 | 13 |
| 16 » /* |x| */ | 14 /* |x| */ |
| 17 » u.i.se = e; | 15 u.i.se = e; |
| 18 » x = u.f; | 16 x = u.f; |
| 19 | 17 |
| 20 » if (e < 0x3ff - 1) { | 18 if (e < 0x3ff - 1) { |
| 21 » » if (e < 0x3ff - LDBL_MANT_DIG/2) { | 19 if (e < 0x3ff - LDBL_MANT_DIG / 2) { |
| 22 » » » /* handle underflow */ | 20 /* handle underflow */ |
| 23 » » » if (e == 0) | 21 if (e == 0) |
| 24 » » » » FORCE_EVAL((float)x); | 22 FORCE_EVAL((float)x); |
| 25 » » } else { | 23 } else { |
| 26 » » » /* |x| < 0.5, up to 1.7ulp error */ | 24 /* |x| < 0.5, up to 1.7ulp error */ |
| 27 » » » x = 0.5*log1pl(2*x + 2*x*x/(1-x)); | 25 x = 0.5 * log1pl(2 * x + 2 * x * x / (1 - x)); |
| 28 » » } | 26 } |
| 29 » } else { | 27 } else { |
| 30 » » /* avoid overflow */ | 28 /* avoid overflow */ |
| 31 » » x = 0.5*log1pl(2*(x/(1-x))); | 29 x = 0.5 * log1pl(2 * (x / (1 - x))); |
| 32 » } | 30 } |
| 33 » return s ? -x : x; | 31 return s ? -x : x; |
| 34 } | 32 } |
| 35 #endif | 33 #endif |
| OLD | NEW |