1 /* 2 * Double-precision vector exp(x) - 1 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 "v_math.h" 9 #include "test_sig.h" 10 #include "test_defs.h" 11 #include "v_expm1_inline.h" 12 13 static const struct data 14 { 15 struct v_expm1_data d; 16 #if WANT_SIMD_EXCEPT 17 uint64x2_t thresh, tiny_bound; 18 #else 19 float64x2_t oflow_bound; 20 #endif 21 } data = { 22 .d = V_EXPM1_DATA, 23 #if WANT_SIMD_EXCEPT 24 /* asuint64(oflow_bound) - asuint64(0x1p-51), shifted left by 1 for abs 25 compare. */ 26 .thresh = V2 (0x78c56fa6d34b552), 27 /* asuint64(0x1p-51) << 1. */ 28 .tiny_bound = V2 (0x3cc0000000000000 << 1), 29 #else 30 /* Value above which expm1(x) should overflow. Absolute value of the 31 underflow bound is greater than this, so it catches both cases - there is 32 a small window where fallbacks are triggered unnecessarily. */ 33 .oflow_bound = V2 (0x1.62b7d369a5aa9p+9), 34 #endif 35 }; 36 37 static float64x2_t VPCS_ATTR NOINLINE 38 special_case (float64x2_t x, uint64x2_t special, const struct data *d) 39 { 40 return v_call_f64 (expm1, x, expm1_inline (v_zerofy_f64 (x, special), &d->d), 41 special); 42 } 43 44 /* Double-precision vector exp(x) - 1 function. 45 The maximum error observed error is 2.05 ULP: 46 _ZGVnN2v_expm1(0x1.6329669eb8c87p-2) got 0x1.a8897eef87b34p-2 47 want 0x1.a8897eef87b32p-2. */ 48 float64x2_t VPCS_ATTR V_NAME_D1 (expm1) (float64x2_t x) 49 { 50 const struct data *d = ptr_barrier (&data); 51 52 #if WANT_SIMD_EXCEPT 53 uint64x2_t ix = vreinterpretq_u64_f64 (x); 54 /* If fp exceptions are to be triggered correctly, fall back to scalar for 55 |x| < 2^-51, |x| > oflow_bound, Inf & NaN. Add ix to itself for 56 shift-left by 1, and compare with thresh which was left-shifted offline - 57 this is effectively an absolute compare. */ 58 uint64x2_t special 59 = vcgeq_u64 (vsubq_u64 (vaddq_u64 (ix, ix), d->tiny_bound), d->thresh); 60 #else 61 /* Large input, NaNs and Infs. */ 62 uint64x2_t special = vcageq_f64 (x, d->oflow_bound); 63 #endif 64 65 if (unlikely (v_any_u64 (special))) 66 return special_case (x, special, d); 67 68 /* expm1(x) ~= p * t + (t - 1). */ 69 return expm1_inline (x, &d->d); 70 } 71 72 TEST_SIG (V, D, 1, expm1, -9.9, 9.9) 73 TEST_ULP (V_NAME_D1 (expm1), 1.56) 74 TEST_DISABLE_FENV_IF_NOT (V_NAME_D1 (expm1), WANT_SIMD_EXCEPT) 75 TEST_SYM_INTERVAL (V_NAME_D1 (expm1), 0, 0x1p-51, 1000) 76 TEST_SYM_INTERVAL (V_NAME_D1 (expm1), 0x1p-51, 0x1.62b7d369a5aa9p+9, 100000) 77 TEST_SYM_INTERVAL (V_NAME_D1 (expm1), 0x1.62b7d369a5aa9p+9, inf, 100) 78