OLD | NEW |
(Empty) | |
| 1 #include "libm.h" |
| 2 |
| 3 /* cosh(x) = (exp(x) + 1/exp(x))/2 |
| 4 * = 1 + 0.5*(exp(x)-1)*(exp(x)-1)/exp(x) |
| 5 * = 1 + x*x/2 + o(x^4) |
| 6 */ |
| 7 double cosh(double x) |
| 8 { |
| 9 union {double f; uint64_t i;} u = {.f = x}; |
| 10 uint32_t w; |
| 11 double t; |
| 12 |
| 13 /* |x| */ |
| 14 u.i &= (uint64_t)-1/2; |
| 15 x = u.f; |
| 16 w = u.i >> 32; |
| 17 |
| 18 /* |x| < log(2) */ |
| 19 if (w < 0x3fe62e42) { |
| 20 if (w < 0x3ff00000 - (26<<20)) { |
| 21 /* raise inexact if x!=0 */ |
| 22 FORCE_EVAL(x + 0x1p120f); |
| 23 return 1; |
| 24 } |
| 25 t = expm1(x); |
| 26 return 1 + t*t/(2*(1+t)); |
| 27 } |
| 28 |
| 29 /* |x| < log(DBL_MAX) */ |
| 30 if (w < 0x40862e42) { |
| 31 t = exp(x); |
| 32 /* note: if x>log(0x1p26) then the 1/t is not needed */ |
| 33 return 0.5*(t + 1/t); |
| 34 } |
| 35 |
| 36 /* |x| > log(DBL_MAX) or nan */ |
| 37 /* note: the result is stored to handle overflow */ |
| 38 t = __expo2(x); |
| 39 return t; |
| 40 } |
OLD | NEW |