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