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/complex/catanf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/complex/catan.c
6
7const std = @import("../../std.zig");
8const testing = std.testing;
9const math = std.math;
10const cmath = math.complex;
11const Complex = cmath.Complex;
12
13/// Returns the arc-tangent of z.
14pub fn atan(z: anytype) Complex(@TypeOf(z.re, z.im)) {
15 const T = @TypeOf(z.re, z.im);
16 return switch (T) {
17 f32 => atan32(z),
18 f64 => atan64(z),
19 else => @compileError("atan not implemented for " ++ @typeName(z)),
20 };
21}
22
23fn redupif32(x: f32) f32 {
24 const DP1 = 3.140625;
25 const DP2 = 9.67502593994140625e-4;
26 const DP3 = 1.509957990978376432e-7;
27
28 var t = x / math.pi;
29 if (t >= 0.0) {
30 t += 0.5;
31 } else {
32 t -= 0.5;
33 }
34
35 const u: f32 = @trunc(t);
36 return ((x - u * DP1) - u * DP2) - u * DP3;
37}
38
39fn atan32(z: Complex(f32)) Complex(f32) {
40 const x = z.re;
41 const y = z.im;
42
43 const x2 = x * x;
44 var a = 1.0 - x2 - (y * y);
45
46 var t = 0.5 * math.atan2(2.0 * x, a);
47 const w = redupif32(t);
48
49 t = y - 1.0;
50 a = x2 + t * t;
51
52 t = y + 1.0;
53 a = (x2 + (t * t)) / a;
54 return Complex(f32).init(w, 0.25 * @log(a));
55}
56
57fn redupif64(x: f64) f64 {
58 const DP1 = 3.14159265160560607910;
59 const DP2 = 1.98418714791870343106e-9;
60 const DP3 = 1.14423774522196636802e-17;
61
62 var t = x / math.pi;
63 if (t >= 0.0) {
64 t += 0.5;
65 } else {
66 t -= 0.5;
67 }
68
69 const u: f64 = @trunc(t);
70 return ((x - u * DP1) - u * DP2) - u * DP3;
71}
72
73fn atan64(z: Complex(f64)) Complex(f64) {
74 const x = z.re;
75 const y = z.im;
76
77 const x2 = x * x;
78 var a = 1.0 - x2 - (y * y);
79
80 var t = 0.5 * math.atan2(2.0 * x, a);
81 const w = redupif64(t);
82
83 t = y - 1.0;
84 a = x2 + t * t;
85
86 t = y + 1.0;
87 a = (x2 + (t * t)) / a;
88 return Complex(f64).init(w, 0.25 * @log(a));
89}
90
91test atan32 {
92 const epsilon = math.floatEps(f32);
93 const a = Complex(f32).init(5, 3);
94 const c = atan(a);
95
96 try testing.expectApproxEqAbs(1.423679, c.re, epsilon);
97 try testing.expectApproxEqAbs(0.086569, c.im, epsilon);
98}
99
100test atan64 {
101 const epsilon = math.floatEps(f64);
102 const a = Complex(f64).init(5, 3);
103 const c = atan(a);
104
105 try testing.expectApproxEqAbs(1.4236790442393028, c.re, epsilon);
106 try testing.expectApproxEqAbs(0.08656905917945844, c.im, epsilon);
107}