src/math/sinh.c

   1 #include "libm.h"
   2
   3 /* sinh(x) = (exp(x) - 1/exp(x))/2
   4  *         = (exp(x)-1 + (exp(x)-1)/exp(x))/2
   5  *         = x + x^3/6 + o(x^5)
   6  */
   7 double sinh(double x)
   8 {
   9         union {double f; uint64_t i;} u = {.f = x};
  10         uint32_t w;
  11         double t, h, absx;
  12
  13         h = 0.5;
  14         if (u.i >> 63)
  15                 h = -h;
  16         /* |x| */
  17         u.i &= (uint64_t)-1/2;
  18         absx = u.f;
  19         w = u.i >> 32;
  20
  21         /* |x| < log(DBL_MAX) */
  22         if (w < 0x40862e42) {
  23                 t = expm1(absx);
  24                 if (w < 0x3ff00000) {
  25                         if (w < 0x3ff00000 - (26<<20))
  26                                 /* note: inexact is raised by expm1 */
  27                                 /* note: this branch avoids underflow */
  28                                 return x;
  29                         return h*(2*t - t*t/(t+1));
  30                 }
  31                 /* note: |x|>log(0x1p26)+eps could be just h*exp(x) */
  32                 return h*(t + t/(t+1));
  33         }
  34
  35         /* |x| > log(DBL_MAX) or nan */
  36         /* note: the result is stored to handle overflow */
  37         t = 2*h*__expo2(absx);
  38         return t;
  39 }