1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/math/atan2f.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/atan2.c
6
7const std = @import("../std.zig");
8const math = std.math;
9const expect = std.testing.expect;
10
11/// Returns the arc-tangent of y/x.
12///
13/// Special Cases:
14/// | y | x | radians |
15/// |-------|-------|---------|
16/// | fin | nan | nan |
17/// | nan | fin | nan |
18/// | +0 | >=+0 | +0 |
19/// | -0 | >=+0 | -0 |
20/// | +0 | <=-0 | pi |
21/// | -0 | <=-0 | -pi |
22/// | pos | 0 | +pi/2 |
23/// | neg | 0 | -pi/2 |
24/// | +inf | +inf | +pi/4 |
25/// | -inf | +inf | -pi/4 |
26/// | +inf | -inf | 3pi/4 |
27/// | -inf | -inf | -3pi/4 |
28/// | fin | +inf | 0 |
29/// | pos | -inf | +pi |
30/// | neg | -inf | -pi |
31/// | +inf | fin | +pi/2 |
32/// | -inf | fin | -pi/2 |
33pub fn atan2(y: anytype, x: anytype) @TypeOf(x, y) {
34 const T = @TypeOf(x, y);
35 return switch (T) {
36 f32 => atan2_32(y, x),
37 f64 => atan2_64(y, x),
38 else => @compileError("atan2 not implemented for " ++ @typeName(T)),
39 };
40}
41
42fn atan2_32(y: f32, x: f32) f32 {
43 const pi: f32 = 3.1415927410e+00;
44 const pi_lo: f32 = -8.7422776573e-08;
45
46 if (math.isNan(x) or math.isNan(y)) {
47 return x + y;
48 }
49
50 var ix = @as(u32, @bitCast(x));
51 var iy = @as(u32, @bitCast(y));
52
53 // x = 1.0
54 if (ix == 0x3F800000) {
55 return math.atan(y);
56 }
57
58 // 2 * sign(x) + sign(y)
59 const m = ((iy >> 31) & 1) | ((ix >> 30) & 2);
60 ix &= 0x7FFFFFFF;
61 iy &= 0x7FFFFFFF;
62
63 if (iy == 0) {
64 switch (m) {
65 0, 1 => return y, // atan(+-0, +...)
66 2 => return pi, // atan(+0, -...)
67 3 => return -pi, // atan(-0, -...)
68 else => unreachable,
69 }
70 }
71
72 if (ix == 0) {
73 if (m & 1 != 0) {
74 return -pi / 2;
75 } else {
76 return pi / 2;
77 }
78 }
79
80 if (ix == 0x7F800000) {
81 if (iy == 0x7F800000) {
82 switch (m) {
83 0 => return pi / 4, // atan(+inf, +inf)
84 1 => return -pi / 4, // atan(-inf, +inf)
85 2 => return 3 * pi / 4, // atan(+inf, -inf)
86 3 => return -3 * pi / 4, // atan(-inf, -inf)
87 else => unreachable,
88 }
89 } else {
90 switch (m) {
91 0 => return 0.0, // atan(+..., +inf)
92 1 => return -0.0, // atan(-..., +inf)
93 2 => return pi, // atan(+..., -inf)
94 3 => return -pi, // atan(-...f, -inf)
95 else => unreachable,
96 }
97 }
98 }
99
100 // |y / x| > 0x1p26
101 if (ix + (26 << 23) < iy or iy == 0x7F800000) {
102 if (m & 1 != 0) {
103 return -pi / 2;
104 } else {
105 return pi / 2;
106 }
107 }
108
109 // z = atan(|y / x|) with correct underflow
110 const z = z: {
111 if ((m & 2) != 0 and iy + (26 << 23) < ix) {
112 break :z 0.0;
113 } else {
114 break :z math.atan(@abs(y / x));
115 }
116 };
117
118 switch (m) {
119 0 => return z, // atan(+, +)
120 1 => return -z, // atan(-, +)
121 2 => return pi - (z - pi_lo), // atan(+, -)
122 3 => return (z - pi_lo) - pi, // atan(-, -)
123 else => unreachable,
124 }
125}
126
127fn atan2_64(y: f64, x: f64) f64 {
128 const pi: f64 = 3.1415926535897931160E+00;
129 const pi_lo: f64 = 1.2246467991473531772E-16;
130
131 if (math.isNan(x) or math.isNan(y)) {
132 return x + y;
133 }
134
135 const ux: u64 = @bitCast(x);
136 var ix: u32 = @intCast(ux >> 32);
137 const lx: u32 = @intCast(ux & 0xFFFFFFFF);
138
139 const uy: u64 = @bitCast(y);
140 var iy: u32 = @intCast(uy >> 32);
141 const ly: u32 = @intCast(uy & 0xFFFFFFFF);
142
143 // x = 1.0
144 if ((ix -% 0x3FF00000) | lx == 0) {
145 return math.atan(y);
146 }
147
148 // 2 * sign(x) + sign(y)
149 const m = ((iy >> 31) & 1) | ((ix >> 30) & 2);
150 ix &= 0x7FFFFFFF;
151 iy &= 0x7FFFFFFF;
152
153 if (iy | ly == 0) {
154 switch (m) {
155 0, 1 => return y, // atan(+-0, +...)
156 2 => return pi, // atan(+0, -...)
157 3 => return -pi, // atan(-0, -...)
158 else => unreachable,
159 }
160 }
161
162 if (ix | lx == 0) {
163 if (m & 1 != 0) {
164 return -pi / 2;
165 } else {
166 return pi / 2;
167 }
168 }
169
170 if (ix == 0x7FF00000) {
171 if (iy == 0x7FF00000) {
172 switch (m) {
173 0 => return pi / 4, // atan(+inf, +inf)
174 1 => return -pi / 4, // atan(-inf, +inf)
175 2 => return 3 * pi / 4, // atan(+inf, -inf)
176 3 => return -3 * pi / 4, // atan(-inf, -inf)
177 else => unreachable,
178 }
179 } else {
180 switch (m) {
181 0 => return 0.0, // atan(+..., +inf)
182 1 => return -0.0, // atan(-..., +inf)
183 2 => return pi, // atan(+..., -inf)
184 3 => return -pi, // atan(-...f, -inf)
185 else => unreachable,
186 }
187 }
188 }
189
190 // |y / x| > 0x1p64
191 if (ix +% (64 << 20) < iy or iy == 0x7FF00000) {
192 if (m & 1 != 0) {
193 return -pi / 2;
194 } else {
195 return pi / 2;
196 }
197 }
198
199 // z = atan(|y / x|) with correct underflow
200 const z = z: {
201 if ((m & 2) != 0 and iy +% (64 << 20) < ix) {
202 break :z 0.0;
203 } else {
204 break :z math.atan(@abs(y / x));
205 }
206 };
207
208 switch (m) {
209 0 => return z, // atan(+, +)
210 1 => return -z, // atan(-, +)
211 2 => return pi - (z - pi_lo), // atan(+, -)
212 3 => return (z - pi_lo) - pi, // atan(-, -)
213 else => unreachable,
214 }
215}
216
217test atan2 {
218 const y32: f32 = 0.2;
219 const x32: f32 = 0.21;
220 const y64: f64 = 0.2;
221 const x64: f64 = 0.21;
222 try expect(atan2(y32, x32) == atan2_32(0.2, 0.21));
223 try expect(atan2(y64, x64) == atan2_64(0.2, 0.21));
224}
225
226test atan2_32 {
227 const epsilon = 0.000001;
228
229 try expect(math.approxEqAbs(f32, atan2_32(0.0, 0.0), 0.0, epsilon));
230 try expect(math.approxEqAbs(f32, atan2_32(0.2, 0.2), 0.785398, epsilon));
231 try expect(math.approxEqAbs(f32, atan2_32(-0.2, 0.2), -0.785398, epsilon));
232 try expect(math.approxEqAbs(f32, atan2_32(0.2, -0.2), 2.356194, epsilon));
233 try expect(math.approxEqAbs(f32, atan2_32(-0.2, -0.2), -2.356194, epsilon));
234 try expect(math.approxEqAbs(f32, atan2_32(0.34, -0.4), 2.437099, epsilon));
235 try expect(math.approxEqAbs(f32, atan2_32(0.34, 1.243), 0.267001, epsilon));
236}
237
238test atan2_64 {
239 const epsilon = 0.000001;
240
241 try expect(math.approxEqAbs(f64, atan2_64(0.0, 0.0), 0.0, epsilon));
242 try expect(math.approxEqAbs(f64, atan2_64(0.2, 0.2), 0.785398, epsilon));
243 try expect(math.approxEqAbs(f64, atan2_64(-0.2, 0.2), -0.785398, epsilon));
244 try expect(math.approxEqAbs(f64, atan2_64(0.2, -0.2), 2.356194, epsilon));
245 try expect(math.approxEqAbs(f64, atan2_64(-0.2, -0.2), -2.356194, epsilon));
246 try expect(math.approxEqAbs(f64, atan2_64(0.34, -0.4), 2.437099, epsilon));
247 try expect(math.approxEqAbs(f64, atan2_64(0.34, 1.243), 0.267001, epsilon));
248}
249
250test "atan2_32.special" {
251 const epsilon = 0.000001;
252
253 try expect(math.isNan(atan2_32(1.0, math.nan(f32))));
254 try expect(math.isNan(atan2_32(math.nan(f32), 1.0)));
255 try expect(math.isPositiveZero(atan2_32(0.0, 5.0)));
256 try expect(math.isNegativeZero(atan2_32(-0.0, 5.0)));
257 try expect(math.approxEqAbs(f32, atan2_32(0.0, -5.0), math.pi, epsilon));
258 //expect(math.approxEqAbs(f32, atan2_32(-0.0, -5.0), -math.pi, .{.rel=0,.abs=epsilon})); TODO support negative zero?
259 try expect(math.approxEqAbs(f32, atan2_32(1.0, 0.0), math.pi / 2.0, epsilon));
260 try expect(math.approxEqAbs(f32, atan2_32(1.0, -0.0), math.pi / 2.0, epsilon));
261 try expect(math.approxEqAbs(f32, atan2_32(-1.0, 0.0), -math.pi / 2.0, epsilon));
262 try expect(math.approxEqAbs(f32, atan2_32(-1.0, -0.0), -math.pi / 2.0, epsilon));
263 try expect(math.approxEqAbs(f32, atan2_32(math.inf(f32), math.inf(f32)), math.pi / 4.0, epsilon));
264 try expect(math.approxEqAbs(f32, atan2_32(-math.inf(f32), math.inf(f32)), -math.pi / 4.0, epsilon));
265 try expect(math.approxEqAbs(f32, atan2_32(math.inf(f32), -math.inf(f32)), 3.0 * math.pi / 4.0, epsilon));
266 try expect(math.approxEqAbs(f32, atan2_32(-math.inf(f32), -math.inf(f32)), -3.0 * math.pi / 4.0, epsilon));
267 try expect(atan2_32(1.0, math.inf(f32)) == 0.0);
268 try expect(math.approxEqAbs(f32, atan2_32(1.0, -math.inf(f32)), math.pi, epsilon));
269 try expect(math.approxEqAbs(f32, atan2_32(-1.0, -math.inf(f32)), -math.pi, epsilon));
270 try expect(math.approxEqAbs(f32, atan2_32(math.inf(f32), 1.0), math.pi / 2.0, epsilon));
271 try expect(math.approxEqAbs(f32, atan2_32(-math.inf(f32), 1.0), -math.pi / 2.0, epsilon));
272}
273
274test "atan2_64.special" {
275 const epsilon = 0.000001;
276
277 try expect(math.isNan(atan2_64(1.0, math.nan(f64))));
278 try expect(math.isNan(atan2_64(math.nan(f64), 1.0)));
279 try expect(math.isPositiveZero(atan2_64(0.0, 5.0)));
280 try expect(math.isNegativeZero(atan2_64(-0.0, 5.0)));
281 try expect(math.approxEqAbs(f64, atan2_64(0.0, -5.0), math.pi, epsilon));
282 //expect(math.approxEqAbs(f64, atan2_64(-0.0, -5.0), -math.pi, .{.rel=0,.abs=epsilon})); TODO support negative zero?
283 try expect(math.approxEqAbs(f64, atan2_64(1.0, 0.0), math.pi / 2.0, epsilon));
284 try expect(math.approxEqAbs(f64, atan2_64(1.0, -0.0), math.pi / 2.0, epsilon));
285 try expect(math.approxEqAbs(f64, atan2_64(-1.0, 0.0), -math.pi / 2.0, epsilon));
286 try expect(math.approxEqAbs(f64, atan2_64(-1.0, -0.0), -math.pi / 2.0, epsilon));
287 try expect(math.approxEqAbs(f64, atan2_64(math.inf(f64), math.inf(f64)), math.pi / 4.0, epsilon));
288 try expect(math.approxEqAbs(f64, atan2_64(-math.inf(f64), math.inf(f64)), -math.pi / 4.0, epsilon));
289 try expect(math.approxEqAbs(f64, atan2_64(math.inf(f64), -math.inf(f64)), 3.0 * math.pi / 4.0, epsilon));
290 try expect(math.approxEqAbs(f64, atan2_64(-math.inf(f64), -math.inf(f64)), -3.0 * math.pi / 4.0, epsilon));
291 try expect(atan2_64(1.0, math.inf(f64)) == 0.0);
292 try expect(math.approxEqAbs(f64, atan2_64(1.0, -math.inf(f64)), math.pi, epsilon));
293 try expect(math.approxEqAbs(f64, atan2_64(-1.0, -math.inf(f64)), -math.pi, epsilon));
294 try expect(math.approxEqAbs(f64, atan2_64(math.inf(f64), 1.0), math.pi / 2.0, epsilon));
295 try expect(math.approxEqAbs(f64, atan2_64(-math.inf(f64), 1.0), -math.pi / 2.0, epsilon));
296}