1 /* 2 * Double-precision vector asinh(x) function. 3 * 4 * Copyright (c) 2022-2023, Arm Limited. 5 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 6 */ 7 8 #include "v_math.h" 9 #include "poly_advsimd_f64.h" 10 #include "pl_sig.h" 11 #include "pl_test.h" 12 13 #define A(i) v_f64 (__v_log_data.poly[i]) 14 #define N (1 << V_LOG_TABLE_BITS) 15 16 const static struct data 17 { 18 float64x2_t poly[18]; 19 uint64x2_t off, huge_bound, abs_mask; 20 float64x2_t ln2, tiny_bound; 21 } data = { 22 .off = V2 (0x3fe6900900000000), 23 .ln2 = V2 (0x1.62e42fefa39efp-1), 24 .huge_bound = V2 (0x5fe0000000000000), 25 .tiny_bound = V2 (0x1p-26), 26 .abs_mask = V2 (0x7fffffffffffffff), 27 /* Even terms of polynomial s.t. asinh(x) is approximated by 28 asinh(x) ~= x + x^3 * (C0 + C1 * x + C2 * x^2 + C3 * x^3 + ...). 29 Generated using Remez, f = (asinh(sqrt(x)) - sqrt(x))/x^(3/2). */ 30 .poly = { V2 (-0x1.55555555554a7p-3), V2 (0x1.3333333326c7p-4), 31 V2 (-0x1.6db6db68332e6p-5), V2 (0x1.f1c71b26fb40dp-6), 32 V2 (-0x1.6e8b8b654a621p-6), V2 (0x1.1c4daa9e67871p-6), 33 V2 (-0x1.c9871d10885afp-7), V2 (0x1.7a16e8d9d2ecfp-7), 34 V2 (-0x1.3ddca533e9f54p-7), V2 (0x1.0becef748dafcp-7), 35 V2 (-0x1.b90c7099dd397p-8), V2 (0x1.541f2bb1ffe51p-8), 36 V2 (-0x1.d217026a669ecp-9), V2 (0x1.0b5c7977aaf7p-9), 37 V2 (-0x1.e0f37daef9127p-11), V2 (0x1.388b5fe542a6p-12), 38 V2 (-0x1.021a48685e287p-14), V2 (0x1.93d4ba83d34dap-18) }, 39 }; 40 41 static float64x2_t NOINLINE VPCS_ATTR 42 special_case (float64x2_t x, float64x2_t y, uint64x2_t special) 43 { 44 return v_call_f64 (asinh, x, y, special); 45 } 46 47 struct entry 48 { 49 float64x2_t invc; 50 float64x2_t logc; 51 }; 52 53 static inline struct entry 54 lookup (uint64x2_t i) 55 { 56 float64x2_t e0 = vld1q_f64 ( 57 &__v_log_data.table[(i[0] >> (52 - V_LOG_TABLE_BITS)) & (N - 1)].invc); 58 float64x2_t e1 = vld1q_f64 ( 59 &__v_log_data.table[(i[1] >> (52 - V_LOG_TABLE_BITS)) & (N - 1)].invc); 60 return (struct entry){ vuzp1q_f64 (e0, e1), vuzp2q_f64 (e0, e1) }; 61 } 62 63 static inline float64x2_t 64 log_inline (float64x2_t x, const struct data *d) 65 { 66 /* Double-precision vector log, copied from ordinary vector log with some 67 cosmetic modification and special-cases removed. */ 68 uint64x2_t ix = vreinterpretq_u64_f64 (x); 69 uint64x2_t tmp = vsubq_u64 (ix, d->off); 70 int64x2_t k = vshrq_n_s64 (vreinterpretq_s64_u64 (tmp), 52); 71 uint64x2_t iz 72 = vsubq_u64 (ix, vandq_u64 (tmp, vdupq_n_u64 (0xfffULL << 52))); 73 float64x2_t z = vreinterpretq_f64_u64 (iz); 74 struct entry e = lookup (tmp); 75 float64x2_t r = vfmaq_f64 (v_f64 (-1.0), z, e.invc); 76 float64x2_t kd = vcvtq_f64_s64 (k); 77 float64x2_t hi = vfmaq_f64 (vaddq_f64 (e.logc, r), kd, d->ln2); 78 float64x2_t r2 = vmulq_f64 (r, r); 79 float64x2_t y = vfmaq_f64 (A (2), A (3), r); 80 float64x2_t p = vfmaq_f64 (A (0), A (1), r); 81 y = vfmaq_f64 (y, A (4), r2); 82 y = vfmaq_f64 (p, y, r2); 83 y = vfmaq_f64 (hi, y, r2); 84 return y; 85 } 86 87 /* Double-precision implementation of vector asinh(x). 88 asinh is very sensitive around 1, so it is impractical to devise a single 89 low-cost algorithm which is sufficiently accurate on a wide range of input. 90 Instead we use two different algorithms: 91 asinh(x) = sign(x) * log(|x| + sqrt(x^2 + 1) if |x| >= 1 92 = sign(x) * (|x| + |x|^3 * P(x^2)) otherwise 93 where log(x) is an optimized log approximation, and P(x) is a polynomial 94 shared with the scalar routine. The greatest observed error 3.29 ULP, in 95 |x| >= 1: 96 __v_asinh(0x1.2cd9d717e2c9bp+0) got 0x1.ffffcfd0e234fp-1 97 want 0x1.ffffcfd0e2352p-1. */ 98 VPCS_ATTR float64x2_t V_NAME_D1 (asinh) (float64x2_t x) 99 { 100 const struct data *d = ptr_barrier (&data); 101 102 float64x2_t ax = vabsq_f64 (x); 103 uint64x2_t iax = vreinterpretq_u64_f64 (ax); 104 105 uint64x2_t gt1 = vcgeq_f64 (ax, v_f64 (1)); 106 uint64x2_t special = vcgeq_u64 (iax, d->huge_bound); 107 108 #if WANT_SIMD_EXCEPT 109 uint64x2_t tiny = vcltq_f64 (ax, d->tiny_bound); 110 special = vorrq_u64 (special, tiny); 111 #endif 112 113 /* Option 1: |x| >= 1. 114 Compute asinh(x) according by asinh(x) = log(x + sqrt(x^2 + 1)). 115 If WANT_SIMD_EXCEPT is enabled, sidestep special values, which will 116 overflow, by setting special lanes to 1. These will be fixed later. */ 117 float64x2_t option_1 = v_f64 (0); 118 if (likely (v_any_u64 (gt1))) 119 { 120 #if WANT_SIMD_EXCEPT 121 float64x2_t xm = v_zerofy_f64 (ax, special); 122 #else 123 float64x2_t xm = ax; 124 #endif 125 option_1 = log_inline ( 126 vaddq_f64 (xm, vsqrtq_f64 (vfmaq_f64 (v_f64 (1), xm, xm))), d); 127 } 128 129 /* Option 2: |x| < 1. 130 Compute asinh(x) using a polynomial. 131 If WANT_SIMD_EXCEPT is enabled, sidestep special lanes, which will 132 overflow, and tiny lanes, which will underflow, by setting them to 0. They 133 will be fixed later, either by selecting x or falling back to the scalar 134 special-case. The largest observed error in this region is 1.47 ULPs: 135 __v_asinh(0x1.fdfcd00cc1e6ap-1) got 0x1.c1d6bf874019bp-1 136 want 0x1.c1d6bf874019cp-1. */ 137 float64x2_t option_2 = v_f64 (0); 138 if (likely (v_any_u64 (vceqzq_u64 (gt1)))) 139 { 140 #if WANT_SIMD_EXCEPT 141 ax = v_zerofy_f64 (ax, vorrq_u64 (tiny, gt1)); 142 #endif 143 float64x2_t x2 = vmulq_f64 (ax, ax), x3 = vmulq_f64 (ax, x2), 144 z2 = vmulq_f64 (x2, x2), z4 = vmulq_f64 (z2, z2), 145 z8 = vmulq_f64 (z4, z4), z16 = vmulq_f64 (z8, z8); 146 float64x2_t p = v_estrin_17_f64 (x2, z2, z4, z8, z16, d->poly); 147 option_2 = vfmaq_f64 (ax, p, x3); 148 #if WANT_SIMD_EXCEPT 149 option_2 = vbslq_f64 (tiny, x, option_2); 150 #endif 151 } 152 153 /* Choose the right option for each lane. */ 154 float64x2_t y = vbslq_f64 (gt1, option_1, option_2); 155 /* Copy sign. */ 156 y = vbslq_f64 (d->abs_mask, y, x); 157 158 if (unlikely (v_any_u64 (special))) 159 return special_case (x, y, special); 160 return y; 161 } 162 163 PL_SIG (V, D, 1, asinh, -10.0, 10.0) 164 PL_TEST_ULP (V_NAME_D1 (asinh), 2.80) 165 PL_TEST_EXPECT_FENV (V_NAME_D1 (asinh), WANT_SIMD_EXCEPT) 166 /* Test vector asinh 3 times, with control lane < 1, > 1 and special. 167 Ensures the v_sel is choosing the right option in all cases. */ 168 #define V_ASINH_INTERVAL(lo, hi, n) \ 169 PL_TEST_SYM_INTERVAL_C (V_NAME_D1 (asinh), lo, hi, n, 0.5) \ 170 PL_TEST_SYM_INTERVAL_C (V_NAME_D1 (asinh), lo, hi, n, 2) \ 171 PL_TEST_SYM_INTERVAL_C (V_NAME_D1 (asinh), lo, hi, n, 0x1p600) 172 V_ASINH_INTERVAL (0, 0x1p-26, 50000) 173 V_ASINH_INTERVAL (0x1p-26, 1, 50000) 174 V_ASINH_INTERVAL (1, 0x1p511, 50000) 175 V_ASINH_INTERVAL (0x1p511, inf, 40000) 176