From 6c0114e04436404b1c59b66125733350b4be2f5b Mon Sep 17 00:00:00 2001
From: Cody Tapscott <topolarity@tapscott.me>
Date: Mon, 25 Apr 2022 13:25:15 -0700
Subject: [PATCH 1/5] compiler_rt: Implement fmodx for f80

---
 lib/std/special/compiler_rt.zig               |  14 +--
 lib/std/special/compiler_rt/divdf3.zig        |   7 +-
 .../compiler_rt/{floatfmodq.zig => fmodq.zig} |  10 +-
 .../{floatfmodq_test.zig => fmodq_test.zig}   |  16 ++-
 lib/std/special/compiler_rt/fmodx.zig         | 108 ++++++++++++++++++
 lib/std/special/compiler_rt/fmodx_test.zig    |  51 +++++++++
 6 files changed, 185 insertions(+), 21 deletions(-)
 rename lib/std/special/compiler_rt/{floatfmodq.zig => fmodq.zig} (93%)
 rename lib/std/special/compiler_rt/{floatfmodq_test.zig => fmodq_test.zig} (67%)
 create mode 100644 lib/std/special/compiler_rt/fmodx.zig
 create mode 100644 lib/std/special/compiler_rt/fmodx_test.zig

diff --git a/lib/std/special/compiler_rt.zig b/lib/std/special/compiler_rt.zig
index 375de6fced4dae5038020cb26c2dad1b05c4a06e..0f5d6c4d23a0707ae16b040baf6f9cd82f52e1c4 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 2148902de25e216888b6e22fc46ee24e8c933f39..137f5c02f93a1df74cdb4d0ab8e97f8c4d56e291 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 b8da727c90568e7f1e476c6d46edb470d57d0c5b..3f56c4979609837cd112b7f4da89381b769655ea 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 a272b797e3a8af26062f1cf75c06e00f066a8ebb..b8baf8ae9be86006c24a7c487d8f7479afbb56c2 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 0000000000000000000000000000000000000000..efe16f9f160d88a0c6e54a4114698c4e2c8d6993
--- /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 0000000000000000000000000000000000000000..a5d0887ea4a115eef6b241e64bb7dbd498740bb5
--- /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();
+}
-- 
2.54.0


From d930e015c75d60e22d0a116a1f4fa0df4d68b6b6 Mon Sep 17 00:00:00 2001
From: Cody Tapscott <topolarity@tapscott.me>
Date: Mon, 25 Apr 2022 16:31:43 -0700
Subject: [PATCH 2/5] compiler_rt: Implement __divxf3 for f80

---
 lib/std/special/compiler_rt.zig             |  16 +-
 lib/std/special/compiler_rt/divxf3.zig      | 202 ++++++++++++++++++++
 lib/std/special/compiler_rt/divxf3_test.zig |  65 +++++++
 3 files changed, 274 insertions(+), 9 deletions(-)
 create mode 100644 lib/std/special/compiler_rt/divxf3.zig
 create mode 100644 lib/std/special/compiler_rt/divxf3_test.zig

diff --git a/lib/std/special/compiler_rt.zig b/lib/std/special/compiler_rt.zig
index 0f5d6c4d23a0707ae16b040baf6f9cd82f52e1c4..6a44d859c92c6e2317032b727d3c2d2300499fe3 100644
--- a/lib/std/special/compiler_rt.zig
+++ b/lib/std/special/compiler_rt.zig
@@ -253,6 +253,8 @@ comptime {
     @export(__divsf3, .{ .name = "__divsf3", .linkage = linkage });
     const __divdf3 = @import("compiler_rt/divdf3.zig").__divdf3;
     @export(__divdf3, .{ .name = "__divdf3", .linkage = linkage });
+    const __divxf3 = @import("compiler_rt/divxf3.zig").__divxf3;
+    @export(__divxf3, .{ .name = "__divxf3", .linkage = linkage });
     const __divtf3 = @import("compiler_rt/divtf3.zig").__divtf3;
     @export(__divtf3, .{ .name = "__divtf3", .linkage = linkage });
 
@@ -725,17 +727,13 @@ 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 {
-            @export(fmodq, .{ .name = "fmodq", .linkage = linkage });
+            @export(fmodx, .{ .name = "fmodl", .linkage = linkage });
+        } else if (long_double_is_f128) {
+            @export(fmodq, .{ .name = "fmodl", .linkage = linkage });
         }
+        @export(fmodx, .{ .name = "fmodx", .linkage = linkage });
+        @export(fmodq, .{ .name = "fmodq", .linkage = linkage });
 
         @export(floorf, .{ .name = "floorf", .linkage = linkage });
         @export(floor, .{ .name = "floor", .linkage = linkage });
diff --git a/lib/std/special/compiler_rt/divxf3.zig b/lib/std/special/compiler_rt/divxf3.zig
new file mode 100644
index 0000000000000000000000000000000000000000..5f5ed667ec38d222b8ae67078d73fe0400608420
--- /dev/null
+++ b/lib/std/special/compiler_rt/divxf3.zig
@@ -0,0 +1,202 @@
+const std = @import("std");
+const builtin = @import("builtin");
+const normalize = @import("divdf3.zig").normalize;
+const wideMultiply = @import("divdf3.zig").wideMultiply;
+
+pub fn __divxf3(a: f80, b: f80) callconv(.C) f80 {
+    @setRuntimeSafety(builtin.is_test);
+    const T = f80;
+    const Z = std.meta.Int(.unsigned, @bitSizeOf(T));
+
+    const significandBits = std.math.floatMantissaBits(T);
+    const fractionalBits = std.math.floatFractionalBits(T);
+    const exponentBits = std.math.floatExponentBits(T);
+
+    const signBit = (@as(Z, 1) << (significandBits + exponentBits));
+    const maxExponent = ((1 << exponentBits) - 1);
+    const exponentBias = (maxExponent >> 1);
+
+    const integerBit = (@as(Z, 1) << fractionalBits);
+    const quietBit = integerBit >> 1;
+    const significandMask = (@as(Z, 1) << significandBits) - 1;
+
+    const absMask = signBit - 1;
+    const qnanRep = @bitCast(Z, std.math.nan(T)) | 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 quotientSign: 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 / infinity = NaN
+            if (bAbs == infRep) {
+                return @bitCast(T, qnanRep);
+            }
+            // infinity / anything else = +/- infinity
+            else {
+                return @bitCast(T, aAbs | quotientSign);
+            }
+        }
+
+        // anything else / infinity = +/- 0
+        if (bAbs == infRep) return @bitCast(T, quotientSign);
+
+        if (aAbs == 0) {
+            // zero / zero = NaN
+            if (bAbs == 0) {
+                return @bitCast(T, qnanRep);
+            }
+            // zero / anything else = +/- zero
+            else {
+                return @bitCast(T, quotientSign);
+            }
+        }
+        // anything else / zero = +/- infinity
+        if (bAbs == 0) return @bitCast(T, infRep | quotientSign);
+
+        // 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 < integerBit) scale +%= normalize(T, &aSignificand);
+        if (bAbs < integerBit) scale -%= normalize(T, &bSignificand);
+    }
+    var quotientExponent: i32 = @bitCast(i32, aExponent -% bExponent) +% scale;
+
+    // Align the significand of b as a Q63 fixed-point number in the range
+    // [1, 2.0) and get a Q64 approximate reciprocal using a small minimax
+    // polynomial approximation: reciprocal = 3/4 + 1/sqrt(2) - b/2.  This
+    // is accurate to about 3.5 binary digits.
+    const q63b = @intCast(u64, bSignificand);
+    var recip64 = @as(u64, 0x7504f333F9DE6484) -% q63b;
+    // 0x7504f333F9DE6484 / 2^64 + 1 = 3/4 + 1/sqrt(2)
+
+    // Now refine the reciprocal estimate using a Newton-Raphson iteration:
+    //
+    //     x1 = x0 * (2 - x0 * b)
+    //
+    // This doubles the number of correct binary digits in the approximation
+    // with each iteration.
+    var correction64: u64 = undefined;
+    correction64 = @truncate(u64, ~(@as(u128, recip64) *% q63b >> 64) +% 1);
+    recip64 = @truncate(u64, @as(u128, recip64) *% correction64 >> 63);
+    correction64 = @truncate(u64, ~(@as(u128, recip64) *% q63b >> 64) +% 1);
+    recip64 = @truncate(u64, @as(u128, recip64) *% correction64 >> 63);
+    correction64 = @truncate(u64, ~(@as(u128, recip64) *% q63b >> 64) +% 1);
+    recip64 = @truncate(u64, @as(u128, recip64) *% correction64 >> 63);
+    correction64 = @truncate(u64, ~(@as(u128, recip64) *% q63b >> 64) +% 1);
+    recip64 = @truncate(u64, @as(u128, recip64) *% correction64 >> 63);
+    correction64 = @truncate(u64, ~(@as(u128, recip64) *% q63b >> 64) +% 1);
+    recip64 = @truncate(u64, @as(u128, recip64) *% correction64 >> 63);
+
+    // The reciprocal may have overflowed to zero if the upper half of b is
+    // exactly 1.0.  This would sabatoge the full-width final stage of the
+    // computation that follows, so we adjust the reciprocal down by one bit.
+    recip64 -%= 1;
+
+    // We need to perform one more iteration to get us to 112 binary digits;
+    // The last iteration needs to happen with extra precision.
+
+    // NOTE: This operation is equivalent to __multi3, which is not implemented
+    //       in some architechures
+    var reciprocal: u128 = undefined;
+    var correction: u128 = undefined;
+    var dummy: u128 = undefined;
+    wideMultiply(u128, recip64, q63b, &dummy, &correction);
+
+    correction = -%correction;
+
+    const cHi = @truncate(u64, correction >> 64);
+    const cLo = @truncate(u64, correction);
+
+    var r64cH: u128 = undefined;
+    var r64cL: u128 = undefined;
+    wideMultiply(u128, recip64, cHi, &dummy, &r64cH);
+    wideMultiply(u128, recip64, cLo, &dummy, &r64cL);
+
+    reciprocal = r64cH + (r64cL >> 64);
+
+    // Adjust the final 128-bit reciprocal estimate downward to ensure that it
+    // is strictly smaller than the infinitely precise exact reciprocal. Because
+    // the computation of the Newton-Raphson step is truncating at every step,
+    // this adjustment is small; most of the work is already done.
+    reciprocal -%= 2;
+
+    // The numerical reciprocal is accurate to within 2^-112, lies in the
+    // interval [0.5, 1.0), and is strictly smaller than the true reciprocal
+    // of b.  Multiplying a by this reciprocal thus gives a numerical q = a/b
+    // in Q127 with the following properties:
+    //
+    //    1. q < a/b
+    //    2. q is in the interval [0.5, 2.0)
+    //    3. The error in q is bounded away from 2^-63 (actually, we have
+    //       many bits to spare, but this is all we need).
+
+    // We need a 128 x 128 multiply high to compute q.
+    var quotient128: u128 = undefined;
+    var quotientLo: u128 = undefined;
+    wideMultiply(u128, aSignificand << 2, reciprocal, &quotient128, &quotientLo);
+
+    // Two cases: quotient is in [0.5, 1.0) or quotient is in [1.0, 2.0).
+    // Right shift the quotient if it falls in the [1,2) range and adjust the
+    // exponent accordingly.
+    var quotient: u64 = if (quotient128 < (integerBit << 1)) b: {
+        quotientExponent -= 1;
+        break :b @intCast(u64, quotient128);
+    } else @intCast(u64, quotient128 >> 1);
+
+    // We are going to compute a residual of the form
+    //
+    //     r = a - q*b
+    //
+    // We know from the construction of q that r satisfies:
+    //
+    //     0 <= r < ulp(q)*b
+    //
+    // If r is greater than 1/2 ulp(q)*b, then q rounds up.  Otherwise, we
+    // already have the correct result.  The exact halfway case cannot occur.
+    var residual: u64 = -%(quotient *% q63b);
+
+    const writtenExponent = quotientExponent + exponentBias;
+    if (writtenExponent >= maxExponent) {
+        // If we have overflowed the exponent, return infinity.
+        return @bitCast(T, infRep | quotientSign);
+    } else if (writtenExponent < 1) {
+        if (writtenExponent == 0) {
+            // Check whether the rounded result is normal.
+            if (residual > (bSignificand >> 1)) { // round
+                if (quotient == (integerBit - 1)) // If the rounded result is normal, return it
+                    return @bitCast(T, @bitCast(Z, std.math.floatMin(T)) | quotientSign);
+            }
+        }
+        // Flush denormals to zero.  In the future, it would be nice to add
+        // code to round them correctly.
+        return @bitCast(T, quotientSign);
+    } else {
+        const round = @boolToInt(residual > (bSignificand >> 1));
+        // Insert the exponent
+        var absResult = quotient | (@intCast(Z, writtenExponent) << significandBits);
+        // Round
+        absResult +%= round;
+        // Insert the sign and return
+        return @bitCast(T, absResult | quotientSign | integerBit);
+    }
+}
+
+test {
+    _ = @import("divxf3_test.zig");
+}
diff --git a/lib/std/special/compiler_rt/divxf3_test.zig b/lib/std/special/compiler_rt/divxf3_test.zig
new file mode 100644
index 0000000000000000000000000000000000000000..b79b90c6fb05cd7686bdd031e909a212f4c911f5
--- /dev/null
+++ b/lib/std/special/compiler_rt/divxf3_test.zig
@@ -0,0 +1,65 @@
+const std = @import("std");
+const math = std.math;
+const testing = std.testing;
+
+const __divxf3 = @import("divxf3.zig").__divxf3;
+
+fn compareResult(result: f80, expected: u80) bool {
+    const rep = @bitCast(u80, result);
+
+    if (rep == expected) return true;
+    // test other possible NaN representations (signal NaN)
+    if (math.isNan(result) and math.isNan(@bitCast(f80, expected))) return true;
+
+    return false;
+}
+
+fn expect__divxf3_result(a: f80, b: f80, expected: u80) !void {
+    const x = __divxf3(a, b);
+    const ret = compareResult(x, expected);
+    try testing.expect(ret == true);
+}
+
+fn test__divxf3(a: f80, b: f80) !void {
+    const integerBit = 1 << math.floatFractionalBits(f80);
+    const x = __divxf3(a, b);
+
+    // Next float (assuming normal, non-zero result)
+    const x_plus_eps = @bitCast(f80, (@bitCast(u80, x) + 1) | integerBit);
+    // Prev float (assuming normal, non-zero result)
+    const x_minus_eps = @bitCast(f80, (@bitCast(u80, x) - 1) | integerBit);
+
+    // Make sure result is more accurate than the adjacent floats
+    const err_x = std.math.fabs(@mulAdd(f80, x, b, -a));
+    const err_x_plus_eps = std.math.fabs(@mulAdd(f80, x_plus_eps, b, -a));
+    const err_x_minus_eps = std.math.fabs(@mulAdd(f80, x_minus_eps, b, -a));
+
+    try testing.expect(err_x_minus_eps > err_x);
+    try testing.expect(err_x_plus_eps > err_x);
+}
+
+test "divxf3" {
+    // qNaN / any = qNaN
+    try expect__divxf3_result(math.qnan_f80, 0x1.23456789abcdefp+5, 0x7fffC000000000000000);
+    // NaN / any = NaN
+    try expect__divxf3_result(math.nan_f80, 0x1.23456789abcdefp+5, 0x7fffC000000000000000);
+    // inf / any(except inf and nan) = inf
+    try expect__divxf3_result(math.inf(f80), 0x1.23456789abcdefp+5, 0x7fff8000000000000000);
+    // inf / inf = nan
+    try expect__divxf3_result(math.inf(f80), math.inf(f80), 0x7fffC000000000000000);
+    // inf / nan = nan
+    try expect__divxf3_result(math.inf(f80), math.nan(f80), 0x7fffC000000000000000);
+
+    try test__divxf3(0x1.a23b45362464523375893ab4cdefp+5, 0x1.eedcbaba3a94546558237654321fp-1);
+    try test__divxf3(0x1.a2b34c56d745382f9abf2c3dfeffp-50, 0x1.ed2c3ba15935332532287654321fp-9);
+    try test__divxf3(0x1.2345f6aaaa786555f42432abcdefp+456, 0x1.edacbba9874f765463544dd3621fp+6400);
+    try test__divxf3(0x1.2d3456f789ba6322bc665544edefp-234, 0x1.eddcdba39f3c8b7a36564354321fp-4455);
+    try test__divxf3(0x1.2345f6b77b7a8953365433abcdefp+234, 0x1.edcba987d6bb3aa467754354321fp-4055);
+    try test__divxf3(0x1.a23b45362464523375893ab4cdefp+5, 0x1.a2b34c56d745382f9abf2c3dfeffp-50);
+    try test__divxf3(0x1.a23b45362464523375893ab4cdefp+5, 0x1.1234567890abcdef987654321123p0);
+    try test__divxf3(0x1.a23b45362464523375893ab4cdefp+5, 0x1.12394205810257120adae8929f23p+16);
+    try test__divxf3(0x1.a23b45362464523375893ab4cdefp+5, 0x1.febdcefa1231245f9abf2c3dfeffp-50);
+
+    // Result rounds down to zero
+    try expect__divxf3_result(6.72420628622418701252535563464350521E-4932, 2.0, 0x0);
+}
-- 
2.54.0


From f5540778b56636cffd89f794c6b1a9266d7d36f8 Mon Sep 17 00:00:00 2001
From: Cody Tapscott <topolarity@tapscott.me>
Date: Mon, 25 Apr 2022 16:32:42 -0700
Subject: [PATCH 3/5] stdlib: Fix hex-float printing for f80

---
 lib/std/fmt.zig          | 7 ++++---
 lib/std/math/signbit.zig | 6 ++++++
 2 files changed, 10 insertions(+), 3 deletions(-)

diff --git a/lib/std/fmt.zig b/lib/std/fmt.zig
index 80a2e87c5c8c395049f48b65da205fea6937a114..994350ffa719da38440b2dd0094704e6073cfe7a 100644
--- a/lib/std/fmt.zig
+++ b/lib/std/fmt.zig
@@ -1124,6 +1124,7 @@ pub fn formatFloatHexadecimal(
     const TU = std.meta.Int(.unsigned, std.meta.bitCount(T));
 
     const mantissa_bits = math.floatMantissaBits(T);
+    const fractional_bits = math.floatFractionalBits(T);
     const exponent_bits = math.floatExponentBits(T);
     const mantissa_mask = (1 << mantissa_bits) - 1;
     const exponent_mask = (1 << exponent_bits) - 1;
@@ -1155,14 +1156,14 @@ pub fn formatFloatHexadecimal(
         // Adjust the exponent for printing.
         exponent += 1;
     } else {
-        // Add the implicit 1.
-        mantissa |= 1 << mantissa_bits;
+        if (fractional_bits == mantissa_bits)
+            mantissa |= 1 << fractional_bits; // Add the implicit integer bit.
     }
 
     // Fill in zeroes to round the mantissa width to a multiple of 4.
     if (T == f16) mantissa <<= 2 else if (T == f32) mantissa <<= 1;
 
-    const mantissa_digits = (mantissa_bits + 3) / 4;
+    const mantissa_digits = (fractional_bits + 3) / 4;
 
     if (options.precision) |precision| {
         // Round if needed.
diff --git a/lib/std/math/signbit.zig b/lib/std/math/signbit.zig
index 056add5ae6d8d8eebece47a93c70f69dd6cd03d6..a72d6ab11997ce640b91c395c4d98d8a0a37024a 100644
--- a/lib/std/math/signbit.zig
+++ b/lib/std/math/signbit.zig
@@ -9,6 +9,7 @@ pub fn signbit(x: anytype) bool {
         f16 => signbit16(x),
         f32 => signbit32(x),
         f64 => signbit64(x),
+        f80 => signbit80(x),
         f128 => signbit128(x),
         else => @compileError("signbit not implemented for " ++ @typeName(T)),
     };
@@ -29,6 +30,11 @@ fn signbit64(x: f64) bool {
     return bits >> 63 != 0;
 }
 
+fn signbit80(x: f80) bool {
+    const bits = @bitCast(u80, x);
+    return bits >> 79 != 0;
+}
+
 fn signbit128(x: f128) bool {
     const bits = @bitCast(u128, x);
     return bits >> 127 != 0;
-- 
2.54.0


From 96d86e346517b1917165b6f817b2147524a9a0b5 Mon Sep 17 00:00:00 2001
From: Cody Tapscott <topolarity@tapscott.me>
Date: Mon, 25 Apr 2022 16:42:13 -0700
Subject: [PATCH 4/5] compiler_rt: Fix rounding edge case for __mulxf3

---
 lib/std/special/compiler_rt/README.md       | 34 ++++++++++-----------
 lib/std/special/compiler_rt/addXf3.zig      | 12 +++-----
 lib/std/special/compiler_rt/addXf3_test.zig |  2 ++
 lib/std/special/compiler_rt/mulXf3.zig      | 10 +++++-
 lib/std/special/compiler_rt/mulXf3_test.zig |  4 +++
 5 files changed, 37 insertions(+), 25 deletions(-)

diff --git a/lib/std/special/compiler_rt/README.md b/lib/std/special/compiler_rt/README.md
index 7115a5756f9796a0b2f1e603a83509ab068d7aa5..cf1bf29522467f895c181f86df3d2d7cefa1f013 100644
--- a/lib/std/special/compiler_rt/README.md
+++ b/lib/std/special/compiler_rt/README.md
@@ -152,53 +152,53 @@ Bugs should be solved by trying to duplicate the bug upstream, if possible.
 - todo todo __fixsfsi      // convert a to i32, rounding towards zero
 - todo todo __fixdfsi      //
 - todo todo __fixtfsi      //
-- none none __fixxfsi      // missing
+- todo todo __fixxfsi      //
 - todo todo __fixsfdi      // convert a to i64, rounding towards zero
 - todo todo __fixdfdi      //
 - todo todo __fixtfdi      //
-- none none __fixxfdi      // missing
+- todo todo __fixxfdi      //
 - todo todo __fixsfti      // convert a to i128, rounding towards zero
 - todo todo __fixdfti      //
 - todo todo __fixtfdi      //
-- none none __fixxfti      // missing
+- todo todo __fixxfti      //
 
 - __fixunssfsi   // convert to u32, rounding towards zero. negative values become 0.
 - __fixunsdfsi   //
 - __fixunstfsi   //
-- __fixunsxfsi   // missing
+- __fixunsxfsi   //
 - __fixunssfdi   // convert to u64, rounding towards zero. negative values become 0.
 - __fixunsdfdi   //
 - __fixunstfdi   //
-- __fixunsxfdi   // missing
+- __fixunsxfdi   //
 - __fixunssfti   // convert to u128, rounding towards zero. negative values become 0.
 - __fixunsdfti   //
 - __fixunstfdi   //
-- __fixunsxfti   // missing
+- __fixunsxfti   //
 
 - __floatsisf    // convert i32 to floating point
 - __floatsidf    //
 - __floatsitf    //
-- __floatsixf    // missing
+- __floatsixf    //
 - __floatdisf    // convert i64 to floating point
 - __floatdidf    //
 - __floatditf    //
-- __floatdixf    // missing
+- __floatdixf    //
 - __floattisf    // convert i128 to floating point
 - __floattidf    //
-- __floattixf    // missing
+- __floattixf    //
 
-- __floatunsisf  // convert i32 to floating point
+- __floatunsisf  // convert u32 to floating point
 - __floatunsidf  //
 - __floatunsitf  //
-- __floatunsixf  // missing
-- __floatundisf  // convert i64 to floating point
+- __floatunsixf  //
+- __floatundisf  // convert u64 to floating point
 - __floatundidf  //
 - __floatunditf  //
-- __floatundixf  // missing
-- __floatuntisf  // convert i128 to floating point
+- __floatundixf  //
+- __floatuntisf  // convert u128 to floating point
 - __floatuntidf  //
 - __floatuntitf  //
-- __floatuntixf  // missing
+- __floatuntixf  //
 
 #### Float Comparison
 - __cmpsf2       // return (a<b)=>-1,(a==b)=>0,(a>b)=>1,Nan=>1 dont rely on this
@@ -242,11 +242,11 @@ Bugs should be solved by trying to duplicate the bug upstream, if possible.
 - __mulsf3       // a * b
 - __muldf3       // a * b
 - __multf3       // a * b
-- __mulxf3       // a * b missing
+- __mulxf3       // a * b
 - __divsf3       // a / b
 - __divdf3       // a / b
 - __divtf3       // a / b
-- __divxf3       // a / b missing
+- __divxf3       // a / b
 - __negsf2       // -a symbol-level compatibility: libgcc uses this for the rl78
 - __negdf2       // -a unnecessary: can be lowered directly to a xor
 - __negtf2       // -a
diff --git a/lib/std/special/compiler_rt/addXf3.zig b/lib/std/special/compiler_rt/addXf3.zig
index c53d60600d7924490792c9f6e493adb3136eaf99..1a9de0fb748788c135afe4d025483fcdd426e39b 100644
--- a/lib/std/special/compiler_rt/addXf3.zig
+++ b/lib/std/special/compiler_rt/addXf3.zig
@@ -74,7 +74,7 @@ fn normalize(comptime T: type, significand: *std.meta.Int(.unsigned, @typeInfo(T
 
     const shift = @clz(std.meta.Int(.unsigned, bits), significand.*) - @clz(Z, integerBit);
     significand.* <<= @intCast(S, shift);
-    return 1 - shift;
+    return @as(i32, 1) - shift;
 }
 
 // TODO: restore inline keyword, see: https://github.com/ziglang/zig/issues/2154
@@ -210,12 +210,10 @@ fn addXf3(comptime T: type, a: T, b: T) T {
     if (aExponent >= maxExponent) return @bitCast(T, infRep | resultSign);
 
     if (aExponent <= 0) {
-        // Result is denormal before rounding; the exponent is zero and we
-        // need to shift the significand.
-        const shift = @intCast(Z, 1 - aExponent);
-        const sticky = if (aSignificand << @intCast(S, typeWidth - shift) != 0) @as(Z, 1) else 0;
-        aSignificand = aSignificand >> @intCast(S, shift | sticky);
-        aExponent = 0;
+        // Result is denormal; the exponent and round/sticky bits are zero.
+        // All we need to do is shift the significand and apply the correct sign.
+        aSignificand >>= @intCast(S, 4 - aExponent);
+        return @bitCast(T, resultSign | aSignificand);
     }
 
     // Low three bits are round, guard, and sticky.
diff --git a/lib/std/special/compiler_rt/addXf3_test.zig b/lib/std/special/compiler_rt/addXf3_test.zig
index a4a18b41e0d7106e4fbb77063dc5e6e99228ed79..f38c9d60185730e85da31e140b5f80f83d09d546 100644
--- a/lib/std/special/compiler_rt/addXf3_test.zig
+++ b/lib/std/special/compiler_rt/addXf3_test.zig
@@ -152,4 +152,6 @@ test "addxf3" {
     try test__addxf3(0x1.0fff_ffff_ffff_fffep+0, 0x1.8p-63, 0x3FFF_8800000000000000); // round down to even
     try test__addxf3(0x1.0fff_ffff_ffff_fffep+0, 0x1.9p-63, 0x3FFF_8800000000000001); // round up
     try test__addxf3(0x1.0fff_ffff_ffff_fffep+0, 0x2.0p-63, 0x3FFF_8800000000000001); // exact
+    try test__addxf3(0x0.ffff_ffff_ffff_fffcp-16382, 0x0.0000_0000_0000_0002p-16382, 0x0000_7FFFFFFFFFFFFFFF); // exact
+    try test__addxf3(0x0.1fff_ffff_ffff_fffcp-16382, 0x0.0000_0000_0000_0002p-16382, 0x0000_0FFFFFFFFFFFFFFF); // exact
 }
diff --git a/lib/std/special/compiler_rt/mulXf3.zig b/lib/std/special/compiler_rt/mulXf3.zig
index befd7b7b43b225ddde8fd3b5b47329236aba283e..147efd0ab510e0ce296749c9fe3f5c8732c151aa 100644
--- a/lib/std/special/compiler_rt/mulXf3.zig
+++ b/lib/std/special/compiler_rt/mulXf3.zig
@@ -152,6 +152,10 @@ fn mulXf3(comptime T: type, a: T, b: T) T {
         const sticky = wideShrWithTruncation(ZSignificand, &productHi, &productLo, shift);
         productLo |= @boolToInt(sticky);
         result = productHi;
+
+        // We include the integer bit so that rounding will carry to the exponent,
+        // but it will be removed later if the result is still denormal
+        if (significandBits != fractionalBits) result |= integerBit;
     } else {
         // Result is normal before rounding; insert the exponent.
         result = productHi & significandMask;
@@ -166,7 +170,11 @@ fn mulXf3(comptime T: type, a: T, b: T) T {
 
     // Restore any explicit integer bit, if it was rounded off
     if (significandBits != fractionalBits) {
-        if ((result >> significandBits) != 0) result |= integerBit;
+        if ((result >> significandBits) != 0) {
+            result |= integerBit;
+        } else {
+            result &= ~integerBit;
+        }
     }
 
     // Insert the sign of the result:
diff --git a/lib/std/special/compiler_rt/mulXf3_test.zig b/lib/std/special/compiler_rt/mulXf3_test.zig
index 2dd345e69ebd0693277faaef94510e31ecf52fc7..6b4f8ca95363f2afcc72a7a1f1207a65b75f56d7 100644
--- a/lib/std/special/compiler_rt/mulXf3_test.zig
+++ b/lib/std/special/compiler_rt/mulXf3_test.zig
@@ -105,6 +105,7 @@ test "multf3" {
 
     try test__multf3(0x1.0000_0000_0000_0000_0000_0000_0001p+0, 0x1.8p+5, 0x4004_8000_0000_0000, 0x0000_0000_0000_0002);
     try test__multf3(0x1.0000_0000_0000_0000_0000_0000_0002p+0, 0x1.8p+5, 0x4004_8000_0000_0000, 0x0000_0000_0000_0003);
+    try test__multf3(2.0, math.floatTrueMin(f128), 0x0000_0000_0000_0000, 0x0000_0000_0000_0002);
 }
 
 const qnan80 = @bitCast(f80, @bitCast(u80, math.nan(f80)) | (1 << (math.floatFractionalBits(f80) - 1)));
@@ -164,4 +165,7 @@ test "mulxf3" {
 
     try test__mulxf3(0x1.0000_0001p+0, 0x1.0000_0001p+0, 0x3FFF_8000_0001_0000_0000); // round down to even
     try test__mulxf3(0x1.0000_0001p+0, 0x1.0000_0001_0002p+0, 0x3FFF_8000_0001_0001_0001); // round up
+    try test__mulxf3(0x0.8000_0000_0000_0000p-16382, 2.0, 0x0001_8000_0000_0000_0000); // denormal -> normal
+    try test__mulxf3(0x0.7fff_ffff_ffff_fffep-16382, 0x2.0000_0000_0000_0008p0, 0x0001_8000_0000_0000_0000); // denormal -> normal
+    try test__mulxf3(0x0.7fff_ffff_ffff_fffep-16382, 0x1.0000_0000_0000_0000p0, 0x0000_3FFF_FFFF_FFFF_FFFF); // denormal -> denormal
 }
-- 
2.54.0


From 28c05b0a53e451cbe06d321f4e2f92cf3b1a7be4 Mon Sep 17 00:00:00 2001
From: Cody Tapscott <topolarity@tapscott.me>
Date: Mon, 25 Apr 2022 17:49:56 -0700
Subject: [PATCH 5/5] compiler_rt: Bypass `fmodx` impl. on stage2

This function is codegen'd incorrectly in stage2, since it fails to
generate the correct soft-float operations. This will be fixed once
issue #11161 is implemented
---
 lib/std/special/compiler_rt.zig | 13 ++++++++++++-
 1 file changed, 12 insertions(+), 1 deletion(-)

diff --git a/lib/std/special/compiler_rt.zig b/lib/std/special/compiler_rt.zig
index 6a44d859c92c6e2317032b727d3c2d2300499fe3..93e0ffbe1a178099bebdcfbf433600755b434c21 100644
--- a/lib/std/special/compiler_rt.zig
+++ b/lib/std/special/compiler_rt.zig
@@ -731,8 +731,13 @@ comptime {
             @export(fmodx, .{ .name = "fmodl", .linkage = linkage });
         } else if (long_double_is_f128) {
             @export(fmodq, .{ .name = "fmodl", .linkage = linkage });
+        } else {
+            @export(fmodl, .{ .name = "fmodl", .linkage = linkage });
+        }
+        if (long_double_is_f80 or builtin.zig_backend == .stage1) {
+            // TODO: https://github.com/ziglang/zig/issues/11161
+            @export(fmodx, .{ .name = "fmodx", .linkage = linkage });
         }
-        @export(fmodx, .{ .name = "fmodx", .linkage = linkage });
         @export(fmodq, .{ .name = "fmodq", .linkage = linkage });
 
         @export(floorf, .{ .name = "floorf", .linkage = linkage });
@@ -889,6 +894,12 @@ fn ceill(x: c_longdouble) callconv(.C) c_longdouble {
 
 const fmodq = @import("compiler_rt/fmodq.zig").fmodq;
 const fmodx = @import("compiler_rt/fmodx.zig").fmodx;
+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));
+}
 
 // Avoid dragging in the runtime safety mechanisms into this .o file,
 // unless we're trying to test this file.
-- 
2.54.0