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 |