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