| OLD | NEW |
| 1 #include "libm.h" | 1 #include "libm.h" |
| 2 | 2 |
| 3 #if FLT_EVAL_METHOD==2 | 3 #if FLT_EVAL_METHOD == 2 |
| 4 #undef sqrt | 4 #undef sqrt |
| 5 #define sqrt sqrtl | 5 #define sqrt sqrtl |
| 6 #endif | 6 #endif |
| 7 | 7 |
| 8 /* acosh(x) = log(x + sqrt(x*x-1)) */ | 8 /* acosh(x) = log(x + sqrt(x*x-1)) */ |
| 9 double acosh(double x) | 9 double acosh(double x) { |
| 10 { | 10 union { |
| 11 » union {double f; uint64_t i;} u = {.f = x}; | 11 double f; |
| 12 » unsigned e = u.i >> 52 & 0x7ff; | 12 uint64_t i; |
| 13 } u = {.f = x}; |
| 14 unsigned e = u.i >> 52 & 0x7ff; |
| 13 | 15 |
| 14 » /* x < 1 domain error is handled in the called functions */ | 16 /* x < 1 domain error is handled in the called functions */ |
| 15 | 17 |
| 16 » if (e < 0x3ff + 1) | 18 if (e < 0x3ff + 1) |
| 17 » » /* |x| < 2, up to 2ulp error in [1,1.125] */ | 19 /* |x| < 2, up to 2ulp error in [1,1.125] */ |
| 18 » » return log1p(x-1 + sqrt((x-1)*(x-1)+2*(x-1))); | 20 return log1p(x - 1 + sqrt((x - 1) * (x - 1) + 2 * (x - 1))); |
| 19 » if (e < 0x3ff + 26) | 21 if (e < 0x3ff + 26) |
| 20 » » /* |x| < 0x1p26 */ | 22 /* |x| < 0x1p26 */ |
| 21 » » return log(2*x - 1/(x+sqrt(x*x-1))); | 23 return log(2 * x - 1 / (x + sqrt(x * x - 1))); |
| 22 » /* |x| >= 0x1p26 or nan */ | 24 /* |x| >= 0x1p26 or nan */ |
| 23 » return log(x) + 0.693147180559945309417232121458176568; | 25 return log(x) + 0.693147180559945309417232121458176568; |
| 24 } | 26 } |
| OLD | NEW |