1 | /* Single-precision AdvSIMD atan2 |
2 | |
3 | Copyright (C) 2023-2024 Free Software Foundation, Inc. |
4 | This file is part of the GNU C Library. |
5 | |
6 | The GNU C Library is free software; you can redistribute it and/or |
7 | modify it under the terms of the GNU Lesser General Public |
8 | License as published by the Free Software Foundation; either |
9 | version 2.1 of the License, or (at your option) any later version. |
10 | |
11 | The GNU C Library is distributed in the hope that it will be useful, |
12 | but WITHOUT ANY WARRANTY; without even the implied warranty of |
13 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
14 | Lesser General Public License for more details. |
15 | |
16 | You should have received a copy of the GNU Lesser General Public |
17 | License along with the GNU C Library; if not, see |
18 | <https://www.gnu.org/licenses/>. */ |
19 | |
20 | #include "v_math.h" |
21 | #include "poly_advsimd_f32.h" |
22 | |
23 | static const struct data |
24 | { |
25 | float32x4_t poly[8]; |
26 | float32x4_t pi_over_2; |
27 | } data = { |
28 | /* Coefficients of polynomial P such that atan(x)~x+x*P(x^2) on |
29 | [2**-128, 1.0]. |
30 | Generated using fpminimax between FLT_MIN and 1. */ |
31 | .poly = { V4 (-0x1.55555p-2f), V4 (0x1.99935ep-3f), V4 (-0x1.24051ep-3f), |
32 | V4 (0x1.bd7368p-4f), V4 (-0x1.491f0ep-4f), V4 (0x1.93a2c0p-5f), |
33 | V4 (-0x1.4c3c60p-6f), V4 (0x1.01fd88p-8f) }, |
34 | .pi_over_2 = V4 (0x1.921fb6p+0f), |
35 | }; |
36 | |
37 | #define SignMask v_u32 (0x80000000) |
38 | |
39 | /* Special cases i.e. 0, infinity and nan (fall back to scalar calls). */ |
40 | static float32x4_t VPCS_ATTR NOINLINE |
41 | special_case (float32x4_t y, float32x4_t x, float32x4_t ret, uint32x4_t cmp) |
42 | { |
43 | return v_call2_f32 (atan2f, y, x, ret, cmp); |
44 | } |
45 | |
46 | /* Returns 1 if input is the bit representation of 0, infinity or nan. */ |
47 | static inline uint32x4_t |
48 | zeroinfnan (uint32x4_t i) |
49 | { |
50 | /* 2 * i - 1 >= 2 * 0x7f800000lu - 1. */ |
51 | return vcgeq_u32 (vsubq_u32 (vmulq_n_u32 (i, 2), v_u32 (x: 1)), |
52 | v_u32 (x: 2 * 0x7f800000lu - 1)); |
53 | } |
54 | |
55 | /* Fast implementation of vector atan2f. Maximum observed error is |
56 | 2.95 ULP in [0x1.9300d6p+6 0x1.93c0c6p+6] x [0x1.8c2dbp+6 0x1.8cea6p+6]: |
57 | _ZGVnN4vv_atan2f (0x1.93836cp+6, 0x1.8cae1p+6) got 0x1.967f06p-1 |
58 | want 0x1.967f00p-1. */ |
59 | float32x4_t VPCS_ATTR NOINLINE V_NAME_F2 (atan2) (float32x4_t y, float32x4_t x) |
60 | { |
61 | const struct data *data_ptr = ptr_barrier (&data); |
62 | |
63 | uint32x4_t ix = vreinterpretq_u32_f32 (x); |
64 | uint32x4_t iy = vreinterpretq_u32_f32 (y); |
65 | |
66 | uint32x4_t special_cases = vorrq_u32 (zeroinfnan (ix), zeroinfnan (iy)); |
67 | |
68 | uint32x4_t sign_x = vandq_u32 (ix, SignMask); |
69 | uint32x4_t sign_y = vandq_u32 (iy, SignMask); |
70 | uint32x4_t sign_xy = veorq_u32 (sign_x, sign_y); |
71 | |
72 | float32x4_t ax = vabsq_f32 (x); |
73 | float32x4_t ay = vabsq_f32 (y); |
74 | |
75 | uint32x4_t pred_xlt0 = vcltzq_f32 (x); |
76 | uint32x4_t pred_aygtax = vcgtq_f32 (ay, ax); |
77 | |
78 | /* Set up z for call to atanf. */ |
79 | float32x4_t n = vbslq_f32 (pred_aygtax, vnegq_f32 (ax), ay); |
80 | float32x4_t d = vbslq_f32 (pred_aygtax, ay, ax); |
81 | float32x4_t z = vdivq_f32 (n, d); |
82 | |
83 | /* Work out the correct shift. */ |
84 | float32x4_t shift = vreinterpretq_f32_u32 ( |
85 | vandq_u32 (pred_xlt0, vreinterpretq_u32_f32 (v_f32 (-2.0f)))); |
86 | shift = vbslq_f32 (pred_aygtax, vaddq_f32 (shift, v_f32 (1.0f)), shift); |
87 | shift = vmulq_f32 (shift, data_ptr->pi_over_2); |
88 | |
89 | /* Calculate the polynomial approximation. |
90 | Use 2-level Estrin scheme for P(z^2) with deg(P)=7. However, |
91 | a standard implementation using z8 creates spurious underflow |
92 | in the very last fma (when z^8 is small enough). |
93 | Therefore, we split the last fma into a mul and an fma. |
94 | Horner and single-level Estrin have higher errors that exceed |
95 | threshold. */ |
96 | float32x4_t z2 = vmulq_f32 (z, z); |
97 | float32x4_t z4 = vmulq_f32 (z2, z2); |
98 | |
99 | float32x4_t ret = vfmaq_f32 ( |
100 | v_pairwise_poly_3_f32 (z2, z4, data_ptr->poly), z4, |
101 | vmulq_f32 (z4, v_pairwise_poly_3_f32 (z2, z4, data_ptr->poly + 4))); |
102 | |
103 | /* y = shift + z * P(z^2). */ |
104 | ret = vaddq_f32 (vfmaq_f32 (z, ret, vmulq_f32 (z2, z)), shift); |
105 | |
106 | /* Account for the sign of y. */ |
107 | ret = vreinterpretq_f32_u32 ( |
108 | veorq_u32 (vreinterpretq_u32_f32 (ret), sign_xy)); |
109 | |
110 | if (__glibc_unlikely (v_any_u32 (special_cases))) |
111 | { |
112 | return special_case (y, x, ret, special_cases); |
113 | } |
114 | |
115 | return ret; |
116 | } |
117 | libmvec_hidden_def (V_NAME_F2 (atan2)) |
118 | HALF_WIDTH_ALIAS_F2(atan2) |
119 | |