From a9cf481470c324a04f2254c5745eee26c45cb309 Mon Sep 17 00:00:00 2001 From: Aina Niemetz Date: Thu, 9 Nov 2017 04:47:02 -0800 Subject: [PATCH] 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. --- src/util/integer_cln_imp.cpp | 32 +++++++++++++ src/util/integer_cln_imp.h | 26 +++++++++++ src/util/integer_gmp_imp.cpp | 27 +++++++++++ src/util/integer_gmp_imp.h | 25 ++++++++++ test/unit/util/integer_black.h | 85 ++++++++++++++++++++++++++++++++++ 5 files changed, 195 insertions(+) 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 #include #include +#include #include #include #include @@ -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 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. @@ -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 ac86051dd..c4c551363 100644 --- 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); + } + } + } }; -- 2.30.2