feat(numerics): add finite field Z/pZ

Signed-off-by: Leonardo de Moura <leonardo@microsoft.com>

src/tests/util/numerics/CMakeLists.txt

@@ -28,3 +28,6 @@ add_test(numeric_traits ${CMAKE_CURRENT_BINARY_DIR}/numeric_traits)
 add_executable(gcd gcd.cpp)
 target_link_libraries(gcd ${EXTRA_LIBS})
 add_test(gcd ${CMAKE_CURRENT_BINARY_DIR}/gcd)
+add_executable(zpz zpz.cpp)
+target_link_libraries(zpz ${EXTRA_LIBS})
+add_test(zpz ${CMAKE_CURRENT_BINARY_DIR}/zpz)

src/tests/util/numerics/zpz.cpp

@@ -0,0 +1,116 @@
+/*
+Copyright (c) 2013 Microsoft Corporation. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+
+Author: Leonardo de Moura
+*/
+#include <iostream>
+#include <sstream>
+#include "util/test.h"
+#include "util/numerics/zpz.h"
+using namespace lean;
+
+static void tst1() {
+    zpz z(1, 7);
+    lean_assert(z.p() == 7);
+    lean_assert(z.get_unsigned_int() == 1);
+    lean_assert(z.hash() == 1);
+    z += 2;
+    lean_assert(z == 3);
+    z -= 1;
+    lean_assert(z == zpz(2, 7));
+    z -= zpz(3, 7);
+    lean_assert(z == 6);
+    z.neg();
+    lean_assert(z == 1);
+    z += zpz(1, 7);
+    z *= zpz(2, 7);
+    lean_assert(z == 4);
+    z *= 2;
+    lean_assert(z == 1);
+    z.inv();
+    lean_assert(z == 1);
+    z /= 3;
+    lean_assert(z == 5);
+    z /= zpz(2, 7);
+    lean_assert(z == 6);
+    z.neg();
+    lean_assert(z == 1);
+    z++;
+    lean_assert(z == 2);
+    z--;
+    --z;
+    z--;
+    lean_assert(z == 6);
+    ++z;
+    lean_assert(z == 0);
+    z = 4;
+    lean_assert(z == 4);
+    z.set_p(3);
+    lean_assert(z == 1);
+    lean_assert(z.p() == 3);
+    z.set_p(7);
+    lean_assert(z.p() == 7);
+    zpz w(3, 13);
+    swap(z, w);
+    lean_assert(z == 3 && z.p() == 13);
+    lean_assert(w == 1 && w.p() == 7);
+    w = z;
+    lean_assert(w == 3 && w.p() == 13);
+    lean_assert(zpz(3, 7) == zpz(3, 13));
+    lean_assert(zpz(3, 7) != zpz(4, 7));
+    lean_assert(zpz(3, 7) < zpz(5, 13));
+    lean_assert(zpz(5, 7) > zpz(3, 13));
+    lean_assert(zpz(5, 7) >= zpz(3, 13));
+    lean_assert(zpz(5, 7) >= zpz(5, 13));
+    lean_assert(zpz(3, 7) < zpz(5, 13));
+    lean_assert(zpz(3, 7) <= zpz(5, 13));
+    lean_assert(zpz(5, 7) <= zpz(5, 13));
+
+    lean_assert(zpz(3, 7) == 3);
+    lean_assert(zpz(3, 7) != 4);
+    lean_assert(zpz(3, 7) < 5);
+    lean_assert(zpz(5, 7) > 3);
+    lean_assert(zpz(5, 7) >= 3);
+    lean_assert(zpz(5, 7) >= 5);
+    lean_assert(zpz(3, 7) < 5);
+    lean_assert(zpz(3, 7) <= 5);
+    lean_assert(zpz(5, 7) <= 5);
+
+    lean_assert(3 == zpz(3, 13));
+    lean_assert(3 != zpz(4, 7));
+    lean_assert(3 < zpz(5, 13));
+    lean_assert(5 > zpz(3, 13));
+    lean_assert(5 >= zpz(3, 13));
+    lean_assert(5 >= zpz(5, 13));
+    lean_assert(3 < zpz(5, 13));
+    lean_assert(3 <= zpz(5, 13));
+    lean_assert(5 <= zpz(5, 13));
+
+    lean_assert(zpz(1, 7) + zpz(4, 7) == 5);
+    lean_assert(zpz(1, 7) - zpz(3, 7) == 5);
+    lean_assert(zpz(2, 7) * zpz(3, 7) == 6);
+    lean_assert(zpz(2, 7) * zpz(4, 7) == 1);
+    lean_assert(zpz(1, 7) / zpz(3, 7) == 5);
+
+    lean_assert(zpz(1, 7) + 4 == 5);
+    lean_assert(zpz(1, 7) - 3 == 5);
+    lean_assert(zpz(2, 7) * 3 == 6);
+    lean_assert(zpz(2, 7) * 4 == 1);
+    lean_assert(zpz(1, 7) / 3 == 5);
+
+    lean_assert(1 + zpz(4, 7) == 5);
+    lean_assert(1 - zpz(3, 7) == 5);
+    lean_assert(2 * zpz(3, 7) == 6);
+    lean_assert(2 * zpz(4, 7) == 1);
+    lean_assert(1 / zpz(3, 7) == 5);
+    std::ostringstream out;
+    z = 3;
+    out << z;
+    lean_assert(out.str() == "3");
+}
+
+int main() {
+    tst1();
+    return has_violations() ? 1 : 0;
+}

src/util/int64.h

@@ -5,6 +5,8 @@ Released under Apache 2.0 license as described in the file LICENSE.
 Author: Leonardo de Moura
 */
 namespace lean {
+typedef long long int64;
 typedef unsigned long long uint64;
+static_assert(sizeof(int64 ) == 8,  "unexpected int64 size");
 static_assert(sizeof(uint64 ) == 8, "unexpected uint64 size");
 }

src/util/numerics/CMakeLists.txt

@@ -1,2 +1,4 @@
-add_library(numerics gmp_init.cpp mpz.cpp mpq.cpp mpbq.cpp mpfp.cpp float.cpp double.cpp numeric_traits.cpp primes.cpp)
+add_library(numerics gmp_init.cpp mpz.cpp mpq.cpp mpbq.cpp mpfp.cpp
+float.cpp double.cpp numeric_traits.cpp primes.cpp zpz.cpp)
+
 target_link_libraries(numerics ${LEAN_LIBS} ${EXTRA_LIBS})

src/util/numerics/power.h

@@ -0,0 +1,25 @@
+/*
+Copyright (c) 2013 Microsoft Corporation. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+
+Author: Leonardo de Moura
+*/
+
+namespace lean {
+/**
+   \brief Template for computing <tt>a^k</tt>.
+*/
+template<typename T>
+T power(T const & a, unsigned k) {
+    unsigned mask = 1;
+    T power = a;
+    T b     = 1;
+    while (mask <= k) {
+        if (mask & k)
+            b *= power;
+        power *= power;
+        mask = mask << 1;
+    }
+    return b;
+}
+}

src/util/numerics/primes.cpp

@@ -6,7 +6,7 @@ Author: Leonardo de Moura
 */
 #include <vector>
 #include <mutex>
-#include "util/uint64.h"
+#include "util/int64.h"
 #include "util/debug.h"
 #include "util/exception.h"
 #include "util/numerics/primes.h"

src/util/numerics/primes.h

@@ -4,7 +4,7 @@ Released under Apache 2.0 license as described in the file LICENSE.
 
 Author: Leonardo de Moura
 */
-#include "util/uint64.h"
+#include "util/int64.h"
 
 namespace lean {
 /** \brief Prime number iterator. It can be used to enumerate the first LEAN_PRIME_LIST_MAX_SIZE primes. */

src/util/numerics/remainder.h

@@ -0,0 +1,20 @@
+/*
+Copyright (c) 2013 Microsoft Corporation. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+
+Author: Leonardo de Moura
+*/
+#include "util/debug.h"
+
+namespace lean {
+template<class T>
+T remainder(T a, T b) {
+    lean_assert(b != 0);
+    if (a > 0)
+        return a % b;
+    else if (b > 0)
+        return a % b + b;
+    else
+        return a % b - b;
+}
+}

src/util/numerics/zpz.cpp

@@ -0,0 +1,21 @@
+/*
+Copyright (c) 2013 Microsoft Corporation. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+
+Author: Leonardo de Moura
+*/
+#include "util/debug.h"
+#include "util/numerics/gcd.h"
+#include "util/numerics/zpz.h"
+
+namespace lean {
+void zpz::inv() {
+    lean_assert(m_value != 0);
+    lean_assert(is_normalized());
+    // eulers theorem a^(p - 2), but gcd could be more efficient
+    // a*t1 + p*t2 = 1 => a*t1 = 1 (mod p) => t1 is the inverse (t3 == 1)
+    int64 g, a, b;
+    gcdext(g, a, b, static_cast<int64>(m_value), static_cast<int64>(m_p));
+    m_value = remainder(a, static_cast<int64>(m_p));
+}
+}

src/util/numerics/zpz.h

@@ -0,0 +1,101 @@
+/*
+Copyright (c) 2013 Microsoft Corporation. All rights reserved.
+Released under Apache 2.0 license as described in the file LICENSE.
+
+Author: Leonardo de Moura
+*/
+#include <algorithm>
+#include "util/debug.h"
+#include "util/int64.h"
+#include "util/numerics/remainder.h"
+#include "util/numerics/primes.h"
+
+namespace lean {
+/**
+   \brief The Z/pZ field (the set of integers modulo a prime p).
+   We use machine integers to represent the values. That is, we only
+   consider primes < 2^32 - 1.
+
+   The values are encoded as a pair (value, p). We want to be able to
+   dynamically change the prime p. This feature is needed when
+   implementing some algorithms based on modular arithmetic.
+*/
+class zpz {
+    unsigned m_value;
+    unsigned m_p;
+    bool is_normalized() const { return m_value < m_p; }
+    void normalize() { m_value %= m_p; }
+public:
+    zpz():m_value(0), m_p(2) {}
+    zpz(unsigned v, unsigned p):m_value(v), m_p(p) { lean_assert(is_prime(p)); }
+
+    unsigned p() { return m_p; }
+    unsigned hash() const { return m_value; }
+    unsigned get_unsigned_int() const { return m_value; }
+    void set_p(unsigned p) { lean_assert(is_prime(p)); m_p = p; normalize(); }
+
+    friend void swap(zpz & a, zpz & b) { std::swap(a.m_value, b.m_value); std::swap(a.m_p, b.m_p); }
+
+    friend bool operator==(zpz const & a, zpz const & b) { return a.m_value == b.m_value; }
+    friend bool operator!=(zpz const & a, zpz const & b) { return !(a == b); }
+    friend bool operator<(zpz const & a, zpz const & b)  { return a.m_value < b.m_value; }
+    friend bool operator>(zpz const & a, zpz const & b)  { return a.m_value > b.m_value; }
+    friend bool operator<=(zpz const & a, zpz const & b) { return a.m_value <= b.m_value; }
+    friend bool operator>=(zpz const & a, zpz const & b) { return a.m_value >= b.m_value; }
+
+    friend bool operator==(zpz const & a, unsigned b) { return a.m_value == b; }
+    friend bool operator!=(zpz const & a, unsigned b) { return !(a == b); }
+    friend bool operator<(zpz const & a, unsigned b)  { return a.m_value < b; }
+    friend bool operator>(zpz const & a, unsigned b)  { return a.m_value > b; }
+    friend bool operator<=(zpz const & a, unsigned b) { return a.m_value <= b; }
+    friend bool operator>=(zpz const & a, unsigned b) { return a.m_value >= b; }
+
+    friend bool operator==(unsigned a, zpz const & b) { return a == b.m_value; }
+    friend bool operator!=(unsigned a, zpz const & b) { return !(a == b); }
+    friend bool operator<(unsigned a, zpz const & b)  { return a < b.m_value; }
+    friend bool operator>(unsigned a, zpz const & b)  { return a > b.m_value; }
+    friend bool operator<=(unsigned a, zpz const & b) { return a <= b.m_value; }
+    friend bool operator>=(unsigned a, zpz const & b) { return a >= b.m_value; }
+
+    zpz & operator=(zpz const & v) { m_value = v.m_value; m_p = v.m_p; lean_assert(is_normalized()); return *this; }
+    zpz & operator=(unsigned v)    { m_value = v; normalize(); return *this; }
+
+    zpz & operator+=(unsigned v)    { m_value = (static_cast<uint64>(m_value) + static_cast<uint64>(v)) % m_p; return *this; }
+    zpz & operator+=(zpz const & v) { return operator+=(v.m_value); }
+
+    zpz & operator*=(unsigned v)    { m_value = (static_cast<uint64>(m_value) * static_cast<uint64>(v)) % m_p; return *this; }
+    zpz & operator*=(zpz const & v) { return operator*=(v.m_value); }
+
+    zpz & operator-=(unsigned v)    { m_value = remainder(static_cast<int64>(m_value) - static_cast<int64>(v), static_cast<int64>(m_p)); return *this; }
+    zpz & operator-=(zpz const & v) { return operator-=(v.m_value); }
+
+    zpz & operator++() { m_value++; if (m_value == m_p) m_value = 0; return *this; }
+    zpz & operator--() { if (m_value == 0) m_value = m_p - 1; else m_value--; return *this; }
+    zpz operator++(int) { zpz tmp(*this); operator++(); return tmp; }
+    zpz operator--(int) { zpz tmp(*this); operator--(); return tmp; }
+
+    void inv();
+    void neg() { m_value = remainder(-static_cast<int64>(m_value), static_cast<int64>(m_p)); }
+
+    zpz & operator/=(zpz v) { v.inv(); return operator*=(v); return *this; }
+    zpz & operator/=(unsigned v) { return operator/=(zpz(v, m_p)); }
+
+    friend zpz operator+(zpz a, zpz const & b) { return a += b; }
+    friend zpz operator+(zpz a, unsigned b)    { return a += b; }
+    friend zpz operator+(unsigned a, zpz b)    { return b += a; }
+
+    friend zpz operator-(zpz a, zpz const & b) { return a -= b; }
+    friend zpz operator-(zpz a, unsigned b)    { return a -= b; }
+    friend zpz operator-(unsigned a, zpz b)    { b.neg(); return b += a; }
+
+    friend zpz operator*(zpz a, zpz const & b) { return a *= b; }
+    friend zpz operator*(zpz a, unsigned b)    { return a *= b; }
+    friend zpz operator*(unsigned a, zpz b)    { return b *= a; }
+
+    friend zpz operator/(zpz a, zpz const & b) { return a /= b; }
+    friend zpz operator/(zpz a, unsigned b)    { return a /= b; }
+    friend zpz operator/(unsigned a, zpz b)    { b.inv(); return b *= a; }
+
+    friend std::ostream & operator<<(std::ostream & out, zpz const & z) { out << z.m_value; return out; }
+};
+}