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 = 2*h*__expo2(absx); 38 return t; 39 } 40