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