cosh.c raw

   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, 1.0);
  39  	return t;
  40  }
  41