diff --git a/std/math/exp.zig b/std/math/exp.zig index 53cac3facd..2efbbeccc9 100644 --- a/std/math/exp.zig +++ b/std/math/exp.zig @@ -11,7 +11,7 @@ pub fn exp(x: var) -> @typeOf(x) { } fn exp32(x_: f32) -> f32 { - const half = []const f32 { 0.5, -0.5 }; + const half = []f32 { 0.5, -0.5 }; const ln2hi = 6.9314575195e-1; const ln2lo = 1.4286067653e-6; const invln2 = 1.4426950216e+0; diff --git a/std/math/ln.zig b/std/math/ln.zig index aff84303b7..5e13be8118 100644 --- a/std/math/ln.zig +++ b/std/math/ln.zig @@ -120,12 +120,12 @@ fn lnd(x_: f64) -> f64 { s * (hfsq + R) + dk * ln2_lo - hfsq + f + dk * ln2_hi } -test "log" { +test "math.ln" { assert(ln(f32(0.2)) == lnf(0.2)); assert(ln(f64(0.2)) == lnd(0.2)); } -test "logf" { +test "math.ln32" { const epsilon = 0.000001; assert(math.approxEq(f32, lnf(0.2), -1.609438, epsilon)); @@ -136,7 +136,7 @@ test "logf" { assert(math.approxEq(f32, lnf(123123.234375), 11.720941, epsilon)); } -test "logd" { +test "math.ln64" { const epsilon = 0.000001; assert(math.approxEq(f64, lnd(0.2), -1.609438, epsilon)); diff --git a/std/math/pow.zig b/std/math/pow.zig index ab281b67a0..69dd458afd 100644 --- a/std/math/pow.zig +++ b/std/math/pow.zig @@ -1,21 +1,12 @@ const math = @import("index.zig"); const assert = @import("../debug.zig").assert; -pub fn pow(comptime T: type, x: T, y: T) -> T { - switch (T) { - f32 => @inlineCall(pow32, x, y), - f64 => @inlineCall(pow64, x, y), - else => @compileError("pow not implemented for " ++ @typeName(T)), - } -} - -fn isOddInteger(x: f64) -> bool { - const r = math.modf(x); - r.fpart == 0.0 and i64(r.ipart) & 1 == 1 -} - // This implementation is taken from the go stlib, musl is a bit more complex. -fn pow32(x: f32, y: f32) -> f32 { +pub fn pow(comptime T: type, x: T, y: T) -> T { + 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) { @@ -25,7 +16,7 @@ fn pow32(x: f32, y: f32) -> f32 { // 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(f32); + return math.nan(T); } // pow(x, 1) = x for all x @@ -46,11 +37,11 @@ fn pow32(x: f32, y: f32) -> f32 { if (y < 0) { // pow(+-0, y) = +- 0 for y an odd integer if (isOddInteger(y)) { - return math.copysign(f32, math.inf(f32), x); + return math.copysign(T, math.inf(T), x); } // pow(+-0, y) = +inf for y an even integer else { - return math.inf(f32); + return math.inf(T); } } else { if (isOddInteger(y)) { @@ -74,13 +65,13 @@ fn pow32(x: f32, y: f32) -> f32 { // pow(x, -inf) = +inf for |x| < 1 // pow(x, +inf) = +inf for |x| > 1 else { - return math.inf(f32); + return math.inf(T); } } if (math.isInf(x)) { if (math.isNegativeInf(x)) { - return pow32(1 / x, -y); + return pow(T, 1 / x, -y); } // pow(+inf, y) = +0 for y < 0 else if (y < 0) { @@ -88,7 +79,7 @@ fn pow32(x: f32, y: f32) -> f32 { } // pow(+inf, y) = +0 for y > 0 else if (y > 0) { - return math.inf(f32); + return math.inf(T); } } @@ -104,14 +95,14 @@ fn pow32(x: f32, y: f32) -> f32 { var yf = r1.fpart; if (yf != 0 and x < 0) { - return math.nan(f32); + return math.nan(T); } - if (yi >= 1 << 31) { + if (yi >= 1 << (T.bit_count - 1)) { return math.exp(y * math.ln(x)); } // a = a1 * 2^ae - var a1: f32 = 1.0; + var a1: T = 1.0; var ae: i32 = 0; // a *= x^yf @@ -151,166 +142,26 @@ fn pow32(x: f32, y: f32) -> f32 { math.scalbn(a1, ae) } -// This implementation is taken from the go stlib, musl is a bit more complex. -fn pow64(x: f64, y: f64) -> f64 { - // 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(f64); - } - - // 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(f64, math.inf(f64), x); - } - // pow(+-0, y) = +inf for y an even integer - else { - return math.inf(f64); - } - } 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; - } - // pow(x, +inf) = +0 for |x| < 1 - // pow(x, -inf) = +0 for |x| > 1 - else if ((math.fabs(x) < 1) == math.isInf(y)) { - return 0; - } - // pow(x, -inf) = +inf for |x| < 1 - // pow(x, +inf) = +inf for |x| > 1 - else { - return math.inf(f64); - } - } - - if (math.isInf(x)) { - if (math.isInf(x)) { - return pow64(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(f64); - } - } - - 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(f64); - } - if (yi >= 1 << 63) { - return math.exp(y * math.ln(x)); - } - - // a = a1 * 2^ae - var a1: f64 = 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 = i64(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 "pow" { - assert(pow(f32, 0.2, 3.3) == pow32(0.2, 3.3)); - assert(pow(f64, 0.2, 3.3) == pow64(0.2, 3.3)); -} - -test "pow32" { +test "math.pow" { const epsilon = 0.000001; - // assert(math.approxEq(f32, pow32(0.0, 3.3), 0.0, epsilon)); // TODO: Handle div zero - assert(math.approxEq(f32, pow32(0.8923, 3.3), 0.686572, epsilon)); - assert(math.approxEq(f32, pow32(0.2, 3.3), 0.004936, epsilon)); - assert(math.approxEq(f32, pow32(1.5, 3.3), 3.811546, epsilon)); - assert(math.approxEq(f32, pow32(37.45, 3.3), 155736.703125, epsilon)); - assert(math.approxEq(f32, pow32(89.123, 3.3), 2722489.5, epsilon)); -} + // assert(math.approxEq(f32, pow(f32, 0.0, 3.3), 0.0, epsilon)); // TODO: Handle div zero + 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)); -test "pow64" { - const epsilon = 0.000001; - // assert(math.approxEq(f32, pow32(0.0, 3.3), 0.0, epsilon)); // TODO: Handle div zero - assert(math.approxEq(f64, pow64(0.8923, 3.3), 0.686572, epsilon)); - assert(math.approxEq(f64, pow64(0.2, 3.3), 0.004936, epsilon)); - assert(math.approxEq(f64, pow64(1.5, 3.3), 3.811546, epsilon)); - assert(math.approxEq(f64, pow64(37.45, 3.3), 155736.7160616, epsilon)); - assert(math.approxEq(f64, pow64(89.123, 3.3), 2722490.231436, epsilon)); + // assert(math.approxEq(f32, pow(f64, 0.0, 3.3), 0.0, epsilon)); // TODO: Handle div zero + 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)); } diff --git a/std/special/builtin.zig b/std/special/builtin.zig index a21705d82e..0c209267e4 100644 --- a/std/special/builtin.zig +++ b/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); -}