sinh.c raw

   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 and underflow are raised by expm1 */
  27  				/* note: this branch avoids spurious 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 = __expo2(absx, 2*h);
  38  	return t;
  39  }
  40