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 } 41