diff options
| -rw-r--r-- | std/math/cos.zig | 146 | ||||
| -rw-r--r-- | std/math/fma.zig | 10 | ||||
| -rw-r--r-- | std/math/pow.zig | 96 | ||||
| -rw-r--r-- | std/math/sin.zig | 160 | ||||
| -rw-r--r-- | std/math/tan.zig | 154 |
5 files changed, 214 insertions, 352 deletions
diff --git a/std/math/cos.zig b/std/math/cos.zig index 9479482894..6d0cc80c31 100644 --- a/std/math/cos.zig +++ b/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)))); } diff --git a/std/math/fma.zig b/std/math/fma.zig index a317bc96de..19c306fa2a 100644 --- a/std/math/fma.zig +++ b/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; diff --git a/std/math/pow.zig b/std/math/pow.zig index 81bb2d95d5..18c9f80634 100644 --- a/std/math/pow.zig +++ b/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)))); +} diff --git a/std/math/sin.zig b/std/math/sin.zig index e25b8a292b..f21db4054e 100644 --- a/std/math/sin.zig +++ b/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)))); } diff --git a/std/math/tan.zig b/std/math/tan.zig index fc11ebdef7..e8259ee7ad 100644 --- a/std/math/tan.zig +++ b/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)))); } |