diff --git a/src/ir.cpp b/src/ir.cpp
index ba2a9b3fe2..950d051492 100644
--- a/src/ir.cpp
+++ b/src/ir.cpp
@@ -7710,6 +7710,14 @@ static void float_rem(ConstExprValue *out_val, ConstExprValue *op1, ConstExprVal
     }
 }
 
+// c = a - b * trunc(a / b)
+static void zig_f128M_mod(const float128_t* a, const float128_t* b, float128_t* c) {
+    f128M_div(a, b, c);
+    f128M_roundToInt(c, softfloat_round_min, true, c);
+    f128M_mul(b, c, c);
+    f128M_sub(a, c, c);
+}
+
 static void float_mod(ConstExprValue *out_val, ConstExprValue *op1, ConstExprValue *op2) {
     assert(op1->type == op2->type);
     out_val->type = op1->type;
@@ -7724,9 +7732,7 @@ static void float_mod(ConstExprValue *out_val, ConstExprValue *op1, ConstExprVal
                 out_val->data.x_f64 = fmod(fmod(op1->data.x_f64, op2->data.x_f64) + op2->data.x_f64, op2->data.x_f64);
                 return;
             case 128:
-                f128M_rem(&op1->data.x_f128, &op2->data.x_f128, &out_val->data.x_f128);
-                f128M_add(&out_val->data.x_f128, &op2->data.x_f128, &out_val->data.x_f128);
-                f128M_rem(&out_val->data.x_f128, &op2->data.x_f128, &out_val->data.x_f128);
+                zig_f128M_mod(&op1->data.x_f128, &op2->data.x_f128, &out_val->data.x_f128);
                 return;
             default:
                 zig_unreachable();
diff --git a/test/cases/math.zig b/test/cases/math.zig
index 0bf99cff0e..08388d3df8 100644
--- a/test/cases/math.zig
+++ b/test/cases/math.zig
@@ -434,6 +434,20 @@ test "comptime float rem int" {
     }
 }
 
+test "remainder division" {
+    comptime remdiv(f32);
+    comptime remdiv(f64);
+    comptime remdiv(f128);
+    remdiv(f32);
+    remdiv(f64);
+    remdiv(f128);
+}
+
+fn remdiv(comptime T: type) void {
+    assert(T(1) == T(1) % T(2));
+    assert(T(1) == T(7) % T(3));
+}
+
 test "@sqrt" {
     testSqrt(f64, 12.0);
     comptime testSqrt(f64, 12.0);