Add modular arithmetic operators. (#1321)

This adds functions on Integers to compute modular addition, multiplication and inverse.

This is required for the Gaussian Elimination preprocessing pass for BV.

diff --git a/src/util/integer_cln_imp.cpp b/src/util/integer_cln_imp.cpp
index b09d2429cece1616765eabf1ad9bca6437d76d1a..8f98d908ff79070178d34eedb3a839520f875628 100644 (file)
--- a/src/util/integer_cln_imp.cpp
+++ b/src/util/integer_cln_imp.cpp
@@ -155,4 +155,36 @@ Integer Integer::pow(unsigned long int exp) const {
   }
 }
 
+Integer Integer::modAdd(const Integer& y, const Integer& m) const
+{
+  cln::cl_modint_ring ry = cln::find_modint_ring(m.d_value);
+  cln::cl_MI xm = ry->canonhom(d_value);
+  cln::cl_MI ym = ry->canonhom(y.d_value);
+  cln::cl_MI res = xm + ym;
+  return Integer(ry->retract(res));
+}
+
+Integer Integer::modMultiply(const Integer& y, const Integer& m) const
+{
+  cln::cl_modint_ring ry = cln::find_modint_ring(m.d_value);
+  cln::cl_MI xm = ry->canonhom(d_value);
+  cln::cl_MI ym = ry->canonhom(y.d_value);
+  cln::cl_MI res = xm * ym;
+  return Integer(ry->retract(res));
+}
+
+Integer Integer::modInverse(const Integer& m) const
+{
+  PrettyCheckArgument(m > 0, m, "m must be greater than zero");
+  cln::cl_modint_ring ry = cln::find_modint_ring(m.d_value);
+  cln::cl_MI xm = ry->canonhom(d_value);
+  /* normalize to modulo m for coprime check */
+  cln::cl_I x = ry->retract(xm);
+  if (x == 0 || cln::gcd(x, m.d_value) != 1)
+  {
+    return Integer(-1);
+  }
+  cln::cl_MI res = cln::recip(xm);
+  return Integer(ry->retract(res));
+}
 } /* namespace CVC4 */

diff --git a/src/util/integer_cln_imp.h b/src/util/integer_cln_imp.h
index c2791af522f0c238b31b512b15e4fc45a851dfa1..0433494cc5bad8553fafbfecf8d29b4aef3e4747 100644 (file)
--- a/src/util/integer_cln_imp.h
+++ b/src/util/integer_cln_imp.h
@@ -23,6 +23,7 @@
 #include <cln/input.h>
 #include <cln/integer.h>
 #include <cln/integer_io.h>
+#include <cln/modinteger.h>
 #include <iostream>
 #include <limits>
 #include <sstream>
@@ -165,6 +166,7 @@ public:
   Integer operator*(const Integer& y) const {
     return Integer( d_value * y.d_value );
   }
+
   Integer& operator*=(const Integer& y) {
     d_value *= y.d_value;
     return *this;
@@ -347,6 +349,30 @@ public:
     return Integer(result);
   }
 
+  /**
+   * Compute addition of this Integer x + y modulo m.
+   */
+  Integer modAdd(const Integer& y, const Integer& m) const;
+
+  /**
+   * Compute multiplication of this Integer x * y modulo m.
+   */
+  Integer modMultiply(const Integer& y, const Integer& m) const;
+
+  /**
+   * Compute modular inverse x^-1 of this Integer x modulo m with m > 0.
+   * Returns a value x^-1 with 0 <= x^-1 < m such that x * x^-1 = 1 modulo m
+   * if such an inverse exists, and -1 otherwise.
+   *
+   * Such an inverse only exists if
+   *   - x is non-zero
+   *   - x and m are coprime, i.e., if gcd (x, m) = 1
+   *
+   * Note that if x and m are coprime, then x^-1 > 0 if m > 1 and x^-1 = 0
+   * if m = 1 (the zero ring).
+   */
+  Integer modInverse(const Integer& m) const;
+
   /**
    * Return true if *this exactly divides y.
    */

diff --git a/src/util/integer_gmp_imp.cpp b/src/util/integer_gmp_imp.cpp
index e24e8bad105f18f8e2556cc6d783d644d0b92988..4be1da4fea5e01f824e7f493b0ff24cf2bad9777 100644 (file)
--- a/src/util/integer_gmp_imp.cpp
+++ b/src/util/integer_gmp_imp.cpp
@@ -100,4 +100,31 @@ Integer Integer::exactQuotient(const Integer& y) const {
   return Integer( q );
 }
 
+Integer Integer::modAdd(const Integer& y, const Integer& m) const
+{
+  mpz_class res;
+  mpz_add(res.get_mpz_t(), d_value.get_mpz_t(), y.d_value.get_mpz_t());
+  mpz_mod(res.get_mpz_t(), res.get_mpz_t(), m.d_value.get_mpz_t());
+  return Integer(res);
+}
+
+Integer Integer::modMultiply(const Integer& y, const Integer& m) const
+{
+  mpz_class res;
+  mpz_mul(res.get_mpz_t(), d_value.get_mpz_t(), y.d_value.get_mpz_t());
+  mpz_mod(res.get_mpz_t(), res.get_mpz_t(), m.d_value.get_mpz_t());
+  return Integer(res);
+}
+
+Integer Integer::modInverse(const Integer& m) const
+{
+  PrettyCheckArgument(m > 0, m, "m must be greater than zero");
+  mpz_class res;
+  if (mpz_invert(res.get_mpz_t(), d_value.get_mpz_t(), m.d_value.get_mpz_t())
+      == 0)
+  {
+    return Integer(-1);
+  }
+  return Integer(res);
+}
 } /* namespace CVC4 */

diff --git a/src/util/integer_gmp_imp.h b/src/util/integer_gmp_imp.h
index 5f676dbc5237c1e29e044452b090c6f7b159f09b..9d63ea7f0fa2ad54c795e82ebcc64b01258df89f 100644 (file)
--- a/src/util/integer_gmp_imp.h
+++ b/src/util/integer_gmp_imp.h
@@ -293,6 +293,7 @@ public:
       }
     }
   }
+
   /**
    * Returns the quotient according to Boute's Euclidean definition.
    * See the documentation for euclidianQR.
@@ -391,6 +392,30 @@ public:
     return Integer(result);
   }
 
+  /**
+   * Compute addition of this Integer x + y modulo m.
+   */
+  Integer modAdd(const Integer& y, const Integer& m) const;
+
+  /**
+   * Compute multiplication of this Integer x * y modulo m.
+   */
+  Integer modMultiply(const Integer& y, const Integer& m) const;
+
+  /**
+   * Compute modular inverse x^-1 of this Integer x modulo m with m > 0.
+   * Returns a value x^-1 with 0 <= x^-1 < m such that x * x^-1 = 1 modulo m
+   * if such an inverse exists, and -1 otherwise.
+   *
+   * Such an inverse only exists if
+   *   - x is non-zero
+   *   - x and m are coprime, i.e., if gcd (x, m) = 1
+   *
+   * Note that if x and m are coprime, then x^-1 > 0 if m > 1 and x^-1 = 0
+   * if m = 1 (the zero ring).
+   */
+  Integer modInverse(const Integer& m) const;
+
   /**
    * All non-zero integers z, z.divide(0)
    * ! zero.divides(zero)

diff --git a/test/unit/util/integer_black.h b/test/unit/util/integer_black.h
index ac86051dd4c890cd0c6142c15584f2cb95f7e72d..c4c5513638cccad1545901d96247986758f71595 100644 (file)
--- a/test/unit/util/integer_black.h
+++ b/test/unit/util/integer_black.h
@@ -471,4 +471,89 @@ public:
     Integer one_from_string(leadingZeroes,2);
     TS_ASSERT_EQUALS(one, one_from_string);
   }
+
+  void testModAdd()
+  {
+    for (unsigned i = 0; i <= 10; ++i)
+    {
+      for (unsigned j = 0; j <= 10; ++j)
+      {
+        Integer yy;
+        Integer x(i);
+        Integer y = x + j;
+        Integer yp = x.modAdd(j, 3);
+        for (yy = y; yy >= 3; yy -= 3)
+          ;
+        TS_ASSERT(yp == yy);
+        yp = x.modAdd(j, 7);
+        for (yy = y; yy >= 7; yy -= 7)
+          ;
+        TS_ASSERT(yp == yy);
+        yp = x.modAdd(j, 11);
+        for (yy = y; yy >= 11; yy -= 11)
+          ;
+        TS_ASSERT(yp == yy);
+      }
+    }
+  }
+
+  void testModMultiply()
+  {
+    for (unsigned i = 0; i <= 10; ++i)
+    {
+      for (unsigned j = 0; j <= 10; ++j)
+      {
+        Integer yy;
+        Integer x(i);
+        Integer y = x * j;
+        Integer yp = x.modMultiply(j, 3);
+        for (yy = y; yy >= 3; yy -= 3)
+          ;
+        TS_ASSERT(yp == yy);
+        yp = x.modMultiply(j, 7);
+        for (yy = y; yy >= 7; yy -= 7)
+          ;
+        TS_ASSERT(yp == yy);
+        yp = x.modMultiply(j, 11);
+        for (yy = y; yy >= 11; yy -= 11)
+          ;
+        TS_ASSERT(yp == yy);
+      }
+    }
+  }
+
+  void testModInverse()
+  {
+    for (unsigned i = 0; i <= 10; ++i)
+    {
+      Integer x(i);
+      Integer inv = x.modInverse(3);
+      if (i == 0 || i == 3 || i == 6 || i == 9)
+      {
+        TS_ASSERT(inv == -1); /* no inverse */
+      }
+      else
+      {
+        TS_ASSERT(x.modMultiply(inv, 3) == 1);
+      }
+      inv = x.modInverse(7);
+      if (i == 0 || i == 7)
+      {
+        TS_ASSERT(inv == -1); /* no inverse */
+      }
+      else
+      {
+        TS_ASSERT(x.modMultiply(inv, 7) == 1);
+      }
+      inv = x.modInverse(11);
+      if (i == 0)
+      {
+        TS_ASSERT(inv == -1); /* no inverse */
+      }
+      else
+      {
+        TS_ASSERT(x.modMultiply(inv, 11) == 1);
+      }
+    }
+  }
 };