From: Wojciech Matyjewicz <wmatyjewicz@fastmail.fm>
Date: Mon, 23 Jun 2008 19:39:50 +0000 (+0000)
Subject: First step to fix PR2088. Implement routine to compute the
X-Git-Url: http://demsky.eecs.uci.edu/git/?a=commitdiff_plain;h=300c6c5167d2869d1568d783d0e3e48bf4b03a6c;p=oota-llvm.git

First step to fix PR2088. Implement routine to compute the
multiplicative inverse of a given number. Modify udivrem to allow input and
output pairs of arguments to overlap. Patch is based on the work by Chandler
Carruth.


git-svn-id: https://llvm.org/svn/llvm-project/llvm/trunk@52638 91177308-0d34-0410-b5e6-96231b3b80d8
---

diff --git a/include/llvm/ADT/APInt.h b/include/llvm/ADT/APInt.h
index d2a1baff244..c98ff18d101 100644
--- a/include/llvm/ADT/APInt.h
+++ b/include/llvm/ADT/APInt.h
@@ -668,9 +668,11 @@ public:
     return this->urem(RHS);
   }
 
-  /// Sometimes it is convenient to divide two APInt values and obtain both
-  /// the quotient and remainder. This function does both operations in the
-  /// same computation making it a little more efficient.
+  /// Sometimes it is convenient to divide two APInt values and obtain both the
+  /// quotient and remainder. This function does both operations in the same
+  /// computation making it a little more efficient. The pair of input arguments
+  /// may overlap with the pair of output arguments. It is safe to call
+  /// udivrem(X, Y, X, Y), for example.
   /// @brief Dual division/remainder interface.
   static void udivrem(const APInt &LHS, const APInt &RHS, 
                       APInt &Quotient, APInt &Remainder);
@@ -1107,6 +1109,9 @@ public:
     return *this;
   }
 
+  /// @returns the multiplicative inverse for a given modulo.
+  APInt multiplicativeInverse(const APInt& modulo) const;
+
   /// @}
   /// @name Building-block Operations for APInt and APFloat
   /// @{
@@ -1301,7 +1306,7 @@ inline uint32_t logBase2(const APInt& APIVal) {
 }
 
 /// GreatestCommonDivisor - This function returns the greatest common
-/// divisor of the two APInt values using Enclid's algorithm.
+/// divisor of the two APInt values using Euclid's algorithm.
 /// @returns the greatest common divisor of Val1 and Val2
 /// @brief Compute GCD of two APInt values.
 APInt GreatestCommonDivisor(const APInt& Val1, const APInt& Val2);
diff --git a/lib/Support/APInt.cpp b/lib/Support/APInt.cpp
index 38b379005af..b0577b3680d 100644
--- a/lib/Support/APInt.cpp
+++ b/lib/Support/APInt.cpp
@@ -1426,6 +1426,50 @@ APInt APInt::sqrt() const {
   return x_old + 1;
 }
 
+/// Computes the multiplicative inverse of this APInt for a given modulo. The
+/// iterative extended Euclidean algorithm is used to solve for this value,
+/// however we simplify it to speed up calculating only the inverse, and take
+/// advantage of div+rem calculations. We also use some tricks to avoid copying
+/// (potentially large) APInts around.
+APInt APInt::multiplicativeInverse(const APInt& modulo) const {
+  assert(ult(modulo) && "This APInt must be smaller than the modulo");
+
+  // Using the properties listed at the following web page (accessed 06/21/08):
+  //   http://www.numbertheory.org/php/euclid.html
+  // (especially the properties numbered 3, 4 and 9) it can be proved that
+  // BitWidth bits suffice for all the computations in the algorithm implemented
+  // below. More precisely, this number of bits suffice if the multiplicative
+  // inverse exists, but may not suffice for the general extended Euclidean
+  // algorithm.
+
+  APInt r[2] = { modulo, *this };
+  APInt t[2] = { APInt(BitWidth, 0), APInt(BitWidth, 1) };
+  APInt q(BitWidth, 0);
+  
+  unsigned i;
+  for (i = 0; r[i^1] != 0; i ^= 1) {
+    // An overview of the math without the confusing bit-flipping:
+    // q = r[i-2] / r[i-1]
+    // r[i] = r[i-2] % r[i-1]
+    // t[i] = t[i-2] - t[i-1] * q
+    udivrem(r[i], r[i^1], q, r[i]);
+    t[i] -= t[i^1] * q;
+  }
+
+  // If this APInt and the modulo are not coprime, there is no multiplicative
+  // inverse, so return 0. We check this by looking at the next-to-last
+  // remainder, which is the gcd(*this,modulo) as calculated by the Euclidean
+  // algorithm.
+  if (r[i] != 1)
+    return APInt(BitWidth, 0);
+
+  // The next-to-last t is the multiplicative inverse.  However, we are
+  // interested in a positive inverse. Calcuate a positive one from a negative
+  // one if necessary. A simple addition of the modulo suffices because
+  // abs(t[i]) is known to less than *this/2 (see the link above).
+  return t[i].isNegative() ? t[i] + modulo : t[i];
+}
+
 /// Implementation of Knuth's Algorithm D (Division of nonnegative integers)
 /// from "Art of Computer Programming, Volume 2", section 4.3.1, p. 272. The
 /// variables here have the same names as in the algorithm. Comments explain
@@ -1882,13 +1926,10 @@ void APInt::udivrem(const APInt &LHS, const APInt &RHS,
   
   if (lhsWords == 1 && rhsWords == 1) {
     // There is only one word to consider so use the native versions.
-    if (LHS.isSingleWord()) {
-      Quotient = APInt(LHS.getBitWidth(), LHS.VAL / RHS.VAL);
-      Remainder = APInt(LHS.getBitWidth(), LHS.VAL % RHS.VAL);
-    } else {
-      Quotient = APInt(LHS.getBitWidth(), LHS.pVal[0] / RHS.pVal[0]);
-      Remainder = APInt(LHS.getBitWidth(), LHS.pVal[0] % RHS.pVal[0]);
-    }
+    uint64_t lhsValue = LHS.isSingleWord() ? LHS.VAL : LHS.pVal[0];
+    uint64_t rhsValue = RHS.isSingleWord() ? RHS.VAL : RHS.pVal[0];
+    Quotient = APInt(LHS.getBitWidth(), lhsValue / rhsValue);
+    Remainder = APInt(LHS.getBitWidth(), lhsValue % rhsValue);
     return;
   }