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}