1
2 /* @(#)e_sinh.c 5.1 93/09/24 */
3 /*
4 * ====================================================
5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 *
7 * Developed at SunPro, a Sun Microsystems, Inc. business.
8 * Permission to use, copy, modify, and distribute this
9 * software is freely granted, provided that this notice
10 * is preserved.
11 * ====================================================
12 */
13
14 /* sinh(x)
15 * Method :
16 * mathematically sinh(x) if defined to be (exp(x)-exp(-x))/2
17 * 1. Replace x by |x| (sinh(-x) = -sinh(x)).
18 * 2.
19 * E + E/(E+1)
20 * 0 <= x <= 22 : sinh(x) := --------------, E=expm1(x)
21 * 2
22 *
23 * 22 <= x <= lnovft : sinh(x) := exp(x)/2
24 * lnovft <= x <= ln2ovft: sinh(x) := exp(x/2)/2 * exp(x/2)
25 * ln2ovft < x : sinh(x) := x*shuge (overflow)
26 *
27 * Special cases:
28 * sinh(x) is |x| if x is +INF, -INF, or NaN.
29 * only sinh(0)=0 is exact for finite x.
30 */
31
32 #include "fdlibm.h"
33
34 #ifdef _NEED_FLOAT64
35
36 static const __float64 one = _F_64(1.0), shuge = _F_64(1.0e307);
37
38 __float64
sinh64(__float64 x)39 sinh64(__float64 x)
40 {
41 __float64 t, w, h;
42 __int32_t ix, jx;
43 __uint32_t lx;
44
45 /* High word of |x|. */
46 GET_HIGH_WORD(jx, x);
47 ix = jx & 0x7fffffff;
48
49 /* x is INF or NaN */
50 if (ix >= 0x7ff00000)
51 return x + x;
52
53 h = _F_64(0.5);
54 if (jx < 0)
55 h = -h;
56 /* |x| in [0,22], return sign(x)*0.5*(E+E/(E+1))) */
57 if (ix < 0x40360000) { /* |x|<22 */
58 if (ix < 0x3e300000) /* |x|<2**-28 */
59 if (shuge + x > one)
60 return x; /* sinh(tiny) = tiny with inexact */
61 t = expm164(fabs64(x));
62 if (ix < 0x3ff00000)
63 return h * (_F_64(2.0) * t - t * t / (t + one));
64 return h * (t + t / (t + one));
65 }
66
67 /* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
68 if (ix < 0x40862E42)
69 return h * exp64(fabs64(x));
70
71 /* |x| in [log(maxdouble), overflowthresold] */
72 GET_LOW_WORD(lx, x);
73 if (ix < 0x408633CE || (ix == 0x408633ce && lx <= (__uint32_t)0x8fb9f87d)) {
74 w = exp64(_F_64(0.5) * fabs64(x));
75 t = h * w;
76 return t * w;
77 }
78
79 /* |x| > overflowthresold, sinh(x) overflow */
80 return __math_oflow(x < 0);
81 }
82
83 _MATH_ALIAS_d_d(sinh)
84
85 #endif /* _NEED_FLOAT64 */
86
