From 4e505b2b85a815d76f6bf99c544009745abf80b3 Mon Sep 17 00:00:00 2001
From: "Yukihiro \"Matz\" Matsumoto" <matz@ruby.or.jp>
Date: Wed, 23 Jul 2025 09:31:36 +0900
Subject: [PATCH] mruby-bigint: fix multiplication commutativity bug

Fixed non-commutative multiplication bug where operands with different
limb counts would produce different results based on order (a*b \!= b*a).

Root cause was asymmetric carry propagation in the multiplication algorithm.
The fix ensures consistent operand ordering by always processing the smaller
operand first in the nested loops, making multiplication truly commutative.

Also fixed division algorithm quotient allocation and qhat refinement.

Co-authored-by: Claude <noreply@anthropic.com>
---
 mrbgems/mruby-bigint/core/bigint.c  | 74 ++++++++++++++++++++++-------
 mrbgems/mruby-bigint/test/bigint.rb |  7 +++
 2 files changed, 64 insertions(+), 17 deletions(-)

diff --git a/mrbgems/mruby-bigint/core/bigint.c b/mrbgems/mruby-bigint/core/bigint.c
index 5354cd7fd..d0d460f84 100644
--- a/mrbgems/mruby-bigint/core/bigint.c
+++ b/mrbgems/mruby-bigint/core/bigint.c
@@ -181,8 +181,8 @@ digits(mpz_t *x)
   size_t i;
 
   if (x->sz == 0) return 0;
-  for (i = x->sz - 1; x->p[i] == 0; i--)
-    if (i == 0) break;
+  for (i = x->sz - 1; x->p[i] == 0 && i > 0; i--)
+    ;
   return i+1;
 }
 
@@ -415,24 +415,41 @@ mpz_mul(mrb_state *mrb, mpz_t *ww, mpz_t *u, mpz_t *v)
     return;
   }
 
+  // Ensure consistent operand ordering: smaller operand as second argument
+  mpz_t *first, *second;
+  if (u->sz <= v->sz) {
+    first = u;
+    second = v;
+  }
+  else {
+    first = v;
+    second = u;
+  }
+
   mpz_t w;
   mpz_init(mrb, &w);
-  mpz_realloc(mrb, &w, u->sz + v->sz);
+  mpz_realloc(mrb, &w, first->sz + second->sz + 1);
 
-  for (size_t j = 0; j < u->sz; j++) {
-    size_t i;
-    mp_dbl_limb cc = (mp_limb)0;
-    mp_limb u0 = u->p[j];
+  // Standard multiplication algorithm with consistent operand order
+  for (size_t j = 0; j < first->sz; j++) {
+    mp_limb u0 = first->p[j];
     if (u0 == 0) continue;
-    for (i = 0; i < v->sz; i++) {
-      mp_limb v0 = v->p[i];
-      if (v0 == 0) continue;
+
+    mp_dbl_limb cc = 0;
+    size_t i;
+    for (i = 0; i < second->sz; i++) {
+      mp_limb v0 = second->p[i];
       cc += (mp_dbl_limb)w.p[i + j] + (mp_dbl_limb)u0 * (mp_dbl_limb)v0;
       w.p[i + j] = LOW(cc);
       cc = HIGH(cc);
     }
-    if (cc) {
-      w.p[i + j] = (mp_limb)cc;
+
+    // Propagate carries
+    while (cc && (i + j) < w.sz) {
+      cc += (mp_dbl_limb)w.p[i + j];
+      w.p[i + j] = LOW(cc);
+      cc = HIGH(cc);
+      i++;
     }
   }
 
@@ -666,17 +683,40 @@ udiv(mrb_state *mrb, mpz_t *qq, mpz_t *rr, mpz_t *xx, mpz_t *yy)
   ulshift(mrb, &x, xx, ns);
   ulshift(mrb, &y, yy, ns);
   size_t xd = digits(&x);
-  mpz_realloc(mrb, &q, xd);
+  mpz_realloc(mrb, &q, xd-yd+1);  // Quotient has xd-yd+1 digits maximum
   mp_dbl_limb z = y.p[yd-1];
   if (xd>=yd) {
     for (size_t j=xd-yd;; j--) {
       mp_dbl_limb_signed b=0;
       mp_dbl_limb qhat;
 
-      if (j+yd == xd)
-        qhat = x.p[j+yd-1] / z;
-      else
-        qhat = (((mp_dbl_limb)x.p[j+yd] << DIG_SIZE) + x.p[j+yd-1]) / z;
+
+
+      mp_dbl_limb rhat;
+      if (j+yd == xd) {
+        // Treat missing high limb as 0 and use the same two-limb formula
+        qhat = (((mp_dbl_limb)0 << DIG_SIZE) + x.p[j+yd-1]) / z;
+        rhat = (((mp_dbl_limb)0 << DIG_SIZE) + x.p[j+yd-1]) % z;
+      }
+      else {
+        mp_dbl_limb dividend = ((mp_dbl_limb)x.p[j+yd] << DIG_SIZE) + x.p[j+yd-1];
+        qhat = dividend / z;
+        rhat = dividend % z;
+      }
+
+      // Knuth's qhat refinement step - essential to prevent overestimation
+      if (yd > 1) {  // Apply refinement for all iterations, including j=0
+        mp_dbl_limb left_side = qhat * y.p[yd-2];
+        mp_dbl_limb right_side = (rhat << DIG_SIZE) + (j+yd-2 < x.sz ? x.p[j+yd-2] : 0);
+
+        while (qhat >= ((mp_dbl_limb)1 << DIG_SIZE) || (left_side > right_side)) {
+          qhat--;
+          rhat += z;
+          if (rhat >= ((mp_dbl_limb)1 << DIG_SIZE)) break;
+          left_side = qhat * y.p[yd-2];
+          right_side = (rhat << DIG_SIZE) + (j+yd-2 < x.sz ? x.p[j+yd-2] : 0);
+        }
+      }
       if (qhat) {
         size_t i;
 
diff --git a/mrbgems/mruby-bigint/test/bigint.rb b/mrbgems/mruby-bigint/test/bigint.rb
index 3280c5004..07f6c783e 100644
--- a/mrbgems/mruby-bigint/test/bigint.rb
+++ b/mrbgems/mruby-bigint/test/bigint.rb
@@ -42,6 +42,13 @@ assert 'Bigint *' do
   assert_equal(-1361129467683753853853498429727072845824, -n * n)
   assert_equal(-1361129467683753853853498429727072845824, n * -n)
   assert_equal 1361129467683753853853498429727072845824, -n * -n
+
+  # Test multiplication commutativity for large numbers with different limb counts
+  # This test specifically targets the bug where operands with different
+  # limb counts would produce different results based on order
+  a = (2**512) - 1      # 16 limbs
+  b = 26815615859885194199148049996411692254958731641184786755447122887443528060147093953603748596333806855380063716372972101707507765623893139892867298012168194  # 17 limbs
+  assert_equal(a * b, b * a)
 end
 
 assert 'Bigint /' do