Merge pull request #1 from zig-lang/master

Sync with zig-lang/zig master

28 files changed, 1461 insertions, 283 deletions

diff --git a/std/math/complex/abs.zig b/std/math/complex/abs.zig
new file mode 100644
index 0000000000..4cd095c46b
--- /dev/null
+++ b/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));
+}
diff --git a/std/math/complex/acos.zig b/std/math/complex/acos.zig
new file mode 100644
index 0000000000..a5760b4ace
--- /dev/null
+++ b/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));
+}
diff --git a/std/math/complex/acosh.zig b/std/math/complex/acosh.zig
new file mode 100644
index 0000000000..8dd91b2836
--- /dev/null
+++ b/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));
+}
diff --git a/std/math/complex/arg.zig b/std/math/complex/arg.zig
new file mode 100644
index 0000000000..f24512ac73
--- /dev/null
+++ b/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));
+}
diff --git a/std/math/complex/asin.zig b/std/math/complex/asin.zig
new file mode 100644
index 0000000000..584a3a1a9b
--- /dev/null
+++ b/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));
+}
diff --git a/std/math/complex/asinh.zig b/std/math/complex/asinh.zig
new file mode 100644
index 0000000000..0c4dc2b6e4
--- /dev/null
+++ b/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));
+}
diff --git a/std/math/complex/atan.zig b/std/math/complex/atan.zig
new file mode 100644
index 0000000000..b7bbf930eb
--- /dev/null
+++ b/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));
+}
diff --git a/std/math/complex/atanh.zig b/std/math/complex/atanh.zig
new file mode 100644
index 0000000000..f70c741765
--- /dev/null
+++ b/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));
+}
diff --git a/std/math/complex/conj.zig b/std/math/complex/conj.zig
new file mode 100644
index 0000000000..ad3e8b5036
--- /dev/null
+++ b/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);
+}
diff --git a/std/math/complex/cos.zig b/std/math/complex/cos.zig
new file mode 100644
index 0000000000..96e4ffcdb0
--- /dev/null
+++ b/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));
+}
diff --git a/std/math/complex/cosh.zig b/std/math/complex/cosh.zig
new file mode 100644
index 0000000000..96eac68556
--- /dev/null
+++ b/std/math/complex/cosh.zig
@@ -0,0 +1,165 @@
+const std = @import("../../index.zig");
+const debug = std.debug;
+const math = std.math;
+const cmath = math.complex;
+const Complex = cmath.Complex;
+
+const ldexp_cexp = @import("ldexp.zig").ldexp_cexp;
+
+pub fn cosh(z: var) Complex(@typeOf(z.re)) {
+    const T = @typeOf(z.re);
+    return switch (T) {
+        f32 => cosh32(z),
+        f64 => cosh64(z),
+        else => @compileError("cosh not implemented for " ++ @typeName(z)),
+    };
+}
+
+fn cosh32(z: &const Complex(f32)) Complex(f32) {
+    const x = z.re;
+    const y = z.im;
+
+    const hx = @bitCast(u32, x);
+    const ix = hx & 0x7fffffff;
+
+    const hy = @bitCast(u32, y);
+    const iy = hy & 0x7fffffff;
+
+    if (ix < 0x7f800000 and iy < 0x7f800000) {
+        if (iy == 0) {
+            return Complex(f32).new(math.cosh(x), y);
+        }
+        // small x: normal case
+        if (ix < 0x41100000) {
+            return Complex(f32).new(math.cosh(x) * math.cos(y), math.sinh(x) * math.sin(y));
+        }
+
+        // |x|>= 9, so cosh(x) ~= exp(|x|)
+        if (ix < 0x42b17218) {
+            // x < 88.7: exp(|x|) won't overflow
+            const h = math.exp(math.fabs(x)) * 0.5;
+            return Complex(f32).new(math.copysign(f32, h, x) * math.cos(y), h * math.sin(y));
+        }
+        // x < 192.7: scale to avoid overflow
+        else if (ix < 0x4340b1e7) {
+            const v = Complex(f32).new(math.fabs(x), y);
+            const r = ldexp_cexp(v, -1);
+            return Complex(f32).new(x, y * math.copysign(f32, 1, x));
+        }
+        // x >= 192.7: result always overflows
+        else {
+            const h = 0x1p127 * x;
+            return Complex(f32).new(h * h * math.cos(y), h * math.sin(y));
+        }
+    }
+
+    if (ix == 0 and iy >= 0x7f800000) {
+        return Complex(f32).new(y - y, math.copysign(f32, 0, x * (y - y)));
+    }
+
+    if (iy == 0 and ix >= 0x7f800000) {
+        if (hx & 0x7fffff == 0) {
+            return Complex(f32).new(x * x, math.copysign(f32, 0, x) * y);
+        }
+        return Complex(f32).new(x, math.copysign(f32, 0, (x + x) * y));
+    }
+
+    if (ix < 0x7f800000 and iy >= 0x7f800000) {
+        return Complex(f32).new(y - y, x * (y - y));
+    }
+
+    if (ix >= 0x7f800000 and (hx & 0x7fffff) == 0) {
+        if (iy >= 0x7f800000) {
+            return Complex(f32).new(x * x, x * (y - y));
+        }
+        return Complex(f32).new((x * x) * math.cos(y), x * math.sin(y));
+    }
+
+    return Complex(f32).new((x * x) * (y - y), (x + x) * (y - y));
+}
+
+fn cosh64(z: &const Complex(f64)) Complex(f64) {
+    const x = z.re;
+    const y = z.im;
+
+    const fx = @bitCast(u64, x);
+    const hx = u32(fx >> 32);
+    const lx = @truncate(u32, fx);
+    const ix = hx & 0x7fffffff;
+
+    const fy = @bitCast(u64, y);
+    const hy = u32(fy >> 32);
+    const ly = @truncate(u32, fy);
+    const iy = hy & 0x7fffffff;
+
+    // nearly non-exceptional case where x, y are finite
+    if (ix < 0x7ff00000 and iy < 0x7ff00000) {
+        if (iy | ly == 0) {
+            return Complex(f64).new(math.cosh(x), x * y);
+        }
+        // small x: normal case
+        if (ix < 0x40360000) {
+            return Complex(f64).new(math.cosh(x) * math.cos(y), math.sinh(x) * math.sin(y));
+        }
+
+        // |x|>= 22, so cosh(x) ~= exp(|x|)
+        if (ix < 0x40862e42) {
+            // x < 710: exp(|x|) won't overflow
+            const h = math.exp(math.fabs(x)) * 0.5;
+            return Complex(f64).new(h * math.cos(y), math.copysign(f64, h, x) * math.sin(y));
+        }
+        // x < 1455: scale to avoid overflow
+        else if (ix < 0x4096bbaa) {
+            const v = Complex(f64).new(math.fabs(x), y);
+            const r = ldexp_cexp(v, -1);
+            return Complex(f64).new(x, y * math.copysign(f64, 1, x));
+        }
+        // x >= 1455: result always overflows
+        else {
+            const h = 0x1p1023;
+            return Complex(f64).new(h * h * math.cos(y), h * math.sin(y));
+        }
+    }
+
+    if (ix | lx == 0 and iy >= 0x7ff00000) {
+        return Complex(f64).new(y - y, math.copysign(f64, 0, x * (y - y)));
+    }
+
+    if (iy | ly == 0 and ix >= 0x7ff00000) {
+        if ((hx & 0xfffff) | lx == 0) {
+            return Complex(f64).new(x * x, math.copysign(f64, 0, x) * y);
+        }
+        return Complex(f64).new(x * x, math.copysign(f64, 0, (x + x) * y));
+    }
+
+    if (ix < 0x7ff00000 and iy >= 0x7ff00000) {
+        return Complex(f64).new(y - y, x * (y - y));
+    }
+
+    if (ix >= 0x7ff00000 and (hx & 0xfffff) | lx == 0) {
+        if (iy >= 0x7ff00000) {
+            return Complex(f64).new(x * x, x * (y - y));
+        }
+        return Complex(f64).new(x * x * math.cos(y), x * math.sin(y));
+    }
+
+    return Complex(f64).new((x * x) * (y - y), (x + x) * (y - y));
+}
+
+const epsilon = 0.0001;
+
+test "complex.ccosh32" {
+    const a = Complex(f32).new(5, 3);
+    const c = cosh(a);
+
+    debug.assert(math.approxEq(f32, c.re, -73.467300, epsilon));
+    debug.assert(math.approxEq(f32, c.im, 10.471557, epsilon));
+}
+
+test "complex.ccosh64" {
+    const a = Complex(f64).new(5, 3);
+    const c = cosh(a);
+
+    debug.assert(math.approxEq(f64, c.re, -73.467300, epsilon));
+    debug.assert(math.approxEq(f64, c.im, 10.471557, epsilon));
+}
diff --git a/std/math/complex/exp.zig b/std/math/complex/exp.zig
new file mode 100644
index 0000000000..03f7f9e41b
--- /dev/null
+++ b/std/math/complex/exp.zig
@@ -0,0 +1,140 @@
+const std = @import("../../index.zig");
+const debug = std.debug;
+const math = std.math;
+const cmath = math.complex;
+const Complex = cmath.Complex;
+
+const ldexp_cexp = @import("ldexp.zig").ldexp_cexp;
+
+pub fn exp(z: var) Complex(@typeOf(z.re)) {
+    const T = @typeOf(z.re);
+
+    return switch (T) {
+        f32 => exp32(z),
+        f64 => exp64(z),
+        else => @compileError("exp not implemented for " ++ @typeName(z)),
+    };
+}
+
+fn exp32(z: &const Complex(f32)) Complex(f32) {
+    @setFloatMode(this, @import("builtin").FloatMode.Strict);
+
+    const exp_overflow  = 0x42b17218;    // max_exp * ln2 ~= 88.72283955
+    const cexp_overflow = 0x43400074;   // (max_exp - min_denom_exp) * ln2
+
+    const x = z.re;
+    const y = z.im;
+
+    const hy = @bitCast(u32, y) & 0x7fffffff;
+    // cexp(x + i0) = exp(x) + i0
+    if (hy == 0) {
+        return Complex(f32).new(math.exp(x), y);
+    }
+
+    const hx = @bitCast(u32, x);
+    // cexp(0 + iy) = cos(y) + isin(y)
+    if ((hx & 0x7fffffff) == 0) {
+        return Complex(f32).new(math.cos(y), math.sin(y));
+    }
+
+    if (hy >= 0x7f800000) {
+        // cexp(finite|nan +- i inf|nan) = nan + i nan
+        if ((hx & 0x7fffffff) != 0x7f800000) {
+            return Complex(f32).new(y - y, y - y);
+        }
+        // cexp(-inf +- i inf|nan) = 0 + i0
+        else if (hx & 0x80000000 != 0) {
+            return Complex(f32).new(0, 0);
+        }
+        // cexp(+inf +- i inf|nan) = inf + i nan
+        else {
+            return Complex(f32).new(x, y - y);
+        }
+    }
+
+    // 88.7 <= x <= 192 so must scale
+    if (hx >= exp_overflow and hx <= cexp_overflow) {
+        return ldexp_cexp(z, 0);
+    }
+    // - x < exp_overflow => exp(x) won't overflow (common)
+    // - x > cexp_overflow, so exp(x) * s overflows for s > 0
+    // - x = +-inf
+    // - x = nan
+    else {
+        const exp_x = math.exp(x);
+        return Complex(f32).new(exp_x * math.cos(y), exp_x * math.sin(y));
+    }
+}
+
+fn exp64(z: &const Complex(f64)) Complex(f64) {
+    const exp_overflow  = 0x40862e42;    // high bits of max_exp * ln2 ~= 710
+    const cexp_overflow = 0x4096b8e4;   // (max_exp - min_denorm_exp) * ln2
+
+    const x = z.re;
+    const y = z.im;
+
+    const fy = @bitCast(u64, y);
+    const hy = u32(fy >> 32) & 0x7fffffff;
+    const ly = @truncate(u32, fy);
+
+    // cexp(x + i0) = exp(x) + i0
+    if (hy | ly == 0) {
+        return Complex(f64).new(math.exp(x), y);
+    }
+
+    const fx = @bitCast(u64, x);
+    const hx = u32(fx >> 32);
+    const lx = @truncate(u32, fx);
+
+    // cexp(0 + iy) = cos(y) + isin(y)
+    if ((hx & 0x7fffffff) | lx == 0) {
+        return Complex(f64).new(math.cos(y), math.sin(y));
+    }
+
+    if (hy >= 0x7ff00000) {
+        // cexp(finite|nan +- i inf|nan) = nan + i nan
+        if (lx != 0 or (hx & 0x7fffffff) != 0x7ff00000) {
+            return Complex(f64).new(y - y, y - y);
+        }
+        // cexp(-inf +- i inf|nan) = 0 + i0
+        else if (hx & 0x80000000 != 0) {
+            return Complex(f64).new(0, 0);
+        }
+        // cexp(+inf +- i inf|nan) = inf + i nan
+        else {
+            return Complex(f64).new(x, y - y);
+        }
+    }
+
+    // 709.7 <= x <= 1454.3 so must scale
+    if (hx >= exp_overflow and hx <= cexp_overflow) {
+        const r = ldexp_cexp(z, 0);
+        return *r;
+    }
+    // - x < exp_overflow => exp(x) won't overflow (common)
+    // - x > cexp_overflow, so exp(x) * s overflows for s > 0
+    // - x = +-inf
+    // - x = nan
+    else {
+        const exp_x = math.exp(x);
+        return Complex(f64).new(exp_x * math.cos(y), exp_x * math.sin(y));
+    }
+}
+
+const epsilon = 0.0001;
+
+test "complex.cexp32" {
+    const a = Complex(f32).new(5, 3);
+    const c = exp(a);
+
+    debug.assert(math.approxEq(f32, c.re, -146.927917, epsilon));
+    debug.assert(math.approxEq(f32, c.im, 20.944065, epsilon));
+}
+
+test "complex.cexp64" {
+    const a = Complex(f32).new(5, 3);
+    const c = exp(a);
+
+    debug.assert(math.approxEq(f64, c.re, -146.927917, epsilon));
+    debug.assert(math.approxEq(f64, c.im, 20.944065, epsilon));
+}
diff --git a/std/math/complex/index.zig b/std/math/complex/index.zig
new file mode 100644
index 0000000000..a4d493307e
--- /dev/null
+++ b/std/math/complex/index.zig
@@ -0,0 +1,171 @@
+const std = @import("../../index.zig");
+const debug = std.debug;
+const math = std.math;
+
+pub const abs = @import("abs.zig").abs;
+pub const acosh = @import("acosh.zig").acosh;
+pub const acos = @import("acos.zig").acos;
+pub const arg = @import("arg.zig").arg;
+pub const asinh = @import("asinh.zig").asinh;
+pub const asin = @import("asin.zig").asin;
+pub const atanh = @import("atanh.zig").atanh;
+pub const atan = @import("atan.zig").atan;
+pub const conj = @import("conj.zig").conj;
+pub const cosh = @import("cosh.zig").cosh;
+pub const cos = @import("cos.zig").cos;
+pub const exp = @import("exp.zig").exp;
+pub const log = @import("log.zig").log;
+pub const pow = @import("pow.zig").pow;
+pub const proj = @import("proj.zig").proj;
+pub const sinh = @import("sinh.zig").sinh;
+pub const sin = @import("sin.zig").sin;
+pub const sqrt = @import("sqrt.zig").sqrt;
+pub const tanh = @import("tanh.zig").tanh;
+pub const tan = @import("tan.zig").tan;
+
+pub fn Complex(comptime T: type) type {
+    return struct {
+        const Self = this;
+
+        re: T,
+        im: T,
+
+        pub fn new(re: T, im: T) Self {
+            return Self {
+                .re = re,
+                .im = im,
+            };
+        }
+
+        pub fn add(self: &const Self, other: &const Self) Self {
+            return Self {
+                .re = self.re + other.re,
+                .im = self.im + other.im,
+            };
+        }
+
+        pub fn sub(self: &const Self, other: &const Self) Self {
+            return Self {
+                .re = self.re - other.re,
+                .im = self.im - other.im,
+            };
+        }
+
+        pub fn mul(self: &const Self, other: &const Self) Self {
+            return Self {
+                .re = self.re * other.re - self.im * other.im,
+                .im = self.im * other.re + self.re * other.im,
+            };
+        }
+
+        pub fn div(self: &const Self, other: &const Self) Self {
+            const re_num = self.re * other.re + self.im * other.im;
+            const im_num = self.im * other.re - self.re * other.im;
+            const den = other.re * other.re + other.im * other.im;
+
+            return Self {
+                .re = re_num / den,
+                .im = im_num / den,
+            };
+        }
+
+        pub fn conjugate(self: &const Self) Self {
+            return Self {
+                .re = self.re,
+                .im = -self.im,
+            };
+        }
+
+        pub fn reciprocal(self: &const Self) Self {
+            const m = self.re * self.re + self.im * self.im;
+            return Self {
+                .re = self.re / m,
+                .im = -self.im / m,
+            };
+        }
+
+        pub fn magnitude(self: &const Self) T {
+            return math.sqrt(self.re * self.re + self.im * self.im);
+        }
+    };
+}
+
+const epsilon = 0.0001;
+
+test "complex.add" {
+    const a = Complex(f32).new(5, 3);
+    const b = Complex(f32).new(2, 7);
+    const c = a.add(b);
+
+    debug.assert(c.re == 7 and c.im == 10);
+}
+
+test "complex.sub" {
+    const a = Complex(f32).new(5, 3);
+    const b = Complex(f32).new(2, 7);
+    const c = a.sub(b);
+
+    debug.assert(c.re == 3 and c.im == -4);
+}
+
+test "complex.mul" {
+    const a = Complex(f32).new(5, 3);
+    const b = Complex(f32).new(2, 7);
+    const c = a.mul(b);
+
+    debug.assert(c.re == -11 and c.im == 41);
+}
+
+test "complex.div" {
+    const a = Complex(f32).new(5, 3);
+    const b = Complex(f32).new(2, 7);
+    const c = a.div(b);
+
+    debug.assert(math.approxEq(f32, c.re, f32(31)/53, epsilon) and
+                 math.approxEq(f32, c.im, f32(-29)/53, epsilon));
+}
+
+test "complex.conjugate" {
+    const a = Complex(f32).new(5, 3);
+    const c = a.conjugate();
+
+    debug.assert(c.re == 5 and c.im == -3);
+}
+
+test "complex.reciprocal" {
+    const a = Complex(f32).new(5, 3);
+    const c = a.reciprocal();
+
+    debug.assert(math.approxEq(f32, c.re, f32(5)/34, epsilon) and
+                 math.approxEq(f32, c.im, f32(-3)/34, epsilon));
+}
+
+test "complex.magnitude" {
+    const a = Complex(f32).new(5, 3);
+    const c = a.magnitude();
+
+    debug.assert(math.approxEq(f32, c, 5.83095, epsilon));
+}
+
+test "complex.cmath" {
+    _ = @import("abs.zig");
+    _ = @import("acosh.zig");
+    _ = @import("acos.zig");
+    _ = @import("arg.zig");
+    _ = @import("asinh.zig");
+    _ = @import("asin.zig");
+    _ = @import("atanh.zig");
+    _ = @import("atan.zig");
+    _ = @import("conj.zig");
+    _ = @import("cosh.zig");
+    _ = @import("cos.zig");
+    _ = @import("exp.zig");
+    _ = @import("log.zig");
+    _ = @import("pow.zig");
+    _ = @import("proj.zig");
+    _ = @import("sinh.zig");
+    _ = @import("sin.zig");
+    _ = @import("sqrt.zig");
+    _ = @import("tanh.zig");
+    _ = @import("tan.zig");
+}
diff --git a/std/math/complex/ldexp.zig b/std/math/complex/ldexp.zig
new file mode 100644
index 0000000000..4fb5a6815f
--- /dev/null
+++ b/std/math/complex/ldexp.zig
@@ -0,0 +1,75 @@
+const std = @import("../../index.zig");
+const debug = std.debug;
+const math = std.math;
+const cmath = math.complex;
+const Complex = cmath.Complex;
+
+pub fn ldexp_cexp(z: var, expt: i32) Complex(@typeOf(z.re)) {
+    const T = @typeOf(z.re);
+
+    return switch (T) {
+        f32 => ldexp_cexp32(z, expt),
+        f64 => ldexp_cexp64(z, expt),
+        else => unreachable,
+    };
+}
+
+fn frexp_exp32(x: f32, expt: &i32) f32 {
+    const k = 235;              // reduction constant
+    const kln2 = 162.88958740;  // k * ln2
+
+    const exp_x = math.exp(x - kln2);
+    const hx = @bitCast(u32, exp_x);
+    *expt = i32(hx >> 23) - (0x7f + 127) + k;
+    return @bitCast(f32, (hx & 0x7fffff) | ((0x7f + 127) << 23));
+}
+
+fn ldexp_cexp32(z: &const Complex(f32), expt: i32) Complex(f32) {
+    var ex_expt: i32 = undefined;
+    const exp_x = frexp_exp32(z.re, &ex_expt);
+    const exptf = expt + ex_expt;
+
+    const half_expt1 = @divTrunc(exptf, 2);
+    const scale1 = @bitCast(f32, (0x7f + half_expt1) << 23);
+
+    const half_expt2 = exptf - half_expt1;
+    const scale2 = @bitCast(f32, (0x7f + half_expt2) << 23);
+
+    return Complex(f32).new(
+        math.cos(z.im) * exp_x * scale1 * scale2,
+        math.sin(z.im) * exp_x * scale1 * scale2,
+    );
+}
+
+fn frexp_exp64(x: f64, expt: &i32) f64 {
+    const k = 1799;                         // reduction constant
+    const kln2 = 1246.97177782734161156;    // k * ln2
+
+    const exp_x = math.exp(x - kln2);
+
+    const fx = @bitCast(u64, x);
+    const hx = u32(fx >> 32);
+    const lx = @truncate(u32, fx);
+
+    *expt = i32(hx >> 20) - (0x3ff + 1023) + k;
+
+    const high_word = (hx & 0xfffff) | ((0x3ff + 1023) << 20);
+    return @bitCast(f64, (u64(high_word) << 32) | lx);
+}
+
+fn ldexp_cexp64(z: &const Complex(f64), expt: i32) Complex(f64) {
+    var ex_expt: i32 = undefined;
+    const exp_x = frexp_exp64(z.re, &ex_expt);
+    const exptf = i64(expt + ex_expt);
+
+    const half_expt1 = @divTrunc(exptf, 2);
+    const scale1 = @bitCast(f64, (0x3ff + half_expt1) << 20);
+
+    const half_expt2 = exptf - half_expt1;
+    const scale2 = @bitCast(f64, (0x3ff + half_expt2) << 20);
+
+    return Complex(f64).new(
+        math.cos(z.im) * exp_x * scale1 * scale2,
+        math.sin(z.im) * exp_x * scale1 * scale2,
+    );
+}
diff --git a/std/math/complex/log.zig b/std/math/complex/log.zig
new file mode 100644
index 0000000000..a4a1d1664f
--- /dev/null
+++ b/std/math/complex/log.zig
@@ -0,0 +1,23 @@
+const std = @import("../../index.zig");
+const debug = std.debug;
+const math = std.math;
+const cmath = math.complex;
+const Complex = cmath.Complex;
+
+pub fn log(z: var) Complex(@typeOf(z.re)) {
+    const T = @typeOf(z.re);
+    const r = cmath.abs(z);
+    const phi = cmath.arg(z);
+
+    return Complex(T).new(math.ln(r), phi);
+}
+
+const epsilon = 0.0001;
+
+test "complex.clog" {
+    const a = Complex(f32).new(5, 3);
+    const c = log(a);
+
+    debug.assert(math.approxEq(f32, c.re, 1.763180, epsilon));
+    debug.assert(math.approxEq(f32, c.im, 0.540419, epsilon));
+}
diff --git a/std/math/complex/pow.zig b/std/math/complex/pow.zig
new file mode 100644
index 0000000000..bef9fde542
--- /dev/null
+++ b/std/math/complex/pow.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 pow(comptime T: type, z: &const T, c: &const T) T {
+    const p = cmath.log(z);
+    const q = c.mul(p);
+    return cmath.exp(q);
+}
+
+const epsilon = 0.0001;
+
+test "complex.cpow" {
+    const a = Complex(f32).new(5, 3);
+    const b = Complex(f32).new(2.3, -1.3);
+    const c = pow(Complex(f32), a, b);
+
+    debug.assert(math.approxEq(f32, c.re, 58.049110, epsilon));
+    debug.assert(math.approxEq(f32, c.im, -101.003433, epsilon));
+}
diff --git a/std/math/complex/proj.zig b/std/math/complex/proj.zig
new file mode 100644
index 0000000000..b6c4cc046e
--- /dev/null
+++ b/std/math/complex/proj.zig
@@ -0,0 +1,24 @@
+const std = @import("../../index.zig");
+const debug = std.debug;
+const math = std.math;
+const cmath = math.complex;
+const Complex = cmath.Complex;
+
+pub fn proj(z: var) Complex(@typeOf(z.re)) {
+    const T = @typeOf(z.re);
+
+    if (math.isInf(z.re) or math.isInf(z.im)) {
+        return Complex(T).new(math.inf(T), math.copysign(T, 0, z.re));
+    }
+
+    return Complex(T).new(z.re, z.im);
+}
+
+const epsilon = 0.0001;
+
+test "complex.cproj" {
+    const a = Complex(f32).new(5, 3);
+    const c = proj(a);
+
+    debug.assert(c.re == 5 and c.im == 3);
+}
diff --git a/std/math/complex/sin.zig b/std/math/complex/sin.zig
new file mode 100644
index 0000000000..d32b771d3b
--- /dev/null
+++ b/std/math/complex/sin.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 sin(z: var) Complex(@typeOf(z.re)) {
+    const T = @typeOf(z.re);
+    const p = Complex(T).new(-z.im, z.re);
+    const q = cmath.sinh(p);
+    return Complex(T).new(q.im, -q.re);
+}
+
+const epsilon = 0.0001;
+
+test "complex.csin" {
+    const a = Complex(f32).new(5, 3);
+    const c = sin(a);
+
+    debug.assert(math.approxEq(f32, c.re, -9.654126, epsilon));
+    debug.assert(math.approxEq(f32, c.im, 2.841692, epsilon));
+}
diff --git a/std/math/complex/sinh.zig b/std/math/complex/sinh.zig
new file mode 100644
index 0000000000..09a62ca058
--- /dev/null
+++ b/std/math/complex/sinh.zig
@@ -0,0 +1,164 @@
+const std = @import("../../index.zig");
+const debug = std.debug;
+const math = std.math;
+const cmath = math.complex;
+const Complex = cmath.Complex;
+
+const ldexp_cexp = @import("ldexp.zig").ldexp_cexp;
+
+pub fn sinh(z: var) Complex(@typeOf(z.re)) {
+    const T = @typeOf(z.re);
+    return switch (T) {
+        f32 => sinh32(z),
+        f64 => sinh64(z),
+        else => @compileError("tan not implemented for " ++ @typeName(z)),
+    };
+}
+
+fn sinh32(z: &const Complex(f32)) Complex(f32) {
+    const x = z.re;
+    const y = z.im;
+
+    const hx = @bitCast(u32, x);
+    const ix = hx & 0x7fffffff;
+
+    const hy = @bitCast(u32, y);
+    const iy = hy & 0x7fffffff;
+
+    if (ix < 0x7f800000 and iy < 0x7f800000) {
+        if (iy == 0) {
+            return Complex(f32).new(math.sinh(x), y);
+        }
+        // small x: normal case
+        if (ix < 0x41100000) {
+            return Complex(f32).new(math.sinh(x) * math.cos(y), math.cosh(x) * math.sin(y));
+        }
+
+        // |x|>= 9, so cosh(x) ~= exp(|x|)
+        if (ix < 0x42b17218) {
+            // x < 88.7: exp(|x|) won't overflow
+            const h = math.exp(math.fabs(x)) * 0.5;
+            return Complex(f32).new(math.copysign(f32, h, x) * math.cos(y), h * math.sin(y));
+        }
+        // x < 192.7: scale to avoid overflow
+        else if (ix < 0x4340b1e7) {
+            const v = Complex(f32).new(math.fabs(x), y);
+            const r = ldexp_cexp(v, -1);
+            return Complex(f32).new(x * math.copysign(f32, 1, x), y);
+        }
+        // x >= 192.7: result always overflows
+        else {
+            const h = 0x1p127 * x;
+            return Complex(f32).new(h * math.cos(y), h * h * math.sin(y));
+        }
+    }
+
+    if (ix == 0 and iy >= 0x7f800000) {
+        return Complex(f32).new(math.copysign(f32, 0, x * (y - y)), y - y);
+    }
+
+    if (iy == 0 and ix >= 0x7f800000) {
+        if (hx & 0x7fffff == 0) {
+            return Complex(f32).new(x, y);
+        }
+        return Complex(f32).new(x, math.copysign(f32, 0, y));
+    }
+
+    if (ix < 0x7f800000 and iy >= 0x7f800000) {
+        return Complex(f32).new(y - y, x * (y - y));
+    }
+
+    if (ix >= 0x7f800000 and (hx & 0x7fffff) == 0) {
+        if (iy >= 0x7f800000) {
+            return Complex(f32).new(x * x, x * (y - y));
+        }
+        return Complex(f32).new(x * math.cos(y), math.inf_f32 * math.sin(y));
+    }
+
+    return Complex(f32).new((x * x) * (y - y), (x + x) * (y - y));
+}
+
+fn sinh64(z: &const Complex(f64)) Complex(f64) {
+    const x = z.re;
+    const y = z.im;
+
+    const fx = @bitCast(u64, x);
+    const hx = u32(fx >> 32);
+    const lx = @truncate(u32, fx);
+    const ix = hx & 0x7fffffff;
+
+    const fy = @bitCast(u64, y);
+    const hy = u32(fy >> 32);
+    const ly = @truncate(u32, fy);
+    const iy = hy & 0x7fffffff;
+
+    if (ix < 0x7ff00000 and iy < 0x7ff00000) {
+        if (iy | ly == 0) {
+            return Complex(f64).new(math.sinh(x), y);
+        }
+        // small x: normal case
+        if (ix < 0x40360000) {
+            return Complex(f64).new(math.sinh(x) * math.cos(y), math.cosh(x) * math.sin(y));
+        }
+
+        // |x|>= 22, so cosh(x) ~= exp(|x|)
+        if (ix < 0x40862e42) {
+            // x < 710: exp(|x|) won't overflow
+            const h = math.exp(math.fabs(x)) * 0.5;
+            return Complex(f64).new(math.copysign(f64, h, x) * math.cos(y), h * math.sin(y));
+        }
+        // x < 1455: scale to avoid overflow
+        else if (ix < 0x4096bbaa) {
+            const v = Complex(f64).new(math.fabs(x), y);
+            const r = ldexp_cexp(v, -1);
+            return Complex(f64).new(x * math.copysign(f64, 1, x), y);
+        }
+        // x >= 1455: result always overflows
+        else {
+            const h = 0x1p1023 * x;
+            return Complex(f64).new(h * math.cos(y), h * h * math.sin(y));
+        }
+    }
+
+    if (ix | lx == 0 and iy >= 0x7ff00000) {
+        return Complex(f64).new(math.copysign(f64, 0, x * (y - y)), y - y);
+    }
+
+    if (iy | ly == 0 and ix >= 0x7ff00000) {
+        if ((hx & 0xfffff) | lx == 0) {
+            return Complex(f64).new(x, y);
+        }
+        return Complex(f64).new(x, math.copysign(f64, 0, y));
+    }
+
+    if (ix < 0x7ff00000 and iy >= 0x7ff00000) {
+        return Complex(f64).new(y - y, x * (y - y));
+    }
+
+    if (ix >= 0x7ff00000 and (hx & 0xfffff) | lx == 0) {
+        if (iy >= 0x7ff00000) {
+            return Complex(f64).new(x * x, x * (y - y));
+        }
+        return Complex(f64).new(x * math.cos(y), math.inf_f64 * math.sin(y));
+    }
+
+    return Complex(f64).new((x * x) * (y - y), (x + x) * (y - y));
+}
+
+const epsilon = 0.0001;
+
+test "complex.csinh32" {
+    const a = Complex(f32).new(5, 3);
+    const c = sinh(a);
+
+    debug.assert(math.approxEq(f32, c.re, -73.460617, epsilon));
+    debug.assert(math.approxEq(f32, c.im, 10.472508, epsilon));
+}
+
+test "complex.csinh64" {
+    const a = Complex(f64).new(5, 3);
+    const c = sinh(a);
+
+    debug.assert(math.approxEq(f64, c.re, -73.460617, epsilon));
+    debug.assert(math.approxEq(f64, c.im, 10.472508, epsilon));
+}
diff --git a/std/math/complex/sqrt.zig b/std/math/complex/sqrt.zig
new file mode 100644
index 0000000000..afda69f7c9
--- /dev/null
+++ b/std/math/complex/sqrt.zig
@@ -0,0 +1,133 @@
+const std = @import("../../index.zig");
+const debug = std.debug;
+const math = std.math;
+const cmath = math.complex;
+const Complex = cmath.Complex;
+
+// TODO when #733 is solved this can be @typeOf(z) instead of Complex(@typeOf(z.re))
+pub fn sqrt(z: var) Complex(@typeOf(z.re)) {
+    const T = @typeOf(z.re);
+
+    return switch (T) {
+        f32 => sqrt32(z),
+        f64 => sqrt64(z),
+        else => @compileError("sqrt not implemented for " ++ @typeName(z)),
+    };
+}
+
+fn sqrt32(z: &const Complex(f32)) Complex(f32) {
+    const x = z.re;
+    const y = z.im;
+
+    if (x == 0 and y == 0) {
+        return Complex(f32).new(0, y);
+    }
+    if (math.isInf(y)) {
+        return Complex(f32).new(math.inf(f32), y);
+    }
+    if (math.isNan(x)) {
+        // raise invalid if y is not nan
+        const t = (y - y) / (y - y);
+        return Complex(f32).new(x, t);
+    }
+    if (math.isInf(x)) {
+        // sqrt(inf + i nan)    = inf + nan i
+        // sqrt(inf + iy)       = inf + i0
+        // sqrt(-inf + i nan)   = nan +- inf i
+        // sqrt(-inf + iy)      = 0 + inf i
+        if (math.signbit(x)) {
+            return Complex(f32).new(math.fabs(x - y), math.copysign(f32, x, y));
+        } else {
+            return Complex(f32).new(x, math.copysign(f32, y - y, y));
+        }
+    }
+
+    // y = nan special case is handled fine below
+
+    // double-precision avoids overflow with correct rounding.
+    const dx = f64(x);
+    const dy = f64(y);
+
+    if (dx >= 0) {
+        const t = math.sqrt((dx + math.hypot(f64, dx, dy)) * 0.5);
+        return Complex(f32).new(f32(t), f32(dy / (2.0 * t)));
+    } else {
+        const t = math.sqrt((-dx + math.hypot(f64, dx, dy)) * 0.5);
+        return Complex(f32).new(f32(math.fabs(y) / (2.0 * t)), f32(math.copysign(f64, t, y)));
+    }
+}
+
+fn sqrt64(z: &const Complex(f64)) Complex(f64) {
+    // may encounter overflow for im,re >= DBL_MAX / (1 + sqrt(2))
+    const threshold = 0x1.a827999fcef32p+1022;
+
+    var x = z.re;
+    var y = z.im;
+
+    if (x == 0 and y == 0) {
+        return Complex(f64).new(0, y);
+    }
+    if (math.isInf(y)) {
+        return Complex(f64).new(math.inf(f64), y);
+    }
+    if (math.isNan(x)) {
+        // raise invalid if y is not nan
+        const t = (y - y) / (y - y);
+        return Complex(f64).new(x, t);
+    }
+    if (math.isInf(x)) {
+        // sqrt(inf + i nan)    = inf + nan i
+        // sqrt(inf + iy)       = inf + i0
+        // sqrt(-inf + i nan)   = nan +- inf i
+        // sqrt(-inf + iy)      = 0 + inf i
+        if (math.signbit(x)) {
+            return Complex(f64).new(math.fabs(x - y), math.copysign(f64, x, y));
+        } else {
+            return Complex(f64).new(x, math.copysign(f64, y - y, y));
+        }
+    }
+
+    // y = nan special case is handled fine below
+
+    // scale to avoid overflow
+    var scale = false;
+    if (math.fabs(x) >= threshold or math.fabs(y) >= threshold) {
+        x *= 0.25;
+        y *= 0.25;
+        scale = true;
+    }
+
+    var result: Complex(f64) = undefined;
+    if (x >= 0) {
+        const t = math.sqrt((x + math.hypot(f64, x, y)) * 0.5);
+        result = Complex(f64).new(t, y / (2.0 * t));
+    } else {
+        const t = math.sqrt((-x + math.hypot(f64, x, y)) * 0.5);
+        result = Complex(f64).new(math.fabs(y) / (2.0 * t), math.copysign(f64, t, y));
+    }
+
+    if (scale) {
+        result.re *= 2;
+        result.im *= 2;
+    }
+
+    return result;
+}
+
+const epsilon = 0.0001;
+
+test "complex.csqrt32" {
+    const a = Complex(f32).new(5, 3);
+    const c = sqrt(a);
+
+    debug.assert(math.approxEq(f32, c.re, 2.327117, epsilon));
+    debug.assert(math.approxEq(f32, c.im, 0.644574, epsilon));
+}
+
+test "complex.csqrt64" {
+    const a = Complex(f64).new(5, 3);
+    const c = sqrt(a);
+
+    debug.assert(math.approxEq(f64, c.re, 2.3271175190399496, epsilon));
+    debug.assert(math.approxEq(f64, c.im, 0.6445742373246469, epsilon));
+}
diff --git a/std/math/complex/tan.zig b/std/math/complex/tan.zig
new file mode 100644
index 0000000000..4ea5182fa7
--- /dev/null
+++ b/std/math/complex/tan.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 tan(z: var) Complex(@typeOf(z.re)) {
+    const T = @typeOf(z.re);
+    const q = Complex(T).new(-z.im, z.re);
+    const r = cmath.tanh(q);
+    return Complex(T).new(r.im, -r.re);
+}
+
+const epsilon = 0.0001;
+
+test "complex.ctan" {
+    const a = Complex(f32).new(5, 3);
+    const c = tan(a);
+
+    debug.assert(math.approxEq(f32, c.re, -0.002708233, epsilon));
+    debug.assert(math.approxEq(f32, c.im, 1.004165, epsilon));
+}
diff --git a/std/math/complex/tanh.zig b/std/math/complex/tanh.zig
new file mode 100644
index 0000000000..6af62f48ae
--- /dev/null
+++ b/std/math/complex/tanh.zig
@@ -0,0 +1,111 @@
+const std = @import("../../index.zig");
+const debug = std.debug;
+const math = std.math;
+const cmath = math.complex;
+const Complex = cmath.Complex;
+
+pub fn tanh(z: var) Complex(@typeOf(z.re)) {
+    const T = @typeOf(z.re);
+    return switch (T) {
+        f32 => tanh32(z),
+        f64 => tanh64(z),
+        else => @compileError("tan not implemented for " ++ @typeName(z)),
+    };
+}
+
+fn tanh32(z: &const Complex(f32)) Complex(f32) {
+    const x = z.re;
+    const y = z.im;
+
+    const hx = @bitCast(u32, x);
+    const ix = hx & 0x7fffffff;
+
+    if (ix >= 0x7f800000) {
+        if (ix & 0x7fffff != 0) {
+            const r = if (y == 0) y else x * y;
+            return Complex(f32).new(x, r);
+        }
+        const xx = @bitCast(f32, hx - 0x40000000);
+        const r = if (math.isInf(y)) y else math.sin(y) * math.cos(y);
+        return Complex(f32).new(xx, math.copysign(f32, 0, r));
+    }
+
+    if (!math.isFinite(y)) {
+        const r = if (ix != 0) y - y else x;
+        return Complex(f32).new(r, y - y);
+    }
+
+    // x >= 11
+    if (ix >= 0x41300000) {
+        const exp_mx = math.exp(-math.fabs(x));
+        return Complex(f32).new(math.copysign(f32, 1, x), 4 * math.sin(y) * math.cos(y) * exp_mx * exp_mx);
+    }
+
+    // Kahan's algorithm
+    const t = math.tan(y);
+    const beta = 1.0 + t * t;
+    const s = math.sinh(x);
+    const rho = math.sqrt(1 + s * s);
+    const den = 1 + beta * s * s;
+
+    return Complex(f32).new((beta * rho * s) / den, t / den);
+}
+
+fn tanh64(z: &const Complex(f64)) Complex(f64) {
+    const x = z.re;
+    const y = z.im;
+
+    const fx = @bitCast(u64, x);
+    const hx = u32(fx >> 32);
+    const lx = @truncate(u32, fx);
+    const ix = hx & 0x7fffffff;
+
+    if (ix >= 0x7ff00000) {
+        if ((ix & 0x7fffff) | lx != 0) {
+            const r = if (y == 0) y else x * y;
+            return Complex(f64).new(x, r);
+        }
+
+        const xx = @bitCast(f64, (u64(hx - 0x40000000) << 32) | lx);
+        const r = if (math.isInf(y)) y else math.sin(y) * math.cos(y);
+        return Complex(f64).new(xx, math.copysign(f64, 0, r));
+    }
+
+    if (!math.isFinite(y)) {
+        const r = if (ix != 0) y - y else x;
+        return Complex(f64).new(r, y - y);
+    }
+
+    // x >= 22
+    if (ix >= 0x40360000) {
+        const exp_mx = math.exp(-math.fabs(x));
+        return Complex(f64).new(math.copysign(f64, 1, x), 4 * math.sin(y) * math.cos(y) * exp_mx * exp_mx);
+    }
+
+    // Kahan's algorithm
+    const t = math.tan(y);
+    const beta = 1.0 + t * t;
+    const s = math.sinh(x);
+    const rho = math.sqrt(1 + s * s);
+    const den = 1 + beta * s * s;
+
+    return Complex(f64).new((beta * rho * s) / den, t / den);
+}
+
+const epsilon = 0.0001;
+
+test "complex.ctanh32" {
+    const a = Complex(f32).new(5, 3);
+    const c = tanh(a);
+
+    debug.assert(math.approxEq(f32, c.re,  0.999913, epsilon));
+    debug.assert(math.approxEq(f32, c.im, -0.000025, epsilon));
+}
+
+test "complex.ctanh64" {
+    const a = Complex(f64).new(5, 3);
+    const c = tanh(a);
+
+    debug.assert(math.approxEq(f64, c.re,  0.999913, epsilon));
+    debug.assert(math.approxEq(f64, c.im, -0.000025, epsilon));
+}
diff --git a/std/math/exp.zig b/std/math/exp.zig
index 4032930a43..21aa558c57 100644
--- a/std/math/exp.zig
+++ b/std/math/exp.zig
@@ -6,6 +6,7 @@
 const std = @import("../index.zig");
 const math = std.math;
 const assert = std.debug.assert;
+const builtin = @import("builtin");
 
 pub fn exp(x: var) @typeOf(x) {
     const T = @typeOf(x);
@@ -17,6 +18,8 @@ pub fn exp(x: var) @typeOf(x) {
 }
 
 fn exp32(x_: f32) f32 {
+    @setFloatMode(this, builtin.FloatMode.Strict);
+
     const half = []f32 { 0.5, -0.5 };
     const ln2hi = 6.9314575195e-1;
     const ln2lo = 1.4286067653e-6;
@@ -94,6 +97,8 @@ fn exp32(x_: f32) f32 {
 }
 
 fn exp64(x_: f64) f64 {
+    @setFloatMode(this, builtin.FloatMode.Strict);
+
     const half = []const f64 { 0.5, -0.5 };
     const ln2hi: f64 = 6.93147180369123816490e-01;
     const ln2lo: f64 = 1.90821492927058770002e-10;
diff --git a/std/math/index.zig b/std/math/index.zig
index 477dafcbcc..a549a6bb61 100644
--- a/std/math/index.zig
+++ b/std/math/index.zig
@@ -129,6 +129,9 @@ pub const cos = @import("cos.zig").cos;
 pub const sin = @import("sin.zig").sin;
 pub const tan = @import("tan.zig").tan;
 
+pub const complex = @import("complex/index.zig");
+pub const Complex = complex.Complex;
+
 test "math" {
     _ = @import("nan.zig");
     _ = @import("isnan.zig");
@@ -172,6 +175,8 @@ test "math" {
     _ = @import("sin.zig");
     _ = @import("cos.zig");
     _ = @import("tan.zig");
+
+    _ = @import("complex/index.zig");
 }
 
 
@@ -553,6 +558,32 @@ test "math.floorPowerOfTwo" {
     comptime testFloorPowerOfTwo();
 }
 
+pub fn log2_int(comptime T: type, x: T) Log2Int(T) {
+    assert(x != 0);
+    return Log2Int(T)(T.bit_count - 1 - @clz(x));
+}
+
+pub fn log2_int_ceil(comptime T: type, x: T) Log2Int(T) {
+    assert(x != 0);
+    const log2_val = log2_int(T, x);
+    if (T(1) << log2_val == x)
+        return log2_val;
+    return log2_val + 1;
+}
+
+test "std.math.log2_int_ceil" {
+    assert(log2_int_ceil(u32, 1) == 0);
+    assert(log2_int_ceil(u32, 2) == 1);
+    assert(log2_int_ceil(u32, 3) == 2);
+    assert(log2_int_ceil(u32, 4) == 2);
+    assert(log2_int_ceil(u32, 5) == 3);
+    assert(log2_int_ceil(u32, 6) == 3);
+    assert(log2_int_ceil(u32, 7) == 3);
+    assert(log2_int_ceil(u32, 8) == 3);
+    assert(log2_int_ceil(u32, 9) == 4);
+    assert(log2_int_ceil(u32, 10) == 4);
+}
+
 fn testFloorPowerOfTwo() void {
     assert(floorPowerOfTwo(u32, 63) == 32);
     assert(floorPowerOfTwo(u32, 64) == 64);
diff --git a/std/math/ln.zig b/std/math/ln.zig
index c349ed7c6f..d09494b998 100644
--- a/std/math/ln.zig
+++ b/std/math/ln.zig
@@ -89,6 +89,8 @@ pub fn ln_32(x_: f32) f32 {
 }
 
 pub fn ln_64(x_: f64) f64 {
+    @setFloatMode(this, @import("builtin").FloatMode.Strict);
+
     const ln2_hi: f64 = 6.93147180369123816490e-01;
     const ln2_lo: f64 = 1.90821492927058770002e-10;
     const Lg1: f64 = 6.666666666666735130e-01;
diff --git a/std/math/log2.zig b/std/math/log2.zig
index 998d6d6c5e..d5bbe385c2 100644
--- a/std/math/log2.zig
+++ b/std/math/log2.zig
@@ -31,17 +31,12 @@ pub fn log2(x: var) @typeOf(x) {
             return result;
         },
         TypeId.Int => {
-            return log2_int(T, x);
+            return math.log2_int(T, x);
         },
         else => @compileError("log2 not implemented for " ++ @typeName(T)),
     }
 }
 
-pub fn log2_int(comptime T: type, x: T) T {
-    assert(x != 0);
-    return T.bit_count - 1 - T(@clz(x));
-}
-
 pub fn log2_32(x_: f32) f32 {
     const ivln2hi: f32 =  1.4428710938e+00;
     const ivln2lo: f32 = -1.7605285393e-04;
diff --git a/std/math/sqrt.zig b/std/math/sqrt.zig
index 690f8b6901..982bd28b72 100644
--- a/std/math/sqrt.zig
+++ b/std/math/sqrt.zig
@@ -14,26 +14,8 @@ const TypeId = builtin.TypeId;
 pub fn sqrt(x: var) (if (@typeId(@typeOf(x)) == TypeId.Int) @IntType(false, @typeOf(x).bit_count / 2) else @typeOf(x)) {
     const T = @typeOf(x);
     switch (@typeId(T)) {
-        TypeId.FloatLiteral => {
-            return T(sqrt64(x));
-        },
-        TypeId.Float => {
-            switch (T) {
-                f32 => {
-                    switch (builtin.arch) {
-                        builtin.Arch.x86_64 => return @import("x86_64/sqrt.zig").sqrt32(x),
-                        else => return sqrt32(x),
-                    }
-                },
-                f64 => {
-                    switch (builtin.arch) {
-                        builtin.Arch.x86_64 => return @import("x86_64/sqrt.zig").sqrt64(x),
-                        else => return sqrt64(x),
-                    }
-                },
-                else => @compileError("sqrt not implemented for " ++ @typeName(T)),
-            }
-        },
+        TypeId.FloatLiteral => return T(@sqrt(f64, x)), // TODO upgrade to f128
+        TypeId.Float => return @sqrt(T, x),
         TypeId.IntLiteral => comptime {
             if (x > @maxValue(u128)) {
                 @compileError("sqrt not implemented for comptime_int greater than 128 bits");
@@ -43,269 +25,58 @@ pub fn sqrt(x: var) (if (@typeId(@typeOf(x)) == TypeId.Int) @IntType(false, @typ
             }
             return T(sqrt_int(u128, x));
         },
-        TypeId.Int => {
-            return sqrt_int(T, x);
-        },
+        TypeId.Int => return sqrt_int(T, x),
         else => @compileError("sqrt not implemented for " ++ @typeName(T)),
     }
 }
 
-fn sqrt32(x: f32) f32 {
-    const tiny: f32 = 1.0e-30;
-    const sign: i32 = @bitCast(i32, u32(0x80000000));
-    var ix: i32 = @bitCast(i32, x);
-
-    if ((ix & 0x7F800000) == 0x7F800000) {
-        return x * x + x;   // sqrt(nan) = nan, sqrt(+inf) = +inf, sqrt(-inf) = snan
-    }
-
-    // zero
-    if (ix <= 0) {
-        if (ix & ~sign == 0) {
-            return x;       // sqrt (+-0) = +-0
-        }
-        if (ix < 0) {
-            return math.snan(f32);
-        }
-    }
-
-    // normalize
-    var m = ix >> 23;
-    if (m == 0) {
-        // subnormal
-        var i: i32 = 0;
-        while (ix & 0x00800000 == 0) : (i += 1) {
-            ix <<= 1;
-        }
-        m -= i - 1;
-    }
-
-    m -= 127;               // unbias exponent
-    ix = (ix & 0x007FFFFF) | 0x00800000;
-
-    if (m & 1 != 0) {       // odd m, double x to even
-        ix += ix;
-    }
-
-    m >>= 1;                // m = [m / 2]
-
-    // sqrt(x) bit by bit
-    ix += ix;
-    var q: i32 = 0;              // q = sqrt(x)
-    var s: i32 = 0;
-    var r: i32 = 0x01000000;     // r = moving bit right -> left
-
-    while (r != 0) {
-        const t = s + r;
-        if (t <= ix) {
-            s = t + r;
-            ix -= t;
-            q += r;
-        }
-        ix += ix;
-        r >>= 1;
-    }
-
-    // floating add to find rounding direction
-    if (ix != 0) {
-        var z = 1.0 - tiny;     // inexact
-        if (z >= 1.0) {
-            z = 1.0 + tiny;
-            if (z > 1.0) {
-                q += 2;
-            } else {
-                if (q & 1 != 0) {
-                    q += 1;
-                }
-            }
-        }
-    }
-
-    ix = (q >> 1) + 0x3f000000;
-    ix += m << 23;
-    return @bitCast(f32, ix);
-}
-
-// NOTE: The original code is full of implicit signed -> unsigned assumptions and u32 wraparound
-// behaviour. Most intermediate i32 values are changed to u32 where appropriate but there are
-// potentially some edge cases remaining that are not handled in the same way.
-fn sqrt64(x: f64) f64 {
-    const tiny: f64 = 1.0e-300;
-    const sign: u32 = 0x80000000;
-    const u = @bitCast(u64, x);
-
-    var ix0 = u32(u >> 32);
-    var ix1 = u32(u & 0xFFFFFFFF);
-
-    // sqrt(nan) = nan, sqrt(+inf) = +inf, sqrt(-inf) = nan
-    if (ix0 & 0x7FF00000 == 0x7FF00000) {
-        return x * x + x;
-    }
-
-    // sqrt(+-0) = +-0
-    if (x == 0.0) {
-        return x;
-    }
-    // sqrt(-ve) = snan
-    if (ix0 & sign != 0) {
-        return math.snan(f64);
-    }
-
-    // normalize x
-    var m = i32(ix0 >> 20);
-    if (m == 0) {
-        // subnormal
-        while (ix0 == 0) {
-            m -= 21;
-            ix0 |= ix1 >> 11;
-            ix1 <<= 21;
-        }
-
-        // subnormal
-        var i: u32 = 0;
-        while (ix0 & 0x00100000 == 0) : (i += 1) {
-            ix0 <<= 1;
-        }
-        m -= i32(i) - 1;
-        ix0 |= ix1 >> u5(32 - i);
-        ix1 <<= u5(i);
-    }
-
-    // unbias exponent
-    m -= 1023;
-    ix0 = (ix0 & 0x000FFFFF) | 0x00100000;
-    if (m & 1 != 0) {
-        ix0 += ix0 + (ix1 >> 31);
-        ix1 = ix1 +% ix1;
-    }
-    m >>= 1;
-
-    // sqrt(x) bit by bit
-    ix0 += ix0 + (ix1 >> 31);
-    ix1 = ix1 +% ix1;
-
-    var q: u32 = 0;
-    var q1: u32 = 0;
-    var s0: u32 = 0;
-    var s1: u32 = 0;
-    var r: u32 = 0x00200000;
-    var t: u32 = undefined;
-    var t1: u32 = undefined;
-
-    while (r != 0) {
-        t = s0 +% r;
-        if (t <= ix0) {
-            s0 = t + r;
-            ix0 -= t;
-            q += r;
-        }
-        ix0 = ix0 +% ix0 +% (ix1 >> 31);
-        ix1 = ix1 +% ix1;
-        r >>= 1;
-    }
-
-    r = sign;
-    while (r != 0) {
-        t = s1 +% r;
-        t = s0;
-        if (t < ix0 or (t == ix0 and t1 <= ix1)) {
-            s1 = t1 +% r;
-            if (t1 & sign == sign and s1 & sign == 0) {
-                s0 += 1;
-            }
-            ix0 -= t;
-            if (ix1 < t1) {
-                ix0 -= 1;
-            }
-            ix1 = ix1 -% t1;
-            q1 += r;
-        }
-        ix0 = ix0 +% ix0 +% (ix1 >> 31);
-        ix1 = ix1 +% ix1;
-        r >>= 1;
-    }
-
-    // rounding direction
-    if (ix0 | ix1 != 0) {
-        var z = 1.0 - tiny;   // raise inexact
-        if (z >= 1.0) {
-            z = 1.0 + tiny;
-            if (q1 == 0xFFFFFFFF) {
-                q1 = 0;
-                q += 1;
-            } else if (z > 1.0) {
-                if (q1 == 0xFFFFFFFE) {
-                    q += 1;
-                }
-                q1 += 2;
-            } else {
-                q1 += q1 & 1;
-            }
-        }
-    }
-
-    ix0 = (q >> 1) + 0x3FE00000;
-    ix1 = q1 >> 1;
-    if (q & 1 != 0) {
-        ix1 |= 0x80000000;
-    }
-
-    // NOTE: musl here appears to rely on signed twos-complement wraparound. +% has the same
-    // behaviour at least.
-    var iix0 = i32(ix0);
-    iix0 = iix0 +% (m << 20);
-
-    const uz = (u64(iix0) << 32) | ix1;
-    return @bitCast(f64, uz);
-}
-
 test "math.sqrt" {
-    assert(sqrt(f32(0.0)) == sqrt32(0.0));
-    assert(sqrt(f64(0.0)) == sqrt64(0.0));
+    assert(sqrt(f32(0.0)) == @sqrt(f32, 0.0));
+    assert(sqrt(f64(0.0)) == @sqrt(f64, 0.0));
 }
 
 test "math.sqrt32" {
     const epsilon = 0.000001;
 
-    assert(sqrt32(0.0) == 0.0);
-    assert(math.approxEq(f32, sqrt32(2.0), 1.414214, epsilon));
-    assert(math.approxEq(f32, sqrt32(3.6), 1.897367, epsilon));
-    assert(sqrt32(4.0) == 2.0);
-    assert(math.approxEq(f32, sqrt32(7.539840), 2.745877, epsilon));
-    assert(math.approxEq(f32, sqrt32(19.230934), 4.385309, epsilon));
-    assert(sqrt32(64.0) == 8.0);
-    assert(math.approxEq(f32, sqrt32(64.1), 8.006248, epsilon));
-    assert(math.approxEq(f32, sqrt32(8942.230469), 94.563370, epsilon));
+    assert(@sqrt(f32, 0.0) == 0.0);
+    assert(math.approxEq(f32, @sqrt(f32, 2.0), 1.414214, epsilon));
+    assert(math.approxEq(f32, @sqrt(f32, 3.6), 1.897367, epsilon));
+    assert(@sqrt(f32, 4.0) == 2.0);
+    assert(math.approxEq(f32, @sqrt(f32, 7.539840), 2.745877, epsilon));
+    assert(math.approxEq(f32, @sqrt(f32, 19.230934), 4.385309, epsilon));
+    assert(@sqrt(f32, 64.0) == 8.0);
+    assert(math.approxEq(f32, @sqrt(f32, 64.1), 8.006248, epsilon));
+    assert(math.approxEq(f32, @sqrt(f32, 8942.230469), 94.563370, epsilon));
 }
 
 test "math.sqrt64" {
     const epsilon = 0.000001;
 
-    assert(sqrt64(0.0) == 0.0);
-    assert(math.approxEq(f64, sqrt64(2.0), 1.414214, epsilon));
-    assert(math.approxEq(f64, sqrt64(3.6), 1.897367, epsilon));
-    assert(sqrt64(4.0) == 2.0);
-    assert(math.approxEq(f64, sqrt64(7.539840), 2.745877, epsilon));
-    assert(math.approxEq(f64, sqrt64(19.230934), 4.385309, epsilon));
-    assert(sqrt64(64.0) == 8.0);
-    assert(math.approxEq(f64, sqrt64(64.1), 8.006248, epsilon));
-    assert(math.approxEq(f64, sqrt64(8942.230469), 94.563367, epsilon));
+    assert(@sqrt(f64, 0.0) == 0.0);
+    assert(math.approxEq(f64, @sqrt(f64, 2.0), 1.414214, epsilon));
+    assert(math.approxEq(f64, @sqrt(f64, 3.6), 1.897367, epsilon));
+    assert(@sqrt(f64, 4.0) == 2.0);
+    assert(math.approxEq(f64, @sqrt(f64, 7.539840), 2.745877, epsilon));
+    assert(math.approxEq(f64, @sqrt(f64, 19.230934), 4.385309, epsilon));
+    assert(@sqrt(f64, 64.0) == 8.0);
+    assert(math.approxEq(f64, @sqrt(f64, 64.1), 8.006248, epsilon));
+    assert(math.approxEq(f64, @sqrt(f64, 8942.230469), 94.563367, epsilon));
 }
 
 test "math.sqrt32.special" {
-    assert(math.isPositiveInf(sqrt32(math.inf(f32))));
-    assert(sqrt32(0.0) == 0.0);
-    assert(sqrt32(-0.0) == -0.0);
-    assert(math.isNan(sqrt32(-1.0)));
-    assert(math.isNan(sqrt32(math.nan(f32))));
+    assert(math.isPositiveInf(@sqrt(f32, math.inf(f32))));
+    assert(@sqrt(f32, 0.0) == 0.0);
+    assert(@sqrt(f32, -0.0) == -0.0);
+    assert(math.isNan(@sqrt(f32, -1.0)));
+    assert(math.isNan(@sqrt(f32, math.nan(f32))));
 }
 
 test "math.sqrt64.special" {
-    assert(math.isPositiveInf(sqrt64(math.inf(f64))));
-    assert(sqrt64(0.0) == 0.0);
-    assert(sqrt64(-0.0) == -0.0);
-    assert(math.isNan(sqrt64(-1.0)));
-    assert(math.isNan(sqrt64(math.nan(f64))));
+    assert(math.isPositiveInf(@sqrt(f64, math.inf(f64))));
+    assert(@sqrt(f64, 0.0) == 0.0);
+    assert(@sqrt(f64, -0.0) == -0.0);
+    assert(math.isNan(@sqrt(f64, -1.0)));
+    assert(math.isNan(@sqrt(f64, math.nan(f64))));
 }
 
 fn sqrt_int(comptime T: type, value: T) @IntType(false, T.bit_count / 2) {
diff --git a/std/math/x86_64/sqrt.zig b/std/math/x86_64/sqrt.zig
deleted file mode 100644
index ad9ce0c96c..0000000000
--- a/std/math/x86_64/sqrt.zig
+++ /dev/null
@@ -1,15 +0,0 @@
-pub fn sqrt32(x: f32) f32 {
-    return asm (
-        \\sqrtss %%xmm0, %%xmm0
-        : [ret] "={xmm0}" (-> f32)
-        : [x] "{xmm0}" (x)
-    );
-}
-
-pub fn sqrt64(x: f64) f64 {
-    return asm (
-        \\sqrtsd %%xmm0, %%xmm0
-        : [ret] "={xmm0}" (-> f64)
-        : [x] "{xmm0}" (x)
-    );
-}