1 /* 2 * Single-precision vector acosh(x) function. 3 * Copyright (c) 2023, Arm Limited. 4 * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception 5 */ 6 7 #include "v_math.h" 8 #include "pl_sig.h" 9 #include "pl_test.h" 10 11 #define WANT_V_LOG1P_K0_SHORTCUT 1 12 #include "v_log1p_inline.h" 13 14 const static struct data 15 { 16 struct v_log1p_data log1p_consts; 17 uint64x2_t one, thresh; 18 } data = { 19 .log1p_consts = V_LOG1P_CONSTANTS_TABLE, 20 .one = V2 (0x3ff0000000000000), 21 .thresh = V2 (0x1ff0000000000000) /* asuint64(0x1p511) - asuint64(1). */ 22 }; 23 24 static float64x2_t NOINLINE VPCS_ATTR 25 special_case (float64x2_t x, float64x2_t y, uint64x2_t special, 26 const struct v_log1p_data *d) 27 { 28 return v_call_f64 (acosh, x, log1p_inline (y, d), special); 29 } 30 31 /* Vector approximation for double-precision acosh, based on log1p. 32 The largest observed error is 3.02 ULP in the region where the 33 argument to log1p falls in the k=0 interval, i.e. x close to 1: 34 _ZGVnN2v_acosh(0x1.00798aaf80739p+0) got 0x1.f2d6d823bc9dfp-5 35 want 0x1.f2d6d823bc9e2p-5. */ 36 VPCS_ATTR float64x2_t V_NAME_D1 (acosh) (float64x2_t x) 37 { 38 const struct data *d = ptr_barrier (&data); 39 uint64x2_t special 40 = vcgeq_u64 (vsubq_u64 (vreinterpretq_u64_f64 (x), d->one), d->thresh); 41 float64x2_t special_arg = x; 42 43 #if WANT_SIMD_EXCEPT 44 if (unlikely (v_any_u64 (special))) 45 x = vbslq_f64 (special, vreinterpretq_f64_u64 (d->one), x); 46 #endif 47 48 float64x2_t xm1 = vsubq_f64 (x, v_f64 (1)); 49 float64x2_t y; 50 y = vaddq_f64 (x, v_f64 (1)); 51 y = vmulq_f64 (y, xm1); 52 y = vsqrtq_f64 (y); 53 y = vaddq_f64 (xm1, y); 54 55 if (unlikely (v_any_u64 (special))) 56 return special_case (special_arg, y, special, &d->log1p_consts); 57 return log1p_inline (y, &d->log1p_consts); 58 } 59 60 PL_SIG (V, D, 1, acosh, 1.0, 10.0) 61 PL_TEST_ULP (V_NAME_D1 (acosh), 2.53) 62 PL_TEST_EXPECT_FENV (V_NAME_D1 (acosh), WANT_SIMD_EXCEPT) 63 PL_TEST_INTERVAL (V_NAME_D1 (acosh), 1, 0x1p511, 90000) 64 PL_TEST_INTERVAL (V_NAME_D1 (acosh), 0x1p511, inf, 10000) 65 PL_TEST_INTERVAL (V_NAME_D1 (acosh), 0, 1, 1000) 66 PL_TEST_INTERVAL (V_NAME_D1 (acosh), -0, -inf, 10000) 67