From d83836825f1d386e0872fd87bacfaa6f545783a9 Mon Sep 17 00:00:00 2001 From: Andrew Kelley Date: Tue, 19 Mar 2019 14:42:29 -0400 Subject: [PATCH] add mulXf3 to compiler-rt this adds the following functions to compiler-rt: * `__mulsf3` * `__muldf3` * `__multf3` See #1290 --- CMakeLists.txt | 1 + std/math.zig | 13 +- std/special/compiler_rt.zig | 4 + std/special/compiler_rt/mulXf3.zig | 285 ++++++++++++++++++++++++ std/special/compiler_rt/mulXf3_test.zig | 86 +++++++ 5 files changed, 388 insertions(+), 1 deletion(-) create mode 100644 std/special/compiler_rt/mulXf3.zig create mode 100644 std/special/compiler_rt/mulXf3_test.zig diff --git a/CMakeLists.txt b/CMakeLists.txt index 4383b086e..f0abc2e59 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -662,6 +662,7 @@ set(ZIG_STD_FILES "special/compiler_rt/floatuntisf.zig" "special/compiler_rt/floatuntitf.zig" "special/compiler_rt/muloti4.zig" + "special/compiler_rt/mulXf3.zig" "special/compiler_rt/multi3.zig" "special/compiler_rt/popcountdi2.zig" "special/compiler_rt/truncXfYf2.zig" diff --git a/std/math.zig b/std/math.zig index cead68577..d42d6df0e 100644 --- a/std/math.zig +++ b/std/math.zig @@ -593,7 +593,16 @@ fn testRem() void { /// Returns the absolute value of the integer parameter. /// Result is an unsigned integer. -pub fn absCast(x: var) @IntType(false, @typeOf(x).bit_count) { +pub fn absCast(x: var) t: { + if (@typeOf(x) == comptime_int) { + break :t comptime_int; + } else { + break :t @IntType(false, @typeOf(x).bit_count); + } +} { + if (@typeOf(x) == comptime_int) { + return if (x < 0) -x else x; + } const uint = @IntType(false, @typeOf(x).bit_count); if (x >= 0) return @intCast(uint, x); @@ -609,6 +618,8 @@ test "math.absCast" { testing.expect(absCast(i32(minInt(i32))) == -minInt(i32)); testing.expect(@typeOf(absCast(i32(minInt(i32)))) == u32); + + testing.expect(absCast(-999) == 999); } /// Returns the negation of the integer parameter. diff --git a/std/special/compiler_rt.zig b/std/special/compiler_rt.zig index d0f29a68d..b62bc6dd6 100644 --- a/std/special/compiler_rt.zig +++ b/std/special/compiler_rt.zig @@ -24,6 +24,10 @@ comptime { @export("__addtf3", @import("compiler_rt/addXf3.zig").__addtf3, linkage); @export("__subtf3", @import("compiler_rt/addXf3.zig").__subtf3, linkage); + @export("__mulsf3", @import("compiler_rt/mulXf3.zig").__mulsf3, linkage); + @export("__muldf3", @import("compiler_rt/mulXf3.zig").__muldf3, linkage); + @export("__multf3", @import("compiler_rt/mulXf3.zig").__multf3, linkage); + @export("__floattitf", @import("compiler_rt/floattitf.zig").__floattitf, linkage); @export("__floattidf", @import("compiler_rt/floattidf.zig").__floattidf, linkage); @export("__floattisf", @import("compiler_rt/floattisf.zig").__floattisf, linkage); diff --git a/std/special/compiler_rt/mulXf3.zig b/std/special/compiler_rt/mulXf3.zig new file mode 100644 index 000000000..400b60a23 --- /dev/null +++ b/std/special/compiler_rt/mulXf3.zig @@ -0,0 +1,285 @@ +// Ported from: +// +// https://github.com/llvm/llvm-project/blob/2ffb1b0413efa9a24eb3c49e710e36f92e2cb50b/compiler-rt/lib/builtins/fp_mul_impl.inc + +const std = @import("std"); +const builtin = @import("builtin"); +const compiler_rt = @import("../compiler_rt.zig"); + +pub extern fn __multf3(a: f128, b: f128) f128 { + return mulXf3(f128, a, b); +} +pub extern fn __muldf3(a: f64, b: f64) f64 { + return mulXf3(f64, a, b); +} +pub extern fn __mulsf3(a: f32, b: f32) f32 { + return mulXf3(f32, a, b); +} + +fn mulXf3(comptime T: type, a: T, b: T) T { + const Z = @IntType(false, T.bit_count); + + const typeWidth = T.bit_count; + const significandBits = std.math.floatMantissaBits(T); + const exponentBits = std.math.floatExponentBits(T); + + const signBit = (Z(1) << (significandBits + exponentBits)); + const maxExponent = ((1 << exponentBits) - 1); + const exponentBias = (maxExponent >> 1); + + const implicitBit = (Z(1) << significandBits); + const quietBit = implicitBit >> 1; + const significandMask = implicitBit - 1; + + const absMask = signBit - 1; + const exponentMask = absMask ^ significandMask; + const qnanRep = exponentMask | quietBit; + const infRep = @bitCast(Z, std.math.inf(T)); + + const aExponent = @truncate(u32, (@bitCast(Z, a) >> significandBits) & maxExponent); + const bExponent = @truncate(u32, (@bitCast(Z, b) >> significandBits) & maxExponent); + const productSign: Z = (@bitCast(Z, a) ^ @bitCast(Z, b)) & signBit; + + var aSignificand: Z = @bitCast(Z, a) & significandMask; + var bSignificand: Z = @bitCast(Z, b) & significandMask; + var scale: i32 = 0; + + // Detect if a or b is zero, denormal, infinity, or NaN. + if (aExponent -% 1 >= maxExponent -% 1 or bExponent -% 1 >= maxExponent -% 1) { + const aAbs: Z = @bitCast(Z, a) & absMask; + const bAbs: Z = @bitCast(Z, b) & absMask; + + // NaN * anything = qNaN + if (aAbs > infRep) return @bitCast(T, @bitCast(Z, a) | quietBit); + // anything * NaN = qNaN + if (bAbs > infRep) return @bitCast(T, @bitCast(Z, b) | quietBit); + + if (aAbs == infRep) { + // infinity * non-zero = +/- infinity + if (bAbs != 0) { + return @bitCast(T, aAbs | productSign); + } else { + // infinity * zero = NaN + return @bitCast(T, qnanRep); + } + } + + if (bAbs == infRep) { + //? non-zero * infinity = +/- infinity + if (aAbs != 0) { + return @bitCast(T, bAbs | productSign); + } else { + // zero * infinity = NaN + return @bitCast(T, qnanRep); + } + } + + // zero * anything = +/- zero + if (aAbs == 0) return @bitCast(T, productSign); + // anything * zero = +/- zero + if (bAbs == 0) return @bitCast(T, productSign); + + // one or both of a or b is denormal, the other (if applicable) is a + // normal number. Renormalize one or both of a and b, and set scale to + // include the necessary exponent adjustment. + if (aAbs < implicitBit) scale +%= normalize(T, &aSignificand); + if (bAbs < implicitBit) scale +%= normalize(T, &bSignificand); + } + + // Or in the implicit significand bit. (If we fell through from the + // denormal path it was already set by normalize( ), but setting it twice + // won't hurt anything.) + aSignificand |= implicitBit; + bSignificand |= implicitBit; + + // Get the significand of a*b. Before multiplying the significands, shift + // one of them left to left-align it in the field. Thus, the product will + // have (exponentBits + 2) integral digits, all but two of which must be + // zero. Normalizing this result is just a conditional left-shift by one + // and bumping the exponent accordingly. + var productHi: Z = undefined; + var productLo: Z = undefined; + wideMultiply(Z, aSignificand, bSignificand << exponentBits, &productHi, &productLo); + + var productExponent: i32 = @bitCast(i32, aExponent +% bExponent) -% exponentBias +% scale; + + // Normalize the significand, adjust exponent if needed. + if ((productHi & implicitBit) != 0) { + productExponent +%= 1; + } else { + productHi = (productHi << 1) | (productLo >> (typeWidth - 1)); + productLo = productLo << 1; + } + + // If we have overflowed the type, return +/- infinity. + if (productExponent >= maxExponent) return @bitCast(T, infRep | productSign); + + if (productExponent <= 0) { + // Result is denormal before rounding + // + // If the result is so small that it just underflows to zero, return + // a zero of the appropriate sign. Mathematically there is no need to + // handle this case separately, but we make it a special case to + // simplify the shift logic. + const shift: u32 = @truncate(u32, Z(1) -% @bitCast(u32, productExponent)); + if (shift >= typeWidth) return @bitCast(T, productSign); + + // Otherwise, shift the significand of the result so that the round + // bit is the high bit of productLo. + wideRightShiftWithSticky(Z, &productHi, &productLo, shift); + } else { + // Result is normal before rounding; insert the exponent. + productHi &= significandMask; + productHi |= Z(@bitCast(u32, productExponent)) << significandBits; + } + + // Insert the sign of the result: + productHi |= productSign; + + // Final rounding. The final result may overflow to infinity, or underflow + // to zero, but those are the correct results in those cases. We use the + // default IEEE-754 round-to-nearest, ties-to-even rounding mode. + if (productLo > signBit) productHi +%= 1; + if (productLo == signBit) productHi +%= productHi & 1; + return @bitCast(T, productHi); +} + +fn wideMultiply(comptime Z: type, a: Z, b: Z, hi: *Z, lo: *Z) void { + switch (Z) { + u32 => { + // 32x32 --> 64 bit multiply + const product = u64(a) * u64(b); + hi.* = @truncate(u32, product >> 32); + lo.* = @truncate(u32, product); + }, + u64 => { + const S = struct { + fn loWord(x: u64) u64 { + return @truncate(u32, x); + } + fn hiWord(x: u64) u64 { + return @truncate(u32, x >> 32); + } + }; + // 64x64 -> 128 wide multiply for platforms that don't have such an operation; + // many 64-bit platforms have this operation, but they tend to have hardware + // floating-point, so we don't bother with a special case for them here. + // Each of the component 32x32 -> 64 products + const plolo: u64 = S.loWord(a) * S.loWord(b); + const plohi: u64 = S.loWord(a) * S.hiWord(b); + const philo: u64 = S.hiWord(a) * S.loWord(b); + const phihi: u64 = S.hiWord(a) * S.hiWord(b); + // Sum terms that contribute to lo in a way that allows us to get the carry + const r0: u64 = S.loWord(plolo); + const r1: u64 = S.hiWord(plolo) +% S.loWord(plohi) +% S.loWord(philo); + lo.* = r0 +% (r1 << 32); + // Sum terms contributing to hi with the carry from lo + hi.* = S.hiWord(plohi) +% S.hiWord(philo) +% S.hiWord(r1) +% phihi; + }, + u128 => { + const Word_LoMask = u64(0x00000000ffffffff); + const Word_HiMask = u64(0xffffffff00000000); + const Word_FullMask = u64(0xffffffffffffffff); + const S = struct { + fn Word_1(x: u128) u64 { + return @truncate(u32, x >> 96); + } + fn Word_2(x: u128) u64 { + return @truncate(u32, x >> 64); + } + fn Word_3(x: u128) u64 { + return @truncate(u32, x >> 32); + } + fn Word_4(x: u128) u64 { + return @truncate(u32, x); + } + }; + // 128x128 -> 256 wide multiply for platforms that don't have such an operation; + // many 64-bit platforms have this operation, but they tend to have hardware + // floating-point, so we don't bother with a special case for them here. + + const product11: u64 = S.Word_1(a) * S.Word_1(b); + const product12: u64 = S.Word_1(a) * S.Word_2(b); + const product13: u64 = S.Word_1(a) * S.Word_3(b); + const product14: u64 = S.Word_1(a) * S.Word_4(b); + const product21: u64 = S.Word_2(a) * S.Word_1(b); + const product22: u64 = S.Word_2(a) * S.Word_2(b); + const product23: u64 = S.Word_2(a) * S.Word_3(b); + const product24: u64 = S.Word_2(a) * S.Word_4(b); + const product31: u64 = S.Word_3(a) * S.Word_1(b); + const product32: u64 = S.Word_3(a) * S.Word_2(b); + const product33: u64 = S.Word_3(a) * S.Word_3(b); + const product34: u64 = S.Word_3(a) * S.Word_4(b); + const product41: u64 = S.Word_4(a) * S.Word_1(b); + const product42: u64 = S.Word_4(a) * S.Word_2(b); + const product43: u64 = S.Word_4(a) * S.Word_3(b); + const product44: u64 = S.Word_4(a) * S.Word_4(b); + + const sum0: u128 = u128(product44); + const sum1: u128 = u128(product34) +% + u128(product43); + const sum2: u128 = u128(product24) +% + u128(product33) +% + u128(product42); + const sum3: u128 = u128(product14) +% + u128(product23) +% + u128(product32) +% + u128(product41); + const sum4: u128 = u128(product13) +% + u128(product22) +% + u128(product31); + const sum5: u128 = u128(product12) +% + u128(product21); + const sum6: u128 = u128(product11); + + const r0: u128 = (sum0 & Word_FullMask) +% + ((sum1 & Word_LoMask) << 32); + const r1: u128 = (sum0 >> 64) +% + ((sum1 >> 32) & Word_FullMask) +% + (sum2 & Word_FullMask) +% + ((sum3 << 32) & Word_HiMask); + + lo.* = r0 +% (r1 << 64); + hi.* = (r1 >> 64) +% + (sum1 >> 96) +% + (sum2 >> 64) +% + (sum3 >> 32) +% + sum4 +% + (sum5 << 32) +% + (sum6 << 64); + }, + else => @compileError("unsupported"), + } +} + +fn normalize(comptime T: type, significand: *@IntType(false, T.bit_count)) i32 { + const Z = @IntType(false, T.bit_count); + const significandBits = std.math.floatMantissaBits(T); + const implicitBit = Z(1) << significandBits; + + const shift = @clz(significand.*) - @clz(implicitBit); + significand.* <<= @intCast(std.math.Log2Int(Z), shift); + return 1 - shift; +} + +fn wideRightShiftWithSticky(comptime Z: type, hi: *Z, lo: *Z, count: u32) void { + const typeWidth = Z.bit_count; + const S = std.math.Log2Int(Z); + if (count < typeWidth) { + const sticky = @truncate(u8, lo.* << @intCast(S, typeWidth -% count)); + lo.* = (hi.* << @intCast(S, typeWidth -% count)) | (lo.* >> @intCast(S, count)) | sticky; + hi.* = hi.* >> @intCast(S, count); + } else if (count < 2 * typeWidth) { + const sticky = @truncate(u8, hi.* << @intCast(S, 2 * typeWidth -% count) | lo.*); + lo.* = hi.* >> @intCast(S, count -% typeWidth) | sticky; + hi.* = 0; + } else { + const sticky = @truncate(u8, hi.* | lo.*); + lo.* = sticky; + hi.* = 0; + } +} + +test "import mulXf3" { + _ = @import("mulXf3_test.zig"); +} diff --git a/std/special/compiler_rt/mulXf3_test.zig b/std/special/compiler_rt/mulXf3_test.zig new file mode 100644 index 000000000..1c0c0fc04 --- /dev/null +++ b/std/special/compiler_rt/mulXf3_test.zig @@ -0,0 +1,86 @@ +// Ported from: +// +// https://github.com/llvm/llvm-project/blob/2ffb1b0413efa9a24eb3c49e710e36f92e2cb50b/compiler-rt/test/builtins/Unit/multf3_test.c + +const qnan128 = @bitCast(f128, u128(0x7fff800000000000) << 64); +const inf128 = @bitCast(f128, u128(0x7fff000000000000) << 64); + +const __multf3 = @import("mulXf3.zig").__multf3; + +// return true if equal +// use two 64-bit integers intead of one 128-bit integer +// because 128-bit integer constant can't be assigned directly +fn compareResultLD(result: f128, expectedHi: u64, expectedLo: u64) bool { + const rep = @bitCast(u128, result); + const hi = @intCast(u64, rep >> 64); + const lo = @truncate(u64, rep); + + if (hi == expectedHi and lo == expectedLo) { + return true; + } + // test other possible NaN representation(signal NaN) + if (expectedHi == 0x7fff800000000000 and expectedLo == 0x0) { + if ((hi & 0x7fff000000000000) == 0x7fff000000000000 and + ((hi & 0xffffffffffff) > 0 or lo > 0)) + { + return true; + } + } + return false; +} + +fn test__multf3(a: f128, b: f128, expected_hi: u64, expected_lo: u64) void { + const x = __multf3(a, b); + + if (compareResultLD(x, expected_hi, expected_lo)) + return; + + @panic("__multf3 test failure"); +} + +fn makeNaN128(rand: u64) f128 { + const int_result = u128(0x7fff000000000000 | (rand & 0xffffffffffff)) << 64; + const float_result = @bitCast(f128, int_result); + return float_result; +} +test "multf3" { + // qNaN * any = qNaN + test__multf3(qnan128, 0x1.23456789abcdefp+5, 0x7fff800000000000, 0x0); + + // NaN * any = NaN + const a = makeNaN128(0x800030000000); + test__multf3(a, 0x1.23456789abcdefp+5, 0x7fff800000000000, 0x0); + // inf * any = inf + test__multf3(inf128, 0x1.23456789abcdefp+5, 0x7fff000000000000, 0x0); + + // any * any + test__multf3( + @bitCast(f128, u128(0x40042eab345678439abcdefea5678234)), + @bitCast(f128, u128(0x3ffeedcb34a235253948765432134675)), + 0x400423e7f9e3c9fc, + 0xd906c2c2a85777c4, + ); + + test__multf3( + @bitCast(f128, u128(0x3fcd353e45674d89abacc3a2ebf3ff50)), + @bitCast(f128, u128(0x3ff6ed8764648369535adf4be3214568)), + 0x3fc52a163c6223fc, + 0xc94c4bf0430768b4, + ); + + test__multf3( + 0x1.234425696abcad34a35eeffefdcbap+456, + 0x451.ed98d76e5d46e5f24323dff21ffp+600, + 0x44293a91de5e0e94, + 0xe8ed17cc2cdf64ac, + ); + + test__multf3( + @bitCast(f128, u128(0x3f154356473c82a9fabf2d22ace345df)), + @bitCast(f128, u128(0x3e38eda98765476743ab21da23d45679)), + 0x3d4f37c1a3137cae, + 0xfc6807048bc2836a, + ); + + test__multf3(0x1.23456734245345p-10000, 0x1.edcba524498724p-6497, 0x0, 0x0); +}