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