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
special_case(float64x2_t x,float64x2_t y,uint64x2_t special)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. */
V_NAME_D1(asin)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