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.

4 files changed, 110 insertions, 0 deletions

diff --git a/src/util/integer_cln_imp.cpp b/src/util/integer_cln_imp.cpp
index b09d2429c..8f98d908f 100644
--- 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 c2791af52..0433494cc 100644
--- 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;
@@ -348,6 +350,30 @@ public:
   }
 
   /**
+   * 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.
    */
   bool divides(const Integer& y) const {
diff --git a/src/util/integer_gmp_imp.cpp b/src/util/integer_gmp_imp.cpp
index e24e8bad1..4be1da4fe 100644
--- 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 5f676dbc5..9d63ea7f0 100644
--- 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.
@@ -392,6 +393,30 @@ public:
   }
 
   /**
+   * 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)
    */