Merge pull request #949 from zig-lang/complex-math

Add initial complex-number support

Author: andrewrk

CMakeLists.txt

@@ -498,6 +498,28 @@ set(ZIG_STD_FILES
     "math/tan.zig"
     "math/tanh.zig"
     "math/trunc.zig"
+    "math/complex/abs.zig"
+    "math/complex/acosh.zig"
+    "math/complex/acos.zig"
+    "math/complex/arg.zig"
+    "math/complex/asinh.zig"
+    "math/complex/asin.zig"
+    "math/complex/atanh.zig"
+    "math/complex/atan.zig"
+    "math/complex/conj.zig"
+    "math/complex/cosh.zig"
+    "math/complex/cos.zig"
+    "math/complex/exp.zig"
+    "math/complex/index.zig"
+    "math/complex/ldexp.zig"
+    "math/complex/log.zig"
+    "math/complex/pow.zig"
+    "math/complex/proj.zig"
+    "math/complex/sinh.zig"
+    "math/complex/sin.zig"
+    "math/complex/sqrt.zig"
+    "math/complex/tanh.zig"
+    "math/complex/tan.zig"
     "mem.zig"
     "net.zig"
     "os/child_process.zig"

std/math/complex/abs.zig

@@ -0,0 +1,18 @@
+const std = @import("../../index.zig");
+const debug = std.debug;
+const math = std.math;
+const cmath = math.complex;
+const Complex = cmath.Complex;
+
+pub fn abs(z: var) @typeOf(z.re) {
+    const T = @typeOf(z.re);
+    return math.hypot(T, z.re, z.im);
+}
+
+const epsilon = 0.0001;
+
+test "complex.cabs" {
+    const a = Complex(f32).new(5, 3);
+    const c = abs(a);
+    debug.assert(math.approxEq(f32, c, 5.83095, epsilon));
+}

std/math/complex/acos.zig

@@ -0,0 +1,21 @@
+const std = @import("../../index.zig");
+const debug = std.debug;
+const math = std.math;
+const cmath = math.complex;
+const Complex = cmath.Complex;
+
+pub fn acos(z: var) Complex(@typeOf(z.re)) {
+    const T = @typeOf(z.re);
+    const q = cmath.asin(z);
+    return Complex(T).new(T(math.pi) / 2 - q.re, -q.im);
+}
+
+const epsilon = 0.0001;
+
+test "complex.cacos" {
+    const a = Complex(f32).new(5, 3);
+    const c = acos(a);
+
+    debug.assert(math.approxEq(f32, c.re, 0.546975, epsilon));
+    debug.assert(math.approxEq(f32, c.im, -2.452914, epsilon));
+}

std/math/complex/acosh.zig

@@ -0,0 +1,21 @@
+const std = @import("../../index.zig");
+const debug = std.debug;
+const math = std.math;
+const cmath = math.complex;
+const Complex = cmath.Complex;
+
+pub fn acosh(z: var) Complex(@typeOf(z.re)) {
+    const T = @typeOf(z.re);
+    const q = cmath.acos(z);
+    return Complex(T).new(-q.im, q.re);
+}
+
+const epsilon = 0.0001;
+
+test "complex.cacosh" {
+    const a = Complex(f32).new(5, 3);
+    const c = acosh(a);
+
+    debug.assert(math.approxEq(f32, c.re, 2.452914, epsilon));
+    debug.assert(math.approxEq(f32, c.im, 0.546975, epsilon));
+}

std/math/complex/arg.zig

@@ -0,0 +1,18 @@
+const std = @import("../../index.zig");
+const debug = std.debug;
+const math = std.math;
+const cmath = math.complex;
+const Complex = cmath.Complex;
+
+pub fn arg(z: var) @typeOf(z.re) {
+    const T = @typeOf(z.re);
+    return math.atan2(T, z.im, z.re);
+}
+
+const epsilon = 0.0001;
+
+test "complex.carg" {
+    const a = Complex(f32).new(5, 3);
+    const c = arg(a);
+    debug.assert(math.approxEq(f32, c, 0.540420, epsilon));
+}

std/math/complex/asin.zig

@@ -0,0 +1,27 @@
+const std = @import("../../index.zig");
+const debug = std.debug;
+const math = std.math;
+const cmath = math.complex;
+const Complex = cmath.Complex;
+
+pub fn asin(z: var) Complex(@typeOf(z.re)) {
+    const T = @typeOf(z.re);
+    const x = z.re;
+    const y = z.im;
+
+    const p = Complex(T).new(1.0 - (x - y) * (x + y), -2.0 * x * y);
+    const q = Complex(T).new(-y, x);
+    const r = cmath.log(q.add(cmath.sqrt(p)));
+
+    return Complex(T).new(r.im, -r.re);
+}
+
+const epsilon = 0.0001;
+
+test "complex.casin" {
+    const a = Complex(f32).new(5, 3);
+    const c = asin(a);
+
+    debug.assert(math.approxEq(f32, c.re, 1.023822, epsilon));
+    debug.assert(math.approxEq(f32, c.im, 2.452914, epsilon));
+}

std/math/complex/asinh.zig

@@ -0,0 +1,22 @@
+const std = @import("../../index.zig");
+const debug = std.debug;
+const math = std.math;
+const cmath = math.complex;
+const Complex = cmath.Complex;
+
+pub fn asinh(z: var) Complex(@typeOf(z.re)) {
+    const T = @typeOf(z.re);
+    const q = Complex(T).new(-z.im, z.re);
+    const r = cmath.asin(q);
+    return Complex(T).new(r.im, -r.re);
+}
+
+const epsilon = 0.0001;
+
+test "complex.casinh" {
+    const a = Complex(f32).new(5, 3);
+    const c = asinh(a);
+
+    debug.assert(math.approxEq(f32, c.re, 2.459831, epsilon));
+    debug.assert(math.approxEq(f32, c.im, 0.533999, epsilon));
+}

std/math/complex/atan.zig

@@ -0,0 +1,130 @@
+const std = @import("../../index.zig");
+const debug = std.debug;
+const math = std.math;
+const cmath = math.complex;
+const Complex = cmath.Complex;
+
+pub fn atan(z: var) Complex(@typeOf(z.re)) {
+    const T = @typeOf(z.re);
+    return switch (T) {
+        f32 => atan32(z),
+        f64 => atan64(z),
+        else => @compileError("atan not implemented for " ++ @typeName(z)),
+    };
+}
+
+fn redupif32(x: f32) f32 {
+    const DP1 = 3.140625;
+    const DP2 = 9.67502593994140625e-4;
+    const DP3 = 1.509957990978376432e-7;
+
+    var t = x / math.pi;
+    if (t >= 0.0) {
+        t += 0.5;
+    } else {
+        t -= 0.5;
+    }
+
+    const u = f32(i32(t));
+    return ((x - u * DP1) - u * DP2) - t * DP3;
+}
+
+fn atan32(z: &const Complex(f32)) Complex(f32) {
+    const maxnum = 1.0e38;
+
+    const x = z.re;
+    const y = z.im;
+
+    if ((x == 0.0) and (y > 1.0)) {
+        // overflow
+        return Complex(f32).new(maxnum, maxnum);
+    }
+
+    const x2 = x * x;
+    var a = 1.0 - x2 - (y * y);
+    if (a == 0.0) {
+        // overflow
+        return Complex(f32).new(maxnum, maxnum);
+    }
+
+    var t = 0.5 * math.atan2(f32, 2.0 * x, a);
+    var w = redupif32(t);
+
+    t = y - 1.0;
+    a = x2 + t * t;
+    if (a == 0.0) {
+        // overflow
+        return Complex(f32).new(maxnum, maxnum);
+    }
+
+    t = y + 1.0;
+    a = (x2 + (t * t)) / a;
+    return Complex(f32).new(w, 0.25 * math.ln(a));
+}
+
+fn redupif64(x: f64) f64 {
+    const DP1 = 3.14159265160560607910;
+    const DP2 = 1.98418714791870343106e-9;
+    const DP3 = 1.14423774522196636802e-17;
+
+    var t = x / math.pi;
+    if (t >= 0.0) {
+        t += 0.5;
+    } else {
+        t -= 0.5;
+    }
+
+    const u = f64(i64(t));
+    return ((x - u * DP1) - u * DP2) - t * DP3;
+}
+
+fn atan64(z: &const Complex(f64)) Complex(f64) {
+    const maxnum = 1.0e308;
+
+    const x = z.re;
+    const y = z.im;
+
+    if ((x == 0.0) and (y > 1.0)) {
+        // overflow
+        return Complex(f64).new(maxnum, maxnum);
+    }
+
+    const x2 = x * x;
+    var a = 1.0 - x2 - (y * y);
+    if (a == 0.0) {
+        // overflow
+        return Complex(f64).new(maxnum, maxnum);
+    }
+
+    var t = 0.5 * math.atan2(f64, 2.0 * x, a);
+    var w = redupif64(t);
+
+    t = y - 1.0;
+    a = x2 + t * t;
+    if (a == 0.0) {
+        // overflow
+        return Complex(f64).new(maxnum, maxnum);
+    }
+
+    t = y + 1.0;
+    a = (x2 + (t * t)) / a;
+    return Complex(f64).new(w, 0.25 * math.ln(a));
+}
+
+const epsilon = 0.0001;
+
+test "complex.catan32" {
+    const a = Complex(f32).new(5, 3);
+    const c = atan(a);
+
+    debug.assert(math.approxEq(f32, c.re, 1.423679, epsilon));
+    debug.assert(math.approxEq(f32, c.im, 0.086569, epsilon));
+}
+
+test "complex.catan64" {
+    const a = Complex(f64).new(5, 3);
+    const c = atan(a);
+
+    debug.assert(math.approxEq(f64, c.re, 1.423679, epsilon));
+    debug.assert(math.approxEq(f64, c.im, 0.086569, epsilon));
+}

std/math/complex/atanh.zig

@@ -0,0 +1,22 @@
+const std = @import("../../index.zig");
+const debug = std.debug;
+const math = std.math;
+const cmath = math.complex;
+const Complex = cmath.Complex;
+
+pub fn atanh(z: var) Complex(@typeOf(z.re)) {
+    const T = @typeOf(z.re);
+    const q = Complex(T).new(-z.im, z.re);
+    const r = cmath.atan(q);
+    return Complex(T).new(r.im, -r.re);
+}
+
+const epsilon = 0.0001;
+
+test "complex.catanh" {
+    const a = Complex(f32).new(5, 3);
+    const c = atanh(a);
+
+    debug.assert(math.approxEq(f32, c.re, 0.146947, epsilon));
+    debug.assert(math.approxEq(f32, c.im, 1.480870, epsilon));
+}

std/math/complex/conj.zig

@@ -0,0 +1,17 @@
+const std = @import("../../index.zig");
+const debug = std.debug;
+const math = std.math;
+const cmath = math.complex;
+const Complex = cmath.Complex;
+
+pub fn conj(z: var) Complex(@typeOf(z.re)) {
+    const T = @typeOf(z.re);
+    return Complex(T).new(z.re, -z.im);
+}
+
+test "complex.conj" {
+    const a = Complex(f32).new(5, 3);
+    const c = a.conjugate();
+
+    debug.assert(c.re == 5 and c.im == -3);
+}

std/math/complex/cos.zig

@@ -0,0 +1,21 @@
+const std = @import("../../index.zig");
+const debug = std.debug;
+const math = std.math;
+const cmath = math.complex;
+const Complex = cmath.Complex;
+
+pub fn cos(z: var) Complex(@typeOf(z.re)) {
+    const T = @typeOf(z.re);
+    const p = Complex(T).new(-z.im, z.re);
+    return cmath.cosh(p);
+}
+
+const epsilon = 0.0001;
+
+test "complex.ccos" {
+    const a = Complex(f32).new(5, 3);
+    const c = cos(a);
+
+    debug.assert(math.approxEq(f32, c.re, 2.855815, epsilon));
+    debug.assert(math.approxEq(f32, c.im, 9.606383, epsilon));
+}