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