diff --git a/lib/std/special/compiler_rt.zig b/lib/std/special/compiler_rt.zig index 375de6fced..0f5d6c4d23 100644 --- a/lib/std/special/compiler_rt.zig +++ b/lib/std/special/compiler_rt.zig @@ -726,6 +726,11 @@ comptime { if (!is_test) { @export(fmodl, .{ .name = "fmodl", .linkage = linkage }); + if (long_double_is_f80) { + @export(fmodl, .{ .name = "fmodx", .linkage = linkage }); + } else { + @export(fmodx, .{ .name = "fmodx", .linkage = linkage }); + } if (long_double_is_f128) { @export(fmodl, .{ .name = "fmodq", .linkage = linkage }); } else { @@ -884,13 +889,8 @@ fn ceill(x: c_longdouble) callconv(.C) c_longdouble { return math.ceil(x); } -const fmodq = @import("compiler_rt/floatfmodq.zig").fmodq; -fn fmodl(x: c_longdouble, y: c_longdouble) callconv(.C) c_longdouble { - if (!long_double_is_f128) { - @panic("TODO implement this"); - } - return @floatCast(c_longdouble, fmodq(x, y)); -} +const fmodq = @import("compiler_rt/fmodq.zig").fmodq; +const fmodx = @import("compiler_rt/fmodx.zig").fmodx; // Avoid dragging in the runtime safety mechanisms into this .o file, // unless we're trying to test this file. diff --git a/lib/std/special/compiler_rt/divdf3.zig b/lib/std/special/compiler_rt/divdf3.zig index 2148902de2..137f5c02f9 100644 --- a/lib/std/special/compiler_rt/divdf3.zig +++ b/lib/std/special/compiler_rt/divdf3.zig @@ -314,12 +314,11 @@ pub fn wideMultiply(comptime Z: type, a: Z, b: Z, hi: *Z, lo: *Z) void { pub fn normalize(comptime T: type, significand: *std.meta.Int(.unsigned, @typeInfo(T).Float.bits)) i32 { @setRuntimeSafety(builtin.is_test); const Z = std.meta.Int(.unsigned, @typeInfo(T).Float.bits); - const significandBits = std.math.floatMantissaBits(T); - const implicitBit = @as(Z, 1) << significandBits; + const integerBit = @as(Z, 1) << std.math.floatFractionalBits(T); - const shift = @clz(Z, significand.*) - @clz(Z, implicitBit); + const shift = @clz(Z, significand.*) - @clz(Z, integerBit); significand.* <<= @intCast(std.math.Log2Int(Z), shift); - return 1 - shift; + return @as(i32, 1) - shift; } pub fn __aeabi_ddiv(a: f64, b: f64) callconv(.AAPCS) f64 { diff --git a/lib/std/special/compiler_rt/floatfmodq.zig b/lib/std/special/compiler_rt/fmodq.zig similarity index 93% rename from lib/std/special/compiler_rt/floatfmodq.zig rename to lib/std/special/compiler_rt/fmodq.zig index b8da727c90..3f56c49796 100644 --- a/lib/std/special/compiler_rt/floatfmodq.zig +++ b/lib/std/special/compiler_rt/fmodq.zig @@ -27,7 +27,7 @@ pub fn fmodq(a: f128, b: f128) callconv(.C) f128 { const signA = aPtr_u16[exp_and_sign_index] & 0x8000; var expA = @intCast(i32, (aPtr_u16[exp_and_sign_index] & 0x7fff)); - var expB = bPtr_u16[exp_and_sign_index] & 0x7fff; + var expB = @intCast(i32, (bPtr_u16[exp_and_sign_index] & 0x7fff)); // There are 3 cases where the answer is undefined, check for: // - fmodq(val, 0) @@ -55,12 +55,12 @@ pub fn fmodq(a: f128, b: f128) callconv(.C) f128 { if (expA == 0) { amod *= 0x1p120; - expA = aPtr_u16[exp_and_sign_index] -% 120; + expA = @as(i32, aPtr_u16[exp_and_sign_index]) - 120; } if (expB == 0) { bmod *= 0x1p120; - expB = bPtr_u16[exp_and_sign_index] -% 120; + expB = @as(i32, bPtr_u16[exp_and_sign_index]) - 120; } // OR in extra non-stored mantissa digit @@ -73,7 +73,7 @@ pub fn fmodq(a: f128, b: f128) callconv(.C) f128 { var high = highA -% highB; var low = lowA -% lowB; if (lowA < lowB) { - high = highA -% 1; + high -%= 1; } if (high >> 63 == 0) { if ((high | low) == 0) { @@ -122,5 +122,5 @@ pub fn fmodq(a: f128, b: f128) callconv(.C) f128 { } test { - _ = @import("floatfmodq_test.zig"); + _ = @import("fmodq_test.zig"); } diff --git a/lib/std/special/compiler_rt/floatfmodq_test.zig b/lib/std/special/compiler_rt/fmodq_test.zig similarity index 67% rename from lib/std/special/compiler_rt/floatfmodq_test.zig rename to lib/std/special/compiler_rt/fmodq_test.zig index a272b797e3..b8baf8ae9b 100644 --- a/lib/std/special/compiler_rt/floatfmodq_test.zig +++ b/lib/std/special/compiler_rt/fmodq_test.zig @@ -1,5 +1,5 @@ const std = @import("std"); -const fmodq = @import("floatfmodq.zig"); +const fmodq = @import("fmodq.zig"); const testing = std.testing; fn test_fmodq(a: f128, b: f128, exp: f128) !void { @@ -34,12 +34,18 @@ test "fmodq" { try test_fmodq(-0.0, 1.0, -0.0); try test_fmodq(7046119.0, 5558362.0, 1487757.0); try test_fmodq(9010357.0, 1957236.0, 1181413.0); + try test_fmodq(5192296858534827628530496329220095, 10.0, 5.0); + try test_fmodq(5192296858534827628530496329220095, 922337203681230954775807, 220474884073715748246157); // Denormals - const a: f128 = 0xedcb34a235253948765432134674p-16494; - const b: f128 = 0x5d2e38791cfbc0737402da5a9518p-16494; - const exp: f128 = 0x336ec3affb2db8618e4e7d5e1c44p-16494; - try test_fmodq(a, b, exp); + const a1: f128 = 0xedcb34a235253948765432134674p-16494; + const b1: f128 = 0x5d2e38791cfbc0737402da5a9518p-16494; + const exp1: f128 = 0x336ec3affb2db8618e4e7d5e1c44p-16494; + try test_fmodq(a1, b1, exp1); + const a2: f128 = 0x0.7654_3210_fdec_ba98_7654_3210_fdecp-16382; + const b2: f128 = 0x0.0012_fdac_bdef_1234_fdec_3222_1111p-16382; + const exp2: f128 = 0x0.0001_aecd_9d66_4a6e_67b7_d7d0_a901p-16382; + try test_fmodq(a2, b2, exp2); try test_fmodq_nans(); try test_fmodq_infs(); diff --git a/lib/std/special/compiler_rt/fmodx.zig b/lib/std/special/compiler_rt/fmodx.zig new file mode 100644 index 0000000000..efe16f9f16 --- /dev/null +++ b/lib/std/special/compiler_rt/fmodx.zig @@ -0,0 +1,108 @@ +const builtin = @import("builtin"); +const std = @import("std"); +const math = std.math; +const normalize = @import("divdf3.zig").normalize; + +// fmodx - floating modulo large, returns the remainder of division for f80 types +// Logic and flow heavily inspired by MUSL fmodl for 113 mantissa digits +pub fn fmodx(a: f80, b: f80) callconv(.C) f80 { + @setRuntimeSafety(builtin.is_test); + + const T = f80; + const Z = std.meta.Int(.unsigned, @bitSizeOf(T)); + + const significandBits = math.floatMantissaBits(T); + const fractionalBits = math.floatFractionalBits(T); + const exponentBits = math.floatExponentBits(T); + + const signBit = (@as(Z, 1) << (significandBits + exponentBits)); + const maxExponent = ((1 << exponentBits) - 1); + + var aRep = @bitCast(Z, a); + var bRep = @bitCast(Z, b); + + const signA = aRep & signBit; + var expA = @intCast(i32, (@bitCast(Z, a) >> significandBits) & maxExponent); + var expB = @intCast(i32, (@bitCast(Z, b) >> significandBits) & maxExponent); + + // There are 3 cases where the answer is undefined, check for: + // - fmodx(val, 0) + // - fmodx(val, NaN) + // - fmodx(inf, val) + // The sign on checked values does not matter. + // Doing (a * b) / (a * b) procudes undefined results + // because the three cases always produce undefined calculations: + // - 0 / 0 + // - val * NaN + // - inf / inf + if (b == 0 or math.isNan(b) or expA == maxExponent) { + return (a * b) / (a * b); + } + + // Remove the sign from both + aRep &= ~signBit; + bRep &= ~signBit; + if (aRep <= bRep) { + if (aRep == bRep) { + return 0 * a; + } + return a; + } + + if (expA == 0) expA = normalize(f80, &aRep); + if (expB == 0) expB = normalize(f80, &bRep); + + var highA: u64 = 0; + var highB: u64 = 0; + var lowA: u64 = @truncate(u64, aRep); + var lowB: u64 = @truncate(u64, bRep); + + while (expA > expB) : (expA -= 1) { + var high = highA -% highB; + var low = lowA -% lowB; + if (lowA < lowB) { + high -%= 1; + } + if (high >> 63 == 0) { + if ((high | low) == 0) { + return 0 * a; + } + highA = 2 *% high + (low >> 63); + lowA = 2 *% low; + } else { + highA = 2 *% highA + (lowA >> 63); + lowA = 2 *% lowA; + } + } + + var high = highA -% highB; + var low = lowA -% lowB; + if (lowA < lowB) { + high -%= 1; + } + if (high >> 63 == 0) { + if ((high | low) == 0) { + return 0 * a; + } + highA = high; + lowA = low; + } + + while ((lowA >> fractionalBits) == 0) { + lowA = 2 *% lowA; + expA = expA - 1; + } + + // Combine the exponent with the sign and significand, normalize if happened to be denormalized + if (expA < -fractionalBits) { + return @bitCast(T, signA); + } else if (expA <= 0) { + return @bitCast(T, (lowA >> @intCast(math.Log2Int(u64), 1 - expA)) | signA); + } else { + return @bitCast(T, lowA | (@as(Z, @intCast(u16, expA)) << significandBits) | signA); + } +} + +test { + _ = @import("fmodx_test.zig"); +} diff --git a/lib/std/special/compiler_rt/fmodx_test.zig b/lib/std/special/compiler_rt/fmodx_test.zig new file mode 100644 index 0000000000..a5d0887ea4 --- /dev/null +++ b/lib/std/special/compiler_rt/fmodx_test.zig @@ -0,0 +1,51 @@ +const std = @import("std"); +const fmodx = @import("fmodx.zig"); +const testing = std.testing; + +fn test_fmodx(a: f80, b: f80, exp: f80) !void { + const res = fmodx.fmodx(a, b); + try testing.expect(exp == res); +} + +fn test_fmodx_nans() !void { + try testing.expect(std.math.isNan(fmodx.fmodx(1.0, std.math.nan(f80)))); + try testing.expect(std.math.isNan(fmodx.fmodx(1.0, -std.math.nan(f80)))); + try testing.expect(std.math.isNan(fmodx.fmodx(std.math.nan(f80), 1.0))); + try testing.expect(std.math.isNan(fmodx.fmodx(-std.math.nan(f80), 1.0))); +} + +fn test_fmodx_infs() !void { + try testing.expect(fmodx.fmodx(1.0, std.math.inf(f80)) == 1.0); + try testing.expect(fmodx.fmodx(1.0, -std.math.inf(f80)) == 1.0); + try testing.expect(std.math.isNan(fmodx.fmodx(std.math.inf(f80), 1.0))); + try testing.expect(std.math.isNan(fmodx.fmodx(-std.math.inf(f80), 1.0))); +} + +test "fmodx" { + try test_fmodx(6.4, 4.0, 2.4); + try test_fmodx(6.4, -4.0, 2.4); + try test_fmodx(-6.4, 4.0, -2.4); + try test_fmodx(-6.4, -4.0, -2.4); + try test_fmodx(3.0, 2.0, 1.0); + try test_fmodx(-5.0, 3.0, -2.0); + try test_fmodx(3.0, 2.0, 1.0); + try test_fmodx(1.0, 2.0, 1.0); + try test_fmodx(0.0, 1.0, 0.0); + try test_fmodx(-0.0, 1.0, -0.0); + try test_fmodx(7046119.0, 5558362.0, 1487757.0); + try test_fmodx(9010357.0, 1957236.0, 1181413.0); + try test_fmodx(9223372036854775807, 10.0, 7.0); + + // Denormals + const a1: f80 = 0x0.76e5_9a51_1a92_9ca4p-16381; + const b1: f80 = 0x0.2e97_1c3c_8e7d_e03ap-16381; + const exp1: f80 = 0x0.19b7_61d7_fd96_dc30p-16381; + try test_fmodx(a1, b1, exp1); + const a2: f80 = 0x0.76e5_9a51_1a92_9ca4p-16381; + const b2: f80 = 0x0.0e97_1c3c_8e7d_e03ap-16381; + const exp2: f80 = 0x0.022c_b86c_a6a3_9ad4p-16381; + try test_fmodx(a2, b2, exp2); + + try test_fmodx_nans(); + try test_fmodx_infs(); +}