1*072a4ba8SAndrew Turner /*
2*072a4ba8SAndrew Turner * Single-precision scalar atan2(x) function.
3*072a4ba8SAndrew Turner *
4*072a4ba8SAndrew Turner * Copyright (c) 2021-2023, Arm Limited.
5*072a4ba8SAndrew Turner * SPDX-License-Identifier: MIT OR Apache-2.0 WITH LLVM-exception
6*072a4ba8SAndrew Turner */
7*072a4ba8SAndrew Turner
8*072a4ba8SAndrew Turner #include <stdbool.h>
9*072a4ba8SAndrew Turner
10*072a4ba8SAndrew Turner #include "atanf_common.h"
11*072a4ba8SAndrew Turner #include "math_config.h"
12*072a4ba8SAndrew Turner #include "pl_sig.h"
13*072a4ba8SAndrew Turner #include "pl_test.h"
14*072a4ba8SAndrew Turner
15*072a4ba8SAndrew Turner #define Pi (0x1.921fb6p+1f)
16*072a4ba8SAndrew Turner #define PiOver2 (0x1.921fb6p+0f)
17*072a4ba8SAndrew Turner #define PiOver4 (0x1.921fb6p-1f)
18*072a4ba8SAndrew Turner #define SignMask (0x80000000)
19*072a4ba8SAndrew Turner
20*072a4ba8SAndrew Turner /* We calculate atan2f by P(n/d), where n and d are similar to the input
21*072a4ba8SAndrew Turner arguments, and P is a polynomial. The polynomial may underflow.
22*072a4ba8SAndrew Turner POLY_UFLOW_BOUND is the lower bound of the difference in exponents of n and d
23*072a4ba8SAndrew Turner for which P underflows, and is used to special-case such inputs. */
24*072a4ba8SAndrew Turner #define POLY_UFLOW_BOUND 24
25*072a4ba8SAndrew Turner
26*072a4ba8SAndrew Turner static inline int32_t
biased_exponent(float f)27*072a4ba8SAndrew Turner biased_exponent (float f)
28*072a4ba8SAndrew Turner {
29*072a4ba8SAndrew Turner uint32_t fi = asuint (f);
30*072a4ba8SAndrew Turner int32_t ex = (int32_t) ((fi & 0x7f800000) >> 23);
31*072a4ba8SAndrew Turner if (unlikely (ex == 0))
32*072a4ba8SAndrew Turner {
33*072a4ba8SAndrew Turner /* Subnormal case - we still need to get the exponent right for subnormal
34*072a4ba8SAndrew Turner numbers as division may take us back inside the normal range. */
35*072a4ba8SAndrew Turner return ex - __builtin_clz (fi << 9);
36*072a4ba8SAndrew Turner }
37*072a4ba8SAndrew Turner return ex;
38*072a4ba8SAndrew Turner }
39*072a4ba8SAndrew Turner
40*072a4ba8SAndrew Turner /* Fast implementation of scalar atan2f. Largest observed error is
41*072a4ba8SAndrew Turner 2.88ulps in [99.0, 101.0] x [99.0, 101.0]:
42*072a4ba8SAndrew Turner atan2f(0x1.9332d8p+6, 0x1.8cb6c4p+6) got 0x1.964646p-1
43*072a4ba8SAndrew Turner want 0x1.964640p-1. */
44*072a4ba8SAndrew Turner float
atan2f(float y,float x)45*072a4ba8SAndrew Turner atan2f (float y, float x)
46*072a4ba8SAndrew Turner {
47*072a4ba8SAndrew Turner uint32_t ix = asuint (x);
48*072a4ba8SAndrew Turner uint32_t iy = asuint (y);
49*072a4ba8SAndrew Turner
50*072a4ba8SAndrew Turner uint32_t sign_x = ix & SignMask;
51*072a4ba8SAndrew Turner uint32_t sign_y = iy & SignMask;
52*072a4ba8SAndrew Turner
53*072a4ba8SAndrew Turner uint32_t iax = ix & ~SignMask;
54*072a4ba8SAndrew Turner uint32_t iay = iy & ~SignMask;
55*072a4ba8SAndrew Turner
56*072a4ba8SAndrew Turner /* x or y is NaN. */
57*072a4ba8SAndrew Turner if ((iax > 0x7f800000) || (iay > 0x7f800000))
58*072a4ba8SAndrew Turner return x + y;
59*072a4ba8SAndrew Turner
60*072a4ba8SAndrew Turner /* m = 2 * sign(x) + sign(y). */
61*072a4ba8SAndrew Turner uint32_t m = ((iy >> 31) & 1) | ((ix >> 30) & 2);
62*072a4ba8SAndrew Turner
63*072a4ba8SAndrew Turner /* The following follows glibc ieee754 implementation, except
64*072a4ba8SAndrew Turner that we do not use +-tiny shifts (non-nearest rounding mode). */
65*072a4ba8SAndrew Turner
66*072a4ba8SAndrew Turner int32_t exp_diff = biased_exponent (x) - biased_exponent (y);
67*072a4ba8SAndrew Turner
68*072a4ba8SAndrew Turner /* Special case for (x, y) either on or very close to the x axis. Either y =
69*072a4ba8SAndrew Turner 0, or y is tiny and x is huge (difference in exponents >=
70*072a4ba8SAndrew Turner POLY_UFLOW_BOUND). In the second case, we only want to use this special
71*072a4ba8SAndrew Turner case when x is negative (i.e. quadrants 2 or 3). */
72*072a4ba8SAndrew Turner if (unlikely (iay == 0 || (exp_diff >= POLY_UFLOW_BOUND && m >= 2)))
73*072a4ba8SAndrew Turner {
74*072a4ba8SAndrew Turner switch (m)
75*072a4ba8SAndrew Turner {
76*072a4ba8SAndrew Turner case 0:
77*072a4ba8SAndrew Turner case 1:
78*072a4ba8SAndrew Turner return y; /* atan(+-0,+anything)=+-0. */
79*072a4ba8SAndrew Turner case 2:
80*072a4ba8SAndrew Turner return Pi; /* atan(+0,-anything) = pi. */
81*072a4ba8SAndrew Turner case 3:
82*072a4ba8SAndrew Turner return -Pi; /* atan(-0,-anything) =-pi. */
83*072a4ba8SAndrew Turner }
84*072a4ba8SAndrew Turner }
85*072a4ba8SAndrew Turner /* Special case for (x, y) either on or very close to the y axis. Either x =
86*072a4ba8SAndrew Turner 0, or x is tiny and y is huge (difference in exponents >=
87*072a4ba8SAndrew Turner POLY_UFLOW_BOUND). */
88*072a4ba8SAndrew Turner if (unlikely (iax == 0 || exp_diff <= -POLY_UFLOW_BOUND))
89*072a4ba8SAndrew Turner return sign_y ? -PiOver2 : PiOver2;
90*072a4ba8SAndrew Turner
91*072a4ba8SAndrew Turner /* x is INF. */
92*072a4ba8SAndrew Turner if (iax == 0x7f800000)
93*072a4ba8SAndrew Turner {
94*072a4ba8SAndrew Turner if (iay == 0x7f800000)
95*072a4ba8SAndrew Turner {
96*072a4ba8SAndrew Turner switch (m)
97*072a4ba8SAndrew Turner {
98*072a4ba8SAndrew Turner case 0:
99*072a4ba8SAndrew Turner return PiOver4; /* atan(+INF,+INF). */
100*072a4ba8SAndrew Turner case 1:
101*072a4ba8SAndrew Turner return -PiOver4; /* atan(-INF,+INF). */
102*072a4ba8SAndrew Turner case 2:
103*072a4ba8SAndrew Turner return 3.0f * PiOver4; /* atan(+INF,-INF). */
104*072a4ba8SAndrew Turner case 3:
105*072a4ba8SAndrew Turner return -3.0f * PiOver4; /* atan(-INF,-INF). */
106*072a4ba8SAndrew Turner }
107*072a4ba8SAndrew Turner }
108*072a4ba8SAndrew Turner else
109*072a4ba8SAndrew Turner {
110*072a4ba8SAndrew Turner switch (m)
111*072a4ba8SAndrew Turner {
112*072a4ba8SAndrew Turner case 0:
113*072a4ba8SAndrew Turner return 0.0f; /* atan(+...,+INF). */
114*072a4ba8SAndrew Turner case 1:
115*072a4ba8SAndrew Turner return -0.0f; /* atan(-...,+INF). */
116*072a4ba8SAndrew Turner case 2:
117*072a4ba8SAndrew Turner return Pi; /* atan(+...,-INF). */
118*072a4ba8SAndrew Turner case 3:
119*072a4ba8SAndrew Turner return -Pi; /* atan(-...,-INF). */
120*072a4ba8SAndrew Turner }
121*072a4ba8SAndrew Turner }
122*072a4ba8SAndrew Turner }
123*072a4ba8SAndrew Turner /* y is INF. */
124*072a4ba8SAndrew Turner if (iay == 0x7f800000)
125*072a4ba8SAndrew Turner return sign_y ? -PiOver2 : PiOver2;
126*072a4ba8SAndrew Turner
127*072a4ba8SAndrew Turner uint32_t sign_xy = sign_x ^ sign_y;
128*072a4ba8SAndrew Turner
129*072a4ba8SAndrew Turner float ax = asfloat (iax);
130*072a4ba8SAndrew Turner float ay = asfloat (iay);
131*072a4ba8SAndrew Turner
132*072a4ba8SAndrew Turner bool pred_aygtax = (ay > ax);
133*072a4ba8SAndrew Turner
134*072a4ba8SAndrew Turner /* Set up z for call to atanf. */
135*072a4ba8SAndrew Turner float n = pred_aygtax ? -ax : ay;
136*072a4ba8SAndrew Turner float d = pred_aygtax ? ay : ax;
137*072a4ba8SAndrew Turner float z = n / d;
138*072a4ba8SAndrew Turner
139*072a4ba8SAndrew Turner float ret;
140*072a4ba8SAndrew Turner if (unlikely (m < 2 && exp_diff >= POLY_UFLOW_BOUND))
141*072a4ba8SAndrew Turner {
142*072a4ba8SAndrew Turner /* If (x, y) is very close to x axis and x is positive, the polynomial
143*072a4ba8SAndrew Turner will underflow and evaluate to z. */
144*072a4ba8SAndrew Turner ret = z;
145*072a4ba8SAndrew Turner }
146*072a4ba8SAndrew Turner else
147*072a4ba8SAndrew Turner {
148*072a4ba8SAndrew Turner /* Work out the correct shift. */
149*072a4ba8SAndrew Turner float shift = sign_x ? -2.0f : 0.0f;
150*072a4ba8SAndrew Turner shift = pred_aygtax ? shift + 1.0f : shift;
151*072a4ba8SAndrew Turner shift *= PiOver2;
152*072a4ba8SAndrew Turner
153*072a4ba8SAndrew Turner ret = eval_poly (z, z, shift);
154*072a4ba8SAndrew Turner }
155*072a4ba8SAndrew Turner
156*072a4ba8SAndrew Turner /* Account for the sign of x and y. */
157*072a4ba8SAndrew Turner return asfloat (asuint (ret) ^ sign_xy);
158*072a4ba8SAndrew Turner }
159*072a4ba8SAndrew Turner
160*072a4ba8SAndrew Turner /* Arity of 2 means no mathbench entry emitted. See test/mathbench_funcs.h. */
161*072a4ba8SAndrew Turner PL_SIG (S, F, 2, atan2)
162*072a4ba8SAndrew Turner PL_TEST_ULP (atan2f, 2.4)
163*072a4ba8SAndrew Turner PL_TEST_INTERVAL (atan2f, -10.0, 10.0, 50000)
164*072a4ba8SAndrew Turner PL_TEST_INTERVAL (atan2f, -1.0, 1.0, 40000)
165*072a4ba8SAndrew Turner PL_TEST_INTERVAL (atan2f, 0.0, 1.0, 40000)
166*072a4ba8SAndrew Turner PL_TEST_INTERVAL (atan2f, 1.0, 100.0, 40000)
167*072a4ba8SAndrew Turner PL_TEST_INTERVAL (atan2f, 1e6, 1e32, 40000)
168