src/math/cosh.c

   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 }