1 #include "libm.h"
2 
3 /* cosh(x) = (exp(x) + 1/exp(x))/2
4  *         = 1 + 0.5*(exp(x)-1)*(exp(x)-1)/exp(x)
5  *         = 1 + x*x/2 + o(x^4)
6  */
cosh(double x)7 double cosh(double x) {
8     union {
9         double f;
10         uint64_t i;
11     } u = {.f = x};
12     uint32_t w;
13     double t;
14 
15     /* |x| */
16     u.i &= (uint64_t)-1 / 2;
17     x = u.f;
18     w = u.i >> 32;
19 
20     /* |x| < log(2) */
21     if (w < 0x3fe62e42) {
22         if (w < 0x3ff00000 - (26 << 20)) {
23             /* raise inexact if x!=0 */
24             FORCE_EVAL(x + 0x1p120f);
25             return 1;
26         }
27         t = expm1(x);
28         return 1 + t * t / (2 * (1 + t));
29     }
30 
31     /* |x| < log(DBL_MAX) */
32     if (w < 0x40862e42) {
33         t = exp(x);
34         /* note: if x>log(0x1p26) then the 1/t is not needed */
35         return 0.5 * (t + 1 / t);
36     }
37 
38     /* |x| > log(DBL_MAX) or nan */
39     /* note: the result is stored to handle overflow */
40     t = __expo2(x);
41     return t;
42 }
43