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 */
7double 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#ifdef __wasilibc_unmodified_upstream // Wasm doesn't have alternate rounding modes
38 t = __expo2(absx, 2*h);
39#else
40 t = 2*h*__expo2(absx);
41#endif
42 return t;
43}