Merge branch 'math'

This merges the standard library math functions that
Marc Tiehuis (@tiehuis) has been working on. Marc has
joined the Zig organization and now has commit access.
Thank you for this huge contribution to Zig.

Closes #374.

CMakeLists.txt

@@ -246,9 +246,51 @@ install(FILES "${CMAKE_SOURCE_DIR}/std/hash_map.zig" DESTINATION "${ZIG_STD_DEST
 install(FILES "${CMAKE_SOURCE_DIR}/std/index.zig" DESTINATION "${ZIG_STD_DEST}")
 install(FILES "${CMAKE_SOURCE_DIR}/std/io.zig" DESTINATION "${ZIG_STD_DEST}")
 install(FILES "${CMAKE_SOURCE_DIR}/std/linked_list.zig" DESTINATION "${ZIG_STD_DEST}")
-install(FILES "${CMAKE_SOURCE_DIR}/std/math/index.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/acos.zig" DESTINATION "${ZIG_STD_DEST}/math")
-install(FILES "${CMAKE_SOURCE_DIR}/std/math/frexp.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/acosh.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/asin.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/asinh.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/atan.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/atan2.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/atanh.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/cbrt.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/ceil.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/copysign.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/cos.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/cosh.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/exp.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/exp2.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/expm1.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/expo2.zig" DESTINATION "${ZIG_STD_DEST}/math")
 install(FILES "${CMAKE_SOURCE_DIR}/std/math/fabs.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/floor.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/fma.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/frexp.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/hypot.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/ilogb.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/index.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/inf.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/isfinite.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/isinf.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/isnan.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/isnormal.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/ln.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/log.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/log10.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/log1p.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/log2.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/modf.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/nan.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/pow.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/round.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/scalbn.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/signbit.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/sin.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/sinh.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/sqrt.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/tan.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/tanh.zig" DESTINATION "${ZIG_STD_DEST}/math")
+install(FILES "${CMAKE_SOURCE_DIR}/std/math/trunc.zig" DESTINATION "${ZIG_STD_DEST}/math")
 install(FILES "${CMAKE_SOURCE_DIR}/std/mem.zig" DESTINATION "${ZIG_STD_DEST}")
 install(FILES "${CMAKE_SOURCE_DIR}/std/net.zig" DESTINATION "${ZIG_STD_DEST}")
 install(FILES "${CMAKE_SOURCE_DIR}/std/os/child_process.zig" DESTINATION "${ZIG_STD_DEST}/os")

src/codegen.cpp

@@ -438,6 +438,7 @@ static LLVMValueRef fn_llvm_value(CodeGen *g, FnTableEntry *fn_table_entry) {
     }
 
     addLLVMFnAttr(fn_table_entry->llvm_value, "nounwind");
+    addLLVMFnAttr(fn_table_entry->llvm_value, "nobuiltin");
     if (g->build_mode == BuildModeDebug && fn_table_entry->fn_inline != FnInlineAlways) {
         ZigLLVMAddFunctionAttr(fn_table_entry->llvm_value, "no-frame-pointer-elim", "true");
         ZigLLVMAddFunctionAttr(fn_table_entry->llvm_value, "no-frame-pointer-elim-non-leaf", nullptr);

std/math/acos.zig

@@ -0,0 +1,182 @@
+// Special Cases:
+//
+// - acos(x)   = nan if x < -1 or x > 1
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+pub const acos = acos_workaround;
+
+// TODO issue #393
+pub fn acos_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(acos32, x),
+        f64 => @inlineCall(acos64, x),
+        else => @compileError("acos not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn r32(z: f32) -> f32 {
+    const pS0 =  1.6666586697e-01;
+    const pS1 = -4.2743422091e-02;
+    const pS2 = -8.6563630030e-03;
+    const qS1 = -7.0662963390e-01;
+
+    const p = z * (pS0 + z * (pS1 + z * pS2));
+    const q = 1.0 + z * qS1;
+    p / q
+}
+
+fn acos32(x: f32) -> f32 {
+    const pio2_hi = 1.5707962513e+00;
+    const pio2_lo = 7.5497894159e-08;
+
+    const hx: u32 = @bitCast(u32, x);
+    const ix: u32 = hx & 0x7FFFFFFF;
+
+    // |x| >= 1 or nan
+    if (ix >= 0x3F800000) {
+        if (ix == 0x3F800000) {
+            if (hx >> 31 != 0) {
+                return 2.0 * pio2_hi + 0x1.0p-120;
+            } else {
+                return 0;
+            }
+        } else {
+            return math.nan(f32);
+        }
+    }
+
+    // |x| < 0.5
+    if (ix < 0x3F000000) {
+        if (ix <= 0x32800000) { // |x| < 2^(-26)
+            return pio2_hi + 0x1.0p-120;
+        } else {
+            return pio2_hi - (x - (pio2_lo - x * r32(x * x)));
+        }
+    }
+
+    // x < -0.5
+    if (hx >> 31 != 0) {
+        const z = (1 + x) * 0.5;
+        const s = math.sqrt(z);
+        const w = r32(z) * s - pio2_lo;
+        return 2 * (pio2_hi - (s + w));
+    }
+
+    // x > 0.5
+    const z = (1.0 - x) * 0.5;
+    const s = math.sqrt(z);
+    const jx = @bitCast(u32, s);
+    const df = @bitCast(f32, jx & 0xFFFFF000);
+    const c = (z - df * df) / (s + df);
+    const w = r32(z) * s + c;
+    2 * (df + w)
+}
+
+fn r64(z: f64) -> f64 {
+    const pS0: f64 =  1.66666666666666657415e-01;
+    const pS1: f64 = -3.25565818622400915405e-01;
+    const pS2: f64 =  2.01212532134862925881e-01;
+    const pS3: f64 = -4.00555345006794114027e-02;
+    const pS4: f64 =  7.91534994289814532176e-04;
+    const pS5: f64 =  3.47933107596021167570e-05;
+    const qS1: f64 = -2.40339491173441421878e+00;
+    const qS2: f64 =  2.02094576023350569471e+00;
+    const qS3: f64 = -6.88283971605453293030e-01;
+    const qS4: f64 =  7.70381505559019352791e-02;
+
+    const p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5)))));
+    const q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4)));
+    p / q
+}
+
+fn acos64(x: f64) -> f64 {
+    const pio2_hi: f64 = 1.57079632679489655800e+00;
+    const pio2_lo: f64 = 6.12323399573676603587e-17;
+
+    const ux = @bitCast(u64, x);
+    const hx = u32(ux >> 32);
+    const ix = hx & 0x7FFFFFFF;
+
+    // |x| >= 1 or nan
+    if (ix >= 0x3FF00000) {
+        const lx = u32(ux & 0xFFFFFFFF);
+
+        // acos(1) = 0, acos(-1) = pi
+        if ((ix - 0x3FF00000) | lx == 0) {
+            if (hx >> 31 != 0) {
+                return 2 * pio2_hi + 0x1.0p-120;
+            } else {
+                return 0;
+            }
+        }
+
+        return math.nan(f32);
+    }
+
+    // |x| < 0.5
+    if (ix < 0x3FE00000) {
+        // |x| < 2^(-57)
+        if (ix <= 0x3C600000) {
+            return pio2_hi + 0x1.0p-120;
+        } else {
+            return pio2_hi - (x - (pio2_lo - x * r64(x * x)));
+        }
+    }
+
+    // x < -0.5
+    if (hx >> 31 != 0) {
+        const z = (1.0 + x) * 0.5;
+        const s = math.sqrt(z);
+        const w = r64(z) * s - pio2_lo;
+        return 2 * (pio2_hi - (s + w));
+    }
+
+    // x > 0.5
+    const z = (1.0 - x) * 0.5;
+    const s = math.sqrt(z);
+    const jx = @bitCast(u64, s);
+    const df = @bitCast(f64, jx & 0xFFFFFFFF00000000);
+    const c = (z - df * df) / (s + df);
+    const w = r64(z) * s + c;
+    2 * (df + w)
+}
+
+test "math.acos" {
+    assert(acos_workaround(f32(0.0)) == acos32(0.0));
+    assert(acos_workaround(f64(0.0)) == acos64(0.0));
+}
+
+test "math.acos32" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, acos32(0.0), 1.570796, epsilon));
+    assert(math.approxEq(f32, acos32(0.2), 1.369438, epsilon));
+    assert(math.approxEq(f32, acos32(0.3434), 1.220262, epsilon));
+    assert(math.approxEq(f32, acos32(0.5), 1.047198, epsilon));
+    assert(math.approxEq(f32, acos32(0.8923), 0.468382, epsilon));
+    assert(math.approxEq(f32, acos32(-0.2), 1.772154, epsilon));
+}
+
+test "math.acos64" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f64, acos64(0.0), 1.570796, epsilon));
+    assert(math.approxEq(f64, acos64(0.2), 1.369438, epsilon));
+    assert(math.approxEq(f64, acos64(0.3434), 1.220262, epsilon));
+    assert(math.approxEq(f64, acos64(0.5), 1.047198, epsilon));
+    assert(math.approxEq(f64, acos64(0.8923), 0.468382, epsilon));
+    assert(math.approxEq(f64, acos64(-0.2), 1.772154, epsilon));
+}
+
+test "math.acos32.special" {
+    assert(math.isNan(acos32(-2)));
+    assert(math.isNan(acos32(1.5)));
+}
+
+test "math.acos64.special" {
+    assert(math.isNan(acos64(-2)));
+    assert(math.isNan(acos64(1.5)));
+}

std/math/acosh.zig

@@ -0,0 +1,89 @@
+// Special Cases:
+//
+// - acosh(x)   = snan if x < 1
+// - acosh(nan) = nan
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+pub const acosh = acosh_workaround;
+
+// TODO issue #393
+pub fn acosh_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(acosh32, x),
+        f64 => @inlineCall(acosh64, x),
+        else => @compileError("acosh not implemented for " ++ @typeName(T)),
+    }
+}
+
+// acosh(x) = log(x + sqrt(x * x - 1))
+fn acosh32(x: f32) -> f32 {
+    const u = @bitCast(u32, x);
+    const i = u & 0x7FFFFFFF;
+
+    // |x| < 2, invalid if x < 1 or nan
+    if (i < 0x3F800000 + (1 << 23)) {
+        math.log1p(x - 1 + math.sqrt((x - 1) * (x - 1) + 2 * (x - 1)))
+    }
+    // |x| < 0x1p12
+    else if (i < 0x3F800000 + (12 << 23)) {
+        math.ln(2 * x - 1 / (x + math.sqrt(x * x - 1)))
+    }
+    // |x| >= 0x1p12
+    else {
+        math.ln(x) + 0.693147180559945309417232121458176568
+    }
+}
+
+fn acosh64(x: f64) -> f64 {
+    const u = @bitCast(u64, x);
+    const e = (u >> 52) & 0x7FF;
+
+    // |x| < 2, invalid if x < 1 or nan
+    if (e < 0x3FF + 1) {
+        math.log1p(x - 1 + math.sqrt((x - 1) * (x - 1) + 2 * (x - 1)))
+    }
+    // |x| < 0x1p26
+    else if (e < 0x3FF + 26) {
+        math.ln(2 * x - 1 / (x + math.sqrt(x * x - 1)))
+    }
+    // |x| >= 0x1p26 or nan
+    else {
+        math.ln(x) + 0.693147180559945309417232121458176568
+    }
+}
+
+test "math.acosh" {
+    assert(acosh_workaround(f32(1.5)) == acosh32(1.5));
+    assert(acosh_workaround(f64(1.5)) == acosh64(1.5));
+}
+
+test "math.acosh32" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, acosh32(1.5), 0.962424, epsilon));
+    assert(math.approxEq(f32, acosh32(37.45), 4.315976, epsilon));
+    assert(math.approxEq(f32, acosh32(89.123), 5.183133, epsilon));
+    assert(math.approxEq(f32, acosh32(123123.234375), 12.414088, epsilon));
+}
+
+test "math.acosh64" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f64, acosh64(1.5), 0.962424, epsilon));
+    assert(math.approxEq(f64, acosh64(37.45), 4.315976, epsilon));
+    assert(math.approxEq(f64, acosh64(89.123), 5.183133, epsilon));
+    assert(math.approxEq(f64, acosh64(123123.234375), 12.414088, epsilon));
+}
+
+test "math.acosh32.special" {
+    assert(math.isNan(acosh32(math.nan(f32))));
+    assert(math.isSignalNan(acosh32(0.5)));
+}
+
+test "math.acosh64.special" {
+    assert(math.isNan(acosh64(math.nan(f64))));
+    assert(math.isSignalNan(acosh64(0.5)));
+}

std/math/asin.zig

@@ -0,0 +1,179 @@
+// Special Cases:
+//
+// - asin(+-0) = +-0
+// - asin(x)   = nan if x < -1 or x > 1
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+pub const asin = asin_workaround;
+
+// TODO issue #393
+pub fn asin_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(asin32, x),
+        f64 => @inlineCall(asin64, x),
+        else => @compileError("asin not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn r32(z: f32) -> f32 {
+    const pS0 =  1.6666586697e-01;
+    const pS1 = -4.2743422091e-02;
+    const pS2 = -8.6563630030e-03;
+    const qS1 = -7.0662963390e-01;
+
+    const p = z * (pS0 + z * (pS1 + z * pS2));
+    const q = 1.0 + z * qS1;
+    p / q
+}
+
+fn asin32(x: f32) -> f32 {
+    const pio2 = 1.570796326794896558e+00;
+
+    const hx: u32 = @bitCast(u32, x);
+    const ix: u32 = hx & 0x7FFFFFFF;
+
+    // |x| >= 1
+    if (ix >= 0x3F800000) {
+        // |x| >= 1
+        if (ix == 0x3F800000) {
+            return x * pio2 + 0x1.0p-120;   // asin(+-1) = +-pi/2 with inexact
+        } else {
+            return math.nan(f32);           // asin(|x| > 1) is nan
+        }
+    }
+
+    // |x| < 0.5
+    if (ix < 0x3F000000) {
+        // 0x1p-126 <= |x| < 0x1p-12
+        if (ix < 0x39800000 and ix >= 0x00800000) {
+            return x;
+        } else {
+            return x + x * r32(x * x);
+        }
+    }
+
+    // 1 > |x| >= 0.5
+    const z = (1 - math.fabs(x)) * 0.5;
+    const s = math.sqrt(z);
+    const fx = pio2 - 2 * (s + s * r32(z));
+
+    if (hx >> 31 != 0) {
+        -fx
+    } else {
+        fx
+    }
+}
+
+fn r64(z: f64) -> f64 {
+    const pS0: f64 =  1.66666666666666657415e-01;
+    const pS1: f64 = -3.25565818622400915405e-01;
+    const pS2: f64 =  2.01212532134862925881e-01;
+    const pS3: f64 = -4.00555345006794114027e-02;
+    const pS4: f64 =  7.91534994289814532176e-04;
+    const pS5: f64 =  3.47933107596021167570e-05;
+    const qS1: f64 = -2.40339491173441421878e+00;
+    const qS2: f64 =  2.02094576023350569471e+00;
+    const qS3: f64 = -6.88283971605453293030e-01;
+    const qS4: f64 =  7.70381505559019352791e-02;
+
+    const p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5)))));
+    const q = 1.0 + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4)));
+    p / q
+}
+
+fn asin64(x: f64) -> f64 {
+    const pio2_hi: f64 = 1.57079632679489655800e+00;
+    const pio2_lo: f64 = 6.12323399573676603587e-17;
+
+    const ux = @bitCast(u64, x);
+    const hx = u32(ux >> 32);
+    const ix = hx & 0x7FFFFFFF;
+
+    // |x| >= 1 or nan
+    if (ix >= 0x3FF00000) {
+        const lx = u32(ux & 0xFFFFFFFF);
+
+        // asin(1) = +-pi/2 with inexact
+        if ((ix - 0x3FF00000) | lx == 0) {
+            return x * pio2_hi + 0x1.0p-120;
+        } else {
+            return math.nan(f64);
+        }
+    }
+
+    // |x| < 0.5
+    if (ix < 0x3FE00000) {
+        // if 0x1p-1022 <= |x| < 0x1p-26 avoid raising overflow
+        if (ix < 0x3E500000 and ix >= 0x00100000) {
+            return x;
+        } else {
+            return x + x * r64(x * x);
+        }
+    }
+
+    // 1 > |x| >= 0.5
+    const z = (1 - math.fabs(x)) * 0.5;
+    const s = math.sqrt(z);
+    const r = r64(z);
+    var fx: f64 = undefined;
+
+    // |x| > 0.975
+    if (ix >= 0x3FEF3333) {
+        fx = pio2_hi - 2 * (s + s * r)
+    } else {
+        const jx = @bitCast(u64, s);
+        const df = @bitCast(f64, jx & 0xFFFFFFFF00000000);
+        const c = (z - df * df) / (s + df);
+        fx = 0.5 * pio2_hi - (2 * s * r - (pio2_lo - 2 * c) - (0.5 * pio2_hi - 2 * df));
+    }
+
+    if (hx >> 31 != 0) {
+        -fx
+    } else {
+        fx
+    }
+}
+
+test "math.asin" {
+    assert(asin_workaround(f32(0.0)) == asin32(0.0));
+    assert(asin_workaround(f64(0.0)) == asin64(0.0));
+}
+
+test "math.asin32" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, asin32(0.0), 0.0, epsilon));
+    assert(math.approxEq(f32, asin32(0.2), 0.201358, epsilon));
+    assert(math.approxEq(f32, asin32(-0.2), -0.201358, epsilon));
+    assert(math.approxEq(f32, asin32(0.3434), 0.350535, epsilon));
+    assert(math.approxEq(f32, asin32(0.5), 0.523599, epsilon));
+    assert(math.approxEq(f32, asin32(0.8923), 1.102415, epsilon));
+}
+
+test "math.asin64" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f64, asin64(0.0), 0.0, epsilon));
+    assert(math.approxEq(f64, asin64(0.2), 0.201358, epsilon));
+    assert(math.approxEq(f64, asin64(-0.2), -0.201358, epsilon));
+    assert(math.approxEq(f64, asin64(0.3434), 0.350535, epsilon));
+    assert(math.approxEq(f64, asin64(0.5), 0.523599, epsilon));
+    assert(math.approxEq(f64, asin64(0.8923), 1.102415, epsilon));
+}
+
+test "math.asin32.special" {
+    assert(asin32(0.0) == 0.0);
+    assert(asin32(-0.0) == -0.0);
+    assert(math.isNan(asin32(-2)));
+    assert(math.isNan(asin32(1.5)));
+}
+
+test "math.asin64.special" {
+    assert(asin64(0.0) == 0.0);
+    assert(asin64(-0.0) == -0.0);
+    assert(math.isNan(asin64(-2)));
+    assert(math.isNan(asin64(1.5)));
+}

std/math/asinh.zig

@@ -0,0 +1,129 @@
+// Special Cases:
+//
+// - asinh(+-0)   = +-0
+// - asinh(+-inf) = +-inf
+// - asinh(nan)   = nan
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+pub const asinh = asinh_workaround;
+
+// TODO issue #393
+pub fn asinh_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(asinh32, x),
+        f64 => @inlineCall(asinh64, x),
+        else => @compileError("asinh not implemented for " ++ @typeName(T)),
+    }
+}
+
+// asinh(x) = sign(x) * log(|x| + sqrt(x * x + 1)) ~= x - x^3/6 + o(x^5)
+fn asinh32(x: f32) -> f32 {
+    const u = @bitCast(u32, x);
+    const i = u & 0x7FFFFFFF;
+    const s = i >> 31;
+
+    var rx = @bitCast(f32, i); // |x|
+
+    // TODO: Shouldn't need this explicit check.
+    if (math.isNegativeInf(x)) {
+        return x;
+    }
+
+    // |x| >= 0x1p12 or inf or nan
+    if (i >= 0x3F800000 + (12 << 23)) {
+        rx = math.ln(rx) + 0.69314718055994530941723212145817656;
+    }
+    // |x| >= 2
+    else if (i >= 0x3F800000 + (1 << 23)) {
+        rx = math.ln(2 * x + 1 / (math.sqrt(x * x + 1) + x));
+    }
+    // |x| >= 0x1p-12, up to 1.6ulp error
+    else if (i >= 0x3F800000 - (12 << 23)) {
+        rx = math.log1p(x + x * x / (math.sqrt(x * x + 1) + 1));
+    }
+    // |x| < 0x1p-12, inexact if x != 0
+    else {
+        math.forceEval(x + 0x1.0p120);
+    }
+
+    if (s != 0) -rx else rx
+}
+
+fn asinh64(x: f64) -> f64 {
+    const u = @bitCast(u64, x);
+    const e = (u >> 52) & 0x7FF;
+    const s = u >> 63;
+
+    var rx = @bitCast(f64, u & (@maxValue(u64) >> 1)); // |x|
+
+    if (math.isNegativeInf(x)) {
+        return x;
+    }
+
+    // |x| >= 0x1p26 or inf or nan
+    if (e >= 0x3FF + 26) {
+        rx = math.ln(rx) + 0.693147180559945309417232121458176568;
+    }
+    // |x| >= 2
+    else if (e >= 0x3FF + 1) {
+        rx = math.ln(2 * x + 1 / (math.sqrt(x * x + 1) + x));
+    }
+    // |x| >= 0x1p-12, up to 1.6ulp error
+    else if (e >= 0x3FF - 26) {
+        rx = math.log1p(x + x * x / (math.sqrt(x * x + 1) + 1));
+    }
+    // |x| < 0x1p-12, inexact if x != 0
+    else {
+        math.forceEval(x + 0x1.0p120);
+    }
+
+    if (s != 0) -rx else rx
+}
+
+test "math.asinh" {
+    assert(asinh_workaround(f32(0.0)) == asinh32(0.0));
+    assert(asinh_workaround(f64(0.0)) == asinh64(0.0));
+}
+
+test "math.asinh32" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, asinh32(0.0), 0.0, epsilon));
+    assert(math.approxEq(f32, asinh32(0.2), 0.198690, epsilon));
+    assert(math.approxEq(f32, asinh32(0.8923), 0.803133, epsilon));
+    assert(math.approxEq(f32, asinh32(1.5), 1.194763, epsilon));
+    assert(math.approxEq(f32, asinh32(37.45), 4.316332, epsilon));
+    assert(math.approxEq(f32, asinh32(89.123), 5.183196, epsilon));
+    assert(math.approxEq(f32, asinh32(123123.234375), 12.414088, epsilon));
+}
+
+test "math.asinh64" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f64, asinh64(0.0), 0.0, epsilon));
+    assert(math.approxEq(f64, asinh64(0.2), 0.198690, epsilon));
+    assert(math.approxEq(f64, asinh64(0.8923), 0.803133, epsilon));
+    assert(math.approxEq(f64, asinh64(1.5), 1.194763, epsilon));
+    assert(math.approxEq(f64, asinh64(37.45), 4.316332, epsilon));
+    assert(math.approxEq(f64, asinh64(89.123), 5.183196, epsilon));
+    assert(math.approxEq(f64, asinh64(123123.234375), 12.414088, epsilon));
+}
+
+test "math.asinh32.special" {
+    assert(asinh32(0.0) == 0.0);
+    assert(asinh32(-0.0) == -0.0);
+    assert(math.isPositiveInf(asinh32(math.inf(f32))));
+    assert(math.isNegativeInf(asinh32(-math.inf(f32))));
+    assert(math.isNan(asinh32(math.nan(f32))));
+}
+
+test "math.asinh64.special" {
+    assert(asinh64(0.0) == 0.0);
+    assert(asinh64(-0.0) == -0.0);
+    assert(math.isPositiveInf(asinh64(math.inf(f64))));
+    assert(math.isNegativeInf(asinh64(-math.inf(f64))));
+    assert(math.isNan(asinh64(math.nan(f64))));
+}

std/math/atan.zig

@@ -0,0 +1,253 @@
+// Special Cases:
+//
+// - atan(+-0)   = +-0
+// - atan(+-inf) = +-pi/2
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const atan = atan_workaround;
+
+pub fn atan_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(atan32, x),
+        f64 => @inlineCall(atan64, x),
+        else => @compileError("atan not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn atan32(x_: f32) -> f32 {
+    const atanhi = []const f32 {
+        4.6364760399e-01, // atan(0.5)hi
+        7.8539812565e-01, // atan(1.0)hi
+        9.8279368877e-01, // atan(1.5)hi
+        1.5707962513e+00, // atan(inf)hi
+    };
+
+    const atanlo = []const f32 {
+        5.0121582440e-09, // atan(0.5)lo
+        3.7748947079e-08, // atan(1.0)lo
+        3.4473217170e-08, // atan(1.5)lo
+        7.5497894159e-08, // atan(inf)lo
+    };
+
+    const aT = []const f32 {
+        3.3333328366e-01,
+       -1.9999158382e-01,
+        1.4253635705e-01,
+       -1.0648017377e-01,
+        6.1687607318e-02,
+    };
+
+    var x = x_;
+    var ix: u32 = @bitCast(u32, x);
+    const sign = ix >> 31;
+    ix &= 0x7FFFFFFF;
+
+    // |x| >= 2^26
+    if (ix >= 0x4C800000) {
+        if (math.isNan(x)) {
+            return x;
+        } else {
+            const z = atanhi[3] + 0x1.0p-120;
+            return if (sign != 0) -z else z;
+        }
+    }
+
+    var id: ?usize = undefined;
+
+    // |x| < 0.4375
+    if (ix < 0x3EE00000) {
+        // |x| < 2^(-12)
+        if (ix < 0x39800000) {
+            if (ix < 0x00800000) {
+                math.forceEval(x * x);
+            }
+            return x;
+        }
+        id = null;
+    } else {
+        x = math.fabs(x);
+        // |x| < 1.1875
+        if (ix < 0x3F980000) {
+            // 7/16 <= |x| < 11/16
+            if (ix < 0x3F300000) {
+                id = 0;
+                x = (2.0 * x - 1.0) / (2.0 + x);
+            }
+            // 11/16 <= |x| < 19/16
+            else {
+                id = 1;
+                x = (x - 1.0) / (x + 1.0);
+            }
+        }
+        else {
+            // |x| < 2.4375
+            if (ix < 0x401C0000) {
+                id = 2;
+                x = (x - 1.5) / (1.0 + 1.5 * x);
+            }
+            // 2.4375 <= |x| < 2^26
+            else {
+                id = 3;
+                x = -1.0 / x;
+            }
+        }
+    }
+
+    const z = x * x;
+    const w = z * z;
+    const s1 = z * (aT[0] + w * (aT[2] + w * aT[4]));
+    const s2 = w * (aT[1] + w * aT[3]);
+
+    if (id == null) {
+        x - x * (s1 + s2)
+    } else {
+        const zz = atanhi[??id] - ((x * (s1 + s2) - atanlo[??id]) - x);
+        if (sign != 0) -zz else zz
+    }
+}
+
+fn atan64(x_: f64) -> f64 {
+    const atanhi = []const f64 {
+        4.63647609000806093515e-01, // atan(0.5)hi
+        7.85398163397448278999e-01, // atan(1.0)hi
+        9.82793723247329054082e-01, // atan(1.5)hi
+        1.57079632679489655800e+00, // atan(inf)hi
+    };
+
+    const atanlo = []const f64 {
+        2.26987774529616870924e-17, // atan(0.5)lo
+        3.06161699786838301793e-17, // atan(1.0)lo
+        1.39033110312309984516e-17, // atan(1.5)lo
+        6.12323399573676603587e-17, // atan(inf)lo
+    };
+
+    const aT = []const f64 {
+        3.33333333333329318027e-01,
+       -1.99999999998764832476e-01,
+        1.42857142725034663711e-01,
+       -1.11111104054623557880e-01,
+        9.09088713343650656196e-02,
+       -7.69187620504482999495e-02,
+        6.66107313738753120669e-02,
+       -5.83357013379057348645e-02,
+        4.97687799461593236017e-02,
+       -3.65315727442169155270e-02,
+        1.62858201153657823623e-02,
+    };
+
+    var x = x_;
+    var ux = @bitCast(u64, x);
+    var ix = u32(ux >> 32);
+    const sign = ix >> 31;
+    ix &= 0x7FFFFFFF;
+
+    // |x| >= 2^66
+    if (ix >= 0x44100000) {
+        if (math.isNan(x)) {
+            return x;
+        } else {
+            const z = atanhi[3] + 0x1.0p-120;
+            return if (sign != 0) -z else z;
+        }
+    }
+
+    var id: ?usize = undefined;
+
+    // |x| < 0.4375
+    if (ix < 0x3DFC0000) {
+        // |x| < 2^(-27)
+        if (ix < 0x3E400000) {
+            if (ix < 0x00100000) {
+                math.forceEval(f32(x));
+            }
+            return x;
+        }
+        id = null;
+    } else {
+        x = math.fabs(x);
+        // |x| < 1.1875
+        if (ix < 0x3FF30000) {
+            // 7/16 <= |x| < 11/16
+            if (ix < 0x3FE60000) {
+                id = 0;
+                x = (2.0 * x - 1.0) / (2.0 + x);
+            }
+            // 11/16 <= |x| < 19/16
+            else {
+                id = 1;
+                x = (x - 1.0) / (x + 1.0);
+            }
+        }
+        else {
+            // |x| < 2.4375
+            if (ix < 0x40038000) {
+                id = 2;
+                x = (x - 1.5) / (1.0 + 1.5 * x);
+            }
+            // 2.4375 <= |x| < 2^66
+            else {
+                id = 3;
+                x = -1.0 / x;
+            }
+        }
+    }
+
+    const z = x * x;
+    const w = z * z;
+    const s1 = z * (aT[0] + w * (aT[2] + w * (aT[4] + w * (aT[6] + w * (aT[8] + w * aT[10])))));
+    const s2 = w * (aT[1] + w * (aT[3] + w * (aT[5] + w * (aT[7] + w * aT[9]))));
+
+    if (id == null) {
+        x - x * (s1 + s2)
+    } else {
+        const zz = atanhi[??id] - ((x * (s1 + s2) - atanlo[??id]) - x);
+        if (sign != 0) -zz else zz
+    }
+}
+
+test "math.atan" {
+    assert(atan_workaround(f32(0.2)) == atan32(0.2));
+    assert(atan_workaround(f64(0.2)) == atan64(0.2));
+}
+
+test "math.atan32" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, atan32(0.2), 0.197396, epsilon));
+    assert(math.approxEq(f32, atan32(-0.2), -0.197396, epsilon));
+    assert(math.approxEq(f32, atan32(0.3434), 0.330783, epsilon));
+    assert(math.approxEq(f32, atan32(0.8923), 0.728545, epsilon));
+    assert(math.approxEq(f32, atan32(1.5), 0.982794, epsilon));
+}
+
+test "math.atan64" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f64, atan64(0.2), 0.197396, epsilon));
+    assert(math.approxEq(f64, atan64(-0.2), -0.197396, epsilon));
+    assert(math.approxEq(f64, atan64(0.3434), 0.330783, epsilon));
+    assert(math.approxEq(f64, atan64(0.8923), 0.728545, epsilon));
+    assert(math.approxEq(f64, atan64(1.5), 0.982794, epsilon));
+}
+
+test "math.atan32.special" {
+    const epsilon = 0.000001;
+
+    assert(atan32(0.0) == 0.0);
+    assert(atan32(-0.0) == -0.0);
+    assert(math.approxEq(f32, atan32(math.inf(f32)), math.pi_2, epsilon));
+    assert(math.approxEq(f32, atan32(-math.inf(f32)), -math.pi_2, epsilon));
+}
+
+test "math.atan64.special" {
+    const epsilon = 0.000001;
+
+    assert(atan64(0.0) == 0.0);
+    assert(atan64(-0.0) == -0.0);
+    assert(math.approxEq(f64, atan64(math.inf(f64)), math.pi_2, epsilon));
+    assert(math.approxEq(f64, atan64(-math.inf(f64)), -math.pi_2, epsilon));
+}

std/math/atan2.zig

@@ -0,0 +1,285 @@
+// Special Cases:
+//
+//  atan2(y, nan)     = nan
+//  atan2(nan, x)     = nan
+//  atan2(+0, x>=0)   = +0
+//  atan2(-0, x>=0)   = -0
+//  atan2(+0, x<=-0)  = +pi
+//  atan2(-0, x<=-0)  = -pi
+//  atan2(y>0, 0)     = +pi/2
+//  atan2(y<0, 0)     = -pi/2
+//  atan2(+inf, +inf) = +pi/4
+//  atan2(-inf, +inf) = -pi/4
+//  atan2(+inf, -inf) = 3pi/4
+//  atan2(-inf, -inf) = -3pi/4
+//  atan2(y, +inf)    = 0
+//  atan2(y>0, -inf)  = +pi
+//  atan2(y<0, -inf)  = -pi
+//  atan2(+inf, x)    = +pi/2
+//  atan2(-inf, x)    = -pi/2
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+pub const atan2 = atan2_workaround;
+
+// TODO issue #393
+pub fn atan2_workaround(comptime T: type, x: T, y: T) -> T {
+    switch (T) {
+        f32 => @inlineCall(atan2_32, x, y),
+        f64 => @inlineCall(atan2_64, x, y),
+        else => @compileError("atan2 not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn atan2_32(y: f32, x: f32) -> f32 {
+    const pi: f32    =  3.1415927410e+00;
+    const pi_lo: f32 = -8.7422776573e-08;
+
+    if (math.isNan(x) or math.isNan(y)) {
+        return x + y;
+    }
+
+    var ix = @bitCast(u32, x);
+    var iy = @bitCast(u32, y);
+
+    // x = 1.0
+    if (ix == 0x3F800000) {
+        return math.atan(y);
+    }
+
+    // 2 * sign(x) + sign(y)
+    const m = ((iy >> 31) & 1) | ((ix >> 30) & 2);
+    ix &= 0x7FFFFFFF;
+    iy &= 0x7FFFFFFF;
+
+    if (iy == 0) {
+        switch (m) {
+            0, 1 => return  y,          // atan(+-0, +...)
+            2    => return  pi,         // atan(+0, -...)
+            3    => return -pi,         // atan(-0, -...)
+            else => unreachable,
+        }
+    }
+
+    if (ix == 0) {
+        if (m & 1 != 0) {
+            return -pi / 2;
+        } else {
+            return pi / 2;
+        }
+    }
+
+    if (ix == 0x7F800000) {
+        if (iy == 0x7F800000) {
+            switch (m) {
+                0    => return    pi / 4,   // atan(+inf, +inf)
+                1    => return   -pi / 4,   // atan(-inf, +inf)
+                2    => return  3*pi / 4,   // atan(+inf, -inf)
+                3    => return -3*pi / 4,   // atan(-inf, -inf)
+                else => unreachable,
+            }
+        } else {
+            switch (m) {
+                0    => return  0.0,   // atan(+..., +inf)
+                1    => return -0.0,   // atan(-..., +inf)
+                2    => return   pi,   // atan(+..., -inf)
+                3    => return  -pi,   // atan(-...f, -inf)
+                else => unreachable,
+            }
+        }
+    }
+
+    // |y / x| > 0x1p26
+    if (ix + (26 << 23) < iy or iy == 0x7F800000) {
+        if (m & 1 != 0) {
+            return -pi / 2;
+        } else {
+            return pi / 2;
+        }
+    }
+
+    // z = atan(|y / x|) with correct underflow
+    var z = {
+        if ((m & 2) != 0 and iy + (26 << 23) < ix) {
+            0.0
+        } else {
+            math.atan(math.fabs(y / x))
+        }
+    };
+
+    switch (m) {
+        0    => return  z,                  // atan(+, +)
+        1    => return -z,                  // atan(-, +)
+        2    => return pi - (z - pi_lo),    // atan(+, -)
+        3    => return (z - pi_lo) - pi,    // atan(-, -)
+        else => unreachable,
+    }
+}
+
+fn atan2_64(y: f64, x: f64) -> f64 {
+    const pi: f64    = 3.1415926535897931160E+00;
+    const pi_lo: f64 = 1.2246467991473531772E-16;
+
+    if (math.isNan(x) or math.isNan(y)) {
+        return x + y;
+    }
+
+    var ux = @bitCast(u64, x);
+    var ix = u32(ux >> 32);
+    var lx = u32(ux & 0xFFFFFFFF);
+
+    var uy = @bitCast(u64, y);
+    var iy = u32(uy >> 32);
+    var ly = u32(uy & 0xFFFFFFFF);
+
+    // x = 1.0
+    if ((ix -% 0x3FF00000) | lx == 0) {
+        return math.atan(y);
+    }
+
+    // 2 * sign(x) + sign(y)
+    const m = ((iy >> 31) & 1) | ((ix >> 30) & 2);
+    ix &= 0x7FFFFFFF;
+    iy &= 0x7FFFFFFF;
+
+    if (iy | ly == 0) {
+        switch (m) {
+            0, 1 => return  y,          // atan(+-0, +...)
+            2    => return  pi,         // atan(+0, -...)
+            3    => return -pi,         // atan(-0, -...)
+            else => unreachable,
+        }
+    }
+
+    if (ix | lx == 0) {
+        if (m & 1 != 0) {
+            return -pi / 2;
+        } else {
+            return pi / 2;
+        }
+    }
+
+    if (ix == 0x7FF00000) {
+        if (iy == 0x7FF00000) {
+            switch (m) {
+                0    => return    pi / 4,   // atan(+inf, +inf)
+                1    => return   -pi / 4,   // atan(-inf, +inf)
+                2    => return  3*pi / 4,   // atan(+inf, -inf)
+                3    => return -3*pi / 4,   // atan(-inf, -inf)
+                else => unreachable,
+            }
+        } else {
+            switch (m) {
+                0    => return  0.0,   // atan(+..., +inf)
+                1    => return -0.0,   // atan(-..., +inf)
+                2    => return   pi,   // atan(+..., -inf)
+                3    => return  -pi,   // atan(-...f, -inf)
+                else => unreachable,
+            }
+        }
+    }
+
+    // |y / x| > 0x1p64
+    if (ix +% (64 << 20) < iy or iy == 0x7FF00000) {
+        if (m & 1 != 0) {
+            return -pi / 2;
+        } else {
+            return pi / 2;
+        }
+    }
+
+    // z = atan(|y / x|) with correct underflow
+    var z = {
+        if ((m & 2) != 0 and iy +% (64 << 20) < ix) {
+            0.0
+        } else {
+            math.atan(math.fabs(y / x))
+        }
+    };
+
+    switch (m) {
+        0    => return  z,                  // atan(+, +)
+        1    => return -z,                  // atan(-, +)
+        2    => return pi - (z - pi_lo),    // atan(+, -)
+        3    => return (z - pi_lo) - pi,    // atan(-, -)
+        else => unreachable,
+    }
+}
+
+test "math.atan2" {
+    assert(atan2_workaround(f32, 0.2, 0.21) == atan2_32(0.2, 0.21));
+    assert(atan2_workaround(f64, 0.2, 0.21) == atan2_64(0.2, 0.21));
+}
+
+test "math.atan2_32" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, atan2_32(0.0, 0.0), 0.0, epsilon));
+    assert(math.approxEq(f32, atan2_32(0.2, 0.2), 0.785398, epsilon));
+    assert(math.approxEq(f32, atan2_32(-0.2, 0.2), -0.785398, epsilon));
+    assert(math.approxEq(f32, atan2_32(0.2, -0.2), 2.356194, epsilon));
+    assert(math.approxEq(f32, atan2_32(-0.2, -0.2), -2.356194, epsilon));
+    assert(math.approxEq(f32, atan2_32(0.34, -0.4), 2.437099, epsilon));
+    assert(math.approxEq(f32, atan2_32(0.34, 1.243), 0.267001, epsilon));
+}
+
+test "math.atan2_64" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f64, atan2_64(0.0, 0.0), 0.0, epsilon));
+    assert(math.approxEq(f64, atan2_64(0.2, 0.2), 0.785398, epsilon));
+    assert(math.approxEq(f64, atan2_64(-0.2, 0.2), -0.785398, epsilon));
+    assert(math.approxEq(f64, atan2_64(0.2, -0.2), 2.356194, epsilon));
+    assert(math.approxEq(f64, atan2_64(-0.2, -0.2), -2.356194, epsilon));
+    assert(math.approxEq(f64, atan2_64(0.34, -0.4), 2.437099, epsilon));
+    assert(math.approxEq(f64, atan2_64(0.34, 1.243), 0.267001, epsilon));
+}
+
+test "math.atan2_32.special" {
+    const epsilon = 0.000001;
+
+    assert(math.isNan(atan2_32(1.0, math.nan(f32))));
+    assert(math.isNan(atan2_32(math.nan(f32), 1.0)));
+    assert(atan2_32(0.0, 5.0) == 0.0);
+    assert(atan2_32(-0.0, 5.0) == -0.0);
+    assert(math.approxEq(f32, atan2_32(0.0, -5.0), math.pi, epsilon));
+    assert(math.approxEq(f32, atan2_32(-0.0, -5.0), -math.pi, epsilon));
+    assert(math.approxEq(f32, atan2_32(1.0, 0.0), math.pi_2, epsilon));
+    assert(math.approxEq(f32, atan2_32(1.0, -0.0), math.pi_2, epsilon));
+    assert(math.approxEq(f32, atan2_32(-1.0, 0.0), -math.pi_2, epsilon));
+    assert(math.approxEq(f32, atan2_32(-1.0, -0.0), -math.pi_2, epsilon));
+    assert(math.approxEq(f32, atan2_32(math.inf(f32), math.inf(f32)), math.pi_4, epsilon));
+    assert(math.approxEq(f32, atan2_32(-math.inf(f32), math.inf(f32)), -math.pi_4, epsilon));
+    assert(math.approxEq(f32, atan2_32(math.inf(f32), -math.inf(f32)), 3.0 * math.pi_4, epsilon));
+    assert(math.approxEq(f32, atan2_32(-math.inf(f32), -math.inf(f32)), -3.0 * math.pi_4, epsilon));
+    assert(atan2_32(1.0, math.inf(f32)) == 0.0);
+    assert(math.approxEq(f32, atan2_32(1.0, -math.inf(f32)), math.pi, epsilon));
+    assert(math.approxEq(f32, atan2_32(-1.0, -math.inf(f32)), -math.pi, epsilon));
+    assert(math.approxEq(f32, atan2_32(math.inf(f32), 1.0), math.pi_2, epsilon));
+    assert(math.approxEq(f32, atan2_32(-math.inf(f32), 1.0), -math.pi_2, epsilon));
+}
+
+test "math.atan2_64.special" {
+    const epsilon = 0.000001;
+
+    assert(math.isNan(atan2_64(1.0, math.nan(f64))));
+    assert(math.isNan(atan2_64(math.nan(f64), 1.0)));
+    assert(atan2_64(0.0, 5.0) == 0.0);
+    assert(atan2_64(-0.0, 5.0) == -0.0);
+    assert(math.approxEq(f64, atan2_64(0.0, -5.0), math.pi, epsilon));
+    assert(math.approxEq(f64, atan2_64(-0.0, -5.0), -math.pi, epsilon));
+    assert(math.approxEq(f64, atan2_64(1.0, 0.0), math.pi_2, epsilon));
+    assert(math.approxEq(f64, atan2_64(1.0, -0.0), math.pi_2, epsilon));
+    assert(math.approxEq(f64, atan2_64(-1.0, 0.0), -math.pi_2, epsilon));
+    assert(math.approxEq(f64, atan2_64(-1.0, -0.0), -math.pi_2, epsilon));
+    assert(math.approxEq(f64, atan2_64(math.inf(f64), math.inf(f64)), math.pi_4, epsilon));
+    assert(math.approxEq(f64, atan2_64(-math.inf(f64), math.inf(f64)), -math.pi_4, epsilon));
+    assert(math.approxEq(f64, atan2_64(math.inf(f64), -math.inf(f64)), 3.0 * math.pi_4, epsilon));
+    assert(math.approxEq(f64, atan2_64(-math.inf(f64), -math.inf(f64)), -3.0 * math.pi_4, epsilon));
+    assert(atan2_64(1.0, math.inf(f64)) == 0.0);
+    assert(math.approxEq(f64, atan2_64(1.0, -math.inf(f64)), math.pi, epsilon));
+    assert(math.approxEq(f64, atan2_64(-1.0, -math.inf(f64)), -math.pi, epsilon));
+    assert(math.approxEq(f64, atan2_64(math.inf(f64), 1.0), math.pi_2, epsilon));
+    assert(math.approxEq(f64, atan2_64(-math.inf(f64), 1.0), -math.pi_2, epsilon));
+}

std/math/atanh.zig

@@ -0,0 +1,116 @@
+// Special Cases:
+//
+// - atanh(+-1) = +-inf with signal
+// - atanh(x)   = nan if |x| > 1 with signal
+// - atanh(nan) = nan
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const atanh = atanh_workaround;
+
+pub fn atanh_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(atanh_32, x),
+        f64 => @inlineCall(atanh_64, x),
+        else => @compileError("atanh not implemented for " ++ @typeName(T)),
+    }
+}
+
+// atanh(x) = log((1 + x) / (1 - x)) / 2 = log1p(2x / (1 - x)) / 2 ~= x + x^3 / 3 + o(x^5)
+fn atanh_32(x: f32) -> f32 {
+    const u = @bitCast(u32, x);
+    const i = u & 0x7FFFFFFF;
+    const s = u >> 31;
+
+    var y = @bitCast(f32, i); // |x|
+
+    if (y == 1.0) {
+        return math.copysign(f32, math.inf(f32), x);
+    }
+
+    if (u < 0x3F800000 - (1 << 23)) {
+        if (u < 0x3F800000 - (32 << 23)) {
+            // underflow
+            if (u < (1 << 23)) {
+                math.forceEval(y * y)
+            }
+        }
+        // |x| < 0.5
+        else {
+            y = 0.5 * math.log1p(2 * y + 2 * y * y / (1 - y));
+        }
+    } else {
+        y = 0.5 * math.log1p(2 * (y / (1 - y)));
+    }
+
+    if (s != 0) -y else y
+}
+
+fn atanh_64(x: f64) -> f64 {
+    const u = @bitCast(u64, x);
+    const e = (u >> 52) & 0x7FF;
+    const s = u >> 63;
+
+    var y = @bitCast(f64, u & (@maxValue(u64) >> 1)); // |x|
+
+    if (y == 1.0) {
+        return math.copysign(f64, math.inf(f64), x);
+    }
+
+    if (e < 0x3FF - 1) {
+        if (e < 0x3FF - 32) {
+            // underflow
+            if (e == 0) {
+                math.forceEval(f32(y));
+            }
+        }
+        // |x| < 0.5
+        else {
+            y = 0.5 * math.log1p(2 * y + 2 * y * y / (1 - y));
+        }
+    } else {
+        y = 0.5 * math.log1p(2 * (y / (1 - y)));
+    }
+
+    if (s != 0) -y else y
+}
+
+test "math.atanh" {
+    assert(atanh(f32(0.0)) == atanh_32(0.0));
+    assert(atanh(f64(0.0)) == atanh_64(0.0));
+}
+
+test "math.atanh_32" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, atanh_32(0.0), 0.0, epsilon));
+    assert(math.approxEq(f32, atanh_32(0.2), 0.202733, epsilon));
+    assert(math.approxEq(f32, atanh_32(0.8923), 1.433099, epsilon));
+}
+
+test "math.atanh_64" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f64, atanh_64(0.0), 0.0, epsilon));
+    assert(math.approxEq(f64, atanh_64(0.2), 0.202733, epsilon));
+    assert(math.approxEq(f64, atanh_64(0.8923), 1.433099, epsilon));
+}
+
+test "math.atanh32.special" {
+    assert(math.isPositiveInf(atanh_32(1)));
+    assert(math.isNegativeInf(atanh_32(-1)));
+    assert(math.isSignalNan(atanh_32(1.5)));
+    assert(math.isSignalNan(atanh_32(-1.5)));
+    assert(math.isNan(atanh_32(math.nan(f32))));
+}
+
+test "math.atanh64.special" {
+    assert(math.isPositiveInf(atanh_64(1)));
+    assert(math.isNegativeInf(atanh_64(-1)));
+    assert(math.isSignalNan(atanh_64(1.5)));
+    assert(math.isSignalNan(atanh_64(-1.5)));
+    assert(math.isNan(atanh_64(math.nan(f64))));
+}

std/math/cbrt.zig

@@ -0,0 +1,159 @@
+// Special Cases:
+//
+// - cbrt(+-0)   = +-0
+// - cbrt(+-inf) = +-inf
+// - cbrt(nan)   = nan
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const cbrt = cbrt_workaround;
+
+pub fn cbrt_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(cbrt32, x),
+        f64 => @inlineCall(cbrt64, x),
+        else => @compileError("cbrt not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn cbrt32(x: f32) -> f32 {
+    const B1: u32 = 709958130; // (127 - 127.0 / 3 - 0.03306235651) * 2^23
+    const B2: u32 = 642849266; // (127 - 127.0 / 3 - 24 / 3 - 0.03306235651) * 2^23
+
+    var u = @bitCast(u32, x);
+    var hx = u & 0x7FFFFFFF;
+
+    // cbrt(nan, inf) = itself
+    if (hx >= 0x7F800000) {
+        return x + x;
+    }
+
+    // cbrt to ~5bits
+    if (hx < 0x00800000) {
+        // cbrt(+-0) = itself
+        if (hx == 0) {
+            return x;
+        }
+        u = @bitCast(u32, x * 0x1.0p24);
+        hx = u & 0x7FFFFFFF;
+        hx = hx / 3 + B2;
+    } else {
+        hx = hx / 3 + B1;
+    }
+
+    u &= 0x80000000;
+    u |= hx;
+
+    // first step newton to 16 bits
+    var t: f64 = @bitCast(f32, u);
+    var r: f64 = t * t * t;
+    t = t * (f64(x) + x + r) / (x + r + r);
+
+    // second step newton to 47 bits
+    r = t * t * t;
+    t = t * (f64(x) + x + r) / (x + r + r);
+
+    f32(t)
+}
+
+fn cbrt64(x: f64) -> f64 {
+    const B1: u32 = 715094163;  // (1023 - 1023 / 3 - 0.03306235651 * 2^20
+    const B2: u32 = 696219795;  // (1023 - 1023 / 3 - 54 / 3 - 0.03306235651 * 2^20
+
+    // |1 / cbrt(x) - p(x)| < 2^(23.5)
+    const P0: f64 =  1.87595182427177009643;
+    const P1: f64 = -1.88497979543377169875;
+    const P2: f64 =  1.621429720105354466140;
+    const P3: f64 = -0.758397934778766047437;
+    const P4: f64 =  0.145996192886612446982;
+
+    var u = @bitCast(u64, x);
+    var hx = u32(u >> 32) & 0x7FFFFFFF;
+
+    // cbrt(nan, inf) = itself
+    if (hx >= 0x7FF00000) {
+        return x + x;
+    }
+
+    // cbrt to ~5bits
+    if (hx < 0x00100000) {
+        u = @bitCast(u64, x * 0x1.0p54);
+        hx = u32(u >> 32) & 0x7FFFFFFF;
+
+        // cbrt(0) is itself
+        if (hx == 0) {
+            return 0;
+        }
+        hx = hx / 3 + B2;
+    } else {
+        hx = hx / 3 + B1;
+    }
+
+    u &= 1 << 63;
+    u |= u64(hx) << 32;
+    var t = @bitCast(f64, u);
+
+    // cbrt to 23 bits
+    // cbrt(x) = t * cbrt(x / t^3) ~= t * P(t^3 / x)
+    var r = (t * t) * (t / x);
+    t = t * ((P0 + r * (P1 + r * P2)) + ((r * r) * r) * (P3 + r * P4));
+
+    // Round t away from 0 to 23 bits
+    u = @bitCast(u64, t);
+    u = (u + 0x80000000) & 0xFFFFFFFFC0000000;
+    t = @bitCast(f64, u);
+
+    // one step newton to 53 bits
+    const s = t * t;
+    var q = x / s;
+    var w = t + t;
+    q = (q - t) / (w + q);
+
+    t + t * q
+}
+
+test "math.cbrt" {
+    assert(cbrt(f32(0.0)) == cbrt32(0.0));
+    assert(cbrt(f64(0.0)) == cbrt64(0.0));
+}
+
+test "math.cbrt32" {
+    const epsilon = 0.000001;
+
+    assert(cbrt32(0.0) == 0.0);
+    assert(math.approxEq(f32, cbrt32(0.2), 0.584804, epsilon));
+    assert(math.approxEq(f32, cbrt32(0.8923), 0.962728, epsilon));
+    assert(math.approxEq(f32, cbrt32(1.5), 1.144714, epsilon));
+    assert(math.approxEq(f32, cbrt32(37.45), 3.345676, epsilon));
+    assert(math.approxEq(f32, cbrt32(123123.234375), 49.748501, epsilon));
+}
+
+test "math.cbrt64" {
+    const epsilon = 0.000001;
+
+    assert(cbrt64(0.0) == 0.0);
+    assert(math.approxEq(f64, cbrt64(0.2), 0.584804, epsilon));
+    assert(math.approxEq(f64, cbrt64(0.8923), 0.962728, epsilon));
+    assert(math.approxEq(f64, cbrt64(1.5), 1.144714, epsilon));
+    assert(math.approxEq(f64, cbrt64(37.45), 3.345676, epsilon));
+    assert(math.approxEq(f64, cbrt64(123123.234375), 49.748501, epsilon));
+}
+
+test "math.cbrt.special" {
+    assert(cbrt32(0.0) == 0.0);
+    assert(cbrt32(-0.0) == -0.0);
+    assert(math.isPositiveInf(cbrt32(math.inf(f32))));
+    assert(math.isNegativeInf(cbrt32(-math.inf(f32))));
+    assert(math.isNan(cbrt32(math.nan(f32))));
+}
+
+test "math.cbrt64.special" {
+    assert(cbrt64(0.0) == 0.0);
+    assert(cbrt64(-0.0) == -0.0);
+    assert(math.isPositiveInf(cbrt64(math.inf(f64))));
+    assert(math.isNegativeInf(cbrt64(-math.inf(f64))));
+    assert(math.isNan(cbrt64(math.nan(f64))));
+}

std/math/ceil.zig

@@ -0,0 +1,119 @@
+// Special Cases:
+//
+// - ceil(+-0)   = +-0
+// - ceil(+-inf) = +-inf
+// - ceil(nan)   = nan
+
+const builtin = @import("builtin");
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const ceil = ceil_workaround;
+
+pub fn ceil_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(ceil32, x),
+        f64 => @inlineCall(ceil64, x),
+        else => @compileError("ceil not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn ceil32(x: f32) -> f32 {
+    var u = @bitCast(u32, x);
+    var e = i32((u >> 23) & 0xFF) - 0x7F;
+    var m: u32 = undefined;
+
+    // TODO: Shouldn't need this explicit check.
+    if (x == 0.0) {
+        return x;
+    }
+
+    if (e >= 23) {
+        return x;
+    }
+    else if (e >= 0) {
+        m = 0x007FFFFF >> u32(e);
+        if (u & m == 0) {
+            return x;
+        }
+        math.forceEval(x + 0x1.0p120);
+        if (u >> 31 == 0) {
+            u += m;
+        }
+        u &= ~m;
+        @bitCast(f32, u)
+    } else {
+        math.forceEval(x + 0x1.0p120);
+        if (u >> 31 != 0) {
+            return -0.0;
+        } else {
+            1.0
+        }
+    }
+}
+
+fn ceil64(x: f64) -> f64 {
+    const u = @bitCast(u64, x);
+    const e = (u >> 52) & 0x7FF;
+    var y: f64 = undefined;
+
+    if (e >= 0x3FF+52 or x == 0) {
+        return x;
+    }
+
+    if (u >> 63 != 0) {
+        @setFloatMode(this, builtin.FloatMode.Strict);
+        y = x - math.f64_toint + math.f64_toint - x;
+    } else {
+        @setFloatMode(this, builtin.FloatMode.Strict);
+        y = x + math.f64_toint - math.f64_toint - x;
+    }
+
+    if (e <= 0x3FF-1) {
+        math.forceEval(y);
+        if (u >> 63 != 0) {
+            return -0.0;    // Compiler requires return.
+        } else {
+            1.0
+        }
+    } else if (y < 0) {
+        x + y + 1
+    } else {
+        x + y
+    }
+}
+
+test "math.ceil" {
+    assert(ceil(f32(0.0)) == ceil32(0.0));
+    assert(ceil(f64(0.0)) == ceil64(0.0));
+}
+
+test "math.ceil32" {
+    assert(ceil32(1.3) == 2.0);
+    assert(ceil32(-1.3) == -1.0);
+    assert(ceil32(0.2) == 1.0);
+}
+
+test "math.ceil64" {
+    assert(ceil64(1.3) == 2.0);
+    assert(ceil64(-1.3) == -1.0);
+    assert(ceil64(0.2) == 1.0);
+}
+
+test "math.ceil32.special" {
+    assert(ceil32(0.0) == 0.0);
+    assert(ceil32(-0.0) == -0.0);
+    assert(math.isPositiveInf(ceil32(math.inf(f32))));
+    assert(math.isNegativeInf(ceil32(-math.inf(f32))));
+    assert(math.isNan(ceil32(math.nan(f32))));
+}
+
+test "math.ceil64.special" {
+    assert(ceil64(0.0) == 0.0);
+    assert(ceil64(-0.0) == -0.0);
+    assert(math.isPositiveInf(ceil64(math.inf(f64))));
+    assert(math.isNegativeInf(ceil64(-math.inf(f64))));
+    assert(math.isNan(ceil64(math.nan(f64))));
+}

std/math/copysign.zig

@@ -0,0 +1,50 @@
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const copysign = copysign_workaround;
+
+pub fn copysign_workaround(comptime T: type, x: T, y: T) -> T {
+    switch (T) {
+        f32 => @inlineCall(copysign32, x, y),
+        f64 => @inlineCall(copysign64, x, y),
+        else => @compileError("copysign not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn copysign32(x: f32, y: f32) -> f32 {
+    const ux = @bitCast(u32, x);
+    const uy = @bitCast(u32, y);
+
+    const h1 = ux & (@maxValue(u32) / 2);
+    const h2 = uy & (u32(1) << 31);
+    @bitCast(f32, h1 | h2)
+}
+
+fn copysign64(x: f64, y: f64) -> f64 {
+    const ux = @bitCast(u64, x);
+    const uy = @bitCast(u64, y);
+
+    const h1 = ux & (@maxValue(u64) / 2);
+    const h2 = uy & (u64(1) << 63);
+    @bitCast(f64, h1 | h2)
+}
+
+test "math.copysign" {
+    assert(copysign(f32, 1.0, 1.0) == copysign32(1.0, 1.0));
+    assert(copysign(f64, 1.0, 1.0) == copysign64(1.0, 1.0));
+}
+
+test "math.copysign32" {
+    assert(copysign32(5.0, 1.0) == 5.0);
+    assert(copysign32(5.0, -1.0) == -5.0);
+    assert(copysign32(-5.0, -1.0) == -5.0);
+    assert(copysign32(-5.0, 1.0) == 5.0);
+}
+
+test "math.copysign64" {
+    assert(copysign64(5.0, 1.0) == 5.0);
+    assert(copysign64(5.0, -1.0) == -5.0);
+    assert(copysign64(-5.0, -1.0) == -5.0);
+    assert(copysign64(-5.0, 1.0) == 5.0);
+}

std/math/cos.zig

@@ -0,0 +1,181 @@
+// Special Cases:
+//
+// - cos(+-inf) = nan
+// - cos(nan)   = nan
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const cos = cos_workaround;
+
+pub fn cos_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(cos32, x),
+        f64 => @inlineCall(cos64, x),
+        else => @compileError("cos not implemented for " ++ @typeName(T)),
+    }
+}
+
+// sin polynomial coefficients
+const S0 =  1.58962301576546568060E-10;
+const S1 = -2.50507477628578072866E-8;
+const S2 =  2.75573136213857245213E-6;
+const S3 = -1.98412698295895385996E-4;
+const S4 =  8.33333333332211858878E-3;
+const S5 = -1.66666666666666307295E-1;
+
+// cos polynomial coeffiecients
+const C0 = -1.13585365213876817300E-11;
+const C1 =  2.08757008419747316778E-9;
+const C2 = -2.75573141792967388112E-7;
+const C3 =  2.48015872888517045348E-5;
+const C4 = -1.38888888888730564116E-3;
+const C5 =  4.16666666666665929218E-2;
+
+// NOTE: This is taken from the go stdlib. The musl implementation is much more complex.
+//
+// This may have slight differences on some edge cases and may need to replaced if so.
+fn cos32(x_: f32) -> f32 {
+    @setFloatMode(this, @import("builtin").FloatMode.Strict);
+
+    const pi4a = 7.85398125648498535156e-1;
+    const pi4b = 3.77489470793079817668E-8;
+    const pi4c = 2.69515142907905952645E-15;
+    const m4pi = 1.273239544735162542821171882678754627704620361328125;
+
+    var x = x_;
+    if (math.isNan(x) or math.isInf(x)) {
+        return math.nan(f32);
+    }
+
+    var sign = false;
+    if (x < 0) {
+        x = -x;
+    }
+
+    var y = math.floor(x * m4pi);
+    var j = i64(y);
+
+    if (j & 1 == 1) {
+        j += 1;
+        y += 1;
+    }
+
+    j &= 7;
+    if (j > 3) {
+        j -= 4;
+        sign = !sign;
+    }
+    if (j > 1) {
+        sign = !sign;
+    }
+
+    const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
+    const w = z * z;
+
+    const r = {
+        if (j == 1 or j == 2) {
+            z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))))
+        } else {
+            1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))))
+        }
+    };
+
+    if (sign) {
+        -r
+    } else {
+        r
+    }
+}
+
+fn cos64(x_: f64) -> f64 {
+    const pi4a = 7.85398125648498535156e-1;
+    const pi4b = 3.77489470793079817668E-8;
+    const pi4c = 2.69515142907905952645E-15;
+    const m4pi = 1.273239544735162542821171882678754627704620361328125;
+
+    var x = x_;
+    if (math.isNan(x) or math.isInf(x)) {
+        return math.nan(f64);
+    }
+
+    var sign = false;
+    if (x < 0) {
+        x = -x;
+    }
+
+    var y = math.floor(x * m4pi);
+    var j = i64(y);
+
+    if (j & 1 == 1) {
+        j += 1;
+        y += 1;
+    }
+
+    j &= 7;
+    if (j > 3) {
+        j -= 4;
+        sign = !sign;
+    }
+    if (j > 1) {
+        sign = !sign;
+    }
+
+    const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
+    const w = z * z;
+
+    const r = {
+        if (j == 1 or j == 2) {
+            z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))))
+        } else {
+            1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))))
+        }
+    };
+
+    if (sign) {
+        -r
+    } else {
+        r
+    }
+}
+
+test "math.cos" {
+    assert(cos(f32(0.0)) == cos32(0.0));
+    assert(cos(f64(0.0)) == cos64(0.0));
+}
+
+test "math.cos32" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, cos32(0.0), 1.0, epsilon));
+    assert(math.approxEq(f32, cos32(0.2), 0.980067, epsilon));
+    assert(math.approxEq(f32, cos32(0.8923), 0.627623, epsilon));
+    assert(math.approxEq(f32, cos32(1.5), 0.070737, epsilon));
+    assert(math.approxEq(f32, cos32(37.45), 0.969132, epsilon));
+    assert(math.approxEq(f32, cos32(89.123), 0.400798, epsilon));
+}
+
+test "math.cos64" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f64, cos64(0.0), 1.0, epsilon));
+    assert(math.approxEq(f64, cos64(0.2), 0.980067, epsilon));
+    assert(math.approxEq(f64, cos64(0.8923), 0.627623, epsilon));
+    assert(math.approxEq(f64, cos64(1.5), 0.070737, epsilon));
+    assert(math.approxEq(f64, cos64(37.45), 0.969132, epsilon));
+    assert(math.approxEq(f64, cos64(89.123), 0.40080, epsilon));
+}
+
+test "math.cos32.special" {
+    assert(math.isNan(cos32(math.inf(f32))));
+    assert(math.isNan(cos32(-math.inf(f32))));
+    assert(math.isNan(cos32(math.nan(f32))));
+}
+
+test "math.cos64.special" {
+    assert(math.isNan(cos64(math.inf(f64))));
+    assert(math.isNan(cos64(-math.inf(f64))));
+    assert(math.isNan(cos64(math.nan(f64))));
+}

std/math/cosh.zig

@@ -0,0 +1,121 @@
+// Special Cases:
+//
+// - cosh(+-0)   = 1
+// - cosh(+-inf) = +inf
+// - cosh(nan)   = nan
+
+const math = @import("index.zig");
+const expo2 = @import("expo2.zig").expo2;
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const cosh = cosh_workaround;
+
+pub fn cosh_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(cosh32, x),
+        f64 => @inlineCall(cosh64, x),
+        else => @compileError("cosh not implemented for " ++ @typeName(T)),
+    }
+}
+
+// cosh(x) = (exp(x) + 1 / exp(x)) / 2
+//         = 1 + 0.5 * (exp(x) - 1) * (exp(x) - 1) / exp(x)
+//         = 1 + (x * x) / 2 + o(x^4)
+fn cosh32(x: f32) -> f32 {
+    const u = @bitCast(u32, x);
+    const ux = u & 0x7FFFFFFF;
+    const ax = @bitCast(f32, ux);
+
+    // |x| < log(2)
+    if (ux < 0x3F317217) {
+        if (ux < 0x3F800000 - (12 << 23)) {
+            math.raiseOverflow();
+            return 1.0;
+        }
+        const t = math.expm1(ax);
+        return 1 + t * t / (2 * (1 + t));
+    }
+
+    // |x| < log(FLT_MAX)
+    if (ux < 0x42B17217) {
+        const t = math.exp(ax);
+        return 0.5 * (t + 1 / t);
+    }
+
+    // |x| > log(FLT_MAX) or nan
+    expo2(ax)
+}
+
+fn cosh64(x: f64) -> f64 {
+    const u = @bitCast(u64, x);
+    const w = u32(u >> 32);
+    const ax = @bitCast(f64, u & (@maxValue(u64) >> 1));
+
+    // TODO: Shouldn't need this explicit check.
+    if (x == 0.0) {
+        return 1.0;
+    }
+
+    // |x| < log(2)
+    if (w < 0x3FE62E42) {
+        if (w < 0x3FF00000 - (26 << 20)) {
+            if (x != 0) {
+                math.raiseInexact();
+            }
+            return 1.0;
+        }
+        const t = math.expm1(ax);
+        return 1 + t * t / (2 * (1 + t));
+    }
+
+    // |x| < log(DBL_MAX)
+    if (w < 0x40862E42) {
+        const t = math.exp(ax);
+        // NOTE: If x > log(0x1p26) then 1/t is not required.
+        return 0.5 * (t + 1 / t);
+    }
+
+    // |x| > log(CBL_MAX) or nan
+    expo2(ax)
+}
+
+test "math.cosh" {
+    assert(cosh(f32(1.5)) == cosh32(1.5));
+    assert(cosh(f64(1.5)) == cosh64(1.5));
+}
+
+test "math.cosh32" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, cosh32(0.0), 1.0, epsilon));
+    assert(math.approxEq(f32, cosh32(0.2), 1.020067, epsilon));
+    assert(math.approxEq(f32, cosh32(0.8923), 1.425225, epsilon));
+    assert(math.approxEq(f32, cosh32(1.5), 2.352410, epsilon));
+}
+
+test "math.cosh64" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f64, cosh64(0.0), 1.0, epsilon));
+    assert(math.approxEq(f64, cosh64(0.2), 1.020067, epsilon));
+    assert(math.approxEq(f64, cosh64(0.8923), 1.425225, epsilon));
+    assert(math.approxEq(f64, cosh64(1.5), 2.352410, epsilon));
+}
+
+test "math.cosh32.special" {
+    assert(cosh32(0.0) == 1.0);
+    assert(cosh32(-0.0) == 1.0);
+    assert(math.isPositiveInf(cosh32(math.inf(f32))));
+    assert(math.isPositiveInf(cosh32(-math.inf(f32))));
+    assert(math.isNan(cosh32(math.nan(f32))));
+}
+
+test "math.cosh64.special" {
+    assert(cosh64(0.0) == 1.0);
+    assert(cosh64(-0.0) == 1.0);
+    assert(math.isPositiveInf(cosh64(math.inf(f64))));
+    assert(math.isPositiveInf(cosh64(-math.inf(f64))));
+    assert(math.isNan(cosh64(math.nan(f64))));
+}

std/math/exp.zig

@@ -0,0 +1,217 @@
+// Special Cases:
+//
+// - exp(+inf) = +inf
+// - exp(nan)  = nan
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const exp = exp_workaround;
+
+pub fn exp_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(exp32, x),
+        f64 => @inlineCall(exp64, x),
+        else => @compileError("exp not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn exp32(x_: f32) -> f32 {
+    const half = []f32 { 0.5, -0.5 };
+    const ln2hi = 6.9314575195e-1;
+    const ln2lo = 1.4286067653e-6;
+    const invln2  = 1.4426950216e+0;
+    const P1 = 1.6666625440e-1;
+    const P2 = -2.7667332906e-3;
+
+    var x = x_;
+    var hx = @bitCast(u32, x);
+    const sign = i32(hx >> 31);
+    hx &= 0x7FFFFFFF;
+
+    if (math.isNan(x)) {
+        return x;
+    }
+
+    // |x| >= -87.33655 or nan
+    if (hx >= 0x42AEAC50) {
+        // nan
+        if (hx > 0x7F800000) {
+            return x;
+        }
+        // x >= 88.722839
+        if (hx >= 0x42b17218 and sign == 0) {
+            return x * 0x1.0p127;
+        }
+        if (sign != 0) {
+            math.forceEval(-0x1.0p-149 / x);   // overflow
+            // x <= -103.972084
+            if (hx >= 0x42CFF1B5) {
+                return 0;
+            }
+        }
+    }
+
+    var k: i32 = undefined;
+    var hi: f32 = undefined;
+    var lo: f32 = undefined;
+
+    // |x| > 0.5 * ln2
+    if (hx > 0x3EB17218) {
+        // |x| > 1.5 * ln2
+        if (hx > 0x3F851592) {
+            k = i32(invln2 * x + half[usize(sign)]);
+        }
+        else {
+            k = 1 - sign - sign;
+        }
+
+        const fk = f32(k);
+        hi = x - fk * ln2hi;
+        lo = fk * ln2lo;
+        x = hi - lo;
+    }
+    // |x| > 2^(-14)
+    else if (hx > 0x39000000) {
+        k = 0;
+        hi = x;
+        lo = 0;
+    }
+    else {
+        math.forceEval(0x1.0p127 + x); // inexact
+        return 1 + x;
+    }
+
+    const xx = x * x;
+    const c = x - xx * (P1 + xx * P2);
+    const y = 1 + (x * c / (2 - c) - lo + hi);
+
+    if (k == 0) {
+        y
+    } else {
+        math.scalbn(y, k)
+    }
+}
+
+fn exp64(x_: f64) -> f64 {
+    const half = []const f64 { 0.5, -0.5 };
+    const ln2hi: f64 = 6.93147180369123816490e-01;
+    const ln2lo: f64 = 1.90821492927058770002e-10;
+    const invln2: f64 = 1.44269504088896338700e+00;
+    const P1: f64   =  1.66666666666666019037e-01;
+    const P2: f64   = -2.77777777770155933842e-03;
+    const P3: f64   =  6.61375632143793436117e-05;
+    const P4: f64   = -1.65339022054652515390e-06;
+    const P5: f64   =  4.13813679705723846039e-08;
+
+    var x = x_;
+    var ux = @bitCast(u64, x);
+    var hx = ux >> 32;
+    const sign = i32(hx >> 31);
+    hx &= 0x7FFFFFFF;
+
+    if (math.isNan(x)) {
+        return x;
+    }
+
+    // |x| >= 708.39 or nan
+    if (hx >= 0x4086232B) {
+        // nan
+        if (hx > 0x7FF00000) {
+            return x;
+        }
+        if (x > 709.782712893383973096) {
+            // overflow if x != inf
+            if (!math.isInf(x)) {
+                math.raiseOverflow();
+            }
+            return math.inf(f64);
+        }
+        if (x < -708.39641853226410622) {
+            // underflow if x != -inf
+            // math.forceEval(f32(-0x1.0p-149 / x));
+            if (x < -745.13321910194110842) {
+                return 0;
+            }
+        }
+    }
+
+    // argument reduction
+    var k: i32 = undefined;
+    var hi: f64 = undefined;
+    var lo: f64 = undefined;
+
+    // |x| > 0.5 * ln2
+    if (hx > 0x3EB17218) {
+        // |x| >= 1.5 * ln2
+        if (hx > 0x3FF0A2B2) {
+            k = i32(invln2 * x + half[usize(sign)]);
+        }
+        else {
+            k = 1 - sign - sign;
+        }
+
+        const dk = f64(k);
+        hi = x - dk * ln2hi;
+        lo = dk * ln2lo;
+        x = hi - lo;
+    }
+    // |x| > 2^(-28)
+    else if (hx > 0x3E300000) {
+        k = 0;
+        hi = x;
+        lo = 0;
+    }
+    else {
+        // inexact if x != 0
+        // math.forceEval(0x1.0p1023 + x);
+        return 1 + x;
+    }
+
+    const xx = x * x;
+    const c = x - xx * (P1 + xx * (P2 + xx * (P3 + xx * (P4 + xx * P5))));
+    const y = 1 + (x * c / (2 - c) - lo + hi);
+
+    if (k == 0) {
+        y
+    } else {
+        math.scalbn(y, k)
+    }
+}
+
+test "math.exp" {
+    assert(exp(f32(0.0)) == exp32(0.0));
+    assert(exp(f64(0.0)) == exp64(0.0));
+}
+
+test "math.exp32" {
+    const epsilon = 0.000001;
+
+    assert(exp32(0.0) == 1.0);
+    assert(math.approxEq(f32, exp32(0.0), 1.0, epsilon));
+    assert(math.approxEq(f32, exp32(0.2), 1.221403, epsilon));
+    assert(math.approxEq(f32, exp32(0.8923), 2.440737, epsilon));
+    assert(math.approxEq(f32, exp32(1.5), 4.481689, epsilon));
+}
+
+test "math.exp64" {
+    const epsilon = 0.000001;
+
+    assert(exp64(0.0) == 1.0);
+    assert(math.approxEq(f64, exp64(0.0), 1.0, epsilon));
+    assert(math.approxEq(f64, exp64(0.2), 1.221403, epsilon));
+    assert(math.approxEq(f64, exp64(0.8923), 2.440737, epsilon));
+    assert(math.approxEq(f64, exp64(1.5), 4.481689, epsilon));
+}
+
+test "math.exp32.special" {
+    assert(math.isPositiveInf(exp32(math.inf(f32))));
+    assert(math.isNan(exp32(math.nan(f32))));
+}
+
+test "math.exp64.special" {
+    assert(math.isPositiveInf(exp64(math.inf(f64))));
+    assert(math.isNan(exp64(math.nan(f64))));
+}

std/math/exp2.zig

@@ -0,0 +1,455 @@
+// Special Cases:
+//
+// - exp2(+inf) = +inf
+// - exp2(nan)  = nan
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const exp2 = exp2_workaround;
+
+pub fn exp2_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(exp2_32, x),
+        f64 => @inlineCall(exp2_64, x),
+        else => @compileError("exp2 not implemented for " ++ @typeName(T)),
+    }
+}
+
+const exp2ft = []const f64 {
+    0x1.6a09e667f3bcdp-1,
+    0x1.7a11473eb0187p-1,
+    0x1.8ace5422aa0dbp-1,
+    0x1.9c49182a3f090p-1,
+    0x1.ae89f995ad3adp-1,
+    0x1.c199bdd85529cp-1,
+    0x1.d5818dcfba487p-1,
+    0x1.ea4afa2a490dap-1,
+    0x1.0000000000000p+0,
+    0x1.0b5586cf9890fp+0,
+    0x1.172b83c7d517bp+0,
+    0x1.2387a6e756238p+0,
+    0x1.306fe0a31b715p+0,
+    0x1.3dea64c123422p+0,
+    0x1.4bfdad5362a27p+0,
+    0x1.5ab07dd485429p+0,
+};
+
+fn exp2_32(x: f32) -> f32 {
+    @setFloatMode(this, @import("builtin").FloatMode.Strict);
+
+    const tblsiz = u32(exp2ft.len);
+    const redux: f32 = 0x1.8p23 / f32(tblsiz);
+    const P1: f32 = 0x1.62e430p-1;
+    const P2: f32 = 0x1.ebfbe0p-3;
+    const P3: f32 = 0x1.c6b348p-5;
+    const P4: f32 = 0x1.3b2c9cp-7;
+
+    var u = @bitCast(u32, x);
+    const ix = u & 0x7FFFFFFF;
+
+    // |x| > 126
+    if (ix > 0x42FC0000) {
+        // nan
+        if (ix > 0x7F800000) {
+            return x;
+        }
+        // x >= 128
+        if (u >= 0x43000000 and u < 0x80000000) {
+            return x * 0x1.0p127;
+        }
+        // x < -126
+        if (u >= 0x80000000) {
+            if (u >= 0xC3160000 or u & 0x000FFFF != 0) {
+                math.forceEval(-0x1.0p-149 / x);
+            }
+            // x <= -150
+            if (u >= 0x3160000) {
+                return 0;
+            }
+        }
+    }
+    // |x| <= 0x1p-25
+    else if (ix <= 0x33000000) {
+        return 1.0 + x;
+    }
+
+    var uf = x + redux;
+    var i0 = @bitCast(u32, uf);
+    i0 += tblsiz / 2;
+
+    const k = i0 / tblsiz;
+    // NOTE: musl relies on undefined overflow shift behaviour. Appears that this produces the
+    // intended result but should confirm how GCC/Clang handle this to ensure.
+    const uk = @bitCast(f64, u64(0x3FF + k) <<% 52);
+    i0 &= tblsiz - 1;
+    uf -= redux;
+
+    const z: f64 = x - uf;
+    var r: f64 = exp2ft[i0];
+    const t: f64 = r * z;
+    r = r + t * (P1 + z * P2) + t * (z * z) * (P3 + z * P4);
+    f32(r * uk)
+}
+
+const exp2dt = []f64 {
+    //  exp2(z + eps)          eps
+    0x1.6a09e667f3d5dp-1,  0x1.9880p-44,
+    0x1.6b052fa751744p-1,  0x1.8000p-50,
+    0x1.6c012750bd9fep-1, -0x1.8780p-45,
+    0x1.6cfdcddd476bfp-1,  0x1.ec00p-46,
+    0x1.6dfb23c651a29p-1, -0x1.8000p-50,
+    0x1.6ef9298593ae3p-1, -0x1.c000p-52,
+    0x1.6ff7df9519386p-1, -0x1.fd80p-45,
+    0x1.70f7466f42da3p-1, -0x1.c880p-45,
+    0x1.71f75e8ec5fc3p-1,  0x1.3c00p-46,
+    0x1.72f8286eacf05p-1, -0x1.8300p-44,
+    0x1.73f9a48a58152p-1, -0x1.0c00p-47,
+    0x1.74fbd35d7ccfcp-1,  0x1.f880p-45,
+    0x1.75feb564267f1p-1,  0x1.3e00p-47,
+    0x1.77024b1ab6d48p-1, -0x1.7d00p-45,
+    0x1.780694fde5d38p-1, -0x1.d000p-50,
+    0x1.790b938ac1d00p-1,  0x1.3000p-49,
+    0x1.7a11473eb0178p-1, -0x1.d000p-49,
+    0x1.7b17b0976d060p-1,  0x1.0400p-45,
+    0x1.7c1ed0130c133p-1,  0x1.0000p-53,
+    0x1.7d26a62ff8636p-1, -0x1.6900p-45,
+    0x1.7e2f336cf4e3bp-1, -0x1.2e00p-47,
+    0x1.7f3878491c3e8p-1, -0x1.4580p-45,
+    0x1.80427543e1b4ep-1,  0x1.3000p-44,
+    0x1.814d2add1071ap-1,  0x1.f000p-47,
+    0x1.82589994ccd7ep-1, -0x1.1c00p-45,
+    0x1.8364c1eb942d0p-1,  0x1.9d00p-45,
+    0x1.8471a4623cab5p-1,  0x1.7100p-43,
+    0x1.857f4179f5bbcp-1,  0x1.2600p-45,
+    0x1.868d99b4491afp-1, -0x1.2c40p-44,
+    0x1.879cad931a395p-1, -0x1.3000p-45,
+    0x1.88ac7d98a65b8p-1, -0x1.a800p-45,
+    0x1.89bd0a4785800p-1, -0x1.d000p-49,
+    0x1.8ace5422aa223p-1,  0x1.3280p-44,
+    0x1.8be05bad619fap-1,  0x1.2b40p-43,
+    0x1.8cf3216b54383p-1, -0x1.ed00p-45,
+    0x1.8e06a5e08664cp-1, -0x1.0500p-45,
+    0x1.8f1ae99157807p-1,  0x1.8280p-45,
+    0x1.902fed0282c0ep-1, -0x1.cb00p-46,
+    0x1.9145b0b91ff96p-1, -0x1.5e00p-47,
+    0x1.925c353aa2ff9p-1,  0x1.5400p-48,
+    0x1.93737b0cdc64ap-1,  0x1.7200p-46,
+    0x1.948b82b5f98aep-1, -0x1.9000p-47,
+    0x1.95a44cbc852cbp-1,  0x1.5680p-45,
+    0x1.96bdd9a766f21p-1, -0x1.6d00p-44,
+    0x1.97d829fde4e2ap-1, -0x1.1000p-47,
+    0x1.98f33e47a23a3p-1,  0x1.d000p-45,
+    0x1.9a0f170ca0604p-1, -0x1.8a40p-44,
+    0x1.9b2bb4d53ff89p-1,  0x1.55c0p-44,
+    0x1.9c49182a3f15bp-1,  0x1.6b80p-45,
+    0x1.9d674194bb8c5p-1, -0x1.c000p-49,
+    0x1.9e86319e3238ep-1,  0x1.7d00p-46,
+    0x1.9fa5e8d07f302p-1,  0x1.6400p-46,
+    0x1.a0c667b5de54dp-1, -0x1.5000p-48,
+    0x1.a1e7aed8eb8f6p-1,  0x1.9e00p-47,
+    0x1.a309bec4a2e27p-1,  0x1.ad80p-45,
+    0x1.a42c980460a5dp-1, -0x1.af00p-46,
+    0x1.a5503b23e259bp-1,  0x1.b600p-47,
+    0x1.a674a8af46213p-1,  0x1.8880p-44,
+    0x1.a799e1330b3a7p-1,  0x1.1200p-46,
+    0x1.a8bfe53c12e8dp-1,  0x1.6c00p-47,
+    0x1.a9e6b5579fcd2p-1, -0x1.9b80p-45,
+    0x1.ab0e521356fb8p-1,  0x1.b700p-45,
+    0x1.ac36bbfd3f381p-1,  0x1.9000p-50,
+    0x1.ad5ff3a3c2780p-1,  0x1.4000p-49,
+    0x1.ae89f995ad2a3p-1, -0x1.c900p-45,
+    0x1.afb4ce622f367p-1,  0x1.6500p-46,
+    0x1.b0e07298db790p-1,  0x1.fd40p-45,
+    0x1.b20ce6c9a89a9p-1,  0x1.2700p-46,
+    0x1.b33a2b84f1a4bp-1,  0x1.d470p-43,
+    0x1.b468415b747e7p-1, -0x1.8380p-44,
+    0x1.b59728de5593ap-1,  0x1.8000p-54,
+    0x1.b6c6e29f1c56ap-1,  0x1.ad00p-47,
+    0x1.b7f76f2fb5e50p-1,  0x1.e800p-50,
+    0x1.b928cf22749b2p-1, -0x1.4c00p-47,
+    0x1.ba5b030a10603p-1, -0x1.d700p-47,
+    0x1.bb8e0b79a6f66p-1,  0x1.d900p-47,
+    0x1.bcc1e904bc1ffp-1,  0x1.2a00p-47,
+    0x1.bdf69c3f3a16fp-1, -0x1.f780p-46,
+    0x1.bf2c25bd71db8p-1, -0x1.0a00p-46,
+    0x1.c06286141b2e9p-1, -0x1.1400p-46,
+    0x1.c199bdd8552e0p-1,  0x1.be00p-47,
+    0x1.c2d1cd9fa64eep-1, -0x1.9400p-47,
+    0x1.c40ab5fffd02fp-1, -0x1.ed00p-47,
+    0x1.c544778fafd15p-1,  0x1.9660p-44,
+    0x1.c67f12e57d0cbp-1, -0x1.a100p-46,
+    0x1.c7ba88988c1b6p-1, -0x1.8458p-42,
+    0x1.c8f6d9406e733p-1, -0x1.a480p-46,
+    0x1.ca3405751c4dfp-1,  0x1.b000p-51,
+    0x1.cb720dcef9094p-1,  0x1.1400p-47,
+    0x1.ccb0f2e6d1689p-1,  0x1.0200p-48,
+    0x1.cdf0b555dc412p-1,  0x1.3600p-48,
+    0x1.cf3155b5bab3bp-1, -0x1.6900p-47,
+    0x1.d072d4a0789bcp-1,  0x1.9a00p-47,
+    0x1.d1b532b08c8fap-1, -0x1.5e00p-46,
+    0x1.d2f87080d8a85p-1,  0x1.d280p-46,
+    0x1.d43c8eacaa203p-1,  0x1.1a00p-47,
+    0x1.d5818dcfba491p-1,  0x1.f000p-50,
+    0x1.d6c76e862e6a1p-1, -0x1.3a00p-47,
+    0x1.d80e316c9834ep-1, -0x1.cd80p-47,
+    0x1.d955d71ff6090p-1,  0x1.4c00p-48,
+    0x1.da9e603db32aep-1,  0x1.f900p-48,
+    0x1.dbe7cd63a8325p-1,  0x1.9800p-49,
+    0x1.dd321f301b445p-1, -0x1.5200p-48,
+    0x1.de7d5641c05bfp-1, -0x1.d700p-46,
+    0x1.dfc97337b9aecp-1, -0x1.6140p-46,
+    0x1.e11676b197d5ep-1,  0x1.b480p-47,
+    0x1.e264614f5a3e7p-1,  0x1.0ce0p-43,
+    0x1.e3b333b16ee5cp-1,  0x1.c680p-47,
+    0x1.e502ee78b3fb4p-1, -0x1.9300p-47,
+    0x1.e653924676d68p-1, -0x1.5000p-49,
+    0x1.e7a51fbc74c44p-1, -0x1.7f80p-47,
+    0x1.e8f7977cdb726p-1, -0x1.3700p-48,
+    0x1.ea4afa2a490e8p-1,  0x1.5d00p-49,
+    0x1.eb9f4867ccae4p-1,  0x1.61a0p-46,
+    0x1.ecf482d8e680dp-1,  0x1.5500p-48,
+    0x1.ee4aaa2188514p-1,  0x1.6400p-51,
+    0x1.efa1bee615a13p-1, -0x1.e800p-49,
+    0x1.f0f9c1cb64106p-1, -0x1.a880p-48,
+    0x1.f252b376bb963p-1, -0x1.c900p-45,
+    0x1.f3ac948dd7275p-1,  0x1.a000p-53,
+    0x1.f50765b6e4524p-1, -0x1.4f00p-48,
+    0x1.f6632798844fdp-1,  0x1.a800p-51,
+    0x1.f7bfdad9cbe38p-1,  0x1.abc0p-48,
+    0x1.f91d802243c82p-1, -0x1.4600p-50,
+    0x1.fa7c1819e908ep-1, -0x1.b0c0p-47,
+    0x1.fbdba3692d511p-1, -0x1.0e00p-51,
+    0x1.fd3c22b8f7194p-1, -0x1.0de8p-46,
+    0x1.fe9d96b2a23eep-1,  0x1.e430p-49,
+    0x1.0000000000000p+0,  0x0.0000p+0,
+    0x1.00b1afa5abcbep+0, -0x1.3400p-52,
+    0x1.0163da9fb3303p+0, -0x1.2170p-46,
+    0x1.02168143b0282p+0,  0x1.a400p-52,
+    0x1.02c9a3e77806cp+0,  0x1.f980p-49,
+    0x1.037d42e11bbcap+0, -0x1.7400p-51,
+    0x1.04315e86e7f89p+0,  0x1.8300p-50,
+    0x1.04e5f72f65467p+0, -0x1.a3f0p-46,
+    0x1.059b0d315855ap+0, -0x1.2840p-47,
+    0x1.0650a0e3c1f95p+0,  0x1.1600p-48,
+    0x1.0706b29ddf71ap+0,  0x1.5240p-46,
+    0x1.07bd42b72a82dp+0, -0x1.9a00p-49,
+    0x1.0874518759bd0p+0,  0x1.6400p-49,
+    0x1.092bdf66607c8p+0, -0x1.0780p-47,
+    0x1.09e3ecac6f383p+0, -0x1.8000p-54,
+    0x1.0a9c79b1f3930p+0,  0x1.fa00p-48,
+    0x1.0b5586cf988fcp+0, -0x1.ac80p-48,
+    0x1.0c0f145e46c8ap+0,  0x1.9c00p-50,
+    0x1.0cc922b724816p+0,  0x1.5200p-47,
+    0x1.0d83b23395dd8p+0, -0x1.ad00p-48,
+    0x1.0e3ec32d3d1f3p+0,  0x1.bac0p-46,
+    0x1.0efa55fdfa9a6p+0, -0x1.4e80p-47,
+    0x1.0fb66affed2f0p+0, -0x1.d300p-47,
+    0x1.1073028d7234bp+0,  0x1.1500p-48,
+    0x1.11301d0125b5bp+0,  0x1.c000p-49,
+    0x1.11edbab5e2af9p+0,  0x1.6bc0p-46,
+    0x1.12abdc06c31d5p+0,  0x1.8400p-49,
+    0x1.136a814f2047dp+0, -0x1.ed00p-47,
+    0x1.1429aaea92de9p+0,  0x1.8e00p-49,
+    0x1.14e95934f3138p+0,  0x1.b400p-49,
+    0x1.15a98c8a58e71p+0,  0x1.5300p-47,
+    0x1.166a45471c3dfp+0,  0x1.3380p-47,
+    0x1.172b83c7d5211p+0,  0x1.8d40p-45,
+    0x1.17ed48695bb9fp+0, -0x1.5d00p-47,
+    0x1.18af9388c8d93p+0, -0x1.c880p-46,
+    0x1.1972658375d66p+0,  0x1.1f00p-46,
+    0x1.1a35beb6fcba7p+0,  0x1.0480p-46,
+    0x1.1af99f81387e3p+0, -0x1.7390p-43,
+    0x1.1bbe084045d54p+0,  0x1.4e40p-45,
+    0x1.1c82f95281c43p+0, -0x1.a200p-47,
+    0x1.1d4873168b9b2p+0,  0x1.3800p-49,
+    0x1.1e0e75eb44031p+0,  0x1.ac00p-49,
+    0x1.1ed5022fcd938p+0,  0x1.1900p-47,
+    0x1.1f9c18438cdf7p+0, -0x1.b780p-46,
+    0x1.2063b88628d8fp+0,  0x1.d940p-45,
+    0x1.212be3578a81ep+0,  0x1.8000p-50,
+    0x1.21f49917ddd41p+0,  0x1.b340p-45,
+    0x1.22bdda2791323p+0,  0x1.9f80p-46,
+    0x1.2387a6e7561e7p+0, -0x1.9c80p-46,
+    0x1.2451ffb821427p+0,  0x1.2300p-47,
+    0x1.251ce4fb2a602p+0, -0x1.3480p-46,
+    0x1.25e85711eceb0p+0,  0x1.2700p-46,
+    0x1.26b4565e27d16p+0,  0x1.1d00p-46,
+    0x1.2780e341de00fp+0,  0x1.1ee0p-44,
+    0x1.284dfe1f5633ep+0, -0x1.4c00p-46,
+    0x1.291ba7591bb30p+0, -0x1.3d80p-46,
+    0x1.29e9df51fdf09p+0,  0x1.8b00p-47,
+    0x1.2ab8a66d10e9bp+0, -0x1.27c0p-45,
+    0x1.2b87fd0dada3ap+0,  0x1.a340p-45,
+    0x1.2c57e39771af9p+0, -0x1.0800p-46,
+    0x1.2d285a6e402d9p+0, -0x1.ed00p-47,
+    0x1.2df961f641579p+0, -0x1.4200p-48,
+    0x1.2ecafa93e2ecfp+0, -0x1.4980p-45,
+    0x1.2f9d24abd8822p+0, -0x1.6300p-46,
+    0x1.306fe0a31b625p+0, -0x1.2360p-44,
+    0x1.31432edeea50bp+0, -0x1.0df8p-40,
+    0x1.32170fc4cd7b8p+0, -0x1.2480p-45,
+    0x1.32eb83ba8e9a2p+0, -0x1.5980p-45,
+    0x1.33c08b2641766p+0,  0x1.ed00p-46,
+    0x1.3496266e3fa27p+0, -0x1.c000p-50,
+    0x1.356c55f929f0fp+0, -0x1.0d80p-44,
+    0x1.36431a2de88b9p+0,  0x1.2c80p-45,
+    0x1.371a7373aaa39p+0,  0x1.0600p-45,
+    0x1.37f26231e74fep+0, -0x1.6600p-46,
+    0x1.38cae6d05d838p+0, -0x1.ae00p-47,
+    0x1.39a401b713ec3p+0, -0x1.4720p-43,
+    0x1.3a7db34e5a020p+0,  0x1.8200p-47,
+    0x1.3b57fbfec6e95p+0,  0x1.e800p-44,
+    0x1.3c32dc313a8f2p+0,  0x1.f800p-49,
+    0x1.3d0e544ede122p+0, -0x1.7a00p-46,
+    0x1.3dea64c1234bbp+0,  0x1.6300p-45,
+    0x1.3ec70df1c4eccp+0, -0x1.8a60p-43,
+    0x1.3fa4504ac7e8cp+0, -0x1.cdc0p-44,
+    0x1.40822c367a0bbp+0,  0x1.5b80p-45,
+    0x1.4160a21f72e95p+0,  0x1.ec00p-46,
+    0x1.423fb27094646p+0, -0x1.3600p-46,
+    0x1.431f5d950a920p+0,  0x1.3980p-45,
+    0x1.43ffa3f84b9ebp+0,  0x1.a000p-48,
+    0x1.44e0860618919p+0, -0x1.6c00p-48,
+    0x1.45c2042a7d201p+0, -0x1.bc00p-47,
+    0x1.46a41ed1d0016p+0, -0x1.2800p-46,
+    0x1.4786d668b3326p+0,  0x1.0e00p-44,
+    0x1.486a2b5c13c00p+0, -0x1.d400p-45,
+    0x1.494e1e192af04p+0,  0x1.c200p-47,
+    0x1.4a32af0d7d372p+0, -0x1.e500p-46,
+    0x1.4b17dea6db801p+0,  0x1.7800p-47,
+    0x1.4bfdad53629e1p+0, -0x1.3800p-46,
+    0x1.4ce41b817c132p+0,  0x1.0800p-47,
+    0x1.4dcb299fddddbp+0,  0x1.c700p-45,
+    0x1.4eb2d81d8ab96p+0, -0x1.ce00p-46,
+    0x1.4f9b2769d2d02p+0,  0x1.9200p-46,
+    0x1.508417f4531c1p+0, -0x1.8c00p-47,
+    0x1.516daa2cf662ap+0, -0x1.a000p-48,
+    0x1.5257de83f51eap+0,  0x1.a080p-43,
+    0x1.5342b569d4edap+0, -0x1.6d80p-45,
+    0x1.542e2f4f6ac1ap+0, -0x1.2440p-44,
+    0x1.551a4ca5d94dbp+0,  0x1.83c0p-43,
+    0x1.56070dde9116bp+0,  0x1.4b00p-45,
+    0x1.56f4736b529dep+0,  0x1.15a0p-43,
+    0x1.57e27dbe2c40ep+0, -0x1.9e00p-45,
+    0x1.58d12d497c76fp+0, -0x1.3080p-45,
+    0x1.59c0827ff0b4cp+0,  0x1.dec0p-43,
+    0x1.5ab07dd485427p+0, -0x1.4000p-51,
+    0x1.5ba11fba87af4p+0,  0x1.0080p-44,
+    0x1.5c9268a59460bp+0, -0x1.6c80p-45,
+    0x1.5d84590998e3fp+0,  0x1.69a0p-43,
+    0x1.5e76f15ad20e1p+0, -0x1.b400p-46,
+    0x1.5f6a320dcebcap+0,  0x1.7700p-46,
+    0x1.605e1b976dcb8p+0,  0x1.6f80p-45,
+    0x1.6152ae6cdf715p+0,  0x1.1000p-47,
+    0x1.6247eb03a5531p+0, -0x1.5d00p-46,
+    0x1.633dd1d1929b5p+0, -0x1.2d00p-46,
+    0x1.6434634ccc313p+0, -0x1.a800p-49,
+    0x1.652b9febc8efap+0, -0x1.8600p-45,
+    0x1.6623882553397p+0,  0x1.1fe0p-40,
+    0x1.671c1c708328ep+0, -0x1.7200p-44,
+    0x1.68155d44ca97ep+0,  0x1.6800p-49,
+    0x1.690f4b19e9471p+0, -0x1.9780p-45,
+};
+
+fn exp2_64(x: f64) -> f64 {
+    @setFloatMode(this, @import("builtin").FloatMode.Strict);
+
+    const tblsiz     = u32(exp2dt.len / 2);
+    const redux: f64 = 0x1.8p52 / f64(tblsiz);
+    const P1: f64    = 0x1.62e42fefa39efp-1;
+    const P2: f64    = 0x1.ebfbdff82c575p-3;
+    const P3: f64    = 0x1.c6b08d704a0a6p-5;
+    const P4: f64    = 0x1.3b2ab88f70400p-7;
+    const P5: f64    = 0x1.5d88003875c74p-10;
+
+    const ux = @bitCast(u64, x);
+    const ix = u32(ux >> 32) & 0x7FFFFFFF;
+
+    // TODO: This should be handled beneath.
+    if (math.isNan(x)) {
+        return math.nan(f64);
+    }
+
+    // |x| >= 1022 or nan
+    if (ix >= 0x408FF000) {
+        // x >= 1024 or nan
+        if (ix >= 0x40900000 and ux >> 63 == 0) {
+            math.raiseOverflow();
+            return math.inf(f64);
+        }
+        // -inf or -nan
+        if (ix >= 0x7FF00000) {
+            return -1 / x;
+        }
+        // x <= -1022
+        if (ux >> 63 != 0) {
+            // underflow
+            if (x <= -1075 or x - 0x1.0p52 + 0x1.0p52 != x) {
+                math.forceEval(f32(-0x1.0p-149 / x));
+            }
+            if (x <= -1075) {
+                return 0;
+            }
+        }
+    }
+    // |x| < 0x1p-54
+    else if (ix < 0x3C900000) {
+        return 1.0 + x;
+    }
+
+    // reduce x
+    var uf = x + redux;
+    // NOTE: musl performs an implicit 64-bit to 32-bit u32 truncation here
+    var i0 = @truncate(u32, @bitCast(u64, uf));
+    i0 += tblsiz / 2;
+
+    const k: u32 = i0 / tblsiz * tblsiz;
+    const ik = @bitCast(i32, k / tblsiz);
+    i0 %= tblsiz;
+    uf -= redux;
+
+    // r = exp2(y) = exp2t[i0] * p(z - eps[i])
+    var z = x - uf;
+    const t = exp2dt[2 * i0];
+    z -= exp2dt[2 * i0 + 1];
+    const r = t + t * z * (P1 + z * (P2 + z * (P3 + z * (P4 + z * P5))));
+
+    math.scalbn(r, ik)
+}
+
+test "math.exp2" {
+    assert(exp2(f32(0.8923)) == exp2_32(0.8923));
+    assert(exp2(f64(0.8923)) == exp2_64(0.8923));
+}
+
+test "math.exp2_32" {
+    const epsilon = 0.000001;
+
+    assert(exp2_32(0.0) == 1.0);
+    assert(math.approxEq(f32, exp2_32(0.2), 1.148698, epsilon));
+    assert(math.approxEq(f32, exp2_32(0.8923), 1.856133, epsilon));
+    assert(math.approxEq(f32, exp2_32(1.5), 2.828427, epsilon));
+    assert(math.approxEq(f32, exp2_32(37.45), 187747237888, epsilon));
+}
+
+test "math.exp2_64" {
+    const epsilon = 0.000001;
+
+    assert(exp2_64(0.0) == 1.0);
+    assert(math.approxEq(f64, exp2_64(0.2), 1.148698, epsilon));
+    assert(math.approxEq(f64, exp2_64(0.8923), 1.856133, epsilon));
+    assert(math.approxEq(f64, exp2_64(1.5), 2.828427, epsilon));
+}
+
+test "math.exp2_32.special" {
+    assert(math.isPositiveInf(exp2_32(math.inf(f32))));
+    assert(math.isNan(exp2_32(math.nan(f32))));
+}
+
+test "math.exp2_64.special" {
+    assert(math.isPositiveInf(exp2_64(math.inf(f64))));
+    assert(math.isNan(exp2_64(math.nan(f64))));
+}

std/math/expm1.zig

@@ -0,0 +1,316 @@
+// Special Cases:
+//
+// - expm1(+inf) = +inf
+// - expm1(-inf) = -1
+// - expm1(nan)  = nan
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const expm1 = expm1_workaround;
+
+pub fn expm1_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(expm1_32, x),
+        f64 => @inlineCall(expm1_64, x),
+        else => @compileError("exp1m not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn expm1_32(x_: f32) -> f32 {
+    const o_threshold: f32 = 8.8721679688e+01;
+    const ln2_hi: f32      = 6.9313812256e-01;
+    const ln2_lo: f32      = 9.0580006145e-06;
+    const invln2: f32      = 1.4426950216e+00;
+    const Q1: f32 = -3.3333212137e-2;
+    const Q2: f32 =  1.5807170421e-3;
+
+    var x = x_;
+    const ux = @bitCast(u32, x);
+    const hx = ux & 0x7FFFFFFF;
+    const sign = hx >> 31;
+
+    // TODO: Shouldn't need this check explicitly.
+    if (math.isNegativeInf(x)) {
+        return -1.0;
+    }
+
+    // |x| >= 27 * ln2
+    if (hx >= 0x4195B844) {
+        // nan
+        if (hx > 0x7F800000) {
+            return x;
+        }
+        if (sign != 0) {
+            return -1;
+        }
+        if (x > o_threshold) {
+            x *= 0x1.0p127;
+            return x;
+        }
+    }
+
+    var hi: f32 = undefined;
+    var lo: f32 = undefined;
+    var c: f32 = undefined;
+    var k: i32 = undefined;
+
+    // |x| > 0.5 * ln2
+    if (hx > 0x3EB17218) {
+        // |x| < 1.5 * ln2
+        if (hx < 0x3F851592) {
+            if (sign == 0) {
+                hi = x - ln2_hi;
+                lo = ln2_lo;
+                k = 1;
+            } else {
+                hi = x + ln2_hi;
+                lo = -ln2_lo;
+                k = -1;
+            }
+        } else {
+            var kf = invln2 * x;
+            if (sign != 0) {
+                kf -= 0.5;
+            } else {
+                kf += 0.5;
+            }
+
+            k = i32(kf);
+            const t = f32(k);
+            hi = x - t * ln2_hi;
+            lo = t * ln2_lo;
+        }
+
+        x = hi - lo;
+        c = (hi - x) - lo;
+    }
+    // |x| < 2^(-25)
+    else if (hx < 0x33000000) {
+        if (hx < 0x00800000) {
+            math.forceEval(x * x);
+        }
+        return x;
+    }
+    else {
+        k = 0;
+    }
+
+    const hfx = 0.5 * x;
+    const hxs = x * hfx;
+    const r1 = 1.0 + hxs * (Q1 + hxs * Q2);
+    const t = 3.0 - r1 * hfx;
+    var e = hxs * ((r1 - t) / (6.0 - x * t));
+
+    // c is 0
+    if (k == 0) {
+        return x - (x * e - hxs);
+    }
+
+    e = x * (e - c) - c;
+    e -= hxs;
+
+    // exp(x) ~ 2^k (x_reduced - e + 1)
+    if (k == -1) {
+        return 0.5 * (x - e) - 0.5;
+    }
+    if (k == 1) {
+        if (x < -0.25) {
+            return -2.0 * (e - (x + 0.5));
+        } else {
+            return 1.0 + 2.0 * (x - e);
+        }
+    }
+
+    const twopk = @bitCast(f32, u32((0x7F + k) <<% 23));
+
+    if (k < 0 or k > 56) {
+        var y = x - e + 1.0;
+        if (k == 128) {
+            y = y * 2.0 * 0x1.0p127;
+        } else {
+            y = y * twopk;
+        }
+
+        return y - 1.0;
+    }
+
+    const uf = @bitCast(f32, u32(0x7F - k) << 23);
+    if (k < 23) {
+        return (x - e + (1 - uf)) * twopk;
+    } else {
+        return (x - (e + uf) + 1) * twopk;
+    }
+}
+
+fn expm1_64(x_: f64) -> f64 {
+    const o_threshold: f64 = 7.09782712893383973096e+02;
+    const ln2_hi: f64      = 6.93147180369123816490e-01;
+    const ln2_lo: f64      = 1.90821492927058770002e-10;
+    const invln2: f64      = 1.44269504088896338700e+00;
+    const Q1: f64 = -3.33333333333331316428e-02;
+    const Q2: f64 =  1.58730158725481460165e-03;
+    const Q3: f64 = -7.93650757867487942473e-05;
+    const Q4: f64 =  4.00821782732936239552e-06;
+    const Q5: f64 = -2.01099218183624371326e-07;
+
+    var x = x_;
+    const ux = @bitCast(u64, x);
+    const hx = u32(ux >> 32) & 0x7FFFFFFF;
+    const sign = hx >> 63;
+
+    if (math.isNegativeInf(x)) {
+        return -1.0;
+    }
+
+    // |x| >= 56 * ln2
+    if (hx >= 0x4043687A) {
+        // exp1md(nan) = nan
+        if (hx > 0x7FF00000) {
+            return x;
+        }
+        // exp1md(-ve) = -1
+        if (sign != 0) {
+            return -1;
+        }
+        if (x > o_threshold) {
+            math.raiseOverflow();
+            return math.inf(f64);
+        }
+    }
+
+    var hi: f64 = undefined;
+    var lo: f64 = undefined;
+    var c: f64 = undefined;
+    var k: i32 = undefined;
+
+    // |x| > 0.5 * ln2
+    if (hx > 0x3FD62E42) {
+        // |x| < 1.5 * ln2
+        if (hx < 0x3FF0A2B2) {
+            if (sign == 0) {
+                hi = x - ln2_hi;
+                lo = ln2_lo;
+                k = 1;
+            } else {
+                hi = x + ln2_hi;
+                lo = -ln2_lo;
+                k = -1;
+            }
+        } else {
+            var kf = invln2 * x;
+            if (sign != 0) {
+                kf -= 0.5;
+            } else {
+                kf += 0.5;
+            }
+
+            k = i32(kf);
+            const t = f64(k);
+            hi = x - t * ln2_hi;
+            lo = t * ln2_lo;
+        }
+
+        x = hi - lo;
+        c = (hi - x) - lo;
+    }
+    // |x| < 2^(-54)
+    else if (hx < 0x3C900000) {
+        if (hx < 0x00100000) {
+            math.forceEval(f32(x));
+        }
+        return x;
+    }
+    else {
+        k = 0;
+    }
+
+    const hfx = 0.5 * x;
+    const hxs = x * hfx;
+    const r1 = 1.0 + hxs * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5))));
+    const t = 3.0 - r1 * hfx;
+    var e = hxs * ((r1 - t) / (6.0 - x * t));
+
+    // c is 0
+    if (k == 0) {
+        return x - (x * e - hxs);
+    }
+
+    e = x * (e - c) - c;
+    e -= hxs;
+
+    // exp(x) ~ 2^k (x_reduced - e + 1)
+    if (k == -1) {
+        return 0.5 * (x - e) - 0.5;
+    }
+    if (k == 1) {
+        if (x < -0.25) {
+            return -2.0 * (e - (x + 0.5));
+        } else {
+            return 1.0 + 2.0 * (x - e);
+        }
+    }
+
+    const twopk = @bitCast(f64, u64(0x3FF + k) <<% 52);
+
+    if (k < 0 or k > 56) {
+        var y = x - e + 1.0;
+        if (k == 1024) {
+            y = y * 2.0 * 0x1.0p1022 * 10;
+        } else {
+            y = y * twopk;
+        }
+
+        return y - 1.0;
+    }
+
+    const uf = @bitCast(f64, u64(0x3FF - k) << 52);
+    if (k < 20) {
+        return (x - e + (1 - uf)) * twopk;
+    } else {
+        return (x - (e + uf) + 1) * twopk;
+    }
+}
+
+test "math.exp1m" {
+    assert(expm1(f32(0.0)) == expm1_32(0.0));
+    assert(expm1(f64(0.0)) == expm1_64(0.0));
+}
+
+test "math.expm1_32" {
+    const epsilon = 0.000001;
+
+    assert(expm1_32(0.0) == 0.0);
+    assert(math.approxEq(f32, expm1_32(0.0), 0.0, epsilon));
+    assert(math.approxEq(f32, expm1_32(0.2), 0.221403, epsilon));
+    assert(math.approxEq(f32, expm1_32(0.8923), 1.440737, epsilon));
+    assert(math.approxEq(f32, expm1_32(1.5), 3.481689, epsilon));
+}
+
+test "math.expm1_64" {
+    const epsilon = 0.000001;
+
+    assert(expm1_64(0.0) == 0.0);
+    assert(math.approxEq(f64, expm1_64(0.0), 0.0, epsilon));
+    assert(math.approxEq(f64, expm1_64(0.2), 0.221403, epsilon));
+    assert(math.approxEq(f64, expm1_64(0.8923), 1.440737, epsilon));
+    assert(math.approxEq(f64, expm1_64(1.5), 3.481689, epsilon));
+}
+
+test "math.expm1_32.special" {
+    const epsilon = 0.000001;
+
+    assert(math.isPositiveInf(expm1_32(math.inf(f32))));
+    assert(expm1_32(-math.inf(f32)) == -1.0);
+    assert(math.isNan(expm1_32(math.nan(f32))));
+}
+
+test "math.expm1_64.special" {
+    const epsilon = 0.000001;
+
+    assert(math.isPositiveInf(expm1_64(math.inf(f64))));
+    assert(expm1_64(-math.inf(f64)) == -1.0);
+    assert(math.isNan(expm1_64(math.nan(f64))));
+}

std/math/expo2.zig

@@ -0,0 +1,28 @@
+const math = @import("index.zig");
+
+pub fn expo2(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => expo2f(x),
+        f64 => expo2d(x),
+        else => @compileError("expo2 not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn expo2f(x: f32) -> f32 {
+    const k: u32 = 235;
+    const kln2 = 0x1.45C778p+7;
+
+    const u = (0x7F + k / 2) << 23;
+    const scale = @bitCast(f32, u);
+    math.exp(x - kln2) * scale * scale
+}
+
+fn expo2d(x: f64) -> f64 {
+    const k: u32 = 2043;
+    const kln2 = 0x1.62066151ADD8BP+10;
+
+    const u = (0x3FF + k / 2) << 20;
+    const scale = @bitCast(f64, u64(u) << 32);
+    math.exp(x - kln2) * scale * scale
+}

std/math/fabs.zig

@@ -1,10 +1,19 @@
+// Special Cases:
+//
+// - fabs(+-inf) = +inf
+// - fabs(nan)   = nan
+
+const math = @import("index.zig");
 const assert = @import("../debug.zig").assert;
 
-pub fn fabs(x: var) -> @typeOf(x) {
+// TODO issue #393
+pub const fabs = fabs_workaround;
+
+pub fn fabs_workaround(x: var) -> @typeOf(x) {
     const T = @typeOf(x);
     switch (T) {
-        f32 => fabs32(x),
+        f32 => @inlineCall(fabs32, x),
-        f64 => fabs64(x),
+        f64 => @inlineCall(fabs64, x),
         else => @compileError("fabs not implemented for " ++ @typeName(T)),
     }
 }
@@ -21,29 +30,29 @@ fn fabs64(x: f64) -> f64 {
     @bitCast(f64, u)
 }
 
-test "fabs" {
+test "math.fabs" {
     assert(fabs(f32(1.0)) == fabs32(1.0));
     assert(fabs(f64(1.0)) == fabs64(1.0));
-    comptime {
-        assert(fabs(f32(1.0)) == fabs32(1.0));
-        assert(fabs(f64(1.0)) == fabs64(1.0));
-    }
 }
 
-test "fabs32" {
+test "math.fabs32" {
     assert(fabs64(1.0) == 1.0);
     assert(fabs64(-1.0) == 1.0);
-    comptime {
-        assert(fabs64(1.0) == 1.0);
-        assert(fabs64(-1.0) == 1.0);
-    }
 }
 
-test "fabs64" {
+test "math.fabs64" {
     assert(fabs64(1.0) == 1.0);
     assert(fabs64(-1.0) == 1.0);
-    comptime {
+}
-        assert(fabs64(1.0) == 1.0);
+
-        assert(fabs64(-1.0) == 1.0);
+test "math.fabs32.special" {
-    }
+    assert(math.isPositiveInf(fabs(math.inf(f32))));
+    assert(math.isPositiveInf(fabs(-math.inf(f32))));
+    assert(math.isNan(fabs(math.nan(f32))));
+}
+
+test "math.fabs64.special" {
+    assert(math.isPositiveInf(fabs(math.inf(f64))));
+    assert(math.isPositiveInf(fabs(-math.inf(f64))));
+    assert(math.isNan(fabs(math.nan(f64))));
 }

std/math/floor.zig

@@ -0,0 +1,119 @@
+// Special Cases:
+//
+// - floor(+-0)   = +-0
+// - floor(+-inf) = +-inf
+// - floor(nan)   = nan
+
+const builtin = @import("builtin");
+const assert = @import("../debug.zig").assert;
+const math = @import("index.zig");
+
+// TODO issue #393
+pub const floor = floor_workaround;
+
+pub fn floor_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(floor32, x),
+        f64 => @inlineCall(floor64, x),
+        else => @compileError("floor not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn floor32(x: f32) -> f32 {
+    var u = @bitCast(u32, x);
+    const e = i32((u >> 23) & 0xFF) - 0x7F;
+    var m: u32 = undefined;
+
+    // TODO: Shouldn't need this explicit check.
+    if (x == 0.0) {
+        return x;
+    }
+
+    if (e >= 23) {
+        return x;
+    }
+
+    if (e >= 0) {
+        m = 0x007FFFFF >> u32(e);
+        if (u & m == 0) {
+            return x;
+        }
+        math.forceEval(x + 0x1.0p120);
+        if (u >> 31 != 0) {
+            u += m;
+        }
+        @bitCast(f32, u & ~m)
+    } else {
+        math.forceEval(x + 0x1.0p120);
+        if (u >> 31 == 0) {
+            return 0.0; // Compiler requires return
+        } else {
+            -1.0
+        }
+    }
+}
+
+fn floor64(x: f64) -> f64 {
+    const u = @bitCast(u64, x);
+    const e = (u >> 52) & 0x7FF;
+    var y: f64 = undefined;
+
+    if (e >= 0x3FF+52 or x == 0) {
+        return x;
+    }
+
+    if (u >> 63 != 0) {
+        @setFloatMode(this, builtin.FloatMode.Strict);
+        y = x - math.f64_toint + math.f64_toint - x;
+    } else {
+        @setFloatMode(this, builtin.FloatMode.Strict);
+        y = x + math.f64_toint - math.f64_toint - x;
+    }
+
+    if (e <= 0x3FF-1) {
+        math.forceEval(y);
+        if (u >> 63 != 0) {
+            return -1.0;    // Compiler requires return.
+        } else {
+            0.0
+        }
+    } else if (y > 0) {
+        x + y - 1
+    } else {
+        x + y
+    }
+}
+
+test "math.floor" {
+    assert(floor(f32(1.3)) == floor32(1.3));
+    assert(floor(f64(1.3)) == floor64(1.3));
+}
+
+test "math.floor32" {
+    assert(floor32(1.3) == 1.0);
+    assert(floor32(-1.3) == -2.0);
+    assert(floor32(0.2) == 0.0);
+}
+
+test "math.floor64" {
+    assert(floor64(1.3) == 1.0);
+    assert(floor64(-1.3) == -2.0);
+    assert(floor64(0.2) == 0.0);
+}
+
+test "math.floor32.special" {
+    assert(floor32(0.0) == 0.0);
+    assert(floor32(-0.0) == -0.0);
+    assert(math.isPositiveInf(floor32(math.inf(f32))));
+    assert(math.isNegativeInf(floor32(-math.inf(f32))));
+    assert(math.isNan(floor32(math.nan(f32))));
+}
+
+test "math.floor64.special" {
+    assert(floor64(0.0) == 0.0);
+    assert(floor64(-0.0) == -0.0);
+    assert(math.isPositiveInf(floor64(math.inf(f64))));
+    assert(math.isNegativeInf(floor64(-math.inf(f64))));
+    assert(math.isNan(floor64(math.nan(f64))));
+}

std/math/fma.zig

@@ -0,0 +1,163 @@
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const fma = fma_workaround;
+
+pub fn fma_workaround(comptime T: type, x: T, y: T, z: T) -> T {
+    switch (T) {
+        f32 => @inlineCall(fma32, x, y, z),
+        f64 => @inlineCall(fma64, x, y ,z),
+        else => @compileError("fma not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn fma32(x: f32, y: f32, z: f32) -> f32 {
+    const xy = f64(x) * y;
+    const xy_z = xy + z;
+    const u = @bitCast(u64, xy_z);
+    const e = (u >> 52) & 0x7FF;
+
+    if ((u & 0x1FFFFFFF) != 0x10000000 or e == 0x7FF or xy_z - xy == z) {
+        f32(xy_z)
+    } else {
+        // TODO: Handle inexact case with double-rounding
+        f32(xy_z)
+    }
+}
+
+fn fma64(x: f64, y: f64, z: f64) -> f64 {
+    if (!math.isFinite(x) or !math.isFinite(y)) {
+        return x * y + z;
+    }
+    if (!math.isFinite(z)) {
+        return z;
+    }
+    if (x == 0.0 or y == 0.0) {
+        return x * y + z;
+    }
+    if (z == 0.0) {
+        return x * y;
+    }
+
+    const x1 = math.frexp(x);
+    var ex = x1.exponent;
+    var xs = x1.significand;
+    const x2 = math.frexp(y);
+    var ey = x2.exponent;
+    var ys = x2.significand;
+    const x3 = math.frexp(z);
+    var ez = x3.exponent;
+    var zs = x3.significand;
+
+    var spread = ex + ey - ez;
+    if (spread <= 53 * 2) {
+        zs = math.scalbn(zs, -spread);
+    } else {
+        zs = math.copysign(f64, math.f64_min, zs);
+    }
+
+    const xy = dd_mul(xs, ys);
+    const r = dd_add(xy.hi, zs);
+    spread = ex + ey;
+
+    if (r.hi == 0.0) {
+        return xy.hi + zs + math.scalbn(xy.lo, spread);
+    }
+
+    const adj = add_adjusted(r.lo, xy.lo);
+    if (spread + math.ilogb(r.hi) > -1023) {
+        math.scalbn(r.hi + adj, spread)
+    } else {
+        add_and_denorm(r.hi, adj, spread)
+    }
+}
+
+const dd = struct { hi: f64, lo: f64, };
+
+fn dd_add(a: f64, b: f64) -> dd {
+    var ret: dd = undefined;
+    ret.hi = a + b;
+    const s = ret.hi - a;
+    ret.lo = (a - (ret.hi - s)) + (b - s);
+    ret
+}
+
+fn dd_mul(a: f64, b: f64) -> dd {
+    var ret: dd = undefined;
+    const split: f64 = 0x1.0p27 + 1.0;
+
+    var p = a * split;
+    var ha = a - p;
+    ha += p;
+    var la = a - ha;
+
+    p = b * split;
+    var hb = b - p;
+    hb += p;
+    var lb = b - hb;
+
+    p = ha * hb;
+    var q = ha * lb + la * hb;
+
+    ret.hi = p + q;
+    ret.lo = p - ret.hi + q + la * lb;
+    ret
+}
+
+fn add_adjusted(a: f64, b: f64) -> f64 {
+    var sum = dd_add(a, b);
+    if (sum.lo != 0) {
+        var uhii = @bitCast(u64, sum.hi);
+        if (uhii & 1 == 0) {
+            // hibits += copysign(1.0, sum.hi, sum.lo)
+            const uloi = @bitCast(u64, sum.lo);
+            uhii += 1 - ((uhii ^ uloi) >> 62);
+            sum.hi = @bitCast(f64, uhii);
+        }
+    }
+    sum.hi
+}
+
+fn add_and_denorm(a: f64, b: f64, scale: i32) -> f64 {
+    var sum = dd_add(a, b);
+    if (sum.lo != 0) {
+        var uhii = @bitCast(u64, sum.hi);
+        const bits_lost = -i32((uhii >> 52) & 0x7FF) - scale + 1;
+        if ((bits_lost != 1) == (uhii & 1 != 0)) {
+            const uloi = @bitCast(u64, sum.lo);
+            uhii += 1 - (((uhii ^ uloi) >> 62) & 2);
+            sum.hi = @bitCast(f64, uhii);
+        }
+    }
+    math.scalbn(sum.hi, scale)
+}
+
+test "math.fma" {
+    assert(fma(f32, 0.0, 1.0, 1.0) == fma32(0.0, 1.0, 1.0));
+    assert(fma(f64, 0.0, 1.0, 1.0) == fma64(0.0, 1.0, 1.0));
+}
+
+test "math.fma32" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, fma32(0.0, 5.0, 9.124), 9.124, epsilon));
+    assert(math.approxEq(f32, fma32(0.2, 5.0, 9.124), 10.124, epsilon));
+    assert(math.approxEq(f32, fma32(0.8923, 5.0, 9.124), 13.5855, epsilon));
+    assert(math.approxEq(f32, fma32(1.5, 5.0, 9.124), 16.624, epsilon));
+    assert(math.approxEq(f32, fma32(37.45, 5.0, 9.124), 196.374004, epsilon));
+    assert(math.approxEq(f32, fma32(89.123, 5.0, 9.124), 454.739005, epsilon));
+    assert(math.approxEq(f32, fma32(123123.234375, 5.0, 9.124), 615625.295875, epsilon));
+}
+
+test "math.fma64" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f64, fma64(0.0, 5.0, 9.124), 9.124, epsilon));
+    assert(math.approxEq(f64, fma64(0.2, 5.0, 9.124), 10.124, epsilon));
+    assert(math.approxEq(f64, fma64(0.8923, 5.0, 9.124), 13.5855, epsilon));
+    assert(math.approxEq(f64, fma64(1.5, 5.0, 9.124), 16.624, epsilon));
+    assert(math.approxEq(f64, fma64(37.45, 5.0, 9.124), 196.374, epsilon));
+    assert(math.approxEq(f64, fma64(89.123, 5.0, 9.124), 454.739, epsilon));
+    assert(math.approxEq(f64, fma64(123123.234375, 5.0, 9.124), 615625.295875, epsilon));
+}

std/math/frexp.zig

@@ -1,89 +1,173 @@
-const assert = @import("../debug.zig").assert;
+// Special Cases:
-const math = @import("index.zig");
+//
+// - frexp(+-0)   = +-0, 0
+// - frexp(+-inf) = +-inf, 0
+// - frexp(nan)   = nan, undefined
 
-pub fn frexp(x: var, e: &i32) -> @typeOf(x) {
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const frexp = frexp_workaround;
+
+fn frexp_result(comptime T: type) -> type {
+    struct {
+        significand: T,
+        exponent: i32,
+    }
+}
+pub const frexp32_result = frexp_result(f32);
+pub const frexp64_result = frexp_result(f64);
+
+pub fn frexp_workaround(x: var) -> frexp_result(@typeOf(x)) {
     const T = @typeOf(x);
     switch (T) {
-        f32 => frexp32(x, e),
+        f32 => @inlineCall(frexp32, x),
-        f64 => frexp64(x, e),
+        f64 => @inlineCall(frexp64, x),
         else => @compileError("frexp not implemented for " ++ @typeName(T)),
     }
 }
 
-fn frexp32(x_: f32, e: &i32) -> f32 {
+fn frexp32(x: f32) -> frexp32_result {
-    var x = x_;
+    var result: frexp32_result = undefined;
-    var y = @bitCast(u32, x);
-    const ee = i32(y >> 23) & 0xFF;
 
-    if (ee == 0) {
+    var y = @bitCast(u32, x);
+    const e = i32(y >> 23) & 0xFF;
+
+    if (e == 0) {
         if (x != 0) {
-            x = frexp32(x * 0x1.0p64, e);
+            // subnormal
-            *e -= 64;
+            result = frexp32(x * 0x1.0p64);
+            result.exponent -= 64;
         } else {
-            *e = 0;
+            // frexp(+-0) = (+-0, 0)
+            result.significand = x;
+            result.exponent = 0;
         }
-        return x;
+        return result;
-    } else if (ee == 0xFF) {
+    } else if (e == 0xFF) {
-        return x;
+        // frexp(nan) = (nan, undefined)
+        result.significand = x;
+        result.exponent = undefined;
+
+        // frexp(+-inf) = (+-inf, 0)
+        if (math.isInf(x)) {
+            result.exponent = 0;
+        }
+
+        return result;
     }
 
-    *e = ee - 0x7E;
+    result.exponent = e - 0x7E;
     y &= 0x807FFFFF;
     y |= 0x3F000000;
-    @bitCast(f32, y)
+    result.significand = @bitCast(f32, y);
+    result
 }
 
-fn frexp64(x_: f64, e: &i32) -> f64 {
+fn frexp64(x: f64) -> frexp64_result {
-    var x = x_;
+    var result: frexp64_result = undefined;
-    var y = @bitCast(u64, x);
-    const ee = i32(y >> 52) & 0x7FF;
 
-    if (ee == 0) {
+    var y = @bitCast(u64, x);
+    const e = i32(y >> 52) & 0x7FF;
+
+    if (e == 0) {
         if (x != 0) {
-            x = frexp64(x * 0x1.0p64, e);
+            // subnormal
-            *e -= 64;
+            result = frexp64(x * 0x1.0p64);
+            result.exponent -= 64;
         } else {
-            *e = 0;
+            // frexp(+-0) = (+-0, 0)
+            result.significand = x;
+            result.exponent = 0;
         }
-        return x;
+        return result;
-    } else if (ee == 0x7FF) {
+    } else if (e == 0x7FF) {
-        return x;
+        // frexp(nan) = (nan, undefined)
+        result.significand = x;
+        result.exponent = undefined;
+
+        // frexp(+-inf) = (+-inf, 0)
+        if (math.isInf(x)) {
+            result.exponent = 0;
+        }
+
+        return result;
     }
 
-    *e = ee - 0x3FE;
+    result.exponent = e - 0x3FE;
     y &= 0x800FFFFFFFFFFFFF;
     y |= 0x3FE0000000000000;
-    @bitCast(f64, y)
+    result.significand = @bitCast(f64, y);
+    result
 }
 
-test "frexp" {
+test "math.frexp" {
-    var i0: i32 = undefined;
+    const a = frexp(f32(1.3));
-    var i1: i32 = undefined;
+    const b = frexp32(1.3);
+    assert(a.significand == b.significand and a.exponent == b.exponent);
 
-    assert(frexp(f32(1.3), &i0) == frexp32(1.3, &i1));
+    const c = frexp(f64(1.3));
-    assert(frexp(f64(1.3), &i0) == frexp64(1.3, &i1));
+    const d = frexp64(1.3);
+    assert(c.significand == d.significand and c.exponent == d.exponent);
 }
 
-test "frexp32" {
+test "math.frexp32" {
     const epsilon = 0.000001;
-    var i: i32 = undefined;
+    var r: frexp32_result = undefined;
-    var d: f32 = undefined;
 
-    d = frexp32(1.3, &i);
+    r = frexp32(1.3);
-    assert(math.approxEq(f32, d, 0.65, epsilon) and i == 1);
+    assert(math.approxEq(f32, r.significand, 0.65, epsilon) and r.exponent == 1);
 
-    d = frexp32(78.0234, &i);
+    r = frexp32(78.0234);
-    assert(math.approxEq(f32, d, 0.609558, epsilon) and i == 7);
+    assert(math.approxEq(f32, r.significand, 0.609558, epsilon) and r.exponent == 7);
 }
 
-test "frexp64" {
+test "math.frexp64" {
     const epsilon = 0.000001;
-    var i: i32 = undefined;
+    var r: frexp64_result = undefined;
-    var d: f64 = undefined;
 
-    d = frexp64(1.3, &i);
+    r = frexp64(1.3);
-    assert(math.approxEq(f64, d, 0.65, epsilon) and i == 1);
+    assert(math.approxEq(f64, r.significand, 0.65, epsilon) and r.exponent == 1);
 
-    d = frexp64(78.0234, &i);
+    r = frexp64(78.0234);
-    assert(math.approxEq(f64, d, 0.609558, epsilon) and i == 7);
+    assert(math.approxEq(f64, r.significand, 0.609558, epsilon) and r.exponent == 7);
+}
+
+test "math.frexp32.special" {
+    var r: frexp32_result = undefined;
+
+    r = frexp32(0.0);
+    assert(r.significand == 0.0 and r.exponent == 0);
+
+    r = frexp32(-0.0);
+    assert(r.significand == -0.0 and r.exponent == 0);
+
+    r = frexp32(math.inf(f32));
+    assert(math.isPositiveInf(r.significand) and r.exponent == 0);
+
+    r = frexp32(-math.inf(f32));
+    assert(math.isNegativeInf(r.significand) and r.exponent == 0);
+
+    r = frexp32(math.nan(f32));
+    assert(math.isNan(r.significand));
+}
+
+test "math.frexp64.special" {
+    var r: frexp64_result = undefined;
+
+    r = frexp64(0.0);
+    assert(r.significand == 0.0 and r.exponent == 0);
+
+    r = frexp64(-0.0);
+    assert(r.significand == -0.0 and r.exponent == 0);
+
+    r = frexp64(math.inf(f64));
+    assert(math.isPositiveInf(r.significand) and r.exponent == 0);
+
+    r = frexp64(-math.inf(f64));
+    assert(math.isNegativeInf(r.significand) and r.exponent == 0);
+
+    r = frexp64(math.nan(f64));
+    assert(math.isNan(r.significand));
 }

std/math/hypot.zig

@@ -0,0 +1,163 @@
+// Special Cases:
+//
+// - hypot(+-inf, y)  = +inf
+// - hypot(x, +-inf)  = +inf
+// - hypot(nan, y)    = nan
+// - hypot(x, nan)    = nan
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const hypot = hypot_workaround;
+
+pub fn hypot_workaround(comptime T: type, x: T, y: T) -> T {
+    switch (T) {
+        f32 => @inlineCall(hypot32, x, y),
+        f64 => @inlineCall(hypot64, x, y),
+        else => @compileError("hypot not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn hypot32(x: f32, y: f32) -> f32 {
+    var ux = @bitCast(u32, x);
+    var uy = @bitCast(u32, y);
+
+    ux &= @maxValue(u32) >> 1;
+    uy &= @maxValue(u32) >> 1;
+    if (ux < uy) {
+        const tmp = ux;
+        ux = uy;
+        uy = tmp;
+    }
+
+    var xx = @bitCast(f32, ux);
+    var yy = @bitCast(f32, uy);
+    if (uy == 0xFF << 23) {
+        return yy;
+    }
+    if (ux >= 0xFF << 23 or uy == 0 or ux - uy >= (25 << 23)) {
+        return xx + yy;
+    }
+
+    var z: f32 = 1.0;
+    if (ux >= (0x7F+60) << 23) {
+        z = 0x1.0p90;
+        xx *= 0x1.0p-90;
+        yy *= 0x1.0p-90;
+    } else if (uy < (0x7F-60) << 23) {
+        z = 0x1.0p-90;
+        xx *= 0x1.0p-90;
+        yy *= 0x1.0p-90;
+    }
+
+    z * math.sqrt(f32(f64(x) * x + f64(y) * y))
+}
+
+fn sq(hi: &f64, lo: &f64, x: f64) {
+    const split: f64 = 0x1.0p27 + 1.0;
+    const xc = x * split;
+    const xh = x - xc + xc;
+    const xl = x - xh;
+    *hi = x * x;
+    *lo = xh * xh - *hi + 2 * xh * xl + xl * xl;
+}
+
+fn hypot64(x: f64, y: f64) -> f64 {
+    var ux = @bitCast(u64, x);
+    var uy = @bitCast(u64, y);
+
+    ux &= @maxValue(u64) >> 1;
+    uy &= @maxValue(u64) >> 1;
+    if (ux < uy) {
+        const tmp = ux;
+        ux = uy;
+        uy = tmp;
+    }
+
+    const ex = ux >> 52;
+    const ey = uy >> 52;
+    var xx = @bitCast(f64, ux);
+    var yy = @bitCast(f64, uy);
+
+    // hypot(inf, nan) == inf
+    if (ey == 0x7FF) {
+        return yy;
+    }
+    if (ex == 0x7FF or uy == 0) {
+        return xx;
+    }
+
+    // hypot(x, y) ~= x + y * y / x / 2 with inexact for small y/x
+    if (ex - ey > 64) {
+        return xx + yy;
+    }
+
+    var z: f64 = 1;
+    if (ex > 0x3FF + 510) {
+        z = 0x1.0p700;
+        xx *= 0x1.0p-700;
+        yy *= 0x1.0p-700;
+    } else if (ey < 0x3FF - 450) {
+        z = 0x1.0p-700;
+        xx *= 0x1.0p700;
+        yy *= 0x1.0p700;
+    }
+
+    var hx: f64 = undefined;
+    var lx: f64 = undefined;
+    var hy: f64 = undefined;
+    var ly: f64 = undefined;
+
+    sq(&hx, &lx, x);
+    sq(&hy, &ly, y);
+
+    z * math.sqrt(ly + lx + hy + hx)
+}
+
+test "math.hypot" {
+    assert(hypot(f32, 0.0, -1.2) == hypot32(0.0, -1.2));
+    assert(hypot(f64, 0.0, -1.2) == hypot64(0.0, -1.2));
+}
+
+test "math.hypot32" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, hypot32(0.0, -1.2), 1.2, epsilon));
+    assert(math.approxEq(f32, hypot32(0.2, -0.34), 0.394462, epsilon));
+    assert(math.approxEq(f32, hypot32(0.8923, 2.636890), 2.783772, epsilon));
+    assert(math.approxEq(f32, hypot32(1.5, 5.25), 5.460083, epsilon));
+    assert(math.approxEq(f32, hypot32(37.45, 159.835), 164.163742, epsilon));
+    assert(math.approxEq(f32, hypot32(89.123, 382.028905), 392.286865, epsilon));
+    assert(math.approxEq(f32, hypot32(123123.234375, 529428.707813), 543556.875, epsilon));
+}
+
+test "math.hypot64" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f64, hypot64(0.0, -1.2), 1.2, epsilon));
+    assert(math.approxEq(f64, hypot64(0.2, -0.34), 0.394462, epsilon));
+    assert(math.approxEq(f64, hypot64(0.8923, 2.636890), 2.783772, epsilon));
+    assert(math.approxEq(f64, hypot64(1.5, 5.25), 5.460082, epsilon));
+    assert(math.approxEq(f64, hypot64(37.45, 159.835), 164.163728, epsilon));
+    assert(math.approxEq(f64, hypot64(89.123, 382.028905), 392.286876, epsilon));
+    assert(math.approxEq(f64, hypot64(123123.234375, 529428.707813), 543556.885247, epsilon));
+}
+
+test "math.hypot32.special" {
+    assert(math.isPositiveInf(hypot32(math.inf(f32), 0.0)));
+    assert(math.isPositiveInf(hypot32(-math.inf(f32), 0.0)));
+    assert(math.isPositiveInf(hypot32(0.0, math.inf(f32))));
+    assert(math.isPositiveInf(hypot32(0.0, -math.inf(f32))));
+    assert(math.isNan(hypot32(math.nan(f32), 0.0)));
+    assert(math.isNan(hypot32(0.0, math.nan(f32))));
+}
+
+test "math.hypot64.special" {
+    assert(math.isPositiveInf(hypot64(math.inf(f64), 0.0)));
+    assert(math.isPositiveInf(hypot64(-math.inf(f64), 0.0)));
+    assert(math.isPositiveInf(hypot64(0.0, math.inf(f64))));
+    assert(math.isPositiveInf(hypot64(0.0, -math.inf(f64))));
+    assert(math.isNan(hypot64(math.nan(f64), 0.0)));
+    assert(math.isNan(hypot64(0.0, math.nan(f64))));
+}

std/math/ilogb.zig

@@ -0,0 +1,132 @@
+// Special Cases:
+//
+// - ilogb(+-inf) = @maxValue(i32)
+// - ilogb(0)     = @maxValue(i32)
+// - ilogb(nan)   = @maxValue(i32)
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const ilogb = ilogb_workaround;
+
+pub fn ilogb_workaround(x: var) -> i32 {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(ilogb32, x),
+        f64 => @inlineCall(ilogb64, x),
+        else => @compileError("ilogb not implemented for " ++ @typeName(T)),
+    }
+}
+
+// NOTE: Should these be exposed publically?
+const fp_ilogbnan = -1 - i32(@maxValue(u32) >> 1);
+const fp_ilogb0 = fp_ilogbnan;
+
+fn ilogb32(x: f32) -> i32 {
+    var u = @bitCast(u32, x);
+    var e = i32((u >> 23) & 0xFF);
+
+    // TODO: We should be able to merge this with the lower check.
+    if (math.isNan(x)) {
+        return @maxValue(i32);
+    }
+
+    if (e == 0) {
+        u <<= 9;
+        if (u == 0) {
+            math.raiseInvalid();
+            return fp_ilogb0;
+        }
+
+        // subnormal
+        e = -0x7F;
+        while (u >> 31 == 0) : (u <<= 1) {
+            e -= 1;
+        }
+        return e;
+    }
+
+    if (e == 0xFF) {
+        math.raiseInvalid();
+        if (u <<% 9 != 0) {
+            return fp_ilogbnan;
+        } else {
+            return @maxValue(i32);
+        }
+    }
+
+    e - 0x7F
+}
+
+fn ilogb64(x: f64) -> i32 {
+    var u = @bitCast(u64, x);
+    var e = i32((u >> 52) & 0x7FF);
+
+    if (math.isNan(x)) {
+        return @maxValue(i32);
+    }
+
+    if (e == 0) {
+        u <<= 12;
+        if (u == 0) {
+            math.raiseInvalid();
+            return fp_ilogb0;
+        }
+
+        // subnormal
+        e = -0x3FF;
+        while (u >> 63 == 0) : (u <<= 1) {
+            e -= 1;
+        }
+        return e;
+    }
+
+    if (e == 0x7FF) {
+        math.raiseInvalid();
+        if (u <<% 12 != 0) {
+            return fp_ilogbnan;
+        } else {
+            return @maxValue(i32);
+        }
+    }
+
+    e - 0x3FF
+}
+
+test "math.ilogb" {
+    assert(ilogb(f32(0.2)) == ilogb32(0.2));
+    assert(ilogb(f64(0.2)) == ilogb64(0.2));
+}
+
+test "math.ilogb32" {
+    assert(ilogb32(0.0) == fp_ilogb0);
+    assert(ilogb32(0.5) == -1);
+    assert(ilogb32(0.8923) == -1);
+    assert(ilogb32(10.0) == 3);
+    assert(ilogb32(-123984) == 16);
+    assert(ilogb32(2398.23) == 11);
+}
+
+test "math.ilogb64" {
+    assert(ilogb64(0.0) == fp_ilogb0);
+    assert(ilogb64(0.5) == -1);
+    assert(ilogb64(0.8923) == -1);
+    assert(ilogb64(10.0) == 3);
+    assert(ilogb64(-123984) == 16);
+    assert(ilogb64(2398.23) == 11);
+}
+
+test "math.ilogb32.special" {
+    assert(ilogb32(math.inf(f32)) == @maxValue(i32));
+    assert(ilogb32(-math.inf(f32)) == @maxValue(i32));
+    assert(ilogb32(0.0) == @minValue(i32));
+    assert(ilogb32(math.nan(f32)) == @maxValue(i32));
+}
+
+test "math.ilogb64.special" {
+    assert(ilogb64(math.inf(f64)) == @maxValue(i32));
+    assert(ilogb64(-math.inf(f64)) == @maxValue(i32));
+    assert(ilogb64(0.0) == @minValue(i32));
+    assert(ilogb64(math.nan(f64)) == @maxValue(i32));
+}

std/math/index.zig

@@ -1,7 +1,189 @@
-const assert = @import("../debug.zig").assert;
 const builtin = @import("builtin");
+const TypeId = builtin.TypeId;
+const assert = @import("../debug.zig").assert;
 
+pub const e         = 2.7182818284590452354;  // e
+pub const log2_e    = 1.4426950408889634074;  // log_2(e)
+pub const log10_e   = 0.43429448190325182765; // log_10(e)
+pub const ln_2      = 0.69314718055994530942; // log_e(2)
+pub const ln_10     = 2.30258509299404568402; // log_e(10)
+pub const pi        = 3.14159265358979323846; // pi
+pub const pi_2      = 1.57079632679489661923; // pi/2
+pub const pi_4      = 0.78539816339744830962; // pi/4
+pub const r1_pi     = 0.31830988618379067154; // 1/pi
+pub const r2_pi     = 0.63661977236758134308; // 2/pi
+pub const r2_sqrtpi = 1.12837916709551257390; // 2/sqrt(pi)
+pub const sqrt2     = 1.41421356237309504880; // sqrt(2)
+pub const r1_sqrt2  = 0.70710678118654752440; // 1/sqrt(2)
+
+// float.h details
+pub const f64_true_min = 4.94065645841246544177e-324;
+pub const f64_min = 2.22507385850720138309e-308;
+pub const f64_max = 1.79769313486231570815e+308;
+pub const f64_epsilon = 2.22044604925031308085e-16;
+pub const f64_toint = 1.0 / f64_epsilon;
+
+pub const f32_true_min = 1.40129846432481707092e-45;
+pub const f32_min = 1.17549435082228750797e-38;
+pub const f32_max = 3.40282346638528859812e+38;
+pub const f32_epsilon = 1.1920928955078125e-07;
+pub const f32_toint = 1.0 / f32_epsilon;
+
+pub const nan_u32 = u32(0x7F800001);
+pub const nan_f32 = @bitCast(f32, nan_u32);
+
+pub const inf_u32 = u32(0x7F800000);
+pub const inf_f32 = @bitCast(f32, inf_u32);
+
+pub const nan_u64 = u64(0x7FF << 52) | 1;
+pub const nan_f64 = @bitCast(f64, nan_u64);
+
+pub const inf_u64 = u64(0x7FF << 52);
+pub const inf_f64 = @bitCast(f64, inf_u64);
+
+pub const nan = @import("nan.zig").nan;
+pub const snan = @import("nan.zig").snan;
+pub const inf = @import("inf.zig").inf;
+
+pub fn approxEq(comptime T: type, x: T, y: T, epsilon: T) -> bool {
+    assert(@typeId(T) == TypeId.Float);
+    fabs(x - y) < epsilon
+}
+
+// TODO: Hide the following in an internal module.
+pub fn forceEval(value: var) {
+    const T = @typeOf(value);
+    switch (T) {
+        f32 => {
+            var x: f32 = undefined;
+            const p = @ptrCast(&volatile f32, &x);
+            *p = x;
+        },
+        f64 => {
+            var x: f64 = undefined;
+            const p = @ptrCast(&volatile f64, &x);
+            *p = x;
+        },
+        else => {
+            @compileError("forceEval not implemented for " ++ @typeName(T));
+        },
+    }
+}
+
+pub fn raiseInvalid() {
+    // Raise INVALID fpu exception
+}
+
+pub fn raiseUnderflow() {
+    // Raise UNDERFLOW fpu exception
+}
+
+pub fn raiseOverflow() {
+    // Raise OVERFLOW fpu exception
+}
+
+pub fn raiseInexact() {
+    // Raise INEXACT fpu exception
+}
+
+pub fn raiseDivByZero() {
+    // Raise INEXACT fpu exception
+}
+
+pub const isNan = @import("isnan.zig").isNan;
+pub const isSignalNan = @import("isnan.zig").isSignalNan;
+pub const fabs = @import("fabs.zig").fabs;
+pub const ceil = @import("ceil.zig").ceil;
+pub const floor = @import("floor.zig").floor;
+pub const trunc = @import("floor.zig").trunc;
+pub const round = @import("round.zig").round;
 pub const frexp = @import("frexp.zig").frexp;
+pub const frexp32_result = @import("frexp.zig").frexp32_result;
+pub const frexp64_result = @import("frexp.zig").frexp64_result;
+pub const modf = @import("modf.zig").modf;
+pub const modf32_result = @import("modf.zig").modf32_result;
+pub const modf64_result = @import("modf.zig").modf64_result;
+pub const copysign = @import("copysign.zig").copysign;
+pub const isFinite = @import("isfinite.zig").isFinite;
+pub const isInf = @import("isinf.zig").isInf;
+pub const isPositiveInf = @import("isinf.zig").isPositiveInf;
+pub const isNegativeInf = @import("isinf.zig").isNegativeInf;
+pub const isNormal = @import("isnormal.zig").isNormal;
+pub const signbit = @import("signbit.zig").signbit;
+pub const scalbn = @import("scalbn.zig").scalbn;
+pub const pow = @import("pow.zig").pow;
+pub const sqrt = @import("sqrt.zig").sqrt;
+pub const cbrt = @import("cbrt.zig").cbrt;
+pub const acos = @import("acos.zig").acos;
+pub const asin = @import("asin.zig").asin;
+pub const atan = @import("atan.zig").atan;
+pub const atan2 = @import("atan2.zig").atan2;
+pub const hypot = @import("hypot.zig").hypot;
+pub const exp = @import("exp.zig").exp;
+pub const exp2 = @import("exp2.zig").exp2;
+pub const expm1 = @import("expm1.zig").expm1;
+pub const ilogb = @import("ilogb.zig").ilogb;
+pub const ln = @import("ln.zig").ln;
+pub const log = @import("log.zig").log;
+pub const log2 = @import("log2.zig").log2;
+pub const log10 = @import("log10.zig").log10;
+pub const log1p = @import("log1p.zig").log1p;
+pub const fma = @import("fma.zig").fma;
+pub const asinh = @import("asinh.zig").asinh;
+pub const acosh = @import("acosh.zig").acosh;
+pub const atanh = @import("atanh.zig").atanh;
+pub const sinh = @import("sinh.zig").sinh;
+pub const cosh = @import("cosh.zig").cosh;
+pub const tanh = @import("tanh.zig").tanh;
+pub const cos = @import("cos.zig").cos;
+pub const sin = @import("sin.zig").sin;
+pub const tan = @import("tan.zig").tan;
+
+test "math" {
+    _ = @import("nan.zig");
+    _ = @import("isnan.zig");
+    _ = @import("fabs.zig");
+    _ = @import("ceil.zig");
+    _ = @import("floor.zig");
+    _ = @import("trunc.zig");
+    _ = @import("round.zig");
+    _ = @import("frexp.zig");
+    _ = @import("modf.zig");
+    _ = @import("copysign.zig");
+    _ = @import("isfinite.zig");
+    _ = @import("isinf.zig");
+    _ = @import("isnormal.zig");
+    _ = @import("signbit.zig");
+    _ = @import("scalbn.zig");
+    _ = @import("pow.zig");
+    _ = @import("sqrt.zig");
+    _ = @import("cbrt.zig");
+    _ = @import("acos.zig");
+    _ = @import("asin.zig");
+    _ = @import("atan.zig");
+    _ = @import("atan2.zig");
+    _ = @import("hypot.zig");
+    _ = @import("exp.zig");
+    _ = @import("exp2.zig");
+    _ = @import("expm1.zig");
+    _ = @import("ilogb.zig");
+    _ = @import("ln.zig");
+    _ = @import("log.zig");
+    _ = @import("log2.zig");
+    _ = @import("log10.zig");
+    _ = @import("log1p.zig");
+    _ = @import("fma.zig");
+    _ = @import("asinh.zig");
+    _ = @import("acosh.zig");
+    _ = @import("atanh.zig");
+    _ = @import("sinh.zig");
+    _ = @import("cosh.zig");
+    _ = @import("tanh.zig");
+    _ = @import("sin.zig");
+    _ = @import("cos.zig");
+    _ = @import("tan.zig");
+}
+
 
 pub const Cmp = enum {
     Less,
@@ -66,25 +248,6 @@ fn testOverflow() {
 }
 
 
-pub fn log(comptime base: usize, value: var) -> @typeOf(value) {
-    const T = @typeOf(value);
-    switch (@typeId(T)) {
-        builtin.TypeId.Int => {
-            if (base == 2) {
-                return T.bit_count - 1 - @clz(value);
-            } else {
-                @compileError("TODO implement log for non base 2 integers");
-            }
-        },
-        builtin.TypeId.Float => {
-            @compileError("TODO implement log for floats");
-        },
-        else => {
-            @compileError("log expects integer or float, found '" ++ @typeName(T) ++ "'");
-        },
-    }
-}
-
 error Overflow;
 pub fn absInt(x: var) -> %@typeOf(x) {
     const T = @typeOf(x);
@@ -244,92 +407,6 @@ fn testRem() {
     if (rem(f32, 10, 0)) |_| unreachable else |err| assert(err == error.DivisionByZero);
 }
 
-fn isNan(comptime T: type, x: T) -> bool {
-    assert(@typeId(T) == builtin.TypeId.Float);
-    if (T == f32) {
-        const bits = @bitCast(u32, x);
-        return (bits & 0x7fffffff) > 0x7f800000;
-    } else if (T == f64) {
-        const bits = @bitCast(u64, x);
-        return (bits & (@maxValue(u64) >> 1)) > (u64(0x7ff) << 52);
-    } else if (T == c_longdouble) {
-        @compileError("TODO support isNan for c_longdouble");
-    } else {
-        unreachable;
-    }
-}
-
-pub fn floor(x: var) -> @typeOf(x) {
-    switch (@typeOf(x)) {
-        f32 => floor_f32(x),
-        f64 => floor_f64(x),
-        c_longdouble => @compileError("TODO support floor for c_longdouble"),
-        else => @compileError("Invalid type for floor: " ++ @typeName(@typeOf(x))),
-    }
-}
-
-fn floor_f32(x: f32) -> f32 {
-    var i = @bitCast(u32, x);
-    const e = i32((i >> 23) & 0xff) -% 0x7f;
-    if (e >= 23)
-        return x;
-    if (e >= 0) {
-        const m = @bitCast(u32, 0x007fffff >> e);
-        if ((i & m) == 0)
-            return x;
-        if (i >> 31 != 0)
-            i +%= m;
-        i &= ~m;
-    } else {
-        if (i >> 31 == 0)
-            return 0;
-        if (i <<% 1 != 0)
-            return -1.0;
-    }
-    return @bitCast(f32, i);
-}
-
-fn floor_f64(x: f64) -> f64 {
-    const DBL_EPSILON = 2.22044604925031308085e-16;
-    const toint = 1.0 / DBL_EPSILON;
-
-    var i = @bitCast(u64, x);
-    const e = (i >> 52) & 0x7ff;
-
-    if (e >= 0x3ff +% 52 or x == 0)
-        return x;
-    // y = int(x) - x, where int(x) is an integer neighbor of x
-    const y = {
-        @setFloatMode(this, builtin.FloatMode.Strict);
-        if (i >> 63 != 0) {
-            x - toint + toint - x
-        } else {
-            x + toint - toint - x
-        }
-    };
-    // special case because of non-nearest rounding modes
-    if (e <= 0x3ff - 1) {
-        if (i >> 63 != 0)
-            return -1.0;
-        return 0.0;
-    }
-    if (y > 0)
-        return x + y - 1;
-    return x + y;
-}
-
-test "math.floor" {
-    assert(floor(f32(1.234)) == 1.0);
-    assert(floor(f32(-1.234)) == -2.0);
-    assert(floor(f32(999.0)) == 999.0);
-    assert(floor(f32(-999.0)) == -999.0);
-
-    assert(floor(f64(1.234)) == 1.0);
-    assert(floor(f64(-1.234)) == -2.0);
-    assert(floor(f64(999.0)) == 999.0);
-    assert(floor(f64(-999.0)) == -999.0);
-}
-
 /// Returns the absolute value of the integer parameter.
 /// Result is an unsigned integer.
 pub fn absCast(x: var) -> @IntType(false, @typeOf(x).bit_count) {
@@ -377,13 +454,3 @@ test "math.negateCast" {
 
     if (negateCast(u32(@maxValue(i32) + 10))) |_| unreachable else |err| assert(err == error.Overflow);
 }
-
-test "math" {
-    _ = @import("frexp.zig");
-}
-
-
-pub fn approxEq(comptime T: type, x: T, y: T, epsilon: T) -> bool {
-    comptime assert(@typeId(T) == builtin.TypeId.Float);
-    absFloat(x - y) < epsilon
-}

std/math/inf.zig

@@ -0,0 +1,10 @@
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+pub fn inf(comptime T: type) -> T {
+    switch (T) {
+        f32 => @bitCast(f32, math.inf_u32),
+        f64 => @bitCast(f64, math.inf_u64),
+        else => @compileError("inf not implemented for " ++ @typeName(T)),
+    }
+}

std/math/isfinite.zig

@@ -0,0 +1,30 @@
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+pub fn isFinite(x: var) -> bool {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => {
+            const bits = @bitCast(u32, x);
+            bits & 0x7FFFFFFF < 0x7F800000
+        },
+        f64 => {
+            const bits = @bitCast(u64, x);
+            bits & (@maxValue(u64) >> 1) < (0x7FF << 52)
+        },
+        else => {
+            @compileError("isFinite not implemented for " ++ @typeName(T));
+        },
+    }
+}
+
+test "math.isFinite" {
+    assert(isFinite(f32(0.0)));
+    assert(isFinite(f32(-0.0)));
+    assert(isFinite(f64(0.0)));
+    assert(isFinite(f64(-0.0)));
+    assert(!isFinite(math.inf(f32)));
+    assert(!isFinite(-math.inf(f32)));
+    assert(!isFinite(math.inf(f64)));
+    assert(!isFinite(-math.inf(f64)));
+}

std/math/isinf.zig

@@ -0,0 +1,82 @@
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+pub fn isInf(x: var) -> bool {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => {
+            const bits = @bitCast(u32, x);
+            bits & 0x7FFFFFFF == 0x7F800000
+        },
+        f64 => {
+            const bits = @bitCast(u64, x);
+            bits & (@maxValue(u64) >> 1) == (0x7FF << 52)
+        },
+        else => {
+            @compileError("isInf not implemented for " ++ @typeName(T));
+        },
+    }
+}
+
+pub fn isPositiveInf(x: var) -> bool {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => {
+            @bitCast(u32, x) == 0x7F800000
+        },
+        f64 => {
+            @bitCast(u64, x) == 0x7FF << 52
+        },
+        else => {
+            @compileError("isPositiveInf not implemented for " ++ @typeName(T));
+        },
+    }
+}
+
+pub fn isNegativeInf(x: var) -> bool {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => {
+            @bitCast(u32, x) == 0xFF800000
+        },
+        f64 => {
+            @bitCast(u64, x) == 0xFFF << 52
+        },
+        else => {
+            @compileError("isNegativeInf not implemented for " ++ @typeName(T));
+        },
+    }
+}
+
+test "math.isInf" {
+    assert(!isInf(f32(0.0)));
+    assert(!isInf(f32(-0.0)));
+    assert(!isInf(f64(0.0)));
+    assert(!isInf(f64(-0.0)));
+    assert(isInf(math.inf(f32)));
+    assert(isInf(-math.inf(f32)));
+    assert(isInf(math.inf(f64)));
+    assert(isInf(-math.inf(f64)));
+}
+
+test "math.isPositiveInf" {
+    assert(!isPositiveInf(f32(0.0)));
+    assert(!isPositiveInf(f32(-0.0)));
+    assert(!isPositiveInf(f64(0.0)));
+    assert(!isPositiveInf(f64(-0.0)));
+    assert(isPositiveInf(math.inf(f32)));
+    assert(!isPositiveInf(-math.inf(f32)));
+    assert(isPositiveInf(math.inf(f64)));
+    assert(!isPositiveInf(-math.inf(f64)));
+}
+
+test "math.isNegativeInf" {
+    assert(!isNegativeInf(f32(0.0)));
+    assert(!isNegativeInf(f32(-0.0)));
+    assert(!isNegativeInf(f64(0.0)));
+    assert(!isNegativeInf(f64(-0.0)));
+    assert(!isNegativeInf(math.inf(f32)));
+    assert(isNegativeInf(-math.inf(f32)));
+    assert(!isNegativeInf(math.inf(f64)));
+    assert(isNegativeInf(-math.inf(f64)));
+}

std/math/isnan.zig

@@ -0,0 +1,32 @@
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+pub fn isNan(x: var) -> bool {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => {
+            const bits = @bitCast(u32, x);
+            bits & 0x7FFFFFFF > 0x7F800000
+        },
+        f64 => {
+            const bits = @bitCast(u64, x);
+            (bits & (@maxValue(u64) >> 1)) > (u64(0x7FF) << 52)
+        },
+        else => {
+            @compileError("isNan not implemented for " ++ @typeName(T));
+        },
+    }
+}
+
+// Note: A signalling nan is identical to a standard right now by may have a different bit
+// representation in the future when required.
+pub fn isSignalNan(x: var) -> bool {
+    isNan(x)
+}
+
+test "math.isNan" {
+    assert(isNan(math.nan(f32)));
+    assert(isNan(math.nan(f64)));
+    assert(!isNan(f32(1.0)));
+    assert(!isNan(f64(1.0)));
+}

std/math/isnormal.zig

@@ -0,0 +1,26 @@
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+pub fn isNormal(x: var) -> bool {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => {
+            const bits = @bitCast(u32, x);
+            (bits + 0x00800000) & 0x7FFFFFFF >= 0x01000000
+        },
+        f64 => {
+            const bits = @bitCast(u64, x);
+            (bits + (1 << 52)) & (@maxValue(u64) >> 1) >= (1 << 53)
+        },
+        else => {
+            @compileError("isNormal not implemented for " ++ @typeName(T));
+        },
+    }
+}
+
+test "math.isNormal" {
+    assert(!isNormal(math.nan(f32)));
+    assert(!isNormal(math.nan(f64)));
+    assert(isNormal(f32(1.0)));
+    assert(isNormal(f64(1.0)));
+}

std/math/ln.zig

@@ -0,0 +1,173 @@
+// Special Cases:
+//
+// - ln(+inf)  = +inf
+// - ln(0)     = -inf
+// - ln(x)     = nan if x < 0
+// - ln(nan)   = nan
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+pub const ln = ln_workaround;
+
+pub fn ln_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(lnf, x),
+        f64 => @inlineCall(lnd, x),
+        else => @compileError("ln not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn lnf(x_: f32) -> f32 {
+    @setFloatMode(this, @import("builtin").FloatMode.Strict);
+
+    const ln2_hi: f32 = 6.9313812256e-01;
+    const ln2_lo: f32 = 9.0580006145e-06;
+    const Lg1: f32 = 0xaaaaaa.0p-24;
+    const Lg2: f32 = 0xccce13.0p-25;
+    const Lg3: f32 = 0x91e9ee.0p-25;
+    const Lg4: f32 = 0xf89e26.0p-26;
+
+    var x = x_;
+    var ix = @bitCast(u32, x);
+    var k: i32 = 0;
+
+    // x < 2^(-126)
+    if (ix < 0x00800000 or ix >> 31 != 0) {
+        // log(+-0) = -inf
+        if (ix <<% 1 == 0) {
+            return -math.inf(f32);
+        }
+        // log(-#) = nan
+        if (ix >> 31 != 0) {
+            return math.nan(f32);
+        }
+
+        // subnormal, scale x
+        k -= 25;
+        x *= 0x1.0p25;
+        ix = @bitCast(u32, x);
+    } else if (ix >= 0x7F800000) {
+        return x;
+    } else if (ix == 0x3F800000) {
+        return 0;
+    }
+
+    // x into [sqrt(2) / 2, sqrt(2)]
+    ix += 0x3F800000 - 0x3F3504F3;
+    k += i32(ix >> 23) - 0x7F;
+    ix = (ix & 0x007FFFFF) + 0x3F3504F3;
+    x = @bitCast(f32, ix);
+
+    const f = x - 1.0;
+    const s = f / (2.0 + f);
+    const z = s * s;
+    const w = z * z;
+    const t1 = w * (Lg2 + w * Lg4);
+    const t2 = z * (Lg1 + w * Lg3);
+    const R = t2 + t1;
+    const hfsq = 0.5 * f * f;
+    const dk = f32(k);
+
+    s * (hfsq + R) + dk * ln2_lo - hfsq + f + dk * ln2_hi
+}
+
+fn lnd(x_: f64) -> f64 {
+    const ln2_hi: f64 = 6.93147180369123816490e-01;
+    const ln2_lo: f64 = 1.90821492927058770002e-10;
+    const Lg1: f64 = 6.666666666666735130e-01;
+    const Lg2: f64 = 3.999999999940941908e-01;
+    const Lg3: f64 = 2.857142874366239149e-01;
+    const Lg4: f64 = 2.222219843214978396e-01;
+    const Lg5: f64 = 1.818357216161805012e-01;
+    const Lg6: f64 = 1.531383769920937332e-01;
+    const Lg7: f64 = 1.479819860511658591e-01;
+
+    var x = x_;
+    var ix = @bitCast(u64, x);
+    var hx = u32(ix >> 32);
+    var k: i32 = 0;
+
+    if (hx < 0x00100000 or hx >> 31 != 0) {
+        // log(+-0) = -inf
+        if (ix <<% 1 == 0) {
+            return -math.inf(f64);
+        }
+        // log(-#) = nan
+        if (hx >> 31 != 0) {
+            return math.nan(f64);
+        }
+
+        // subnormal, scale x
+        k -= 54;
+        x *= 0x1.0p54;
+        hx = u32(@bitCast(u64, ix) >> 32)
+    }
+    else if (hx >= 0x7FF00000) {
+        return x;
+    }
+    else if (hx == 0x3FF00000 and ix << 32 == 0) {
+        return 0;
+    }
+
+    // x into [sqrt(2) / 2, sqrt(2)]
+    hx += 0x3FF00000 - 0x3FE6A09E;
+    k += i32(hx >> 20) - 0x3FF;
+    hx = (hx & 0x000FFFFF) + 0x3FE6A09E;
+    ix = (u64(hx) << 32) | (ix & 0xFFFFFFFF);
+    x = @bitCast(f64, ix);
+
+    const f = x - 1.0;
+    const hfsq = 0.5 * f * f;
+    const s = f / (2.0 + f);
+    const z = s * s;
+    const w = z * z;
+    const t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
+    const t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
+    const R = t2 + t1;
+    const dk = f64(k);
+
+    s * (hfsq + R) + dk * ln2_lo - hfsq + f + dk * ln2_hi
+}
+
+test "math.ln" {
+    assert(ln(f32(0.2)) == lnf(0.2));
+    assert(ln(f64(0.2)) == lnd(0.2));
+}
+
+test "math.ln32" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, lnf(0.2), -1.609438, epsilon));
+    assert(math.approxEq(f32, lnf(0.8923), -0.113953, epsilon));
+    assert(math.approxEq(f32, lnf(1.5), 0.405465, epsilon));
+    assert(math.approxEq(f32, lnf(37.45), 3.623007, epsilon));
+    assert(math.approxEq(f32, lnf(89.123), 4.490017, epsilon));
+    assert(math.approxEq(f32, lnf(123123.234375), 11.720941, epsilon));
+}
+
+test "math.ln64" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f64, lnd(0.2), -1.609438, epsilon));
+    assert(math.approxEq(f64, lnd(0.8923), -0.113953, epsilon));
+    assert(math.approxEq(f64, lnd(1.5), 0.405465, epsilon));
+    assert(math.approxEq(f64, lnd(37.45), 3.623007, epsilon));
+    assert(math.approxEq(f64, lnd(89.123), 4.490017, epsilon));
+    assert(math.approxEq(f64, lnd(123123.234375), 11.720941, epsilon));
+}
+
+test "math.ln32.special" {
+    assert(math.isPositiveInf(lnf(math.inf(f32))));
+    assert(math.isNegativeInf(lnf(0.0)));
+    assert(math.isNan(lnf(-1.0)));
+    assert(math.isNan(lnf(math.nan(f32))));
+}
+
+test "math.ln64.special" {
+    assert(math.isPositiveInf(lnd(math.inf(f64))));
+    assert(math.isNegativeInf(lnd(0.0)));
+    assert(math.isNan(lnd(-1.0)));
+    assert(math.isNan(lnd(math.nan(f64))));
+}

std/math/log.zig

@@ -0,0 +1,64 @@
+const math = @import("index.zig");
+const builtin = @import("builtin");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const log = log_workaround;
+
+pub fn log_workaround(comptime base: usize, x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (@typeId(T)) {
+        builtin.TypeId.Int => {
+            if (base == 2) {
+                return T.bit_count - 1 - @clz(x);
+            } else {
+                @compileError("TODO implement log for non base 2 integers");
+            }
+        },
+
+        builtin.TypeId.Float => switch (T) {
+            f32 => switch (base) {
+                2 => return math.log2(x),
+                10 => return math.log10(x),
+                else => return f32(math.ln(f64(x)) / math.ln(f64(base))),
+            },
+
+            f64 => switch (base) {
+                2 => return math.log2(x),
+                10 => return math.log10(x),
+                // NOTE: This likely is computed with reduced accuracy.
+                else => return math.ln(x) / math.ln(f64(base)),
+            },
+
+            else => @compileError("log not implemented for " ++ @typeName(T)),
+        },
+
+        else => {
+            @compileError("log expects integer or float, found '" ++ @typeName(T) ++ "'");
+        },
+    }
+}
+
+test "math.log integer" {
+    assert(log(2, u8(0x1)) == 0);
+    assert(log(2, u8(0x2)) == 1);
+    assert(log(2, i16(0x72)) == 6);
+    assert(log(2, u32(0xFFFFFF)) == 23);
+    assert(log(2, u64(0x7FF0123456789ABC)) == 62);
+}
+
+test "math.log float" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, log(6, f32(0.23947)), -0.797723, epsilon));
+    assert(math.approxEq(f32, log(89, f32(0.23947)), -0.318432, epsilon));
+    assert(math.approxEq(f64, log(123897, f64(12389216414)), 1.981724596, epsilon));
+}
+
+test "math.log float_special" {
+    assert(log(2, f32(0.2301974)) == math.log2(f32(0.2301974)));
+    assert(log(10, f32(0.2301974)) == math.log10(f32(0.2301974)));
+
+    assert(log(2, f64(213.23019799993)) == math.log2(f64(213.23019799993)));
+    assert(log(10, f64(213.23019799993)) == math.log10(f64(213.23019799993)));
+}

std/math/log10.zig

@@ -0,0 +1,199 @@
+// Special Cases:
+//
+// - log10(+inf)  = +inf
+// - log10(0)     = -inf
+// - log10(x)     = nan if x < 0
+// - log10(nan)   = nan
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const log10 = log10_workaround;
+
+pub fn log10_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(log10_32, x),
+        f64 => @inlineCall(log10_64, x),
+        else => @compileError("log10 not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn log10_32(x_: f32) -> f32 {
+    const ivln10hi: f32  =  4.3432617188e-01;
+    const ivln10lo: f32  = -3.1689971365e-05;
+    const log10_2hi: f32 =  3.0102920532e-01;
+    const log10_2lo: f32 =  7.9034151668e-07;
+    const Lg1: f32 = 0xaaaaaa.0p-24;
+    const Lg2: f32 = 0xccce13.0p-25;
+    const Lg3: f32 = 0x91e9ee.0p-25;
+    const Lg4: f32 = 0xf89e26.0p-26;
+
+    var x = x_;
+    var u = @bitCast(u32, x);
+    var ix = u;
+    var k: i32 = 0;
+
+    // x < 2^(-126)
+    if (ix < 0x00800000 or ix >> 31 != 0) {
+        // log(+-0) = -inf
+        if (ix <<% 1 == 0) {
+            return -math.inf(f32);
+        }
+        // log(-#) = nan
+        if (ix >> 31 != 0) {
+            return math.nan(f32);
+        }
+
+        k -= 25;
+        x *= 0x1.0p25;
+        ix = @bitCast(u32, x);
+    } else if (ix >= 0x7F800000) {
+        return x;
+    } else if (ix == 0x3F800000) {
+        return 0;
+    }
+
+    // x into [sqrt(2) / 2, sqrt(2)]
+    ix += 0x3F800000 - 0x3F3504F3;
+    k += i32(ix >> 23) - 0x7F;
+    ix = (ix & 0x007FFFFF) + 0x3F3504F3;
+    x = @bitCast(f32, ix);
+
+    const f = x - 1.0;
+    const s = f / (2.0 + f);
+    const z = s * s;
+    const w = z * z;
+    const t1 = w * (Lg2 + w * Lg4);
+    const t2 = z * (Lg1 + w * Lg3);
+    const R = t2 + t1;
+    const hfsq = 0.5 * f * f;
+
+    var hi = f - hfsq;
+    u = @bitCast(u32, hi);
+    u &= 0xFFFFF000;
+    hi = @bitCast(f32, u);
+    const lo = f - hi - hfsq + s * (hfsq + R);
+    const dk = f32(k);
+
+    dk * log10_2lo + (lo + hi) * ivln10lo + lo * ivln10hi + hi * ivln10hi + dk * log10_2hi
+}
+
+fn log10_64(x_: f64) -> f64 {
+    const ivln10hi: f64  = 4.34294481878168880939e-01;
+    const ivln10lo: f64  = 2.50829467116452752298e-11;
+    const log10_2hi: f64 = 3.01029995663611771306e-01;
+    const log10_2lo: f64 = 3.69423907715893078616e-13;
+    const Lg1: f64 = 6.666666666666735130e-01;
+    const Lg2: f64 = 3.999999999940941908e-01;
+    const Lg3: f64 = 2.857142874366239149e-01;
+    const Lg4: f64 = 2.222219843214978396e-01;
+    const Lg5: f64 = 1.818357216161805012e-01;
+    const Lg6: f64 = 1.531383769920937332e-01;
+    const Lg7: f64 = 1.479819860511658591e-01;
+
+    var x = x_;
+    var ix = @bitCast(u64, x);
+    var hx = u32(ix >> 32);
+    var k: i32 = 0;
+
+    if (hx < 0x00100000 or hx >> 31 != 0) {
+        // log(+-0) = -inf
+        if (ix <<% 1 == 0) {
+            return -math.inf(f32);
+        }
+        // log(-#) = nan
+        if (hx >> 31 != 0) {
+            return math.nan(f32);
+        }
+
+        // subnormal, scale x
+        k -= 54;
+        x *= 0x1.0p54;
+        hx = u32(@bitCast(u64, x) >> 32)
+    }
+    else if (hx >= 0x7FF00000) {
+        return x;
+    }
+    else if (hx == 0x3FF00000 and ix << 32 == 0) {
+        return 0;
+    }
+
+    // x into [sqrt(2) / 2, sqrt(2)]
+    hx += 0x3FF00000 - 0x3FE6A09E;
+    k += i32(hx >> 20) - 0x3FF;
+    hx = (hx & 0x000FFFFF) + 0x3FE6A09E;
+    ix = (u64(hx) << 32) | (ix & 0xFFFFFFFF);
+    x = @bitCast(f64, ix);
+
+    const f = x - 1.0;
+    const hfsq = 0.5 * f * f;
+    const s = f / (2.0 + f);
+    const z = s * s;
+    const w = z * z;
+    const t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
+    const t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
+    const R = t2 + t1;
+
+    // hi + lo = f - hfsq + s * (hfsq + R) ~ log(1 + f)
+    var hi = f - hfsq;
+    var hii = @bitCast(u64, hi);
+    hii &= @maxValue(u64) << 32;
+    hi = @bitCast(f64, hii);
+    const lo = f - hi - hfsq + s * (hfsq + R);
+
+    // val_hi + val_lo ~ log10(1 + f) + k * log10(2)
+    var val_hi = hi * ivln10hi;
+    const dk = f64(k);
+    const y = dk * log10_2hi;
+    var val_lo = dk * log10_2lo + (lo + hi) * ivln10lo + lo * ivln10hi;
+
+    // Extra precision multiplication
+    const ww = y + val_hi;
+    val_lo += (y - ww) + val_hi;
+    val_hi = ww;
+
+    val_lo + val_hi
+}
+
+test "math.log10" {
+    assert(log10(f32(0.2)) == log10_32(0.2));
+    assert(log10(f64(0.2)) == log10_64(0.2));
+}
+
+test "math.log10_32" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, log10_32(0.2), -0.698970, epsilon));
+    assert(math.approxEq(f32, log10_32(0.8923), -0.049489, epsilon));
+    assert(math.approxEq(f32, log10_32(1.5), 0.176091, epsilon));
+    assert(math.approxEq(f32, log10_32(37.45), 1.573452, epsilon));
+    assert(math.approxEq(f32, log10_32(89.123), 1.94999, epsilon));
+    assert(math.approxEq(f32, log10_32(123123.234375), 5.09034, epsilon));
+}
+
+test "math.log10_64" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f64, log10_64(0.2), -0.698970, epsilon));
+    assert(math.approxEq(f64, log10_64(0.8923), -0.049489, epsilon));
+    assert(math.approxEq(f64, log10_64(1.5), 0.176091, epsilon));
+    assert(math.approxEq(f64, log10_64(37.45), 1.573452, epsilon));
+    assert(math.approxEq(f64, log10_64(89.123), 1.94999, epsilon));
+    assert(math.approxEq(f64, log10_64(123123.234375), 5.09034, epsilon));
+}
+
+test "math.log10_32.special" {
+    assert(math.isPositiveInf(log10_32(math.inf(f32))));
+    assert(math.isNegativeInf(log10_32(0.0)));
+    assert(math.isNan(log10_32(-1.0)));
+    assert(math.isNan(log10_32(math.nan(f32))));
+}
+
+test "math.log10_64.special" {
+    assert(math.isPositiveInf(log10_64(math.inf(f64))));
+    assert(math.isNegativeInf(log10_64(0.0)));
+    assert(math.isNan(log10_64(-1.0)));
+    assert(math.isNan(log10_64(math.nan(f64))));
+}

std/math/log1p.zig

@@ -0,0 +1,226 @@
+// Special Cases:
+//
+// - log1p(+inf)  = +inf
+// - log1p(+-0)   = +-0
+// - log1p(-1)    = -inf
+// - log1p(x)     = nan if x < -1
+// - log1p(nan)   = nan
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const log1p = log1p_workaround;
+
+pub fn log1p_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(log1p_32, x),
+        f64 => @inlineCall(log1p_64, x),
+        else => @compileError("log1p not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn log1p_32(x: f32) -> f32 {
+    const ln2_hi = 6.9313812256e-01;
+    const ln2_lo = 9.0580006145e-06;
+    const Lg1: f32 = 0xaaaaaa.0p-24;
+    const Lg2: f32 = 0xccce13.0p-25;
+    const Lg3: f32 = 0x91e9ee.0p-25;
+    const Lg4: f32 = 0xf89e26.0p-26;
+
+    const u = @bitCast(u32, x);
+    var ix = u;
+    var k: i32 = 1;
+    var f: f32 = undefined;
+    var c: f32 = undefined;
+
+    // 1 + x < sqrt(2)+
+    if (ix < 0x3ED413D0 or ix >> 31 != 0) {
+        // x <= -1.0
+        if (ix >= 0xBF800000) {
+            // log1p(-1) = -inf
+            if (x == -1.0) {
+                return -math.inf(f32);
+            }
+            // log1p(x < -1) = nan
+            else {
+                return math.nan(f32);
+            }
+        }
+        // |x| < 2^(-24)
+        if ((ix <<% 1) < (0x33800000 << 1)) {
+            // underflow if subnormal
+            if (ix & 0x7F800000 == 0) {
+                math.forceEval(x * x);
+            }
+            return x;
+        }
+        // sqrt(2) / 2- <= 1 + x < sqrt(2)+
+        if (ix <= 0xBE95F619) {
+            k = 0;
+            c = 0;
+            f = x;
+        }
+    } else if (ix >= 0x7F800000) {
+        return x;
+    }
+
+    if (k != 0) {
+        const uf = 1 + x;
+        var iu = @bitCast(u32, uf);
+        iu += 0x3F800000 - 0x3F3504F3;
+        k = i32(iu >> 23) - 0x7F;
+
+        // correction to avoid underflow in c / u
+        if (k < 25) {
+            c = if (k >= 2) 1 - (uf - x) else x - (uf - 1);
+            c /= uf;
+        } else {
+            c = 0;
+        }
+
+        // u into [sqrt(2)/2, sqrt(2)]
+        iu = (iu & 0x007FFFFF) + 0x3F3504F3;
+        f = @bitCast(f32, iu) - 1;
+    }
+
+    const s = f / (2.0 + f);
+    const z = s * s;
+    const w = z * z;
+    const t1 = w * (Lg2 + w * Lg4);
+    const t2 = z * (Lg1 + w * Lg3);
+    const R = t2 + t1;
+    const hfsq = 0.5 * f * f;
+    const dk = f32(k);
+
+    s * (hfsq + R) + (dk * ln2_lo + c) - hfsq + f + dk * ln2_hi
+}
+
+fn log1p_64(x: f64) -> f64 {
+    const ln2_hi: f64 = 6.93147180369123816490e-01;
+    const ln2_lo: f64 = 1.90821492927058770002e-10;
+    const Lg1: f64 = 6.666666666666735130e-01;
+    const Lg2: f64 = 3.999999999940941908e-01;
+    const Lg3: f64 = 2.857142874366239149e-01;
+    const Lg4: f64 = 2.222219843214978396e-01;
+    const Lg5: f64 = 1.818357216161805012e-01;
+    const Lg6: f64 = 1.531383769920937332e-01;
+    const Lg7: f64 = 1.479819860511658591e-01;
+
+    var ix = @bitCast(u64, x);
+    var hx = u32(ix >> 32);
+    var k: i32 = 1;
+    var c: f64 = undefined;
+    var f: f64 = undefined;
+
+    // 1 + x < sqrt(2)
+    if (hx < 0x3FDA827A or hx >> 31 != 0) {
+        // x <= -1.0
+        if (hx >= 0xBFF00000) {
+            // log1p(-1) = -inf
+            if (x == -1.0) {
+                return -math.inf(f64);
+            }
+            // log1p(x < -1) = nan
+            else {
+                return math.nan(f64);
+            }
+        }
+        // |x| < 2^(-53)
+        if ((hx <<% 1) < (0x3CA00000 << 1)) {
+            if ((hx & 0x7FF00000) == 0) {
+                math.raiseUnderflow();
+            }
+            return x;
+        }
+        // sqrt(2) / 2- <= 1 + x < sqrt(2)+
+        if (hx <= 0xBFD2BEC4) {
+            k = 0;
+            c = 0;
+            f = x;
+        }
+    }
+    else if (hx >= 0x7FF00000) {
+        return x;
+    }
+
+    if (k != 0) {
+        const uf = 1 + x;
+        const hu = @bitCast(u64, uf);
+        var iu = u32(hu >> 32);
+        iu += 0x3FF00000 - 0x3FE6A09E;
+        k = i32(iu >> 20) - 0x3FF;
+
+        // correction to avoid underflow in c / u
+        if (k < 54) {
+            c = if (k >= 2) 1 - (uf - x) else x - (uf - 1);
+            c /= uf;
+        } else {
+            c = 0;
+        }
+
+        // u into [sqrt(2)/2, sqrt(2)]
+        iu = (iu & 0x000FFFFF) + 0x3FE6A09E;
+        const iq = (u64(iu) << 32) | (hu & 0xFFFFFFFF);
+        f = @bitCast(f64, iq) - 1;
+    }
+
+    const hfsq = 0.5 * f * f;
+    const s = f / (2.0 + f);
+    const z = s * s;
+    const w = z * z;
+    const t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
+    const t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
+    const R = t2 + t1;
+    const dk = f64(k);
+
+    s * (hfsq + R) + (dk * ln2_lo + c) - hfsq + f + dk * ln2_hi
+}
+
+test "math.log1p" {
+    assert(log1p(f32(0.0)) == log1p_32(0.0));
+    assert(log1p(f64(0.0)) == log1p_64(0.0));
+}
+
+test "math.log1p_32" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, log1p_32(0.0), 0.0, epsilon));
+    assert(math.approxEq(f32, log1p_32(0.2), 0.182322, epsilon));
+    assert(math.approxEq(f32, log1p_32(0.8923), 0.637793, epsilon));
+    assert(math.approxEq(f32, log1p_32(1.5), 0.916291, epsilon));
+    assert(math.approxEq(f32, log1p_32(37.45), 3.649359, epsilon));
+    assert(math.approxEq(f32, log1p_32(89.123), 4.501175, epsilon));
+    assert(math.approxEq(f32, log1p_32(123123.234375), 11.720949, epsilon));
+}
+
+test "math.log1p_64" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f64, log1p_64(0.0), 0.0, epsilon));
+    assert(math.approxEq(f64, log1p_64(0.2), 0.182322, epsilon));
+    assert(math.approxEq(f64, log1p_64(0.8923), 0.637793, epsilon));
+    assert(math.approxEq(f64, log1p_64(1.5), 0.916291, epsilon));
+    assert(math.approxEq(f64, log1p_64(37.45), 3.649359, epsilon));
+    assert(math.approxEq(f64, log1p_64(89.123), 4.501175, epsilon));
+    assert(math.approxEq(f64, log1p_64(123123.234375), 11.720949, epsilon));
+}
+
+test "math.log1p_32.special" {
+    assert(math.isPositiveInf(log1p_32(math.inf(f32))));
+    assert(log1p_32(0.0) == 0.0);
+    assert(log1p_32(-0.0) == -0.0);
+    assert(math.isNegativeInf(log1p_32(-1.0)));
+    assert(math.isNan(log1p_32(-2.0)));
+    assert(math.isNan(log1p_32(math.nan(f32))));
+}
+
+test "math.log1p_64.special" {
+    assert(math.isPositiveInf(log1p_64(math.inf(f64))));
+    assert(log1p_64(0.0) == 0.0);
+    assert(log1p_64(-0.0) == -0.0);
+    assert(math.isNegativeInf(log1p_64(-1.0)));
+    assert(math.isNan(log1p_64(-2.0)));
+    assert(math.isNan(log1p_64(math.nan(f64))));
+}

std/math/log2.zig

@@ -0,0 +1,189 @@
+// Special Cases:
+//
+// - log2(+inf)  = +inf
+// - log2(0)     = -inf
+// - log2(x)     = nan if x < 0
+// - log2(nan)   = nan
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const log2 = log2_workaround;
+
+pub fn log2_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(log2_32, x),
+        f64 => @inlineCall(log2_64, x),
+        else => @compileError("log2 not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn log2_32(x_: f32) -> f32 {
+    const ivln2hi: f32 =  1.4428710938e+00;
+    const ivln2lo: f32 = -1.7605285393e-04;
+    const Lg1: f32 = 0xaaaaaa.0p-24;
+    const Lg2: f32 = 0xccce13.0p-25;
+    const Lg3: f32 = 0x91e9ee.0p-25;
+    const Lg4: f32 = 0xf89e26.0p-26;
+
+    var x = x_;
+    var u = @bitCast(u32, x);
+    var ix = u;
+    var k: i32 = 0;
+
+    // x < 2^(-126)
+    if (ix < 0x00800000 or ix >> 31 != 0) {
+        // log(+-0) = -inf
+        if (ix <<% 1 == 0) {
+            return -math.inf(f32);
+        }
+        // log(-#) = nan
+        if (ix >> 31 != 0) {
+            return math.nan(f32);
+        }
+
+        k -= 25;
+        x *= 0x1.0p25;
+        ix = @bitCast(u32, x);
+    } else if (ix >= 0x7F800000) {
+        return x;
+    } else if (ix == 0x3F800000) {
+        return 0;
+    }
+
+    // x into [sqrt(2) / 2, sqrt(2)]
+    ix += 0x3F800000 - 0x3F3504F3;
+    k += i32(ix >> 23) - 0x7F;
+    ix = (ix & 0x007FFFFF) + 0x3F3504F3;
+    x = @bitCast(f32, ix);
+
+    const f = x - 1.0;
+    const s = f / (2.0 + f);
+    const z = s * s;
+    const w = z * z;
+    const t1 = w * (Lg2 + w * Lg4);
+    const t2 = z * (Lg1 + w * Lg3);
+    const R = t2 + t1;
+    const hfsq = 0.5 * f * f;
+
+    var hi = f - hfsq;
+    u = @bitCast(u32, hi);
+    u &= 0xFFFFF000;
+    hi = @bitCast(f32, u);
+    const lo = f - hi - hfsq + s * (hfsq + R);
+    (lo + hi) * ivln2lo + lo * ivln2hi + hi * ivln2hi + f32(k)
+}
+
+fn log2_64(x_: f64) -> f64 {
+    const ivln2hi: f64 = 1.44269504072144627571e+00;
+    const ivln2lo: f64 = 1.67517131648865118353e-10;
+    const Lg1: f64 = 6.666666666666735130e-01;
+    const Lg2: f64 = 3.999999999940941908e-01;
+    const Lg3: f64 = 2.857142874366239149e-01;
+    const Lg4: f64 = 2.222219843214978396e-01;
+    const Lg5: f64 = 1.818357216161805012e-01;
+    const Lg6: f64 = 1.531383769920937332e-01;
+    const Lg7: f64 = 1.479819860511658591e-01;
+
+    var x = x_;
+    var ix = @bitCast(u64, x);
+    var hx = u32(ix >> 32);
+    var k: i32 = 0;
+
+    if (hx < 0x00100000 or hx >> 31 != 0) {
+        // log(+-0) = -inf
+        if (ix <<% 1 == 0) {
+            return -math.inf(f64);
+        }
+        // log(-#) = nan
+        if (hx >> 31 != 0) {
+            return math.nan(f64);
+        }
+
+        // subnormal, scale x
+        k -= 54;
+        x *= 0x1.0p54;
+        hx = u32(@bitCast(u64, x) >> 32);
+    }
+    else if (hx >= 0x7FF00000) {
+        return x;
+    }
+    else if (hx == 0x3FF00000 and ix << 32 == 0) {
+        return 0;
+    }
+
+    // x into [sqrt(2) / 2, sqrt(2)]
+    hx += 0x3FF00000 - 0x3FE6A09E;
+    k += i32(hx >> 20) - 0x3FF;
+    hx = (hx & 0x000FFFFF) + 0x3FE6A09E;
+    ix = (u64(hx) << 32) | (ix & 0xFFFFFFFF);
+    x = @bitCast(f64, ix);
+
+    const f = x - 1.0;
+    const hfsq = 0.5 * f * f;
+    const s = f / (2.0 + f);
+    const z = s * s;
+    const w = z * z;
+    const t1 = w * (Lg2 + w * (Lg4 + w * Lg6));
+    const t2 = z * (Lg1 + w * (Lg3 + w * (Lg5 + w * Lg7)));
+    const R = t2 + t1;
+
+    // hi + lo = f - hfsq + s * (hfsq + R) ~ log(1 + f)
+    var hi = f - hfsq;
+    var hii = @bitCast(u64, hi);
+    hii &= @maxValue(u64) << 32;
+    hi = @bitCast(f64, hii);
+    const lo = f - hi - hfsq + s * (hfsq + R);
+
+    var val_hi = hi * ivln2hi;
+    var val_lo = (lo + hi) * ivln2lo + lo * ivln2hi;
+
+    // spadd(val_hi, val_lo, y)
+    const y = f64(k);
+    const ww = y + val_hi;
+    val_lo += (y - ww) + val_hi;
+    val_hi = ww;
+
+    val_lo + val_hi
+}
+
+test "math.log2" {
+    assert(log2(f32(0.2)) == log2_32(0.2));
+    assert(log2(f64(0.2)) == log2_64(0.2));
+}
+
+test "math.log2_32" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, log2_32(0.2), -2.321928, epsilon));
+    assert(math.approxEq(f32, log2_32(0.8923), -0.164399, epsilon));
+    assert(math.approxEq(f32, log2_32(1.5), 0.584962, epsilon));
+    assert(math.approxEq(f32, log2_32(37.45), 5.226894, epsilon));
+    assert(math.approxEq(f32, log2_32(123123.234375), 16.909744, epsilon));
+}
+
+test "math.log2_64" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f64, log2_64(0.2), -2.321928, epsilon));
+    assert(math.approxEq(f64, log2_64(0.8923), -0.164399, epsilon));
+    assert(math.approxEq(f64, log2_64(1.5), 0.584962, epsilon));
+    assert(math.approxEq(f64, log2_64(37.45), 5.226894, epsilon));
+    assert(math.approxEq(f64, log2_64(123123.234375), 16.909744, epsilon));
+}
+
+test "math.log2_32.special" {
+    assert(math.isPositiveInf(log2_32(math.inf(f32))));
+    assert(math.isNegativeInf(log2_32(0.0)));
+    assert(math.isNan(log2_32(-1.0)));
+    assert(math.isNan(log2_32(math.nan(f32))));
+}
+
+test "math.log2_64.special" {
+    assert(math.isPositiveInf(log2_64(math.inf(f64))));
+    assert(math.isNegativeInf(log2_64(0.0)));
+    assert(math.isNan(log2_64(-1.0)));
+    assert(math.isNan(log2_64(math.nan(f64))));
+}

std/math/modf.zig

@@ -0,0 +1,204 @@
+// Special Cases:
+//
+// - modf(+-inf) = +-inf, nan
+// - modf(nan)   = nan, nan
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const modf = modf_workaround;
+
+fn modf_result(comptime T: type) -> type {
+    struct {
+        fpart: T,
+        ipart: T,
+    }
+}
+pub const modf32_result = modf_result(f32);
+pub const modf64_result = modf_result(f64);
+
+pub fn modf_workaround(x: var) -> modf_result(@typeOf(x)) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(modf32, x),
+        f64 => @inlineCall(modf64, x),
+        else => @compileError("modf not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn modf32(x: f32) -> modf32_result {
+    var result: modf32_result = undefined;
+
+    const u = @bitCast(u32, x);
+    const e = i32((u >> 23) & 0xFF) - 0x7F;
+    const us = u & 0x80000000;
+
+    // TODO: Shouldn't need this.
+    if (math.isInf(x)) {
+        result.ipart = x;
+        result.fpart = math.nan(f32);
+        return result;
+    }
+
+    // no fractional part
+    if (e >= 23) {
+        result.ipart = x;
+        if (e == 0x80 and u <<% 9 != 0) { // nan
+            result.fpart = x;
+        } else {
+            result.fpart = @bitCast(f32, us);
+        }
+        return result;
+    }
+
+    // no integral part
+    if (e < 0) {
+        result.ipart = @bitCast(f32, us);
+        result.fpart = x;
+        return result;
+    }
+
+    const mask = 0x007FFFFF >> u32(e);
+    if (u & mask == 0) {
+        result.ipart = x;
+        result.fpart = @bitCast(f32, us);
+        return result;
+    }
+
+    const uf = @bitCast(f32, u & ~mask);
+    result.ipart = uf;
+    result.fpart = x - uf;
+    result
+}
+
+fn modf64(x: f64) -> modf64_result {
+    var result: modf64_result = undefined;
+
+    const u = @bitCast(u64, x);
+    const e = i32((u >> 52) & 0x7FF) - 0x3FF;
+    const us = u & (1 << 63);
+
+    if (math.isInf(x)) {
+        result.ipart = x;
+        result.fpart = math.nan(f64);
+        return result;
+    }
+
+    // no fractional part
+    if (e >= 52) {
+        result.ipart = x;
+        if (e == 0x400 and u <<% 12 != 0) { // nan
+            result.fpart = x;
+        } else {
+            result.fpart = @bitCast(f64, us);
+        }
+        return result;
+    }
+
+    // no integral part
+    if (e < 0) {
+        result.ipart = @bitCast(f64, us);
+        result.fpart = x;
+        return result;
+    }
+
+    const mask = @maxValue(u64) >> 12 >> u64(e);
+    if (u & mask == 0) {
+        result.ipart = x;
+        result.fpart = @bitCast(f64, us);
+        return result;
+    }
+
+    const uf = @bitCast(f64, u & ~mask);
+    result.ipart = uf;
+    result.fpart = x - uf;
+    result
+}
+
+test "math.modf" {
+    const a = modf(f32(1.0));
+    const b = modf32(1.0);
+    // NOTE: No struct comparison on generic return type function? non-named, makes sense, but still.
+    assert(a.ipart == b.ipart and a.fpart == b.fpart);
+
+    const c = modf(f64(1.0));
+    const d = modf64(1.0);
+    assert(a.ipart == b.ipart and a.fpart == b.fpart);
+}
+
+test "math.modf32" {
+    const epsilon = 0.000001;
+    var r: modf32_result = undefined;
+
+    r = modf32(1.0);
+    assert(math.approxEq(f32, r.ipart, 1.0, epsilon));
+    assert(math.approxEq(f32, r.fpart, 0.0, epsilon));
+
+    r = modf32(2.545);
+    assert(math.approxEq(f32, r.ipart, 2.0, epsilon));
+    assert(math.approxEq(f32, r.fpart, 0.545, epsilon));
+
+    r = modf32(3.978123);
+    assert(math.approxEq(f32, r.ipart, 3.0, epsilon));
+    assert(math.approxEq(f32, r.fpart, 0.978123, epsilon));
+
+    r = modf32(43874.3);
+    assert(math.approxEq(f32, r.ipart, 43874, epsilon));
+    assert(math.approxEq(f32, r.fpart, 0.300781, epsilon));
+
+    r = modf32(1234.340780);
+    assert(math.approxEq(f32, r.ipart, 1234, epsilon));
+    assert(math.approxEq(f32, r.fpart, 0.340820, epsilon));
+}
+
+test "math.modf64" {
+    const epsilon = 0.000001;
+    var r: modf64_result = undefined;
+
+    r = modf64(1.0);
+    assert(math.approxEq(f64, r.ipart, 1.0, epsilon));
+    assert(math.approxEq(f64, r.fpart, 0.0, epsilon));
+
+    r = modf64(2.545);
+    assert(math.approxEq(f64, r.ipart, 2.0, epsilon));
+    assert(math.approxEq(f64, r.fpart, 0.545, epsilon));
+
+    r = modf64(3.978123);
+    assert(math.approxEq(f64, r.ipart, 3.0, epsilon));
+    assert(math.approxEq(f64, r.fpart, 0.978123, epsilon));
+
+    r = modf64(43874.3);
+    assert(math.approxEq(f64, r.ipart, 43874, epsilon));
+    assert(math.approxEq(f64, r.fpart, 0.3, epsilon));
+
+    r = modf64(1234.340780);
+    assert(math.approxEq(f64, r.ipart, 1234, epsilon));
+    assert(math.approxEq(f64, r.fpart, 0.340780, epsilon));
+}
+
+test "math.modf32.special" {
+    var r: modf32_result = undefined;
+
+    r = modf32(math.inf(f32));
+    assert(math.isPositiveInf(r.ipart) and math.isNan(r.fpart));
+
+    r = modf32(-math.inf(f32));
+    assert(math.isNegativeInf(r.ipart) and math.isNan(r.fpart));
+
+    r = modf32(math.nan(f32));
+    assert(math.isNan(r.ipart) and math.isNan(r.fpart));
+}
+
+test "math.modf64.special" {
+    var r: modf64_result = undefined;
+
+    r = modf64(math.inf(f64));
+    assert(math.isPositiveInf(r.ipart) and math.isNan(r.fpart));
+
+    r = modf64(-math.inf(f64));
+    assert(math.isNegativeInf(r.ipart) and math.isNan(r.fpart));
+
+    r = modf64(math.nan(f64));
+    assert(math.isNan(r.ipart) and math.isNan(r.fpart));
+}

std/math/nan.zig

@@ -0,0 +1,23 @@
+const math = @import("index.zig");
+
+pub const nan = nan_workaround;
+
+pub fn nan_workaround(comptime T: type) -> T {
+    switch (T) {
+        f32 => @bitCast(f32, math.nan_u32),
+        f64 => @bitCast(f64, math.nan_u64),
+        else => @compileError("nan not implemented for " ++ @typeName(T)),
+    }
+}
+
+pub const snan = snan_workaround;
+
+// Note: A signalling nan is identical to a standard right now by may have a different bit
+// representation in the future when required.
+pub fn snan_workaround(comptime T: type) -> T {
+    switch (T) {
+        f32 => @bitCast(f32, math.nan_u32),
+        f64 => @bitCast(f64, math.nan_u64),
+        else => @compileError("snan not implemented for " ++ @typeName(T)),
+    }
+}

std/math/pow.zig

@@ -0,0 +1,234 @@
+// Special Cases:
+//
+//  pow(x, +-0)    = 1 for any x
+//  pow(1, y)      = 1 for any y
+//  pow(x, 1)      = x for any x
+//  pow(nan, y)    = nan
+//  pow(x, nan)    = nan
+//  pow(+-0, y)    = +-inf for y an odd integer < 0
+//  pow(+-0, -inf) = +inf
+//  pow(+-0, +inf) = +0
+//  pow(+-0, y)    = +inf for finite y < 0 and not an odd integer
+//  pow(+-0, y)    = +-0 for y an odd integer > 0
+//  pow(+-0, y)    = +0 for finite y > 0 and not an odd integer
+//  pow(-1, +-inf) = 1
+//  pow(x, +inf)   = +inf for |x| > 1
+//  pow(x, -inf)   = +0 for |x| > 1
+//  pow(x, +inf)   = +0 for |x| < 1
+//  pow(x, -inf)   = +inf for |x| < 1
+//  pow(+inf, y)   = +inf for y > 0
+//  pow(+inf, y)   = +0 for y < 0
+//  pow(-inf, y)   = pow(-0, -y)
+//  pow(x, y)      = nan for finite x < 0 and finite non-integer y
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const pow = pow_workaround;
+
+// This implementation is taken from the go stlib, musl is a bit more complex.
+pub fn pow_workaround(comptime T: type, x: T, y: T) -> T {
+
+    @setFloatMode(this, @import("builtin").FloatMode.Strict);
+
+    if (T != f32 and T != f64) {
+        @compileError("pow not implemented for " ++ @typeName(T));
+    }
+
+    // pow(x, +-0) = 1      for all x
+    // pow(1, y) = 1        for all y
+    if (y == 0 or x == 1) {
+        return 1;
+    }
+
+    // pow(nan, y) = nan    for all y
+    // pow(x, nan) = nan    for all x
+    if (math.isNan(x) or math.isNan(y)) {
+        return math.nan(T);
+    }
+
+    // pow(x, 1) = x        for all x
+    if (y == 1) {
+        return x;
+    }
+
+    // special case sqrt
+    if (y == 0.5) {
+        return math.sqrt(x);
+    }
+
+    if (y == -0.5) {
+        return 1 / math.sqrt(x);
+    }
+
+    if (x == 0) {
+        if (y < 0) {
+            // pow(+-0, y) = +- 0   for y an odd integer
+            if (isOddInteger(y)) {
+                return math.copysign(T, math.inf(T), x);
+            }
+            // pow(+-0, y) = +inf   for y an even integer
+            else {
+                return math.inf(T);
+            }
+        } else {
+            if (isOddInteger(y)) {
+                return x;
+            } else {
+                return 0;
+            }
+        }
+    }
+
+    if (math.isInf(y)) {
+        // pow(-1, inf) = 1     for all x
+        if (x == -1) {
+            return 1.0;
+        }
+        // pow(x, +inf) = +0    for |x| < 1
+        // pow(x, -inf) = +0    for |x| > 1
+        else if ((math.fabs(x) < 1) == math.isPositiveInf(y)) {
+            return 0;
+        }
+        // pow(x, -inf) = +inf  for |x| < 1
+        // pow(x, +inf) = +inf  for |x| > 1
+        else {
+            return math.inf(T);
+        }
+    }
+
+    if (math.isInf(x)) {
+        if (math.isNegativeInf(x)) {
+            return pow(T, 1 / x, -y);
+        }
+        // pow(+inf, y) = +0    for y < 0
+        else if (y < 0) {
+            return 0;
+        }
+        // pow(+inf, y) = +0    for y > 0
+        else if (y > 0) {
+            return math.inf(T);
+        }
+    }
+
+    var ay = y;
+    var flip = false;
+    if (ay < 0) {
+        ay = -ay;
+        flip = true;
+    }
+
+    const r1 = math.modf(ay);
+    var yi = r1.ipart;
+    var yf = r1.fpart;
+
+    if (yf != 0 and x < 0) {
+        return math.nan(T);
+    }
+    if (yi >= 1 << (T.bit_count - 1)) {
+        return math.exp(y * math.ln(x));
+    }
+
+    // a = a1 * 2^ae
+    var a1: T = 1.0;
+    var ae: i32 = 0;
+
+    // a *= x^yf
+    if (yf != 0) {
+        if (yf > 0.5) {
+            yf -= 1;
+            yi += 1;
+        }
+        a1 = math.exp(yf * math.ln(x));
+    }
+
+    // a *= x^yi
+    const r2 = math.frexp(x);
+    var xe = r2.exponent;
+    var x1 = r2.significand;
+
+    var i = i32(yi);
+    while (i != 0) : (i >>= 1) {
+        if (i & 1 == 1) {
+            a1 *= x1;
+            ae += xe;
+        }
+        x1 *= x1;
+        xe <<= 1;
+        if (x1 < 0.5) {
+            x1 += x1;
+            xe -= 1;
+        }
+    }
+
+    // a *= a1 * 2^ae
+    if (flip) {
+        a1 = 1 / a1;
+        ae = -ae;
+    }
+
+    math.scalbn(a1, ae)
+}
+
+fn isOddInteger(x: f64) -> bool {
+    const r = math.modf(x);
+    r.fpart == 0.0 and i64(r.ipart) & 1 == 1
+}
+
+test "math.pow" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, pow(f32, 0.0, 3.3), 0.0, epsilon));
+    assert(math.approxEq(f32, pow(f32, 0.8923, 3.3), 0.686572, epsilon));
+    assert(math.approxEq(f32, pow(f32, 0.2, 3.3), 0.004936, epsilon));
+    assert(math.approxEq(f32, pow(f32, 1.5, 3.3), 3.811546, epsilon));
+    assert(math.approxEq(f32, pow(f32, 37.45, 3.3), 155736.703125, epsilon));
+    assert(math.approxEq(f32, pow(f32, 89.123, 3.3), 2722489.5, epsilon));
+
+    assert(math.approxEq(f64, pow(f64, 0.0, 3.3), 0.0, epsilon));
+    assert(math.approxEq(f64, pow(f64, 0.8923, 3.3), 0.686572, epsilon));
+    assert(math.approxEq(f64, pow(f64, 0.2, 3.3), 0.004936, epsilon));
+    assert(math.approxEq(f64, pow(f64, 1.5, 3.3), 3.811546, epsilon));
+    assert(math.approxEq(f64, pow(f64, 37.45, 3.3), 155736.7160616, epsilon));
+    assert(math.approxEq(f64, pow(f64, 89.123, 3.3), 2722490.231436, epsilon));
+}
+
+test "math.pow.special" {
+    const epsilon = 0.000001;
+
+    assert(pow(f32, 4, 0.0) == 1.0);
+    assert(pow(f32, 7, -0.0) == 1.0);
+    assert(pow(f32, 45, 1.0) == 45);
+    assert(pow(f32, -45, 1.0) == -45);
+    assert(math.isNan(pow(f32, math.nan(f32), 5.0)));
+    assert(math.isNan(pow(f32, 5.0, math.nan(f32))));
+    assert(math.isPositiveInf(pow(f32, 0.0, -1.0)));
+    assert(math.isNegativeInf(pow(f32, -0.0, -3.0)));
+    assert(math.isPositiveInf(pow(f32, 0.0, -math.inf(f32))));
+    assert(math.isPositiveInf(pow(f32, -0.0, -math.inf(f32))));
+    assert(pow(f32, 0.0, math.inf(f32)) == 0.0);
+    assert(pow(f32, -0.0, math.inf(f32)) == 0.0);
+    assert(math.isPositiveInf(pow(f32, 0.0, -2.0)));
+    assert(math.isPositiveInf(pow(f32, -0.0, -2.0)));
+    assert(pow(f32, 0.0, 1.0) == 0.0);
+    assert(pow(f32, -0.0, 1.0) == -0.0);
+    assert(pow(f32, 0.0, 2.0) == 0.0);
+    assert(pow(f32, -0.0, 2.0) == 0.0);
+    assert(math.approxEq(f32, pow(f32, -1.0, math.inf(f32)), 1.0, epsilon));
+    assert(math.approxEq(f32, pow(f32, -1.0, -math.inf(f32)), 1.0, epsilon));
+    assert(math.isPositiveInf(pow(f32, 1.2, math.inf(f32))));
+    assert(math.isPositiveInf(pow(f32, -1.2, math.inf(f32))));
+    assert(pow(f32, 1.2, -math.inf(f32)) == 0.0);
+    assert(pow(f32, -1.2, -math.inf(f32)) == 0.0);
+    assert(pow(f32, 0.2, math.inf(f32)) == 0.0);
+    assert(pow(f32, -0.2, math.inf(f32)) == 0.0);
+    assert(math.isPositiveInf(pow(f32, 0.2, -math.inf(f32))));
+    assert(math.isPositiveInf(pow(f32, -0.2, -math.inf(f32))));
+    assert(math.isPositiveInf(pow(f32, math.inf(f32), 1.0)));
+    assert(pow(f32, math.inf(f32), -1.0) == 0.0);
+    assert(pow(f32, -math.inf(f32), 5.0) == pow(f32, -0.0, -5.0));
+    assert(pow(f32, -math.inf(f32), -5.2) == pow(f32, -0.0, 5.2));
+    assert(math.isNan(pow(f32, -1.0, 1.2)));
+    assert(math.isNan(pow(f32, -12.4, 78.5)));
+}

std/math/round.zig

@@ -0,0 +1,130 @@
+// Special Cases:
+//
+// - round(+-0)   = +-0
+// - round(+-inf) = +-inf
+// - round(nan)   = nan
+
+const builtin = @import("builtin");
+const assert = @import("../debug.zig").assert;
+const math = @import("index.zig");
+
+// TODO issue #393
+pub const round = round_workaround;
+
+pub fn round_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(round32, x),
+        f64 => @inlineCall(round64, x),
+        else => @compileError("round not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn round32(x_: f32) -> f32 {
+    var x = x_;
+    const u = @bitCast(u32, x);
+    const e = (u >> 23) & 0xFF;
+    var y: f32 = undefined;
+
+    if (e >= 0x7F+23) {
+        return x;
+    }
+    if (u >> 31 != 0) {
+        x = -x;
+    }
+    if (e < 0x7F-1) {
+        math.forceEval(x + math.f32_toint);
+        return 0 * @bitCast(f32, u);
+    }
+
+    {
+        @setFloatMode(this, builtin.FloatMode.Strict);
+        y = x + math.f32_toint - math.f32_toint - x;
+    }
+
+    if (y > 0.5) {
+        y = y + x - 1;
+    } else if (y <= -0.5) {
+        y = y + x + 1;
+    } else {
+        y = y + x;
+    }
+
+    if (u >> 31 != 0) {
+        -y
+    } else {
+        y
+    }
+}
+
+fn round64(x_: f64) -> f64 {
+    var x = x_;
+    const u = @bitCast(u64, x);
+    const e = (u >> 52) & 0x7FF;
+    var y: f64 = undefined;
+
+    if (e >= 0x3FF+52) {
+        return x;
+    }
+    if (u >> 63 != 0) {
+        x = -x;
+    }
+    if (e < 0x3ff-1) {
+        math.forceEval(x + math.f64_toint);
+        return 0 * @bitCast(f64, u);
+    }
+
+    {
+        @setFloatMode(this, builtin.FloatMode.Strict);
+        y = x + math.f64_toint - math.f64_toint - x;
+    }
+
+    if (y > 0.5) {
+        y = y + x - 1;
+    } else if (y <= -0.5) {
+        y = y + x + 1;
+    } else {
+        y = y + x;
+    }
+
+    if (u >> 63 != 0) {
+        -y
+    } else {
+        y
+    }
+}
+
+test "math.round" {
+    assert(round(f32(1.3)) == round32(1.3));
+    assert(round(f64(1.3)) == round64(1.3));
+}
+
+test "math.round32" {
+    assert(round32(1.3) == 1.0);
+    assert(round32(-1.3) == -1.0);
+    assert(round32(0.2) == 0.0);
+    assert(round32(1.8) == 2.0);
+}
+
+test "math.round64" {
+    assert(round64(1.3) == 1.0);
+    assert(round64(-1.3) == -1.0);
+    assert(round64(0.2) == 0.0);
+    assert(round64(1.8) == 2.0);
+}
+
+test "math.round32.special" {
+    assert(round32(0.0) == 0.0);
+    assert(round32(-0.0) == -0.0);
+    assert(math.isPositiveInf(round32(math.inf(f32))));
+    assert(math.isNegativeInf(round32(-math.inf(f32))));
+    assert(math.isNan(round32(math.nan(f32))));
+}
+
+test "math.round64.special" {
+    assert(round64(0.0) == 0.0);
+    assert(round64(-0.0) == -0.0);
+    assert(math.isPositiveInf(round64(math.inf(f64))));
+    assert(math.isNegativeInf(round64(-math.inf(f64))));
+    assert(math.isNan(round64(math.nan(f64))));
+}

std/math/scalbn.zig

@@ -0,0 +1,87 @@
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const scalbn = scalbn_workaround;
+
+pub fn scalbn_workaround(x: var, n: i32) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(scalbn32, x, n),
+        f64 => @inlineCall(scalbn64, x, n),
+        else => @compileError("scalbn not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn scalbn32(x: f32, n_: i32) -> f32 {
+    var y = x;
+    var n = n_;
+
+    if (n > 127) {
+        y *= 0x1.0p127;
+        n -= 127;
+        if (n > 1023) {
+            y *= 0x1.0p127;
+            n -= 127;
+            if (n > 127) {
+                n = 127;
+            }
+        }
+    } else if (n < -126) {
+        y *= 0x1.0p-126 * 0x1.0p24;
+        n += 126 - 24;
+        if (n < -126) {
+            y *= 0x1.0p-126 * 0x1.0p24;
+            n += 126 - 24;
+            if (n < -126) {
+                n = -126;
+            }
+        }
+    }
+
+    const u = u32(n +% 0x7F) << 23;
+    y * @bitCast(f32, u)
+}
+
+fn scalbn64(x: f64, n_: i32) -> f64 {
+    var y = x;
+    var n = n_;
+
+    if (n > 1023) {
+        y *= 0x1.0p1022 * 2.0;
+        n -= 1023;
+        if (n > 1023) {
+            y *= 0x1.0p1022 * 2.0;
+            n -= 1023;
+            if (n > 1023) {
+                n = 1023;
+            }
+        }
+    } else if (n < -1022) {
+        y *= 0x1.0p-1022 * 0x1.0p53;
+        n += 1022 - 53;
+        if (n < -1022) {
+            y *= 0x1.0p-1022 * 0x1.0p53;
+            n += 1022 - 53;
+            if (n < -1022) {
+                n = -1022;
+            }
+        }
+    }
+
+    const u = u64(n +% 0x3FF) << 52;
+    y * @bitCast(f64, u)
+}
+
+test "math.scalbn" {
+    assert(scalbn(f32(1.5), 4) == scalbn32(1.5, 4));
+    assert(scalbn(f64(1.5), 4) == scalbn64(1.5, 4));
+}
+
+test "math.scalbn32" {
+    assert(scalbn32(1.5, 4) == 24.0);
+}
+
+test "math.scalbn64" {
+    assert(scalbn64(1.5, 4) == 24.0);
+}

std/math/signbit.zig

@@ -0,0 +1,39 @@
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const signbit = signbit_workaround;
+
+pub fn signbit_workaround(x: var) -> bool {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(signbit32, x),
+        f64 => @inlineCall(signbit64, x),
+        else => @compileError("signbit not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn signbit32(x: f32) -> bool {
+    const bits = @bitCast(u32, x);
+    bits >> 31 != 0
+}
+
+fn signbit64(x: f64) -> bool {
+    const bits = @bitCast(u64, x);
+    bits >> 63 != 0
+}
+
+test "math.signbit" {
+    assert(signbit(f32(4.0)) == signbit32(4.0));
+    assert(signbit(f64(4.0)) == signbit64(4.0));
+}
+
+test "math.signbit32" {
+    assert(!signbit32(4.0));
+    assert(signbit32(-3.0));
+}
+
+test "math.signbit64" {
+    assert(!signbit64(4.0));
+    assert(signbit64(-3.0));
+}

std/math/sin.zig

@@ -0,0 +1,188 @@
+// Special Cases:
+//
+// - sin(+-0)   = +-0
+// - sin(+-inf) = nan
+// - sin(nan)   = nan
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const sin = sin_workaround;
+
+pub fn sin_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(sin32, x),
+        f64 => @inlineCall(sin64, x),
+        else => @compileError("sin not implemented for " ++ @typeName(T)),
+    }
+}
+
+// sin polynomial coefficients
+const S0 =  1.58962301576546568060E-10;
+const S1 = -2.50507477628578072866E-8;
+const S2 =  2.75573136213857245213E-6;
+const S3 = -1.98412698295895385996E-4;
+const S4 =  8.33333333332211858878E-3;
+const S5 = -1.66666666666666307295E-1;
+
+// cos polynomial coeffiecients
+const C0 = -1.13585365213876817300E-11;
+const C1 =  2.08757008419747316778E-9;
+const C2 = -2.75573141792967388112E-7;
+const C3 =  2.48015872888517045348E-5;
+const C4 = -1.38888888888730564116E-3;
+const C5 =  4.16666666666665929218E-2;
+
+// NOTE: This is taken from the go stdlib. The musl implementation is much more complex.
+//
+// This may have slight differences on some edge cases and may need to replaced if so.
+fn sin32(x_: f32) -> f32 {
+    @setFloatMode(this, @import("builtin").FloatMode.Strict);
+
+    const pi4a = 7.85398125648498535156e-1;
+    const pi4b = 3.77489470793079817668E-8;
+    const pi4c = 2.69515142907905952645E-15;
+    const m4pi = 1.273239544735162542821171882678754627704620361328125;
+
+    var x = x_;
+    if (x == 0 or math.isNan(x)) {
+        return x;
+    }
+    if (math.isInf(x)) {
+        return math.nan(f32);
+    }
+
+    var sign = false;
+    if (x < 0) {
+        x = -x;
+        sign = true;
+    }
+
+    var y = math.floor(x * m4pi);
+    var j = i64(y);
+
+    if (j & 1 == 1) {
+        j += 1;
+        y += 1;
+    }
+
+    j &= 7;
+    if (j > 3) {
+        j -= 4;
+        sign = !sign;
+    }
+
+    const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
+    const w = z * z;
+
+    const r = {
+        if (j == 1 or j == 2) {
+            1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))))
+        } else {
+            z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))))
+        }
+    };
+
+    if (sign) {
+        -r
+    } else {
+        r
+    }
+}
+
+fn sin64(x_: f64) -> f64 {
+    const pi4a = 7.85398125648498535156e-1;
+    const pi4b = 3.77489470793079817668E-8;
+    const pi4c = 2.69515142907905952645E-15;
+    const m4pi = 1.273239544735162542821171882678754627704620361328125;
+
+    var x = x_;
+    if (x == 0 or math.isNan(x)) {
+        return x;
+    }
+    if (math.isInf(x)) {
+        return math.nan(f64);
+    }
+
+    var sign = false;
+    if (x < 0) {
+        x = -x;
+        sign = true;
+    }
+
+    var y = math.floor(x * m4pi);
+    var j = i64(y);
+
+    if (j & 1 == 1) {
+        j += 1;
+        y += 1;
+    }
+
+    j &= 7;
+    if (j > 3) {
+        j -= 4;
+        sign = !sign;
+    }
+
+    const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
+    const w = z * z;
+
+    const r = {
+        if (j == 1 or j == 2) {
+            1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))))
+        } else {
+            z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))))
+        }
+    };
+
+    if (sign) {
+        -r
+    } else {
+        r
+    }
+}
+
+test "math.sin" {
+    assert(sin(f32(0.0)) == sin32(0.0));
+    assert(sin(f64(0.0)) == sin64(0.0));
+}
+
+test "math.sin32" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, sin32(0.0), 0.0, epsilon));
+    assert(math.approxEq(f32, sin32(0.2), 0.198669, epsilon));
+    assert(math.approxEq(f32, sin32(0.8923), 0.778517, epsilon));
+    assert(math.approxEq(f32, sin32(1.5), 0.997495, epsilon));
+    assert(math.approxEq(f32, sin32(37.45), -0.246544, epsilon));
+    assert(math.approxEq(f32, sin32(89.123), 0.916166, epsilon));
+}
+
+test "math.sin64" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f64, sin64(0.0), 0.0, epsilon));
+    assert(math.approxEq(f64, sin64(0.2), 0.198669, epsilon));
+    assert(math.approxEq(f64, sin64(0.8923), 0.778517, epsilon));
+    assert(math.approxEq(f64, sin64(1.5), 0.997495, epsilon));
+    assert(math.approxEq(f64, sin64(37.45), -0.246543, epsilon));
+    assert(math.approxEq(f64, sin64(89.123), 0.916166, epsilon));
+}
+
+test "math.sin32.special" {
+    assert(sin32(0.0) == 0.0);
+    assert(sin32(-0.0) == -0.0);
+    assert(math.isNan(sin32(math.inf(f32))));
+    assert(math.isNan(sin32(-math.inf(f32))));
+    assert(math.isNan(sin32(math.nan(f32))));
+}
+
+test "math.sin64.special" {
+    assert(sin64(0.0) == 0.0);
+    assert(sin64(-0.0) == -0.0);
+    assert(math.isNan(sin64(math.inf(f64))));
+    assert(math.isNan(sin64(-math.inf(f64))));
+    assert(math.isNan(sin64(math.nan(f64))));
+}

std/math/sinh.zig

@@ -0,0 +1,128 @@
+// Special Cases:
+//
+// - sinh(+-0)   = +-0
+// - sinh(+-inf) = +-inf
+// - sinh(nan)   = nan
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+const expo2 = @import("expo2.zig").expo2;
+
+// TODO issue #393
+pub const sinh = sinh_workaround;
+
+pub fn sinh_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(sinh32, x),
+        f64 => @inlineCall(sinh64, x),
+        else => @compileError("sinh not implemented for " ++ @typeName(T)),
+    }
+}
+
+// sinh(x) = (exp(x) - 1 / exp(x)) / 2
+//         = (exp(x) - 1 + (exp(x) - 1) / exp(x)) / 2
+//         = x + x^3 / 6 + o(x^5)
+fn sinh32(x: f32) -> f32 {
+    const u = @bitCast(u32, x);
+    const ux = u & 0x7FFFFFFF;
+    const ax = @bitCast(f32, ux);
+
+    if (x == 0.0 or math.isNan(x)) {
+        return x;
+    }
+
+    var h: f32 = 0.5;
+    if (u >> 31 != 0) {
+        h = -h;
+    }
+
+    // |x| < log(FLT_MAX)
+    if (ux < 0x42B17217) {
+        const t = math.expm1(ax);
+        if (ux < 0x3F800000) {
+            if (ux < 0x3F800000 - (12 << 23)) {
+                return x;
+            } else {
+                return h * (2 * t - t * t / (t + 1));
+            }
+        }
+        return h * (t + t / (t + 1));
+    }
+
+    // |x| > log(FLT_MAX) or nan
+    2 * h * expo2(ax)
+}
+
+fn sinh64(x: f64) -> f64 {
+    @setFloatMode(this, @import("builtin").FloatMode.Strict);
+
+    const u = @bitCast(u64, x);
+    const w = u32(u >> 32);
+    const ax = @bitCast(f64, u & (@maxValue(u64) >> 1));
+
+    if (x == 0.0 or math.isNan(x)) {
+        return x;
+    }
+
+    var h: f32 = 0.5;
+    if (u >> 63 != 0) {
+        h = -h;
+    }
+
+    // |x| < log(FLT_MAX)
+    if (w < 0x40862E42) {
+        const t = math.expm1(ax);
+        if (w < 0x3FF00000) {
+            if (w < 0x3FF00000 - (26 << 20)) {
+                return x;
+            } else {
+                return h * (2 * t - t * t / (t + 1));
+            }
+        }
+        // NOTE: |x| > log(0x1p26) + eps could be h * exp(x)
+        return h * (t + t / (t + 1));
+    }
+
+    // |x| > log(DBL_MAX) or nan
+    2 * h * expo2(ax)
+}
+
+test "math.sinh" {
+    assert(sinh(f32(1.5)) == sinh32(1.5));
+    assert(sinh(f64(1.5)) == sinh64(1.5));
+}
+
+test "math.sinh32" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, sinh32(0.0), 0.0, epsilon));
+    assert(math.approxEq(f32, sinh32(0.2), 0.201336, epsilon));
+    assert(math.approxEq(f32, sinh32(0.8923), 1.015512, epsilon));
+    assert(math.approxEq(f32, sinh32(1.5), 2.129279, epsilon));
+}
+
+test "math.sinh64" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f64, sinh64(0.0), 0.0, epsilon));
+    assert(math.approxEq(f64, sinh64(0.2), 0.201336, epsilon));
+    assert(math.approxEq(f64, sinh64(0.8923), 1.015512, epsilon));
+    assert(math.approxEq(f64, sinh64(1.5), 2.129279, epsilon));
+}
+
+test "math.sinh32.special" {
+    assert(sinh32(0.0) == 0.0);
+    assert(sinh32(-0.0) == -0.0);
+    assert(math.isPositiveInf(sinh32(math.inf(f32))));
+    assert(math.isNegativeInf(sinh32(-math.inf(f32))));
+    assert(math.isNan(sinh32(math.nan(f32))));
+}
+
+test "math.sinh64.special" {
+    assert(sinh64(0.0) == 0.0);
+    assert(sinh64(-0.0) == -0.0);
+    assert(math.isPositiveInf(sinh64(math.inf(f64))));
+    assert(math.isNegativeInf(sinh64(-math.inf(f64))));
+    assert(math.isNan(sinh64(math.nan(f64))));
+}

std/math/sqrt.zig

@@ -0,0 +1,279 @@
+// Special Cases:
+//
+// - sqrt(+inf)  = +inf
+// - sqrt(+-0)   = +-0
+// - sqrt(x)     = nan if x < 0
+// - sqrt(nan)   = nan
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+// TODO issue #393
+pub const sqrt = sqrt_workaround;
+
+pub fn sqrt_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(sqrt32, x),
+        f64 => @inlineCall(sqrt64, 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;
+    @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 >> (32 - i);
+        ix1 <<= 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;
+    @bitCast(f64, uz)
+}
+
+test "math.sqrt" {
+    assert(sqrt(f32(0.0)) == sqrt32(0.0));
+    assert(sqrt(f64(0.0)) == sqrt64(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));
+}
+
+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));
+}
+
+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))));
+}
+
+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))));
+}

std/math/tan.zig

@@ -0,0 +1,174 @@
+// Special Cases:
+//
+// - tan(+-0)   = +-0
+// - tan(+-inf) = nan
+// - tan(nan)   = nan
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+pub const tan = tan_workaround;
+
+pub fn tan_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(tan32, x),
+        f64 => @inlineCall(tan64, x),
+        else => @compileError("tan not implemented for " ++ @typeName(T)),
+    }
+}
+
+const Tp0 = -1.30936939181383777646E4;
+const Tp1 =  1.15351664838587416140E6;
+const Tp2 = -1.79565251976484877988E7;
+
+const Tq1 =  1.36812963470692954678E4;
+const Tq2 = -1.32089234440210967447E6;
+const Tq3 =  2.50083801823357915839E7;
+const Tq4 = -5.38695755929454629881E7;
+
+// NOTE: This is taken from the go stdlib. The musl implementation is much more complex.
+//
+// This may have slight differences on some edge cases and may need to replaced if so.
+fn tan32(x_: f32) -> f32 {
+    @setFloatMode(this, @import("builtin").FloatMode.Strict);
+
+    const pi4a = 7.85398125648498535156e-1;
+    const pi4b = 3.77489470793079817668E-8;
+    const pi4c = 2.69515142907905952645E-15;
+    const m4pi = 1.273239544735162542821171882678754627704620361328125;
+
+    var x = x_;
+    if (x == 0 or math.isNan(x)) {
+        return x;
+    }
+    if (math.isInf(x)) {
+        return math.nan(f32);
+    }
+
+    var sign = false;
+    if (x < 0) {
+        x = -x;
+        sign = true;
+    }
+
+    var y = math.floor(x * m4pi);
+    var j = i64(y);
+
+    if (j & 1 == 1) {
+        j += 1;
+        y += 1;
+    }
+
+    const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
+    const w = z * z;
+
+    var r = {
+        if (w > 1e-14) {
+            z + z * (w * ((Tp0 * w + Tp1) * w + Tp2) / ((((w + Tq1) * w + Tq2) * w + Tq3) * w + Tq4))
+        } else {
+            z
+        }
+    };
+
+    if (j & 2 == 2) {
+        r = -1 / r;
+    }
+    if (sign) {
+        r = -r;
+    }
+
+    r
+}
+
+fn tan64(x_: f64) -> f64 {
+    const pi4a = 7.85398125648498535156e-1;
+    const pi4b = 3.77489470793079817668E-8;
+    const pi4c = 2.69515142907905952645E-15;
+    const m4pi = 1.273239544735162542821171882678754627704620361328125;
+
+    var x = x_;
+    if (x == 0 or math.isNan(x)) {
+        return x;
+    }
+    if (math.isInf(x)) {
+        return math.nan(f64);
+    }
+
+    var sign = false;
+    if (x < 0) {
+        x = -x;
+        sign = true;
+    }
+
+    var y = math.floor(x * m4pi);
+    var j = i64(y);
+
+    if (j & 1 == 1) {
+        j += 1;
+        y += 1;
+    }
+
+    const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
+    const w = z * z;
+
+    var r = {
+        if (w > 1e-14) {
+            z + z * (w * ((Tp0 * w + Tp1) * w + Tp2) / ((((w + Tq1) * w + Tq2) * w + Tq3) * w + Tq4))
+        } else {
+            z
+        }
+    };
+
+    if (j & 2 == 2) {
+        r = -1 / r;
+    }
+    if (sign) {
+        r = -r;
+    }
+
+    r
+}
+
+test "math.tan" {
+    assert(tan(f32(0.0)) == tan32(0.0));
+    assert(tan(f64(0.0)) == tan64(0.0));
+}
+
+test "math.tan32" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, tan32(0.0), 0.0, epsilon));
+    assert(math.approxEq(f32, tan32(0.2), 0.202710, epsilon));
+    assert(math.approxEq(f32, tan32(0.8923), 1.240422, epsilon));
+    assert(math.approxEq(f32, tan32(1.5), 14.101420, epsilon));
+    assert(math.approxEq(f32, tan32(37.45), -0.254397, epsilon));
+    assert(math.approxEq(f32, tan32(89.123), 2.285852, epsilon));
+}
+
+test "math.tan64" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f64, tan64(0.0), 0.0, epsilon));
+    assert(math.approxEq(f64, tan64(0.2), 0.202710, epsilon));
+    assert(math.approxEq(f64, tan64(0.8923), 1.240422, epsilon));
+    assert(math.approxEq(f64, tan64(1.5), 14.101420, epsilon));
+    assert(math.approxEq(f64, tan64(37.45), -0.254397, epsilon));
+    assert(math.approxEq(f64, tan64(89.123), 2.2858376, epsilon));
+}
+
+test "math.tan32.special" {
+    assert(tan32(0.0) == 0.0);
+    assert(tan32(-0.0) == -0.0);
+    assert(math.isNan(tan32(math.inf(f32))));
+    assert(math.isNan(tan32(-math.inf(f32))));
+    assert(math.isNan(tan32(math.nan(f32))));
+}
+
+test "math.tan64.special" {
+    assert(tan64(0.0) == 0.0);
+    assert(tan64(-0.0) == -0.0);
+    assert(math.isNan(tan64(math.inf(f64))));
+    assert(math.isNan(tan64(-math.inf(f64))));
+    assert(math.isNan(tan64(math.nan(f64))));
+}

std/math/tanh.zig

@@ -0,0 +1,154 @@
+// Special Cases:
+//
+// - sinh(+-0)   = +-0
+// - sinh(+-inf) = +-1
+// - sinh(nan)   = nan
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+const expo2 = @import("expo2.zig").expo2;
+
+// TODO issue #393
+pub const tanh = tanh_workaround;
+
+pub fn tanh_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(tanh32, x),
+        f64 => @inlineCall(tanh64, x),
+        else => @compileError("tanh not implemented for " ++ @typeName(T)),
+    }
+}
+
+// tanh(x) = (exp(x) - exp(-x)) / (exp(x) + exp(-x))
+//         = (exp(2x) - 1) / (exp(2x) - 1 + 2)
+//         = (1 - exp(-2x)) / (exp(-2x) - 1 + 2)
+fn tanh32(x: f32) -> f32 {
+    const u = @bitCast(u32, x);
+    const ux = u & 0x7FFFFFFF;
+    const ax = @bitCast(f32, ux);
+
+    var t: f32 = undefined;
+
+    if (x == 0.0 or math.isNan(x)) {
+        return x;
+    }
+
+    // |x| < log(3) / 2 ~= 0.5493 or nan
+    if (ux > 0x3F0C9F54) {
+        // |x| > 10
+        if (ux > 0x41200000) {
+            t = 1.0;
+        } else {
+            t = math.expm1(2 * x);
+            t = 1 - 2 / (t + 2);
+        }
+    }
+    // |x| > log(5 / 3) / 2 ~= 0.2554
+    else if (ux > 0x3E82C578) {
+        t = math.expm1(2 * x);
+        t = t / (t + 2);
+    }
+    // |x| >= 0x1.0p-126
+    else if (ux >= 0x00800000) {
+        t = math.expm1(-2 * x);
+        t = -t / (t + 2);
+    }
+    // |x| is subnormal
+    else {
+        math.forceEval(x * x);
+        t = x;
+    }
+
+    if (u >> 31 != 0) {
+        -t
+    } else {
+        t
+    }
+}
+
+fn tanh64(x: f64) -> f64 {
+    const u = @bitCast(u64, x);
+    const w = u32(u >> 32);
+    const ax = @bitCast(f64, u & (@maxValue(u64) >> 1));
+
+    var t: f64 = undefined;
+
+    // TODO: Shouldn't need these checks.
+    if (x == 0.0 or math.isNan(x)) {
+        return x;
+    }
+
+    // |x| < log(3) / 2 ~= 0.5493 or nan
+    if (w > 0x3FE193EA) {
+        // |x| > 20 or nan
+        if (w > 0x40340000) {
+            t = 1.0;
+        } else {
+            t = math.expm1(2 * x);
+            t = 1 - 2 / (t + 2);
+        }
+    }
+    // |x| > log(5 / 3) / 2 ~= 0.2554
+    else if (w > 0x3FD058AE) {
+        t = math.expm1(2 * x);
+        t = t / (t + 2);
+    }
+    // |x| >= 0x1.0p-1022
+    else if (w >= 0x00100000) {
+        t = math.expm1(-2 * x);
+        t = -t / (t + 2);
+    }
+    // |x| is subnormal
+    else {
+        math.forceEval(f32(x));
+        t = x;
+    }
+
+    if (u >> 63 != 0) {
+        -t
+    } else {
+        t
+    }
+}
+
+test "math.tanh" {
+    assert(tanh(f32(1.5)) == tanh32(1.5));
+    assert(tanh(f64(1.5)) == tanh64(1.5));
+}
+
+test "math.tanh32" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f32, tanh32(0.0), 0.0, epsilon));
+    assert(math.approxEq(f32, tanh32(0.2), 0.197375, epsilon));
+    assert(math.approxEq(f32, tanh32(0.8923), 0.712528, epsilon));
+    assert(math.approxEq(f32, tanh32(1.5), 0.905148, epsilon));
+    assert(math.approxEq(f32, tanh32(37.45), 1.0, epsilon));
+}
+
+test "math.tanh64" {
+    const epsilon = 0.000001;
+
+    assert(math.approxEq(f64, tanh64(0.0), 0.0, epsilon));
+    assert(math.approxEq(f64, tanh64(0.2), 0.197375, epsilon));
+    assert(math.approxEq(f64, tanh64(0.8923), 0.712528, epsilon));
+    assert(math.approxEq(f64, tanh64(1.5), 0.905148, epsilon));
+    assert(math.approxEq(f64, tanh64(37.45), 1.0, epsilon));
+}
+
+test "math.tanh32.special" {
+    assert(tanh32(0.0) == 0.0);
+    assert(tanh32(-0.0) == -0.0);
+    assert(tanh32(math.inf(f32)) == 1.0);
+    assert(tanh32(-math.inf(f32)) == -1.0);
+    assert(math.isNan(tanh32(math.nan(f32))));
+}
+
+test "math.tanh64.special" {
+    assert(tanh64(0.0) == 0.0);
+    assert(tanh64(-0.0) == -0.0);
+    assert(tanh64(math.inf(f64)) == 1.0);
+    assert(tanh64(-math.inf(f64)) == -1.0);
+    assert(math.isNan(tanh64(math.nan(f64))));
+}

std/math/trunc.zig

@@ -0,0 +1,94 @@
+// Special Cases:
+//
+// - trunc(+-0)   = +-0
+// - trunc(+-inf) = +-inf
+// - trunc(nan)   = nan
+
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+pub const trunc = trunc_workaround;
+
+pub fn trunc_workaround(x: var) -> @typeOf(x) {
+    const T = @typeOf(x);
+    switch (T) {
+        f32 => @inlineCall(trunc32, x),
+        f64 => @inlineCall(trunc64, x),
+        else => @compileError("trunc not implemented for " ++ @typeName(T)),
+    }
+}
+
+fn trunc32(x: f32) -> f32 {
+    const u = @bitCast(u32, x);
+    var e = i32(((u >> 23) & 0xFF)) - 0x7F + 9;
+    var m: u32 = undefined;
+
+    if (e >= 23 + 9) {
+        return x;
+    }
+    if (e < 9) {
+        e = 1;
+    }
+
+    m = @maxValue(u32) >> u32(e);
+    if (u & m == 0) {
+        x
+    } else {
+        math.forceEval(x + 0x1p120);
+        @bitCast(f32, u & ~m)
+    }
+}
+
+fn trunc64(x: f64) -> f64 {
+    const u = @bitCast(u64, x);
+    var e = i32(((u >> 52) & 0x7FF)) - 0x3FF + 12;
+    var m: u64 = undefined;
+
+    if (e >= 52 + 12) {
+        return x;
+    }
+    if (e < 12) {
+        e = 1;
+    }
+
+    m = @maxValue(u64) >> u64(e);
+    if (u & m == 0) {
+        x
+    } else {
+        math.forceEval(x + 0x1p120);
+        @bitCast(f64, u & ~m)
+    }
+}
+
+test "math.trunc" {
+    assert(trunc(f32(1.3)) == trunc32(1.3));
+    assert(trunc(f64(1.3)) == trunc64(1.3));
+}
+
+test "math.trunc32" {
+    assert(trunc32(1.3) == 1.0);
+    assert(trunc32(-1.3) == -1.0);
+    assert(trunc32(0.2) == 0.0);
+}
+
+test "math.trunc64" {
+    assert(trunc64(1.3) == 1.0);
+    assert(trunc64(-1.3) == -1.0);
+    assert(trunc64(0.2) == 0.0);
+}
+
+test "math.trunc32.special" {
+    assert(trunc32(0.0) == 0.0); // 0x3F800000
+    assert(trunc32(-0.0) == -0.0);
+    assert(math.isPositiveInf(trunc32(math.inf(f32))));
+    assert(math.isNegativeInf(trunc32(-math.inf(f32))));
+    assert(math.isNan(trunc32(math.nan(f32))));
+}
+
+test "math.trunc64.special" {
+    assert(trunc64(0.0) == 0.0);
+    assert(trunc64(-0.0) == -0.0);
+    assert(math.isPositiveInf(trunc64(math.inf(f64))));
+    assert(math.isNegativeInf(trunc64(-math.inf(f64))));
+    assert(math.isNan(trunc64(math.nan(f64))));
+}

std/special/builtin.zig

@@ -60,7 +60,7 @@ fn generic_fmod(comptime T: type, x: T, y: T) -> T {
     if (ex == 0) {
         i = ux <<% exp_bits;
         while (i >> bits_minus_1 == 0) : ({ex -= 1; i <<%= 1}) {}
-        ux <<%= twosComplementCast(uint, -ex + 1);
+        ux <<%= @bitCast(u32, -ex + 1);
     } else {
         ux &= @maxValue(uint) >> exp_bits;
         ux |= 1 <<% digits;
@@ -68,7 +68,7 @@ fn generic_fmod(comptime T: type, x: T, y: T) -> T {
     if (ey == 0) {
         i = uy <<% exp_bits;
         while (i >> bits_minus_1 == 0) : ({ey -= 1; i <<%= 1}) {}
-        uy <<= twosComplementCast(uint, -ey + 1);
+        uy <<= @bitCast(u32, -ey + 1);
     } else {
         uy &= @maxValue(uint) >> exp_bits;
         uy |= 1 <<% digits;
@@ -95,9 +95,9 @@ fn generic_fmod(comptime T: type, x: T, y: T) -> T {
     // scale result up
     if (ex > 0) {
         ux -%= 1 <<% digits;
-        ux |= twosComplementCast(uint, ex) <<% digits;
+        ux |= @bitCast(u32, ex) <<% digits;
     } else {
-        ux >>= twosComplementCast(uint, -ex + 1);
+        ux >>= @bitCast(u32, -ex + 1);
     }
     if (T == f32) {
         ux |= sx;
@@ -116,8 +116,3 @@ fn isNan(comptime T: type, bits: T) -> bool {
         unreachable;
     }
 }
-
-// TODO this should be a builtin function and it shouldn't do a ptr cast
-fn twosComplementCast(comptime T: type, src: var) -> T {
-    return *@ptrCast(&const @IntType(T.is_signed, @typeOf(src).bit_count), &src);
-}