1 /* 2 * Double-precision vector asin(x) function. 3 * 4 * Copyright (c) 2023-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 12 static const struct data 13 { 14 float64x2_t c0, c2, c4, c6, c8, c10; 15 float64x2_t pi_over_2; 16 uint64x2_t abs_mask; 17 double c1, c3, c5, c7, c9, c11; 18 } data = { 19 /* Polynomial approximation of (asin(sqrt(x)) - sqrt(x)) / (x * sqrt(x)) 20 on [ 0x1p-106, 0x1p-2 ], relative error: 0x1.c3d8e169p-57. */ 21 .c0 = V2 (0x1.555555555554ep-3), .c1 = 0x1.3333333337233p-4, 22 .c2 = V2 (0x1.6db6db67f6d9fp-5), .c3 = 0x1.f1c71fbd29fbbp-6, 23 .c4 = V2 (0x1.6e8b264d467d6p-6), .c5 = 0x1.1c5997c357e9dp-6, 24 .c6 = V2 (0x1.c86a22cd9389dp-7), .c7 = 0x1.856073c22ebbep-7, 25 .c8 = V2 (0x1.fd1151acb6bedp-8), .c9 = 0x1.087182f799c1dp-6, 26 .c10 = V2 (-0x1.6602748120927p-7), .c11 = 0x1.cfa0dd1f9478p-6, 27 .pi_over_2 = V2 (0x1.921fb54442d18p+0), .abs_mask = V2 (0x7fffffffffffffff), 28 }; 29 30 #define AllMask v_u64 (0xffffffffffffffff) 31 #define One 0x3ff0000000000000 32 #define Small 0x3e50000000000000 /* 2^-12. */ 33 34 #if WANT_SIMD_EXCEPT 35 static float64x2_t VPCS_ATTR NOINLINE 36 special_case (float64x2_t x, float64x2_t y, uint64x2_t special) 37 { 38 return v_call_f64 (asin, x, y, special); 39 } 40 #endif 41 42 /* Double-precision implementation of vector asin(x). 43 44 For |x| < Small, approximate asin(x) by x. Small = 2^-12 for correct 45 rounding. If WANT_SIMD_EXCEPT = 0, Small = 0 and we proceed with the 46 following approximation. 47 48 For |x| in [Small, 0.5], use an order 11 polynomial P such that the final 49 approximation is an odd polynomial: asin(x) ~ x + x^3 P(x^2). 50 51 The largest observed error in this region is 1.01 ulps, 52 _ZGVnN2v_asin (0x1.da9735b5a9277p-2) got 0x1.ed78525a927efp-2 53 want 0x1.ed78525a927eep-2. 54 55 For |x| in [0.5, 1.0], use same approximation with a change of variable 56 57 asin(x) = pi/2 - (y + y * z * P(z)), with z = (1-x)/2 and y = sqrt(z). 58 59 The largest observed error in this region is 2.69 ulps, 60 _ZGVnN2v_asin (0x1.044e8cefee301p-1) got 0x1.1111dd54ddf96p-1 61 want 0x1.1111dd54ddf99p-1. */ 62 float64x2_t VPCS_ATTR V_NAME_D1 (asin) (float64x2_t x) 63 { 64 const struct data *d = ptr_barrier (&data); 65 float64x2_t ax = vabsq_f64 (x); 66 67 #if WANT_SIMD_EXCEPT 68 /* Special values need to be computed with scalar fallbacks so 69 that appropriate exceptions are raised. */ 70 uint64x2_t special 71 = vcgtq_u64 (vsubq_u64 (vreinterpretq_u64_f64 (ax), v_u64 (Small)), 72 v_u64 (One - Small)); 73 if (unlikely (v_any_u64 (special))) 74 return special_case (x, x, AllMask); 75 #endif 76 77 uint64x2_t a_lt_half = vcaltq_f64 (x, v_f64 (0.5)); 78 79 /* Evaluate polynomial Q(x) = y + y * z * P(z) with 80 z = x ^ 2 and y = |x| , if |x| < 0.5 81 z = (1 - |x|) / 2 and y = sqrt(z), if |x| >= 0.5. */ 82 float64x2_t z2 = vbslq_f64 (a_lt_half, vmulq_f64 (x, x), 83 vfmsq_n_f64 (v_f64 (0.5), ax, 0.5)); 84 float64x2_t z = vbslq_f64 (a_lt_half, ax, vsqrtq_f64 (z2)); 85 86 /* Use a single polynomial approximation P for both intervals. */ 87 float64x2_t z4 = vmulq_f64 (z2, z2); 88 float64x2_t z8 = vmulq_f64 (z4, z4); 89 float64x2_t z16 = vmulq_f64 (z8, z8); 90 91 /* order-11 estrin. */ 92 float64x2_t c13 = vld1q_f64 (&d->c1); 93 float64x2_t c57 = vld1q_f64 (&d->c5); 94 float64x2_t c911 = vld1q_f64 (&d->c9); 95 96 float64x2_t p01 = vfmaq_laneq_f64 (d->c0, z2, c13, 0); 97 float64x2_t p23 = vfmaq_laneq_f64 (d->c2, z2, c13, 1); 98 float64x2_t p03 = vfmaq_f64 (p01, z4, p23); 99 100 float64x2_t p45 = vfmaq_laneq_f64 (d->c4, z2, c57, 0); 101 float64x2_t p67 = vfmaq_laneq_f64 (d->c6, z2, c57, 1); 102 float64x2_t p47 = vfmaq_f64 (p45, z4, p67); 103 104 float64x2_t p89 = vfmaq_laneq_f64 (d->c8, z2, c911, 0); 105 float64x2_t p1011 = vfmaq_laneq_f64 (d->c10, z2, c911, 1); 106 float64x2_t p811 = vfmaq_f64 (p89, z4, p1011); 107 108 float64x2_t p07 = vfmaq_f64 (p03, z8, p47); 109 float64x2_t p = vfmaq_f64 (p07, z16, p811); 110 111 /* Finalize polynomial: z + z * z2 * P(z2). */ 112 p = vfmaq_f64 (z, vmulq_f64 (z, z2), p); 113 114 /* asin(|x|) = Q(|x|) , for |x| < 0.5 115 = pi/2 - 2 Q(|x|), for |x| >= 0.5. */ 116 float64x2_t y = vbslq_f64 (a_lt_half, p, vfmsq_n_f64 (d->pi_over_2, p, 2.0)); 117 118 /* Copy sign. */ 119 return vbslq_f64 (d->abs_mask, y, x); 120 } 121 122 TEST_SIG (V, D, 1, asin, -1.0, 1.0) 123 TEST_ULP (V_NAME_D1 (asin), 2.20) 124 TEST_DISABLE_FENV_IF_NOT (V_NAME_D1 (asin), WANT_SIMD_EXCEPT) 125 TEST_INTERVAL (V_NAME_D1 (asin), 0, Small, 5000) 126 TEST_INTERVAL (V_NAME_D1 (asin), Small, 0.5, 50000) 127 TEST_INTERVAL (V_NAME_D1 (asin), 0.5, 1.0, 50000) 128 TEST_INTERVAL (V_NAME_D1 (asin), 1.0, 0x1p11, 50000) 129 TEST_INTERVAL (V_NAME_D1 (asin), 0x1p11, inf, 20000) 130 TEST_INTERVAL (V_NAME_D1 (asin), -0, -inf, 20000) 131