/src/util/mpzzp.h

/*++
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