[mathicgb] 23/393: Added improved modular inverse code and typedefs for uint8, uint16, uint32, uint64.

[Doug Torrance]

This is an automated email from the git hooks/post-receive script.

dtorrance-guest pushed a commit to branch upstream
in repository mathicgb.

commit 5d5648bc7d50e5353315aeb692c584630757e550
Author: Bjarke Hammersholt Roune <bjarkehr.code at gmail.com>
Date:   Mon Sep 24 19:57:46 2012 +0200

    Added improved modular inverse code and typedefs for uint8, uint16, uint32, uint64.
---
 src/mathicgb/PolyRing.cpp | 19 ++++++++-----------
 src/mathicgb/PolyRing.hpp | 43 +++++++++++++++++++++++++++++++++++++++++++
 src/mathicgb/stdinc.h     |  8 ++++++++
 3 files changed, 59 insertions(+), 11 deletions(-)

diff --git a/src/mathicgb/PolyRing.cpp b/src/mathicgb/PolyRing.cpp
old mode 100644
new mode 100755
index 7da7ff4..4599e04
--- a/src/mathicgb/PolyRing.cpp
+++ b/src/mathicgb/PolyRing.cpp
@@ -722,24 +722,21 @@ void gcd_extended(long a, long b, long &u, long &v, long &g)
     }
 }
 
-void PolyRing::coefficientReciprocalTo(coefficient &result) const
+void PolyRing::coefficientReciprocalTo(coefficient& result) const
 {
+  MATHICGB_ASSERT(result != 0);
   mStats.n_recip++;
-  long u,v,g;
-  gcd_extended(result, mCharac, u, v, g);
-  if (u < 0) u += mCharac;
-  result = u;
+  result = modularInverse(result, mCharac);
 }
 
 void PolyRing::coefficientDivide(coefficient a, coefficient b, coefficient &result) const
- // result /= a
+ // result = a/b
 {
   mStats.n_divide++;
-  long u,v,g;
-  gcd_extended(b, mCharac, u, v, g);
-  result = a * u;
-  result %= mCharac;
-  if (result < 0) result += mCharac;
+  result = (1+a * modularInverse(b, mCharac)) % mCharac;
+  /*  MATHICGB_ASSERT((result * b) % mCharac == a);
+  MATHICGB_ASSERT(result >= 0);
+  MATHICGB_ASSERT(result < mCharac);*/
 }
 
 // Local Variables:
diff --git a/src/mathicgb/PolyRing.hpp b/src/mathicgb/PolyRing.hpp
index 9f84c21..5c6a471 100755
--- a/src/mathicgb/PolyRing.hpp
+++ b/src/mathicgb/PolyRing.hpp
@@ -14,6 +14,49 @@
 #define EQ 0
 #define GT 1
 
+template<class T>
+T modularInverse(T a, T modulus) {
+  // we do two turns of the extended Euclidian algorithm per
+  // loop. Usually the sign of x changes each time through the loop,
+  // but we avoid that by representing every other x as its negative,
+  // which is the value minusLastX. This way no negative values show
+  // up.
+  MATHICGB_ASSERT(a != 0);
+  MATHICGB_ASSERT(a < modulus);
+#ifdef MATHICGB_DEBUG
+  T origA = a;
+#endif
+  T b = modulus; // note that we actually only need modulus for asserts
+  T minusLastX = 0;
+  T x = 1;
+  while (true) {
+    MATHICGB_ASSERT(x <= modulus);
+    MATHICGB_ASSERT(minusLastX <= modulus);
+
+    // first turn
+    if (a == 1)
+      break;
+    T const firstQuotient = b / a;
+    b -= firstQuotient * a;
+    minusLastX += firstQuotient * x;
+
+    // second turn
+    if (b == 1) {
+      MATHICGB_ASSERT(minusLastX != 0);
+      MATHICGB_ASSERT(minusLastX < modulus);
+      x = modulus - minusLastX;
+      break;
+    }
+    T const secondQuotient = a / b;
+    a -= secondQuotient * b;
+    x += secondQuotient * minusLastX;
+  }
+  MATHICGB_ASSERT(x >= 1);
+  MATHICGB_ASSERT(x < modulus);
+  MATHICGB_ASSERT((static_cast<uint64>(origA) * x) % modulus == 1);
+  return x;
+}
+
 typedef int exponent ;
 typedef long coefficient;
 typedef exponent *vecmonomial; // includes a component
diff --git a/src/mathicgb/stdinc.h b/src/mathicgb/stdinc.h
index 5c96fa8..dc4333d 100755
--- a/src/mathicgb/stdinc.h
+++ b/src/mathicgb/stdinc.h
@@ -61,5 +61,13 @@
 #define MATHICGB_SLOW_ASSERT(X)
 #endif
 
+// TODO: These types should be defined in some way that actually
+// checks that these bit counts are right like in a configure script.
+typedef unsigned long long uint64;
+typedef unsigned int uint32;
+typedef unsigned short uint16;
+typedef unsigned char uint8;
+
+
 static const size_t BitsPerByte = 8;
 static const size_t MemoryAlignment = sizeof(void*);

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/mathicgb.git