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