Commit f94964cd05
std/math/cos.zig
@@ -1,18 +1,23 @@
-// Special Cases:
+// Ported from go, which is licensed under a BSD-3 license.
+// https://golang.org/LICENSE
//
-// - cos(+-inf) = nan
-// - cos(nan) = nan
+// https://golang.org/src/math/sin.go
const builtin = @import("builtin");
const std = @import("../std.zig");
const math = std.math;
const expect = std.testing.expect;
+/// Returns the cosine of the radian value x.
+///
+/// Special Cases:
+/// - cos(+-inf) = nan
+/// - cos(nan) = nan
pub fn cos(x: var) @typeOf(x) {
const T = @typeOf(x);
return switch (T) {
- f32 => cos32(x),
- f64 => cos64(x),
+ f32 => cos_(f32, x),
+ f64 => cos_(f64, x),
else => @compileError("cos not implemented for " ++ @typeName(T)),
};
}
@@ -33,78 +38,24 @@ 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 {
- 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);
- }
+const pi4a = 7.85398125648498535156e-1;
+const pi4b = 3.77489470793079817668E-8;
+const pi4c = 2.69515142907905952645E-15;
+const m4pi = 1.273239544735162542821171882678754627704620361328125;
- var sign = false;
- if (x < 0) {
- x = -x;
- }
-
- var y = math.floor(x * m4pi);
- var j = @floatToInt(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 = r: {
- if (j == 1 or j == 2) {
- break :r z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))));
- } else {
- break :r 1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))));
- }
- };
-
- if (sign) {
- return -r;
- } else {
- return 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;
+fn cos_(comptime T: type, x_: T) T {
+ const I = @IntType(true, T.bit_count);
var x = x_;
if (math.isNan(x) or math.isInf(x)) {
- return math.nan(f64);
+ return math.nan(f32);
}
var sign = false;
- if (x < 0) {
- x = -x;
- }
+ x = math.fabs(x);
var y = math.floor(x * m4pi);
- var j = @floatToInt(i64, y);
+ var j = @floatToInt(I, y);
if (j & 1 == 1) {
j += 1;
@@ -123,56 +74,51 @@ fn cos64(x_: f64) f64 {
const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
const w = z * z;
- const r = r: {
- if (j == 1 or j == 2) {
- break :r z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))));
- } else {
- break :r 1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))));
- }
- };
+ 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) {
- return -r;
- } else {
- return r;
- }
+ return if (sign) -r else r;
}
test "math.cos" {
- expect(cos(f32(0.0)) == cos32(0.0));
- expect(cos(f64(0.0)) == cos64(0.0));
+ expect(cos(f32(0.0)) == cos_(f32, 0.0));
+ expect(cos(f64(0.0)) == cos_(f64, 0.0));
}
test "math.cos32" {
const epsilon = 0.000001;
- expect(math.approxEq(f32, cos32(0.0), 1.0, epsilon));
- expect(math.approxEq(f32, cos32(0.2), 0.980067, epsilon));
- expect(math.approxEq(f32, cos32(0.8923), 0.627623, epsilon));
- expect(math.approxEq(f32, cos32(1.5), 0.070737, epsilon));
- expect(math.approxEq(f32, cos32(37.45), 0.969132, epsilon));
- expect(math.approxEq(f32, cos32(89.123), 0.400798, epsilon));
+ expect(math.approxEq(f32, cos_(f32, 0.0), 1.0, epsilon));
+ expect(math.approxEq(f32, cos_(f32, 0.2), 0.980067, epsilon));
+ expect(math.approxEq(f32, cos_(f32, 0.8923), 0.627623, epsilon));
+ expect(math.approxEq(f32, cos_(f32, 1.5), 0.070737, epsilon));
+ expect(math.approxEq(f32, cos_(f32, -1.5), 0.070737, epsilon));
+ expect(math.approxEq(f32, cos_(f32, 37.45), 0.969132, epsilon));
+ expect(math.approxEq(f32, cos_(f32, 89.123), 0.400798, epsilon));
}
test "math.cos64" {
const epsilon = 0.000001;
- expect(math.approxEq(f64, cos64(0.0), 1.0, epsilon));
- expect(math.approxEq(f64, cos64(0.2), 0.980067, epsilon));
- expect(math.approxEq(f64, cos64(0.8923), 0.627623, epsilon));
- expect(math.approxEq(f64, cos64(1.5), 0.070737, epsilon));
- expect(math.approxEq(f64, cos64(37.45), 0.969132, epsilon));
- expect(math.approxEq(f64, cos64(89.123), 0.40080, epsilon));
+ expect(math.approxEq(f64, cos_(f64, 0.0), 1.0, epsilon));
+ expect(math.approxEq(f64, cos_(f64, 0.2), 0.980067, epsilon));
+ expect(math.approxEq(f64, cos_(f64, 0.8923), 0.627623, epsilon));
+ expect(math.approxEq(f64, cos_(f64, 1.5), 0.070737, epsilon));
+ expect(math.approxEq(f64, cos_(f64, -1.5), 0.070737, epsilon));
+ expect(math.approxEq(f64, cos_(f64, 37.45), 0.969132, epsilon));
+ expect(math.approxEq(f64, cos_(f64, 89.123), 0.40080, epsilon));
}
test "math.cos32.special" {
- expect(math.isNan(cos32(math.inf(f32))));
- expect(math.isNan(cos32(-math.inf(f32))));
- expect(math.isNan(cos32(math.nan(f32))));
+ expect(math.isNan(cos_(f32, math.inf(f32))));
+ expect(math.isNan(cos_(f32, -math.inf(f32))));
+ expect(math.isNan(cos_(f32, math.nan(f32))));
}
test "math.cos64.special" {
- expect(math.isNan(cos64(math.inf(f64))));
- expect(math.isNan(cos64(-math.inf(f64))));
- expect(math.isNan(cos64(math.nan(f64))));
+ expect(math.isNan(cos_(f64, math.inf(f64))));
+ expect(math.isNan(cos_(f64, -math.inf(f64))));
+ expect(math.isNan(cos_(f64, math.nan(f64))));
}
std/math/fma.zig
@@ -1,7 +1,14 @@
+// Ported from musl, which is licensed under the MIT license:
+// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
+//
+// https://git.musl-libc.org/cgit/musl/tree/src/math/fmaf.c
+// https://git.musl-libc.org/cgit/musl/tree/src/math/fma.c
+
const std = @import("../std.zig");
const math = std.math;
const expect = std.testing.expect;
+/// Returns x * y + z with a single rounding error.
pub fn fma(comptime T: type, x: T, y: T, z: T) T {
return switch (T) {
f32 => fma32(x, y, z),
@@ -16,7 +23,7 @@ fn fma32(x: f32, y: f32, z: f32) f32 {
const u = @bitCast(u64, xy_z);
const e = (u >> 52) & 0x7FF;
- if ((u & 0x1FFFFFFF) != 0x10000000 or e == 0x7FF or xy_z - xy == z) {
+ if ((u & 0x1FFFFFFF) != 0x10000000 or e == 0x7FF or (xy_z - xy == z and xy_z - z == xy)) {
return @floatCast(f32, xy_z);
} else {
// TODO: Handle inexact case with double-rounding
@@ -24,6 +31,7 @@ fn fma32(x: f32, y: f32, z: f32) f32 {
}
}
+// NOTE: Upstream fma.c has been rewritten completely to raise fp exceptions more accurately.
fn fma64(x: f64, y: f64, z: f64) f64 {
if (!math.isFinite(x) or !math.isFinite(y)) {
return x * y + z;
std/math/pow.zig
@@ -1,32 +1,36 @@
-// Special Cases:
+// Ported from go, which is licensed under a BSD-3 license.
+// https://golang.org/LICENSE
//
-// 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
+// https://golang.org/src/math/pow.go
const builtin = @import("builtin");
const std = @import("../std.zig");
const math = std.math;
const expect = std.testing.expect;
-// This implementation is taken from the go stlib, musl is a bit more complex.
+/// Returns x raised to the power of y (x^y).
+///
+/// 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
pub fn pow(comptime T: type, x: T, y: T) T {
if (@typeInfo(T) == builtin.TypeId.Int) {
return math.powi(T, x, y) catch unreachable;
@@ -53,15 +57,6 @@ pub fn pow(comptime T: type, x: T, y: T) T {
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
@@ -112,14 +107,16 @@ pub fn pow(comptime T: type, x: T, y: T) T {
}
}
- var ay = y;
- var flip = false;
- if (ay < 0) {
- ay = -ay;
- flip = true;
+ // special case sqrt
+ if (y == 0.5) {
+ return math.sqrt(x);
+ }
+
+ if (y == -0.5) {
+ return 1 / math.sqrt(x);
}
- const r1 = math.modf(ay);
+ const r1 = math.modf(math.fabs(y));
var yi = r1.ipart;
var yf = r1.fpart;
@@ -148,8 +145,18 @@ pub fn pow(comptime T: type, x: T, y: T) T {
var xe = r2.exponent;
var x1 = r2.significand;
- var i = @floatToInt(i32, yi);
+ var i = @floatToInt(@IntType(true, T.bit_count), yi);
while (i != 0) : (i >>= 1) {
+ const overflow_shift = math.floatExponentBits(T) + 1;
+ if (xe < -(1 << overflow_shift) or (1 << overflow_shift) < xe) {
+ // catch xe before it overflows the left shift below
+ // Since i != 0 it has at least one bit still set, so ae will accumulate xe
+ // on at least one more iteration, ae += xe is a lower bound on ae
+ // the lower bound on ae exceeds the size of a float exp
+ // so the final call to Ldexp will produce under/overflow (0/Inf)
+ ae += xe;
+ break;
+ }
if (i & 1 == 1) {
a1 *= x1;
ae += xe;
@@ -163,7 +170,7 @@ pub fn pow(comptime T: type, x: T, y: T) T {
}
// a *= a1 * 2^ae
- if (flip) {
+ if (y < 0) {
a1 = 1 / a1;
ae = -ae;
}
@@ -202,6 +209,9 @@ test "math.pow.special" {
expect(pow(f32, 45, 1.0) == 45);
expect(pow(f32, -45, 1.0) == -45);
expect(math.isNan(pow(f32, math.nan(f32), 5.0)));
+ expect(math.isPositiveInf(pow(f32, -math.inf(f32), 0.5)));
+ expect(math.isPositiveInf(pow(f32, -0, -0.5)));
+ expect(pow(f32, -0, 0.5) == 0);
expect(math.isNan(pow(f32, 5.0, math.nan(f32))));
expect(math.isPositiveInf(pow(f32, 0.0, -1.0)));
//expect(math.isNegativeInf(pow(f32, -0.0, -3.0))); TODO is this required?
@@ -232,3 +242,11 @@ test "math.pow.special" {
expect(math.isNan(pow(f32, -1.0, 1.2)));
expect(math.isNan(pow(f32, -12.4, 78.5)));
}
+
+test "math.pow.overflow" {
+ expect(math.isPositiveInf(pow(f64, 2, 1 << 32)));
+ expect(pow(f64, 2, -(1 << 32)) == 0);
+ expect(math.isNegativeInf(pow(f64, -2, (1 << 32) + 1)));
+ expect(pow(f64, 0.5, 1 << 45) == 0);
+ expect(math.isPositiveInf(pow(f64, 0.5, -(1 << 45))));
+}
std/math/sin.zig
@@ -1,19 +1,24 @@
-// Special Cases:
+// Ported from go, which is licensed under a BSD-3 license.
+// https://golang.org/LICENSE
//
-// - sin(+-0) = +-0
-// - sin(+-inf) = nan
-// - sin(nan) = nan
+// https://golang.org/src/math/sin.go
const builtin = @import("builtin");
const std = @import("../std.zig");
const math = std.math;
const expect = std.testing.expect;
+/// Returns the sine of the radian value x.
+///
+/// Special Cases:
+/// - sin(+-0) = +-0
+/// - sin(+-inf) = nan
+/// - sin(nan) = nan
pub fn sin(x: var) @typeOf(x) {
const T = @typeOf(x);
return switch (T) {
- f32 => sin32(x),
- f64 => sin64(x),
+ f32 => sin_(T, x),
+ f64 => sin_(T, x),
else => @compileError("sin not implemented for " ++ @typeName(T)),
};
}
@@ -34,83 +39,27 @@ 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 {
- 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 = @floatToInt(i64, y);
-
- if (j & 1 == 1) {
- j += 1;
- y += 1;
- }
+const pi4a = 7.85398125648498535156e-1;
+const pi4b = 3.77489470793079817668E-8;
+const pi4c = 2.69515142907905952645E-15;
+const m4pi = 1.273239544735162542821171882678754627704620361328125;
- j &= 7;
- if (j > 3) {
- j -= 4;
- sign = !sign;
- }
-
- const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
- const w = z * z;
-
- const r = r: {
- if (j == 1 or j == 2) {
- break :r 1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))));
- } else {
- break :r z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))));
- }
- };
-
- if (sign) {
- return -r;
- } else {
- return 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;
+fn sin_(comptime T: type, x_: T) T {
+ const I = @IntType(true, T.bit_count);
var x = x_;
if (x == 0 or math.isNan(x)) {
return x;
}
if (math.isInf(x)) {
- return math.nan(f64);
+ return math.nan(T);
}
- var sign = false;
- if (x < 0) {
- x = -x;
- sign = true;
- }
+ var sign = x < 0;
+ x = math.fabs(x);
var y = math.floor(x * m4pi);
- var j = @floatToInt(i64, y);
+ var j = @floatToInt(I, y);
if (j & 1 == 1) {
j += 1;
@@ -126,61 +75,56 @@ fn sin64(x_: f64) f64 {
const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
const w = z * z;
- const r = r: {
- if (j == 1 or j == 2) {
- break :r 1.0 - 0.5 * w + w * w * (C5 + w * (C4 + w * (C3 + w * (C2 + w * (C1 + w * C0)))));
- } else {
- break :r z + z * w * (S5 + w * (S4 + w * (S3 + w * (S2 + w * (S1 + w * S0)))));
- }
- };
+ 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) {
- return -r;
- } else {
- return r;
- }
+ return if (sign) -r else r;
}
test "math.sin" {
- expect(sin(f32(0.0)) == sin32(0.0));
- expect(sin(f64(0.0)) == sin64(0.0));
+ expect(sin(f32(0.0)) == sin_(f32, 0.0));
+ expect(sin(f64(0.0)) == sin_(f64, 0.0));
expect(comptime (math.sin(f64(2))) == math.sin(f64(2)));
}
test "math.sin32" {
const epsilon = 0.000001;
- expect(math.approxEq(f32, sin32(0.0), 0.0, epsilon));
- expect(math.approxEq(f32, sin32(0.2), 0.198669, epsilon));
- expect(math.approxEq(f32, sin32(0.8923), 0.778517, epsilon));
- expect(math.approxEq(f32, sin32(1.5), 0.997495, epsilon));
- expect(math.approxEq(f32, sin32(37.45), -0.246544, epsilon));
- expect(math.approxEq(f32, sin32(89.123), 0.916166, epsilon));
+ expect(math.approxEq(f32, sin_(f32, 0.0), 0.0, epsilon));
+ expect(math.approxEq(f32, sin_(f32, 0.2), 0.198669, epsilon));
+ expect(math.approxEq(f32, sin_(f32, 0.8923), 0.778517, epsilon));
+ expect(math.approxEq(f32, sin_(f32, 1.5), 0.997495, epsilon));
+ expect(math.approxEq(f32, sin_(f32, -1.5), -0.997495, epsilon));
+ expect(math.approxEq(f32, sin_(f32, 37.45), -0.246544, epsilon));
+ expect(math.approxEq(f32, sin_(f32, 89.123), 0.916166, epsilon));
}
test "math.sin64" {
const epsilon = 0.000001;
- expect(math.approxEq(f64, sin64(0.0), 0.0, epsilon));
- expect(math.approxEq(f64, sin64(0.2), 0.198669, epsilon));
- expect(math.approxEq(f64, sin64(0.8923), 0.778517, epsilon));
- expect(math.approxEq(f64, sin64(1.5), 0.997495, epsilon));
- expect(math.approxEq(f64, sin64(37.45), -0.246543, epsilon));
- expect(math.approxEq(f64, sin64(89.123), 0.916166, epsilon));
+ expect(math.approxEq(f64, sin_(f64, 0.0), 0.0, epsilon));
+ expect(math.approxEq(f64, sin_(f64, 0.2), 0.198669, epsilon));
+ expect(math.approxEq(f64, sin_(f64, 0.8923), 0.778517, epsilon));
+ expect(math.approxEq(f64, sin_(f64, 1.5), 0.997495, epsilon));
+ expect(math.approxEq(f64, sin_(f64, -1.5), -0.997495, epsilon));
+ expect(math.approxEq(f64, sin_(f64, 37.45), -0.246543, epsilon));
+ expect(math.approxEq(f64, sin_(f64, 89.123), 0.916166, epsilon));
}
test "math.sin32.special" {
- expect(sin32(0.0) == 0.0);
- expect(sin32(-0.0) == -0.0);
- expect(math.isNan(sin32(math.inf(f32))));
- expect(math.isNan(sin32(-math.inf(f32))));
- expect(math.isNan(sin32(math.nan(f32))));
+ expect(sin_(f32, 0.0) == 0.0);
+ expect(sin_(f32, -0.0) == -0.0);
+ expect(math.isNan(sin_(f32, math.inf(f32))));
+ expect(math.isNan(sin_(f32, -math.inf(f32))));
+ expect(math.isNan(sin_(f32, math.nan(f32))));
}
test "math.sin64.special" {
- expect(sin64(0.0) == 0.0);
- expect(sin64(-0.0) == -0.0);
- expect(math.isNan(sin64(math.inf(f64))));
- expect(math.isNan(sin64(-math.inf(f64))));
- expect(math.isNan(sin64(math.nan(f64))));
+ expect(sin_(f64, 0.0) == 0.0);
+ expect(sin_(f64, -0.0) == -0.0);
+ expect(math.isNan(sin_(f64, math.inf(f64))));
+ expect(math.isNan(sin_(f64, -math.inf(f64))));
+ expect(math.isNan(sin_(f64, math.nan(f64))));
}
std/math/tan.zig
@@ -1,19 +1,24 @@
-// Special Cases:
+// Ported from go, which is licensed under a BSD-3 license.
+// https://golang.org/LICENSE
//
-// - tan(+-0) = +-0
-// - tan(+-inf) = nan
-// - tan(nan) = nan
+// https://golang.org/src/math/tan.go
const builtin = @import("builtin");
const std = @import("../std.zig");
const math = std.math;
const expect = std.testing.expect;
+/// Returns the tangent of the radian value x.
+///
+/// Special Cases:
+/// - tan(+-0) = +-0
+/// - tan(+-inf) = nan
+/// - tan(nan) = nan
pub fn tan(x: var) @typeOf(x) {
const T = @typeOf(x);
return switch (T) {
- f32 => tan32(x),
- f64 => tan64(x),
+ f32 => tan_(f32, x),
+ f64 => tan_(f64, x),
else => @compileError("tan not implemented for " ++ @typeName(T)),
};
}
@@ -27,80 +32,27 @@ 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 = @floatToInt(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 = r: {
- if (w > 1e-14) {
- break :r z + z * (w * ((Tp0 * w + Tp1) * w + Tp2) / ((((w + Tq1) * w + Tq2) * w + Tq3) * w + Tq4));
- } else {
- break :r z;
- }
- };
-
- if (j & 2 == 2) {
- r = -1 / r;
- }
- if (sign) {
- r = -r;
- }
-
- return r;
-}
+const pi4a = 7.85398125648498535156e-1;
+const pi4b = 3.77489470793079817668E-8;
+const pi4c = 2.69515142907905952645E-15;
+const m4pi = 1.273239544735162542821171882678754627704620361328125;
-fn tan64(x_: f64) f64 {
- const pi4a = 7.85398125648498535156e-1;
- const pi4b = 3.77489470793079817668E-8;
- const pi4c = 2.69515142907905952645E-15;
- const m4pi = 1.273239544735162542821171882678754627704620361328125;
+fn tan_(comptime T: type, x_: T) T {
+ const I = @IntType(true, T.bit_count);
var x = x_;
if (x == 0 or math.isNan(x)) {
return x;
}
if (math.isInf(x)) {
- return math.nan(f64);
+ return math.nan(T);
}
- var sign = false;
- if (x < 0) {
- x = -x;
- sign = true;
- }
+ var sign = x < 0;
+ x = math.fabs(x);
var y = math.floor(x * m4pi);
- var j = @floatToInt(i64, y);
+ var j = @floatToInt(I, y);
if (j & 1 == 1) {
j += 1;
@@ -110,63 +62,57 @@ fn tan64(x_: f64) f64 {
const z = ((x - y * pi4a) - y * pi4b) - y * pi4c;
const w = z * z;
- var r = r: {
- if (w > 1e-14) {
- break :r z + z * (w * ((Tp0 * w + Tp1) * w + Tp2) / ((((w + Tq1) * w + Tq2) * w + Tq3) * w + Tq4));
- } else {
- break :r 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;
- }
- return r;
+ return if (sign) -r else r;
}
test "math.tan" {
- expect(tan(f32(0.0)) == tan32(0.0));
- expect(tan(f64(0.0)) == tan64(0.0));
+ expect(tan(f32(0.0)) == tan_(f32, 0.0));
+ expect(tan(f64(0.0)) == tan_(f64, 0.0));
}
test "math.tan32" {
const epsilon = 0.000001;
- expect(math.approxEq(f32, tan32(0.0), 0.0, epsilon));
- expect(math.approxEq(f32, tan32(0.2), 0.202710, epsilon));
- expect(math.approxEq(f32, tan32(0.8923), 1.240422, epsilon));
- expect(math.approxEq(f32, tan32(1.5), 14.101420, epsilon));
- expect(math.approxEq(f32, tan32(37.45), -0.254397, epsilon));
- expect(math.approxEq(f32, tan32(89.123), 2.285852, epsilon));
+ expect(math.approxEq(f32, tan_(f32, 0.0), 0.0, epsilon));
+ expect(math.approxEq(f32, tan_(f32, 0.2), 0.202710, epsilon));
+ expect(math.approxEq(f32, tan_(f32, 0.8923), 1.240422, epsilon));
+ expect(math.approxEq(f32, tan_(f32, 1.5), 14.101420, epsilon));
+ expect(math.approxEq(f32, tan_(f32, 37.45), -0.254397, epsilon));
+ expect(math.approxEq(f32, tan_(f32, 89.123), 2.285852, epsilon));
}
test "math.tan64" {
const epsilon = 0.000001;
- expect(math.approxEq(f64, tan64(0.0), 0.0, epsilon));
- expect(math.approxEq(f64, tan64(0.2), 0.202710, epsilon));
- expect(math.approxEq(f64, tan64(0.8923), 1.240422, epsilon));
- expect(math.approxEq(f64, tan64(1.5), 14.101420, epsilon));
- expect(math.approxEq(f64, tan64(37.45), -0.254397, epsilon));
- expect(math.approxEq(f64, tan64(89.123), 2.2858376, epsilon));
+ expect(math.approxEq(f64, tan_(f64, 0.0), 0.0, epsilon));
+ expect(math.approxEq(f64, tan_(f64, 0.2), 0.202710, epsilon));
+ expect(math.approxEq(f64, tan_(f64, 0.8923), 1.240422, epsilon));
+ expect(math.approxEq(f64, tan_(f64, 1.5), 14.101420, epsilon));
+ expect(math.approxEq(f64, tan_(f64, 37.45), -0.254397, epsilon));
+ expect(math.approxEq(f64, tan_(f64, 89.123), 2.2858376, epsilon));
}
test "math.tan32.special" {
- expect(tan32(0.0) == 0.0);
- expect(tan32(-0.0) == -0.0);
- expect(math.isNan(tan32(math.inf(f32))));
- expect(math.isNan(tan32(-math.inf(f32))));
- expect(math.isNan(tan32(math.nan(f32))));
+ expect(tan_(f32, 0.0) == 0.0);
+ expect(tan_(f32, -0.0) == -0.0);
+ expect(math.isNan(tan_(f32, math.inf(f32))));
+ expect(math.isNan(tan_(f32, -math.inf(f32))));
+ expect(math.isNan(tan_(f32, math.nan(f32))));
}
test "math.tan64.special" {
- expect(tan64(0.0) == 0.0);
- expect(tan64(-0.0) == -0.0);
- expect(math.isNan(tan64(math.inf(f64))));
- expect(math.isNan(tan64(-math.inf(f64))));
- expect(math.isNan(tan64(math.nan(f64))));
+ expect(tan_(f64, 0.0) == 0.0);
+ expect(tan_(f64, -0.0) == -0.0);
+ expect(math.isNan(tan_(f64, math.inf(f64))));
+ expect(math.isNan(tan_(f64, -math.inf(f64))));
+ expect(math.isNan(tan_(f64, math.nan(f64))));
}