1 /* 2 * Double-precision vector asinh(x) function. 3 * 4 * Copyright (c) 2022-2024, Arm Limited. 5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 6 */ 7 8 #include "test_defs.h" 9 #include "test_sig.h" 10 #include "v_math.h" 11 12 const static struct data 13 { 14 uint64x2_t huge_bound, abs_mask, off, mask; 15 #if WANT_SIMD_EXCEPT 16 float64x2_t tiny_bound; 17 #endif 18 float64x2_t lc0, lc2; 19 double lc1, lc3, ln2, lc4; 20 21 float64x2_t c0, c2, c4, c6, c8, c10, c12, c14, c16, c17; 22 double c1, c3, c5, c7, c9, c11, c13, c15; 23 24 } data = { 25 26 #if WANT_SIMD_EXCEPT 27 .tiny_bound = V2 (0x1p-26), 28 #endif 29 /* Even terms of polynomial s.t. asinh(x) is approximated by 30 asinh(x) ~= x + x^3 * (C0 + C1 * x + C2 * x^2 + C3 * x^3 + ...). 31 Generated using Remez, f = (asinh(sqrt(x)) - sqrt(x))/x^(3/2). */ 32 33 .c0 = V2 (-0x1.55555555554a7p-3), 34 .c1 = 0x1.3333333326c7p-4, 35 .c2 = V2 (-0x1.6db6db68332e6p-5), 36 .c3 = 0x1.f1c71b26fb40dp-6, 37 .c4 = V2 (-0x1.6e8b8b654a621p-6), 38 .c5 = 0x1.1c4daa9e67871p-6, 39 .c6 = V2 (-0x1.c9871d10885afp-7), 40 .c7 = 0x1.7a16e8d9d2ecfp-7, 41 .c8 = V2 (-0x1.3ddca533e9f54p-7), 42 .c9 = 0x1.0becef748dafcp-7, 43 .c10 = V2 (-0x1.b90c7099dd397p-8), 44 .c11 = 0x1.541f2bb1ffe51p-8, 45 .c12 = V2 (-0x1.d217026a669ecp-9), 46 .c13 = 0x1.0b5c7977aaf7p-9, 47 .c14 = V2 (-0x1.e0f37daef9127p-11), 48 .c15 = 0x1.388b5fe542a6p-12, 49 .c16 = V2 (-0x1.021a48685e287p-14), 50 .c17 = V2 (0x1.93d4ba83d34dap-18), 51 52 .lc0 = V2 (-0x1.ffffffffffff7p-2), 53 .lc1 = 0x1.55555555170d4p-2, 54 .lc2 = V2 (-0x1.0000000399c27p-2), 55 .lc3 = 0x1.999b2e90e94cap-3, 56 .lc4 = -0x1.554e550bd501ep-3, 57 .ln2 = 0x1.62e42fefa39efp-1, 58 59 .off = V2 (0x3fe6900900000000), 60 .huge_bound = V2 (0x5fe0000000000000), 61 .abs_mask = V2 (0x7fffffffffffffff), 62 .mask = V2 (0xfffULL << 52), 63 }; 64 65 static float64x2_t NOINLINE VPCS_ATTR 66 special_case (float64x2_t x, float64x2_t y, uint64x2_t abs_mask, 67 uint64x2_t special) 68 { 69 /* Copy sign. */ 70 y = vbslq_f64 (abs_mask, y, x); 71 return v_call_f64 (asinh, x, y, special); 72 } 73 74 #define N (1 << V_LOG_TABLE_BITS) 75 #define IndexMask (N - 1) 76 77 struct entry 78 { 79 float64x2_t invc; 80 float64x2_t logc; 81 }; 82 83 static inline struct entry 84 lookup (uint64x2_t i) 85 { 86 /* Since N is a power of 2, n % N = n & (N - 1). */ 87 struct entry e; 88 uint64_t i0 = (vgetq_lane_u64 (i, 0) >> (52 - V_LOG_TABLE_BITS)) & IndexMask; 89 uint64_t i1 = (vgetq_lane_u64 (i, 1) >> (52 - V_LOG_TABLE_BITS)) & IndexMask; 90 float64x2_t e0 = vld1q_f64 (&__v_log_data.table[i0].invc); 91 float64x2_t e1 = vld1q_f64 (&__v_log_data.table[i1].invc); 92 e.invc = vuzp1q_f64 (e0, e1); 93 e.logc = vuzp2q_f64 (e0, e1); 94 return e; 95 } 96 97 static inline float64x2_t 98 log_inline (float64x2_t xm, const struct data *d) 99 { 100 101 uint64x2_t u = vreinterpretq_u64_f64 (xm); 102 uint64x2_t u_off = vsubq_u64 (u, d->off); 103 104 int64x2_t k = vshrq_n_s64 (vreinterpretq_s64_u64 (u_off), 52); 105 uint64x2_t iz = vsubq_u64 (u, vandq_u64 (u_off, d->mask)); 106 float64x2_t z = vreinterpretq_f64_u64 (iz); 107 108 struct entry e = lookup (u_off); 109 110 /* log(x) = log1p(z/c-1) + log(c) + k*Ln2. */ 111 float64x2_t r = vfmaq_f64 (v_f64 (-1.0), z, e.invc); 112 float64x2_t kd = vcvtq_f64_s64 (k); 113 114 /* hi = r + log(c) + k*Ln2. */ 115 float64x2_t ln2_and_lc4 = vld1q_f64 (&d->ln2); 116 float64x2_t hi = vfmaq_laneq_f64 (vaddq_f64 (e.logc, r), kd, ln2_and_lc4, 0); 117 118 /* y = r2*(A0 + r*A1 + r2*(A2 + r*A3 + r2*A4)) + hi. */ 119 float64x2_t odd_coeffs = vld1q_f64 (&d->lc1); 120 float64x2_t r2 = vmulq_f64 (r, r); 121 float64x2_t y = vfmaq_laneq_f64 (d->lc2, r, odd_coeffs, 1); 122 float64x2_t p = vfmaq_laneq_f64 (d->lc0, r, odd_coeffs, 0); 123 y = vfmaq_laneq_f64 (y, r2, ln2_and_lc4, 1); 124 y = vfmaq_f64 (p, r2, y); 125 return vfmaq_f64 (hi, y, r2); 126 } 127 128 /* Double-precision implementation of vector asinh(x). 129 asinh is very sensitive around 1, so it is impractical to devise a single 130 low-cost algorithm which is sufficiently accurate on a wide range of input. 131 Instead we use two different algorithms: 132 asinh(x) = sign(x) * log(|x| + sqrt(x^2 + 1) if |x| >= 1 133 = sign(x) * (|x| + |x|^3 * P(x^2)) otherwise 134 where log(x) is an optimized log approximation, and P(x) is a polynomial 135 shared with the scalar routine. The greatest observed error 2.79 ULP, in 136 |x| >= 1: 137 _ZGVnN2v_asinh(0x1.2cd9d73ea76a6p+0) got 0x1.ffffd003219dap-1 138 want 0x1.ffffd003219ddp-1. */ 139 VPCS_ATTR float64x2_t V_NAME_D1 (asinh) (float64x2_t x) 140 { 141 const struct data *d = ptr_barrier (&data); 142 float64x2_t ax = vabsq_f64 (x); 143 144 uint64x2_t gt1 = vcgeq_f64 (ax, v_f64 (1)); 145 146 #if WANT_SIMD_EXCEPT 147 uint64x2_t iax = vreinterpretq_u64_f64 (ax); 148 uint64x2_t special = vcgeq_u64 (iax, (d->huge_bound)); 149 uint64x2_t tiny = vcltq_f64 (ax, d->tiny_bound); 150 special = vorrq_u64 (special, tiny); 151 #else 152 uint64x2_t special = vcgeq_f64 (ax, vreinterpretq_f64_u64 (d->huge_bound)); 153 #endif 154 155 /* Option 1: |x| >= 1. 156 Compute asinh(x) according by asinh(x) = log(x + sqrt(x^2 + 1)). 157 If WANT_SIMD_EXCEPT is enabled, sidestep special values, which will 158 overflow, by setting special lanes to 1. These will be fixed later. */ 159 float64x2_t option_1 = v_f64 (0); 160 if (likely (v_any_u64 (gt1))) 161 { 162 #if WANT_SIMD_EXCEPT 163 float64x2_t xm = v_zerofy_f64 (ax, special); 164 #else 165 float64x2_t xm = ax; 166 #endif 167 option_1 = log_inline ( 168 vaddq_f64 (xm, vsqrtq_f64 (vfmaq_f64 (v_f64 (1), xm, xm))), d); 169 } 170 171 /* Option 2: |x| < 1. 172 Compute asinh(x) using a polynomial. 173 If WANT_SIMD_EXCEPT is enabled, sidestep special lanes, which will 174 overflow, and tiny lanes, which will underflow, by setting them to 0. They 175 will be fixed later, either by selecting x or falling back to the scalar 176 special-case. The largest observed error in this region is 1.47 ULPs: 177 _ZGVnN2v_asinh(0x1.fdfcd00cc1e6ap-1) got 0x1.c1d6bf874019bp-1 178 want 0x1.c1d6bf874019cp-1. */ 179 float64x2_t option_2 = v_f64 (0); 180 181 if (likely (v_any_u64 (vceqzq_u64 (gt1)))) 182 { 183 184 #if WANT_SIMD_EXCEPT 185 ax = v_zerofy_f64 (ax, vorrq_u64 (tiny, gt1)); 186 #endif 187 float64x2_t x2 = vmulq_f64 (ax, ax), z2 = vmulq_f64 (x2, x2); 188 /* Order-17 Pairwise Horner scheme. */ 189 float64x2_t c13 = vld1q_f64 (&d->c1); 190 float64x2_t c57 = vld1q_f64 (&d->c5); 191 float64x2_t c911 = vld1q_f64 (&d->c9); 192 float64x2_t c1315 = vld1q_f64 (&d->c13); 193 194 float64x2_t p01 = vfmaq_laneq_f64 (d->c0, x2, c13, 0); 195 float64x2_t p23 = vfmaq_laneq_f64 (d->c2, x2, c13, 1); 196 float64x2_t p45 = vfmaq_laneq_f64 (d->c4, x2, c57, 0); 197 float64x2_t p67 = vfmaq_laneq_f64 (d->c6, x2, c57, 1); 198 float64x2_t p89 = vfmaq_laneq_f64 (d->c8, x2, c911, 0); 199 float64x2_t p1011 = vfmaq_laneq_f64 (d->c10, x2, c911, 1); 200 float64x2_t p1213 = vfmaq_laneq_f64 (d->c12, x2, c1315, 0); 201 float64x2_t p1415 = vfmaq_laneq_f64 (d->c14, x2, c1315, 1); 202 float64x2_t p1617 = vfmaq_f64 (d->c16, x2, d->c17); 203 204 float64x2_t p = vfmaq_f64 (p1415, z2, p1617); 205 p = vfmaq_f64 (p1213, z2, p); 206 p = vfmaq_f64 (p1011, z2, p); 207 p = vfmaq_f64 (p89, z2, p); 208 209 p = vfmaq_f64 (p67, z2, p); 210 p = vfmaq_f64 (p45, z2, p); 211 212 p = vfmaq_f64 (p23, z2, p); 213 214 p = vfmaq_f64 (p01, z2, p); 215 option_2 = vfmaq_f64 (ax, p, vmulq_f64 (ax, x2)); 216 #if WANT_SIMD_EXCEPT 217 option_2 = vbslq_f64 (tiny, x, option_2); 218 #endif 219 } 220 221 /* Choose the right option for each lane. */ 222 float64x2_t y = vbslq_f64 (gt1, option_1, option_2); 223 if (unlikely (v_any_u64 (special))) 224 { 225 return special_case (x, y, d->abs_mask, special); 226 } 227 /* Copy sign. */ 228 return vbslq_f64 (d->abs_mask, y, x); 229 } 230 231 TEST_SIG (V, D, 1, asinh, -10.0, 10.0) 232 TEST_ULP (V_NAME_D1 (asinh), 2.29) 233 TEST_DISABLE_FENV_IF_NOT (V_NAME_D1 (asinh), WANT_SIMD_EXCEPT) 234 TEST_SYM_INTERVAL (V_NAME_D1 (asinh), 0, 0x1p-26, 50000) 235 TEST_SYM_INTERVAL (V_NAME_D1 (asinh), 0x1p-26, 1, 50000) 236 TEST_SYM_INTERVAL (V_NAME_D1 (asinh), 1, 0x1p511, 50000) 237 TEST_SYM_INTERVAL (V_NAME_D1 (asinh), 0x1p511, inf, 40000) 238 /* Test vector asinh 3 times, with control lane < 1, > 1 and special. 239 Ensures the v_sel is choosing the right option in all cases. */ 240 TEST_CONTROL_VALUE (V_NAME_D1 (asinh), 0.5) 241 TEST_CONTROL_VALUE (V_NAME_D1 (asinh), 2) 242 TEST_CONTROL_VALUE (V_NAME_D1 (asinh), 0x1p600) 243