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 */
7double 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