Fix BigInt::div bugs and preclude aliasing via its interface

src/big-int/bigint-func.cc

@@ -94,16 +94,16 @@ modinv (const BigInt &a, const BigInt &m)
 {
   BigInt j (1), i (0);
   BigInt b (m), c (a);
-  BigInt x, y;
+  BigInt y;
   while (!c.is_zero())
     {
-      BigInt::div (b, c, x, y);
+      BigInt::divisiont qr = BigInt::div(b, c);
       b = c;
-      c = y;
+      c = qr.remainder;
       y = j;
 
-      // j = i - j * x; trading clarity for efficiency.
-      j *= x;
+      // j = i - j * qr.quotient; trading clarity for efficiency.
+      j *= qr.quotient;
       j -= i;
       j.negate();
 

src/big-int/bigint.cc

@@ -1033,9 +1033,12 @@ BigInt::operator*= (BigInt const &y)
 
 // Division method returning both quotient and remainder.
 
-void
-BigInt::div (BigInt const &x, BigInt const &y, BigInt &q, BigInt &r)
+BigInt::divisiont BigInt::div(BigInt const &x, BigInt const &y)
 {
+  divisiont result;
+  BigInt &q = result.quotient;
+  BigInt &r = result.remainder;
+
   // Eliminate some trivial cases.
   int cmp = x.ucompare (y);
   if (cmp < 0)
@@ -1044,7 +1047,7 @@ BigInt::div (BigInt const &x, BigInt const &y, BigInt &q, BigInt &r)
       r = x;
       q.length = 0;
       q.positive = true;
-      return;
+      return result;
     }
   if (cmp == 0)
     {
@@ -1055,13 +1058,13 @@ BigInt::div (BigInt const &x, BigInt const &y, BigInt &q, BigInt &r)
       q.digit[0] = 1;
       if (!y.positive)
 	q.negate();
-      return;
+      return result;
     }
   if (y.length == 0)
     {
     zero:
       error ("Division by zero.");
-      return;
+      return result;
     }
   if (x.is_ulong())
     {
@@ -1076,7 +1079,7 @@ BigInt::div (BigInt const &x, BigInt const &y, BigInt &q, BigInt &r)
   else if (y.length == 1)
     {
       // This digit_div() transforms the dividend into the quotient.
-      q = y;
+      q = x;
       r.digit[0] = digit_div (q.digit, q.length, y.digit[0]);
       r.length = r.digit[0] ? 1 : 0;
     }
@@ -1105,9 +1108,11 @@ BigInt::div (BigInt const &x, BigInt const &y, BigInt &q, BigInt &r)
       if (a[al-1] >= b[bl-1])
 	a[al++] = 0;
 
-      // Prepare q for receiving the quotient.
+      // Prepare q for receiving the quotient. Set the length only
+      // after resizing: resize copies `length` digits from the old
+      // buffer, which may be shorter than the new length.
+      q.resize(al - bl);
       q.length = al - bl;
-      q.resize (q.length);
 
       // Divide.
       digit_div (a, b, bl, q.digit, q.length);
@@ -1117,13 +1122,14 @@ BigInt::div (BigInt const &x, BigInt const &y, BigInt &q, BigInt &r)
       if (scale != 1)
 	digit_div (a, al, scale);
       if (al && a[al - 1] == 0) --al;
+      r.resize(al);
       r.length = al;
-      r.resize (r.length);
       memcpy (r.digit, a, al * sizeof (onedig_t));
     }
   q.adjust();
   q.positive = q.length == 0 || x.positive == y.positive;
   r.positive = r.length == 0 || x.positive;
+  return result;
 }
 
 BigInt &

src/big-int/bigint.hh

@@ -270,8 +270,16 @@ public:
   BigInt &operator++ ()	{ return operator+=(1); } // preincrement
   BigInt &operator-- ()	{ return operator-=(1); } // predecrement
 
-  static void div (BigInt const &, BigInt const &,
-		   BigInt &quot, BigInt &rem) _fasta;
+  // Result of the combined quotient/remainder computation. Returning
+  // fresh objects (rather than filling in caller-provided output
+  // arguments) makes aliasing between inputs and outputs impossible by
+  // construction. The struct is defined after this class.
+  struct divisiont;
+
+  // Combined quotient/remainder computation: quotient = x / y,
+  // remainder = x % y (truncated division, the remainder takes the
+  // sign of x).
+  static divisiont div(BigInt const &, BigInt const &) _fasta;
 
   // Returns the largest x such that 2^x <= abs() or 0 if input is 0
   // Not part of original BigInt.
@@ -290,6 +298,11 @@ public:
   }
 };
 
+struct BigInt::divisiont
+{
+  BigInt quotient;
+  BigInt remainder;
+};
 
 // Functions on BigInt. Implementations in bigint-func.cc.
 

unit/big-int/big-int.cpp

@@ -75,6 +75,58 @@ TEST_CASE("arbitrary precision integers", "[core][big-int][bigint]")
     REQUIRE(to_string(i) == "1");
   }
 
+  // =====================================================================
+  // Combined division/remainder.
+  // =====================================================================
+  // BigInt::div used to compute the quotient from the divisor instead of
+  // the dividend when the dividend exceeded 64 bits and the divisor fit a
+  // single digit (32 bits).
+  SECTION("combined division/remainder with single-digit divisor")
+  {
+    const BigInt dividend = pow(BigInt{2}, 89) - 1;
+    const BigInt divisor{97};
+
+    BigInt::divisiont qr = BigInt::div(dividend, divisor);
+    REQUIRE(to_string(qr.quotient) == "6381134223120516880923320");
+    REQUIRE(to_string(qr.remainder) == "71");
+    REQUIRE(qr.quotient * divisor + qr.remainder == dividend);
+
+    // Negative dividend: truncated division, remainder takes the
+    // dividend's sign.
+    qr = BigInt::div(-dividend, divisor);
+    REQUIRE(to_string(qr.quotient) == "-6381134223120516880923320");
+    REQUIRE(to_string(qr.remainder) == "-71");
+    REQUIRE(qr.quotient * divisor + qr.remainder == -dividend);
+
+    // Negative divisor: the quotient flips sign, the remainder keeps
+    // the dividend's sign.
+    qr = BigInt::div(dividend, -divisor);
+    REQUIRE(to_string(qr.quotient) == "-6381134223120516880923320");
+    REQUIRE(to_string(qr.remainder) == "71");
+    REQUIRE(qr.quotient * -divisor + qr.remainder == dividend);
+
+    // Both negative: the quotient is positive, the remainder keeps the
+    // dividend's sign.
+    qr = BigInt::div(-dividend, -divisor);
+    REQUIRE(to_string(qr.quotient) == "6381134223120516880923320");
+    REQUIRE(to_string(qr.remainder) == "-71");
+    REQUIRE(qr.quotient * -divisor + qr.remainder == -dividend);
+
+    // Zero remainder exercises the remainder.length = 0 path.
+    qr = BigInt::div(dividend * divisor, divisor);
+    REQUIRE(qr.quotient == dividend);
+    REQUIRE(qr.remainder.is_zero());
+
+    // A divisor of more than one digit exercises the long-division
+    // branch, which is otherwise only tested via operator/= and
+    // operator%= that have separate copies of that code.
+    const BigInt divisor2 = pow(BigInt{2}, 41) + 3;
+    qr = BigInt::div(dividend, divisor2);
+    REQUIRE(to_string(qr.quotient) == "281474976710272");
+    REQUIRE(to_string(qr.remainder) == "1151");
+    REQUIRE(qr.quotient * divisor2 + qr.remainder == dividend);
+  }
+
   // =====================================================================
   // Test cases from the clisp test suite in number.tst.
   // =====================================================================
@@ -174,15 +226,15 @@ TEST_CASE("arbitrary precision integers", "[core][big-int][bigint]")
         REQUIRE(m == em);
 
         // Also try the method returning both.
-        BigInt::div(a, b, r, m);
+        BigInt::divisiont qr = BigInt::div(a, b);
         // Again, transform to floored divide.
-        if(!m.is_zero() && a.is_positive() != b.is_positive())
+        if(!qr.remainder.is_zero() && a.is_positive() != b.is_positive())
         {
-          r -= 1;
-          m += b;
+          qr.quotient -= 1;
+          qr.remainder += b;
         }
-        REQUIRE(r == er);
-        REQUIRE(m == em);
+        REQUIRE(qr.quotient == er);
+        REQUIRE(qr.remainder == em);
       }
       }
     }