/src/util/mpzzp.h
C Header | 285 lines | 221 code | 33 blank | 31 comment | 32 complexity | 31186fbcf6b65d9eeb28a6aabf9354a7 MD5 | raw file
- /*++
- Copyright (c) 2012 Microsoft Corporation
- Module Name:
- mpzzp.h
- Abstract:
- Combines Z ring, GF(p) finite field, and Z_p ring (when p is not a prime)
- in a single manager;
- That is, the manager may be dynamically configured
- to be Z Ring, GF(p), etc.
- Author:
- Leonardo 2012-01-17.
- Revision History:
- This code is based on mpzp.h.
- In the future, it will replace it.
- --*/
- #ifndef _MPZZP_H_
- #define _MPZZP_H_
- #include "mpz.h"
- class mpzzp_manager {
- typedef unsynch_mpz_manager numeral_manager;
- numeral_manager & m_manager;
-
- bool m_z;
- // instead the usual [0..p) we will keep the numbers in [lower, upper]
- mpz m_p, m_lower, m_upper;
- bool m_p_prime;
- mpz m_inv_tmp1, m_inv_tmp2, m_inv_tmp3;
- mpz m_div_tmp;
- bool is_p_normalized_core(mpz const & x) const {
- return m().ge(x, m_lower) && m().le(x, m_upper);
- }
-
- void setup_p() {
- SASSERT(m().is_pos(m_p) && !m().is_one(m_p));
- bool even = m().is_even(m_p);
- m().div(m_p, 2, m_upper);
- m().set(m_lower, m_upper);
- m().neg(m_lower);
- if (even) {
- m().inc(m_lower);
- }
- TRACE("mpzzp", tout << "lower: " << m_manager.to_string(m_lower) << ", upper: " << m_manager.to_string(m_upper) << "\n";);
- }
- void p_normalize_core(mpz & x) {
- SASSERT(!m_z);
- m().rem(x, m_p, x);
- if (m().gt(x, m_upper)) {
- m().sub(x, m_p, x);
- } else {
- if (m().lt(x, m_lower)) {
- m().add(x, m_p, x);
- }
- }
- SASSERT(is_p_normalized(x));
- }
- public:
- typedef mpz numeral;
- static bool precise() { return true; }
- bool field() { return !m_z && m_p_prime; }
- bool finite() const { return !m_z; }
- bool modular() const { return !m_z; }
- mpzzp_manager(numeral_manager & _m):
- m_manager(_m),
- m_z(true) {
- }
-
- mpzzp_manager(numeral_manager & _m, mpz const & p, bool prime = true):
- m_manager(_m),
- m_z(false) {
- m().set(m_p, p);
- setup_p();
- }
- mpzzp_manager(numeral_manager & _m, uint64 p, bool prime = true):
- m_manager(_m),
- m_z(false) {
- m().set(m_p, p);
- setup_p();
- }
- ~mpzzp_manager() {
- m().del(m_p);
- m().del(m_lower);
- m().del(m_upper);
- m().del(m_inv_tmp1);
- m().del(m_inv_tmp2);
- m().del(m_inv_tmp3);
- m().del(m_div_tmp);
- }
- bool is_p_normalized(mpz const & x) const {
- return m_z || is_p_normalized_core(x);
- }
- void p_normalize(mpz & x) {
- if (!m_z)
- p_normalize_core(x);
- SASSERT(is_p_normalized(x));
- }
-
- numeral_manager & m() const { return m_manager; }
-
- mpz const & p() const { return m_p; }
- void set_z() { m_z = true; }
- void set_zp(mpz const & new_p) { m_z = false; m_p_prime = true; m().set(m_p, new_p); setup_p(); }
- void set_zp(uint64 new_p) { m_z = false; m_p_prime = true; m().set(m_p, new_p); setup_p(); }
- // p = p^2
- void set_p_sq() { SASSERT(!m_z); m_p_prime = false; m().mul(m_p, m_p, m_p); setup_p(); }
- void set_zp_swap(mpz & new_p) { SASSERT(!m_z); m().swap(m_p, new_p); setup_p(); }
- void reset(mpz & a) { m().reset(a); }
- bool is_small(mpz const & a) { return m().is_small(a); }
- void del(mpz & a) { m().del(a); }
- void neg(mpz & a) { m().neg(a); p_normalize(a); }
- void abs(mpz & a) { m().abs(a); p_normalize(a); }
- bool is_zero(mpz const & a) { SASSERT(is_p_normalized(a)); return numeral_manager::is_zero(a); }
- bool is_one(mpz const & a) { SASSERT(is_p_normalized(a)); return numeral_manager::is_one(a); }
- bool is_pos(mpz const & a) { SASSERT(is_p_normalized(a)); return numeral_manager::is_pos(a); }
- bool is_neg(mpz const & a) { SASSERT(is_p_normalized(a)); return numeral_manager::is_neg(a); }
- bool is_nonpos(mpz const & a) { SASSERT(is_p_normalized(a)); return numeral_manager::is_nonpos(a); }
- bool is_nonneg(mpz const & a) { SASSERT(is_p_normalized(a)); return numeral_manager::is_nonneg(a); }
- bool is_minus_one(mpz const & a) { SASSERT(is_p_normalized(a)); return numeral_manager::is_minus_one(a); }
- bool eq(mpz const & a, mpz const & b) { SASSERT(is_p_normalized(a) && is_p_normalized(b)); return m().eq(a, b); }
- bool lt(mpz const & a, mpz const & b) { SASSERT(is_p_normalized(a) && is_p_normalized(b)); return m().lt(a, b); }
- bool le(mpz const & a, mpz const & b) { SASSERT(is_p_normalized(a) && is_p_normalized(b)); return m().le(a, b); }
- bool gt(mpz const & a, mpz const & b) { SASSERT(is_p_normalized(a) && is_p_normalized(b)); return m().gt(a, b); }
- bool ge(mpz const & a, mpz const & b) { SASSERT(is_p_normalized(a) && is_p_normalized(b)); return m().ge(a, b); }
- std::string to_string(mpz const & a) const {
- SASSERT(is_p_normalized(a));
- return m().to_string(a);
- }
- void display(std::ostream & out, mpz const & a) const { m().display(out, a); }
- void add(mpz const & a, mpz const & b, mpz & c) { SASSERT(is_p_normalized(a) && is_p_normalized(b)); m().add(a, b, c); p_normalize(c); }
- void sub(mpz const & a, mpz const & b, mpz & c) { SASSERT(is_p_normalized(a) && is_p_normalized(b)); m().sub(a, b, c); p_normalize(c); }
- void inc(mpz & a) { SASSERT(is_p_normalized(a)); m().inc(a); p_normalize(a); }
- void dec(mpz & a) { SASSERT(is_p_normalized(a)); m().dec(a); p_normalize(a); }
- void mul(mpz const & a, mpz const & b, mpz & c) { SASSERT(is_p_normalized(a) && is_p_normalized(b)); m().mul(a, b, c); p_normalize(c); }
- void addmul(mpz const & a, mpz const & b, mpz const & c, mpz & d) {
- SASSERT(is_p_normalized(a) && is_p_normalized(b) && is_p_normalized(c)); m().addmul(a, b, c, d); p_normalize(d);
- }
- // d <- a - b*c
- void submul(mpz const & a, mpz const & b, mpz const & c, mpz & d) {
- SASSERT(is_p_normalized(a));
- SASSERT(is_p_normalized(b));
- SASSERT(is_p_normalized(c));
- m().submul(a, b, c, d);
- p_normalize(d);
- }
- void inv(mpz & a) {
- if (m_z) {
- UNREACHABLE();
- }
- else {
- SASSERT(!is_zero(a));
- // 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)
- TRACE("mpzp_inv_bug", tout << "a: " << m().to_string(a) << ", p: " << m().to_string(m_p) << "\n";);
- p_normalize(a);
- TRACE("mpzp_inv_bug", tout << "after normalization a: " << m().to_string(a) << "\n";);
- m().gcd(a, m_p, m_inv_tmp1, m_inv_tmp2, m_inv_tmp3);
- TRACE("mpzp_inv_bug", tout << "tmp1: " << m().to_string(m_inv_tmp1) << "\ntmp2: " << m().to_string(m_inv_tmp2)
- << "\ntmp3: " << m().to_string(m_inv_tmp3) << "\n";);
- p_normalize(m_inv_tmp1);
- m().swap(a, m_inv_tmp1);
- SASSERT(m().is_one(m_inv_tmp3)); // otherwise p is not prime and inverse is not defined
- }
- }
-
- void swap(mpz & a, mpz & b) {
- SASSERT(is_p_normalized(a) && is_p_normalized(b));
- m().swap(a, b);
- }
- bool divides(mpz const & a, mpz const & b) { return (field() && !is_zero(a)) || m().divides(a, b); }
- // a/b = a*inv(b)
- void div(mpz const & a, mpz const & b, mpz & c) {
- if (m_z) {
- return m().div(a, b, c);
- }
- else {
- SASSERT(m_p_prime);
- SASSERT(is_p_normalized(a));
- m().set(m_div_tmp, b);
- inv(m_div_tmp);
- mul(a, m_div_tmp, c);
- SASSERT(is_p_normalized(c));
- }
- }
-
- static unsigned hash(mpz const & a) { return numeral_manager::hash(a); }
- void gcd(mpz const & a, mpz const & b, mpz & c) {
- SASSERT(is_p_normalized(a) && is_p_normalized(b));
- m().gcd(a, b, c);
- SASSERT(is_p_normalized(c));
- }
- void gcd(unsigned sz, mpz const * as, mpz & g) {
- m().gcd(sz, as, g);
- SASSERT(is_p_normalized(g));
- }
-
- void gcd(mpz const & r1, mpz const & r2, mpz & a, mpz & b, mpz & g) {
- SASSERT(is_p_normalized(r1) && is_p_normalized(r2));
- m().gcd(r1, r2, a, b, g);
- p_normalize(a);
- p_normalize(b);
- }
- void set(mpz & a, mpz & val) { m().set(a, val); p_normalize(a); }
- void set(mpz & a, int val) { m().set(a, val); p_normalize(a); }
- void set(mpz & a, unsigned val) { m().set(a, val); p_normalize(a); }
- void set(mpz & a, char const * val) { m().set(a, val); p_normalize(a); }
- void set(mpz & a, int64 val) { m().set(a, val); p_normalize(a); }
- void set(mpz & a, uint64 val) { m().set(a, val); p_normalize(a); }
- void set(mpz & a, mpz const & val) { m().set(a, val); p_normalize(a); }
- bool is_uint64(mpz & a) const { const_cast<mpzzp_manager*>(this)->p_normalize(a); return m().is_uint64(a); }
- bool is_int64(mpz & a) const { const_cast<mpzzp_manager*>(this)->p_normalize(a); return m().is_int64(a); }
- uint64 get_uint64(mpz & a) const { const_cast<mpzzp_manager*>(this)->p_normalize(a); return m().get_uint64(a); }
- int64 get_int64(mpz & a) const { const_cast<mpzzp_manager*>(this)->p_normalize(a); return m().get_int64(a); }
- double get_double(mpz & a) const { const_cast<mpzzp_manager*>(this)->p_normalize(a); return m().get_double(a); }
- void power(mpz const & a, unsigned k, mpz & b) {
- SASSERT(is_p_normalized(a));
- unsigned mask = 1;
- mpz power;
- set(power, a);
- set(b, 1);
- while (mask <= k) {
- if (mask & k)
- mul(b, power, b);
- mul(power, power, power);
- mask = mask << 1;
- }
- del(power);
- }
- bool is_perfect_square(mpz const & a, mpz & root) {
- if (m_z) {
- return m().is_perfect_square(a, root);
- }
- else {
- NOT_IMPLEMENTED_YET();
- return false;
- }
- }
- bool is_uint64(mpz const & a) const { return m().is_uint64(a); }
- bool is_int64(mpz const & a) const { return m().is_int64(a); }
- uint64 get_uint64(mpz const & a) const { return m().get_uint64(a); }
- int64 get_int64(mpz const & a) const { return m().get_int64(a); }
- void mul2k(mpz & a, unsigned k) { m().mul2k(a, k); p_normalize(a); }
- void mul2k(mpz const & a, unsigned k, mpz & r) { m().mul2k(a, k, r); p_normalize(r); }
- unsigned power_of_two_multiple(mpz const & n) { return m().power_of_two_multiple(n); }
- unsigned log2(mpz const & n) { return m().log2(n); }
- unsigned mlog2(mpz const & n) { return m().mlog2(n); }
- void machine_div2k(mpz & a, unsigned k) { m().machine_div2k(a, k); SASSERT(is_p_normalized(a)); }
- void machine_div2k(mpz const & a, unsigned k, mpz & r) { m().machine_div2k(a, k, r); SASSERT(is_p_normalized(r)); }
- bool root(mpz & a, unsigned n) { SASSERT(!modular()); return m().root(a, n); }
- bool root(mpz const & a, unsigned n, mpz & r) { SASSERT(!modular()); return m().root(a, n, r); }
- };
- #endif