1*412f47f9SXin Li /*
2*412f47f9SXin Li * Double-precision cosh(x) function.
3*412f47f9SXin Li *
4*412f47f9SXin Li * Copyright (c) 2022-2023, Arm Limited.
5*412f47f9SXin Li * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6*412f47f9SXin Li */
7*412f47f9SXin Li
8*412f47f9SXin Li #include "math_config.h"
9*412f47f9SXin Li #include "pl_sig.h"
10*412f47f9SXin Li #include "pl_test.h"
11*412f47f9SXin Li
12*412f47f9SXin Li #define AbsMask 0x7fffffffffffffff
13*412f47f9SXin Li #define SpecialBound \
14*412f47f9SXin Li 0x40861da04cbafe44 /* 0x1.61da04cbafe44p+9, above which exp overflows. */
15*412f47f9SXin Li
16*412f47f9SXin Li double
17*412f47f9SXin Li __exp_dd (double, double);
18*412f47f9SXin Li
19*412f47f9SXin Li static double
specialcase(double x,uint64_t iax)20*412f47f9SXin Li specialcase (double x, uint64_t iax)
21*412f47f9SXin Li {
22*412f47f9SXin Li if (iax == 0x7ff0000000000000)
23*412f47f9SXin Li return INFINITY;
24*412f47f9SXin Li if (iax > 0x7ff0000000000000)
25*412f47f9SXin Li return __math_invalid (x);
26*412f47f9SXin Li /* exp overflows above SpecialBound. At this magnitude cosh(x) is dominated by
27*412f47f9SXin Li exp(x), so we can approximate cosh(x) by (exp(|x|/2)) ^ 2 / 2. */
28*412f47f9SXin Li double t = __exp_dd (asdouble (iax) / 2, 0);
29*412f47f9SXin Li return (0.5 * t) * t;
30*412f47f9SXin Li }
31*412f47f9SXin Li
32*412f47f9SXin Li /* Approximation for double-precision cosh(x).
33*412f47f9SXin Li cosh(x) = (exp(x) + exp(-x)) / 2.
34*412f47f9SXin Li The greatest observed error is in the special region, 1.93 ULP:
35*412f47f9SXin Li cosh(0x1.628af341989dap+9) got 0x1.fdf28623ef921p+1021
36*412f47f9SXin Li want 0x1.fdf28623ef923p+1021.
37*412f47f9SXin Li
38*412f47f9SXin Li The greatest observed error in the non-special region is 1.03 ULP:
39*412f47f9SXin Li cosh(0x1.502cd8e56ab3bp+0) got 0x1.fe54962842d0ep+0
40*412f47f9SXin Li want 0x1.fe54962842d0fp+0. */
41*412f47f9SXin Li double
cosh(double x)42*412f47f9SXin Li cosh (double x)
43*412f47f9SXin Li {
44*412f47f9SXin Li uint64_t ix = asuint64 (x);
45*412f47f9SXin Li uint64_t iax = ix & AbsMask;
46*412f47f9SXin Li
47*412f47f9SXin Li /* exp overflows a little bit before cosh, so use special-case handler for the
48*412f47f9SXin Li gap, as well as special values. */
49*412f47f9SXin Li if (unlikely (iax >= SpecialBound))
50*412f47f9SXin Li return specialcase (x, iax);
51*412f47f9SXin Li
52*412f47f9SXin Li double ax = asdouble (iax);
53*412f47f9SXin Li /* Use double-precision exp helper to calculate exp(x), then:
54*412f47f9SXin Li cosh(x) = exp(|x|) / 2 + 1 / (exp(|x| * 2). */
55*412f47f9SXin Li double t = __exp_dd (ax, 0);
56*412f47f9SXin Li return 0.5 * t + 0.5 / t;
57*412f47f9SXin Li }
58*412f47f9SXin Li
59*412f47f9SXin Li PL_SIG (S, D, 1, cosh, -10.0, 10.0)
60*412f47f9SXin Li PL_TEST_ULP (cosh, 1.43)
61*412f47f9SXin Li PL_TEST_SYM_INTERVAL (cosh, 0, 0x1.61da04cbafe44p+9, 100000)
62*412f47f9SXin Li PL_TEST_SYM_INTERVAL (cosh, 0x1.61da04cbafe44p+9, 0x1p10, 1000)
63*412f47f9SXin Li PL_TEST_SYM_INTERVAL (cosh, 0x1p10, inf, 100)
64