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