Add math library

This covers the majority of the functions as covered by the C99
specification for a math library.

Code is adapted primarily from musl libc, with the pow and standard
trigonometric functions adapted from the Go stdlib.

Changes:

 - Remove assert expose in index and import as needed.
 - Add float log function and merge with existing base 2 integer
   implementation.

See https://github.com/tiehuis/zig-fmath.
See #374.

1 files changed, 148 insertions, 0 deletions

diff --git a/std/math/tan.zig b/std/math/tan.zig
new file mode 100644
index 0000000000..a8e907a59d
--- /dev/null
+++ b/std/math/tan.zig
@@ -0,0 +1,148 @@
+const math = @import("index.zig");
+const assert = @import("../debug.zig").assert;
+
+pub fn tan(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 {
+    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 "tan" {
+    assert(tan(f32(0.0)) == tan32(0.0));
+    assert(tan(f64(0.0)) == tan64(0.0));
+}
+
+test "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 "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));
+}