1*f3087befSAndrew Turner /*
2*f3087befSAndrew Turner * Double-precision SVE cosh(x) function.
3*f3087befSAndrew Turner *
4*f3087befSAndrew Turner * Copyright (c) 2023-2025, 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 "sv_math.h"
9*f3087befSAndrew Turner #include "test_sig.h"
10*f3087befSAndrew Turner #include "test_defs.h"
11*f3087befSAndrew Turner
12*f3087befSAndrew Turner static const struct data
13*f3087befSAndrew Turner {
14*f3087befSAndrew Turner float64_t poly[3];
15*f3087befSAndrew Turner float64_t inv_ln2, ln2_hi, ln2_lo, shift, thres;
16*f3087befSAndrew Turner uint64_t special_bound;
17*f3087befSAndrew Turner } data = {
18*f3087befSAndrew Turner .poly = { 0x1.fffffffffffd4p-2, 0x1.5555571d6b68cp-3,
19*f3087befSAndrew Turner 0x1.5555576a59599p-5, },
20*f3087befSAndrew Turner
21*f3087befSAndrew Turner .inv_ln2 = 0x1.71547652b82fep8, /* N/ln2. */
22*f3087befSAndrew Turner /* -ln2/N. */
23*f3087befSAndrew Turner .ln2_hi = -0x1.62e42fefa39efp-9,
24*f3087befSAndrew Turner .ln2_lo = -0x1.abc9e3b39803f3p-64,
25*f3087befSAndrew Turner .shift = 0x1.8p+52,
26*f3087befSAndrew Turner .thres = 704.0,
27*f3087befSAndrew Turner
28*f3087befSAndrew Turner /* 0x1.6p9, above which exp overflows. */
29*f3087befSAndrew Turner .special_bound = 0x4086000000000000,
30*f3087befSAndrew Turner };
31*f3087befSAndrew Turner
32*f3087befSAndrew Turner static svfloat64_t NOINLINE
special_case(svfloat64_t x,svbool_t pg,svfloat64_t t,svbool_t special)33*f3087befSAndrew Turner special_case (svfloat64_t x, svbool_t pg, svfloat64_t t, svbool_t special)
34*f3087befSAndrew Turner {
35*f3087befSAndrew Turner svfloat64_t half_t = svmul_x (svptrue_b64 (), t, 0.5);
36*f3087befSAndrew Turner svfloat64_t half_over_t = svdivr_x (pg, t, 0.5);
37*f3087befSAndrew Turner svfloat64_t y = svadd_x (pg, half_t, half_over_t);
38*f3087befSAndrew Turner return sv_call_f64 (cosh, x, y, special);
39*f3087befSAndrew Turner }
40*f3087befSAndrew Turner
41*f3087befSAndrew Turner /* Helper for approximating exp(x). Copied from sv_exp_tail, with no
42*f3087befSAndrew Turner special-case handling or tail. */
43*f3087befSAndrew Turner static inline svfloat64_t
exp_inline(svfloat64_t x,const svbool_t pg,const struct data * d)44*f3087befSAndrew Turner exp_inline (svfloat64_t x, const svbool_t pg, const struct data *d)
45*f3087befSAndrew Turner {
46*f3087befSAndrew Turner /* Calculate exp(x). */
47*f3087befSAndrew Turner svfloat64_t z = svmla_x (pg, sv_f64 (d->shift), x, d->inv_ln2);
48*f3087befSAndrew Turner svfloat64_t n = svsub_x (pg, z, d->shift);
49*f3087befSAndrew Turner
50*f3087befSAndrew Turner svfloat64_t r = svmla_x (pg, x, n, d->ln2_hi);
51*f3087befSAndrew Turner r = svmla_x (pg, r, n, d->ln2_lo);
52*f3087befSAndrew Turner
53*f3087befSAndrew Turner svuint64_t u = svreinterpret_u64 (z);
54*f3087befSAndrew Turner svuint64_t e = svlsl_x (pg, u, 52 - V_EXP_TAIL_TABLE_BITS);
55*f3087befSAndrew Turner svuint64_t i = svand_x (svptrue_b64 (), u, 0xff);
56*f3087befSAndrew Turner
57*f3087befSAndrew Turner svfloat64_t y = svmla_x (pg, sv_f64 (d->poly[1]), r, d->poly[2]);
58*f3087befSAndrew Turner y = svmla_x (pg, sv_f64 (d->poly[0]), r, y);
59*f3087befSAndrew Turner y = svmla_x (pg, sv_f64 (1.0), r, y);
60*f3087befSAndrew Turner y = svmul_x (svptrue_b64 (), r, y);
61*f3087befSAndrew Turner
62*f3087befSAndrew Turner /* s = 2^(n/N). */
63*f3087befSAndrew Turner u = svld1_gather_index (pg, __v_exp_tail_data, i);
64*f3087befSAndrew Turner svfloat64_t s = svreinterpret_f64 (svadd_x (pg, u, e));
65*f3087befSAndrew Turner
66*f3087befSAndrew Turner return svmla_x (pg, s, s, y);
67*f3087befSAndrew Turner }
68*f3087befSAndrew Turner
69*f3087befSAndrew Turner /* Approximation for SVE double-precision cosh(x) using exp_inline.
70*f3087befSAndrew Turner cosh(x) = (exp(x) + exp(-x)) / 2.
71*f3087befSAndrew Turner The greatest observed error is in the scalar fall-back region, so is the
72*f3087befSAndrew Turner same as the scalar routine, 1.93 ULP:
73*f3087befSAndrew Turner _ZGVsMxv_cosh (0x1.628ad45039d2fp+9) got 0x1.fd774e958236dp+1021
74*f3087befSAndrew Turner want 0x1.fd774e958236fp+1021.
75*f3087befSAndrew Turner
76*f3087befSAndrew Turner The greatest observed error in the non-special region is 1.54 ULP:
77*f3087befSAndrew Turner _ZGVsMxv_cosh (0x1.ba5651dd4486bp+2) got 0x1.f5e2bb8d5c98fp+8
78*f3087befSAndrew Turner want 0x1.f5e2bb8d5c991p+8. */
SV_NAME_D1(cosh)79*f3087befSAndrew Turner svfloat64_t SV_NAME_D1 (cosh) (svfloat64_t x, const svbool_t pg)
80*f3087befSAndrew Turner {
81*f3087befSAndrew Turner const struct data *d = ptr_barrier (&data);
82*f3087befSAndrew Turner
83*f3087befSAndrew Turner svfloat64_t ax = svabs_x (pg, x);
84*f3087befSAndrew Turner svbool_t special = svcmpgt (pg, svreinterpret_u64 (ax), d->special_bound);
85*f3087befSAndrew Turner
86*f3087befSAndrew Turner /* Up to the point that exp overflows, we can use it to calculate cosh by
87*f3087befSAndrew Turner exp(|x|) / 2 + 1 / (2 * exp(|x|)). */
88*f3087befSAndrew Turner svfloat64_t t = exp_inline (ax, pg, d);
89*f3087befSAndrew Turner
90*f3087befSAndrew Turner /* Fall back to scalar for any special cases. */
91*f3087befSAndrew Turner if (unlikely (svptest_any (pg, special)))
92*f3087befSAndrew Turner return special_case (x, pg, t, special);
93*f3087befSAndrew Turner
94*f3087befSAndrew Turner svfloat64_t half_t = svmul_x (svptrue_b64 (), t, 0.5);
95*f3087befSAndrew Turner svfloat64_t half_over_t = svdivr_x (pg, t, 0.5);
96*f3087befSAndrew Turner return svadd_x (pg, half_t, half_over_t);
97*f3087befSAndrew Turner }
98*f3087befSAndrew Turner
99*f3087befSAndrew Turner TEST_SIG (SV, D, 1, cosh, -10.0, 10.0)
100*f3087befSAndrew Turner TEST_ULP (SV_NAME_D1 (cosh), 1.43)
101*f3087befSAndrew Turner TEST_DISABLE_FENV (SV_NAME_D1 (cosh))
102*f3087befSAndrew Turner TEST_SYM_INTERVAL (SV_NAME_D1 (cosh), 0, 0x1.6p9, 100000)
103*f3087befSAndrew Turner TEST_SYM_INTERVAL (SV_NAME_D1 (cosh), 0x1.6p9, inf, 1000)
104*f3087befSAndrew Turner CLOSE_SVE_ATTR