z3-src 500.0.0

/*++
Copyright (c) 2006 Microsoft Corporation

Module Name:

    rational.h

Abstract:

    Rational numbers

Author:

    Leonardo de Moura (leonardo) 2006-09-18.

Revision History:

--*/
#pragma once

#include "util/mpq.h"

class rational {
    mpq   m_val;
    static rational                  m_zero;
    static rational                  m_one;
    static rational                  m_minus_one;
    static vector<rational>          m_powers_of_two;
    static synch_mpq_manager *       g_mpq_manager;
    
    static synch_mpq_manager & m() { return *g_mpq_manager; }

    void display_hex(std::ostream & out, unsigned num_bits) const { SASSERT(is_int()); m().display_hex(out, m_val.numerator(), num_bits); }

    void display_bin(std::ostream& out, unsigned num_bits) const { SASSERT(is_int()); m().display_bin(out, m_val.numerator(), num_bits);  }

public:
    static void initialize();
    static void finalize();
    /*
      ADD_INITIALIZER('rational::initialize();')
      ADD_FINALIZER('rational::finalize();')
    */
    rational() = default;

    rational(rational const & r) { m().set(m_val, r.m_val); }
    rational(rational&&) = default;

    explicit rational(int n) { m().set(m_val, n); }

    explicit rational(unsigned n) { m().set(m_val, n); }
      
    rational(int n, int d) { m().set(m_val, n, d); }
    rational(mpq const & q) { m().set(m_val, q); }
    rational(mpq && q) noexcept : m_val(std::move(q)) {}
    rational(mpz const & z) { m().set(m_val, z); }
    rational(mpz && z) noexcept : m_val(std::move(z)) {}

    explicit rational(double  z) { UNREACHABLE(); }
    
    explicit rational(char const * v) { m().set(m_val, v); }
    
    explicit rational(unsigned const * v, unsigned sz) { m().set(m_val, sz, v); }

    struct i64 {};
    rational(int64_t i, i64) { m().set(m_val, i); }

    struct ui64 {};
    rational(uint64_t i, ui64) { m().set(m_val, i); }
    
    ~rational() { synch_mpq_manager::del(g_mpq_manager, m_val); }
    
    mpq const & to_mpq() const { return m_val; }

    unsigned bitsize() const { return m().bitsize(m_val); }

    unsigned storage_size() const { return m().storage_size(m_val); }
    
    void reset() { m().reset(m_val); }

    bool is_int() const { return m().is_int(m_val); }

    bool is_small() const { return m().is_small(m_val); }

    bool is_big() const { return !is_small(); }
    
    unsigned hash() const { return m().hash(m_val); }

    struct hash_proc {  unsigned operator()(rational const& r) const { return r.hash(); }  };

    struct eq_proc { bool operator()(rational const& r1, rational const& r2) const { return r1 == r2; } };
    
    void swap(rational & n) noexcept { m().swap(m_val, n.m_val); }
    
    std::string to_string() const { return m().to_string(m_val); }

    void display(std::ostream & out) const { return m().display(out, m_val); }
    
    void display_decimal(std::ostream & out, unsigned prec, bool truncate = false) const { return m().display_decimal(out, m_val, prec, truncate); }

    void display_smt2(std::ostream & out) const { return m().display_smt2(out, m_val, false); }


    struct as_hex_wrapper {
        rational const& r;
        unsigned bw;
    };

    as_hex_wrapper as_hex(unsigned bw) const { return as_hex_wrapper{*this, bw}; }

    friend inline std::ostream& operator<<(std::ostream& out, as_hex_wrapper const& ab) {
        ab.r.display_hex(out, ab.bw);
        return out;
    }



    struct as_bin_wrapper {
        rational const& r;
        unsigned bw;
    };

    as_bin_wrapper as_bin(unsigned bw) const { return as_bin_wrapper{*this, bw}; }

    friend inline std::ostream& operator<<(std::ostream& out, as_bin_wrapper const& ab) {
        ab.r.display_bin(out, ab.bw);
        return out;
    }


    bool is_uint64() const { return m().is_uint64(m_val); }

    bool is_int64() const { return m().is_int64(m_val); }

    uint64_t get_uint64() const { return m().get_uint64(m_val); }

    int64_t get_int64() const { return m().get_int64(m_val); }
    
    bool is_unsigned() const { return is_uint64() && (get_uint64() < (1ull << 32ull)); }

    unsigned get_unsigned() const {
        SASSERT(is_unsigned());
        return static_cast<unsigned>(get_uint64());
    }

    bool is_int32() const {
        if (is_small() && is_int()) return true; 
        // we don't assume that if it is small, then it is int32.
        if (!is_int64()) return false;
        int64_t v = get_int64();
        return INT_MIN <= v && v <= INT_MAX;
    }

    int get_int32() const {
        SASSERT(is_int32());
        return (int)get_int64();
    }

    double get_double() const { return m().get_double(m_val); }

    rational const & get_rational() const { return *this; }

    rational const & get_infinitesimal() const { return m_zero; }

    rational & operator=(rational&&) = default;

    rational & operator=(rational const & r) {
        m().set(m_val, r.m_val);
        return *this;
    }

    rational & operator=(bool) = delete;
    rational operator*(bool  r1) const = delete;

    rational & operator=(int v) {
        m().set(m_val, v);
        return *this;
    }
    rational & operator=(double v) = delete;

    friend inline rational numerator(rational const & r) { rational result; m().get_numerator(r.m_val, result.m_val); return result; }
    
    friend inline rational denominator(rational const & r) { rational result; m().get_denominator(r.m_val, result.m_val); return result; }

    friend inline rational inv(rational const & r) {
        rational result;
        m().inv(r.m_val, result.m_val);
        return result;
    }
    
    rational & operator+=(rational const & r) { 
        m().add(m_val, r.m_val, m_val);
        return *this; 
    }

    rational & operator+=(int r) {
        (*this) += rational(r);
        return *this;
    }

    rational & operator-=(rational const & r) { 
        m().sub(m_val, r.m_val, m_val);
        return *this; 
    }

    rational& operator-=(int r) {
        (*this) -= rational(r);
        return *this;
    }

    rational & operator*=(rational const & r) {
        m().mul(m_val, r.m_val, m_val);
        return *this; 
    }    

    rational & operator/=(rational const & r) {
        m().div(m_val, r.m_val, m_val);
        return *this; 
    }    
    
    rational & operator%=(rational const & r) {
        m().rem(m_val, r.m_val, m_val);
        return *this; 
    }    

    rational & operator%=(int v) {
        return *this %= rational(v);
    }    

    rational & operator/=(int v) {
        return *this /= rational(v);
    }    

    rational & operator*=(int v) {
        return *this *= rational(v);
    }    

    friend inline rational div(rational const & r1, rational const & r2) {
        rational r;
        rational::m().idiv(r1.m_val, r2.m_val, r.m_val);
        return r;
    }

    friend inline void div(rational const & r1, rational const & r2, rational & r) {
        rational::m().idiv(r1.m_val, r2.m_val, r.m_val);
    }
    
    friend inline rational machine_div(rational const & r1, rational const & r2) {
        rational r;
        rational::m().machine_idiv(r1.m_val, r2.m_val, r.m_val);
        return r;
    }

    friend inline rational machine_div_rem(rational const & r1, rational const & r2, rational & rem) {
        rational r;
        rational::m().machine_idiv_rem(r1.m_val, r2.m_val, r.m_val, rem.m_val);
        return r;
    }

    friend inline rational machine_div2k(rational const & r1, unsigned k) {
        rational r;
        rational::m().machine_idiv2k(r1.m_val, k, r.m_val);
        return r;
    }

    friend inline rational mod(rational const & r1, rational const & r2) {
        rational r;
        rational::m().mod(r1.m_val, r2.m_val, r.m_val);
        return r;
    }
    
    friend inline void mod(rational const & r1, rational const & r2, rational & r) {
        rational::m().mod(r1.m_val, r2.m_val, r.m_val);
    }

    friend inline rational mod2k(rational const & a, unsigned k) {
        if (a.is_int())
            return rational::m().mod2k(a.m_val.numerator(), k);
        return mod(a, power_of_two(k));
    }

    friend inline rational operator%(rational const & r1, rational const & r2) {
        rational r;
        rational::m().rem(r1.m_val, r2.m_val, r.m_val);
        return r;
    }

    friend inline rational mod_hat(rational const & a, rational const & b) {
        SASSERT(b.is_pos());
        rational r = mod(a,b);
        SASSERT(r.is_nonneg());
        rational r2 = r;
        r2 *= rational(2);
        if (operator<(b, r2)) {
            r -= b;
        }
        return r;
    }

    rational & operator++() {
        m().add(m_val, m().mk_q(1), m_val);
        return *this;
    }

    const rational operator++(int) { rational tmp(*this); ++(*this); return tmp; }
    
    rational & operator--() {
        m().sub(m_val, m().mk_q(1), m_val);
        return *this;
    }
    
    const rational operator--(int) { rational tmp(*this); --(*this); return tmp; }

    friend inline bool operator==(rational const & r1, rational const & r2) {
        return rational::m().eq(r1.m_val, r2.m_val);
    }

    friend inline bool operator<(rational const & r1, rational const & r2) { 
        return rational::m().lt(r1.m_val, r2.m_val);
    }
    
    void neg() {
        m().neg(m_val);
    }

    bool is_zero() const {
        return m().is_zero(m_val);
    }

    bool is_one() const {
        return m().is_one(m_val);
    }

    bool is_minus_one() const {
        return m().is_minus_one(m_val);
    }

    bool is_neg() const {
        return m().is_neg(m_val);
    }
    
    bool is_pos() const {
        return m().is_pos(m_val);
    }
    
    bool is_nonneg() const {
        return m().is_nonneg(m_val);
    }
    
    bool is_nonpos() const {
        return m().is_nonpos(m_val);
    }

    bool is_even() const {
        return m().is_even(m_val);
    }

    bool is_odd() const { 
        return !is_even(); 
    }
    
    friend inline rational floor(rational const & r) {
        rational f;
        rational::m().floor(r.m_val, f.m_val);
        return f;
    }

    friend inline rational ceil(rational const & r) {
        rational f;
        rational::m().ceil(r.m_val, f.m_val);
        return f;
    }

    rational expt(int n) const {
        rational result;
        m().power(m_val, n, result.m_val);
        return result;
    }

    static rational power_of_two(unsigned k);

    bool is_power_of_two(unsigned & shift) const {
        return m().is_power_of_two(m_val, shift);
    }
    
    bool is_power_of_two() const {
        unsigned shift = 0;
        return m().is_power_of_two(m_val, shift);
    }

    bool mult_inverse(unsigned num_bits, rational & result) const;
    rational pseudo_inverse(unsigned num_bits) const;

    static rational const & zero() {
        return m_zero;
    }

    static rational const & one() {
        return m_one;
    }

    static rational const & minus_one() {
        return m_minus_one;
    }

    void addmul(rational const & c, rational const & k) {
        if (c.is_one())
            operator+=(k);
        else if (c.is_minus_one())
            operator-=(k);
        else if (k.is_one())
            operator+=(c);
        else if (k.is_minus_one())
            operator-=(c);
        else {
            rational tmp(k);
            tmp *= c;
            operator+=(tmp);
        }
    }

    // Perform:  this -= c * k
    void submul(const rational & c, const rational & k) {
        if (c.is_one())
            operator-=(k);
        else if (c.is_minus_one())
            operator+=(k);
        else {
            rational tmp(k);
            tmp *= c;
            operator-=(tmp);
        }
    }

    bool is_int_perfect_square(rational & root) const {
        return m().is_int_perfect_square(m_val, root.m_val);
    }

    bool is_perfect_square(rational & root) const {
        return m().is_perfect_square(m_val, root.m_val);
    }

    bool root(unsigned n, rational & root) const {
        return m().root(m_val, n, root.m_val);
    }

    friend inline std::ostream & operator<<(std::ostream & target, rational const & r) {
        return target << m().to_string(r.m_val);
    }

    friend inline bool divides(rational const& a, rational const& b) {
        return m().divides(a.to_mpq(), b.to_mpq());
    }

    friend inline rational gcd(rational const & r1, rational const & r2);

    //
    // extended Euclid:
    // r1*a + r2*b = gcd
    //
    friend inline rational gcd(rational const & r1, rational const & r2, rational & a, rational & b);

    friend inline rational lcm(rational const & r1, rational const & r2) {
        rational result;
        m().lcm(r1.m_val, r2.m_val, result.m_val);
        return result;
    }

    friend inline rational bitwise_or(rational const & r1, rational const & r2) {
        rational result;
        m().bitwise_or(r1.m_val, r2.m_val, result.m_val);
        return result;
    }

    friend inline rational bitwise_and(rational const & r1, rational const & r2) {
        rational result;
        m().bitwise_and(r1.m_val, r2.m_val, result.m_val);
        return result;
    }

    friend inline rational bitwise_xor(rational const & r1, rational const & r2) {
        rational result;
        m().bitwise_xor(r1.m_val, r2.m_val, result.m_val);
        return result;
    }

    friend inline rational bitwise_not(unsigned sz, rational const & r1) {
        rational result;
        m().bitwise_not(sz, r1.m_val, result.m_val);
        return result;
    }

    friend inline rational abs(rational const & r);

    rational to_rational() const { return *this; }

    static bool is_rational() { return true; }

    unsigned get_num_digits(rational const& base) const {
        SASSERT(is_int());
        SASSERT(!is_neg());
        rational n(*this);
        unsigned num_digits = 1;
        n = div(n, base);
        while (n.is_pos()) {
            ++num_digits;
            n = div(n, base);
        }
        return num_digits;
    }

    unsigned get_num_bits() const {
        return get_num_digits(rational(2));
    }

    unsigned get_num_decimal() const {
        return get_num_digits(rational(10));
    }

    /**
     * \brief Return the biggest k s.t. 2^k <= a.
     * \remark Return 0 if a is not positive.
     */
    unsigned prev_power_of_two() const { return m().prev_power_of_two(m_val); }

    /**
     * \brief Return the smallest k s.t. a <= 2^k.
     * \remark Return 0 if a is not positive.
     */
    unsigned next_power_of_two() const { return m().next_power_of_two(m_val); }

    bool get_bit(unsigned index) const {
        return m().get_bit(m_val, index);
    }

    unsigned trailing_zeros() const {
        if (is_zero())
            return 0;
        unsigned k = 0;
        for (; !get_bit(k); ++k); 
        return k;
    }

    /** Number of trailing zeros in an N-bit representation */
    unsigned parity(unsigned num_bits) const {
        SASSERT(!is_neg());
        SASSERT(*this < rational::power_of_two(num_bits));
        if (is_zero())
            return num_bits;
        return trailing_zeros();
    }

    static bool limit_denominator(rational &num, rational const& limit);
};

inline bool operator!=(rational const & r1, rational const & r2) { 
    return !operator==(r1, r2); 
}

inline bool operator>(rational const & r1, rational const & r2) { 
    return operator<(r2, r1); 
}

inline bool operator<(int r1, rational const & r2) {
    return rational(r1) < r2;
}

inline bool operator<(rational const & r1, int r2) {
    return r1 < rational(r2);
}

inline bool operator<=(rational const & r1, rational const & r2) { 
    return !operator>(r1, r2); 
}

inline bool operator>=(rational const & r1, rational const & r2) { 
    return !operator<(r1, r2); 
}

inline bool operator>(rational const & a, int b) {
    return a > rational(b);
}

inline bool operator>(int a, rational const & b) {
    return rational(a) > b;
}

inline bool operator>=(rational const& a, int b) {
    return a >= rational(b);
}

inline bool operator>=(int a, rational const& b) {
    return rational(a) >= b;
}

inline bool operator<=(rational const& a, int b) {
    return a <= rational(b);
}

inline bool operator<=(int a, rational const& b) {
    return rational(a) <= b;
}

inline bool operator!=(rational const& a, int b) {
    return !(a == rational(b));
}

inline bool operator==(rational const & a, int b) {
    return a == rational(b);
}

inline rational operator+(rational const & r1, rational const & r2) { 
    return rational(r1) += r2; 
}

inline rational operator+(int r1, rational const & r2) {
    return rational(r1) + r2;
}

inline rational operator+(rational const & r1, int r2) {
    return r1 + rational(r2);
}


inline rational operator-(rational const & r1, rational const & r2) { 
    return rational(r1) -= r2; 
}

inline rational operator-(rational const & r1, int r2) {
    return r1 - rational(r2);
}

inline rational operator-(int r1, rational const & r2) {
    return rational(r1) - r2;
}

inline rational operator-(rational const & r) { 
    rational result(r); 
    result.neg(); 
    return result; 
}

inline rational operator*(rational const & r1, rational const & r2) { 
    return rational(r1) *= r2; 
}

inline rational operator*(rational const & r1, bool r2) {
    UNREACHABLE();
    return r1 * rational(r2);
}
inline rational operator*(rational const & r1, int r2) {
    return r1 * rational(r2);
}
inline rational operator*(bool  r1, rational const & r2) {
    UNREACHABLE();
    return rational(r1) * r2;
}

inline rational operator*(int  r1, rational const & r2) {
    return rational(r1) * r2;
}

inline rational operator/(rational const & r1, rational const & r2) { 
    return rational(r1) /= r2; 
}

inline rational operator/(rational const & r1, int r2) {
    return r1 / rational(r2);
}

inline rational operator/(rational const & r1, bool r2) {
    UNREACHABLE();
    return r1 / rational(r2);
}

inline rational operator/(int r1, rational const &    r2) {
    return rational(r1) / r2;
}

inline rational power(rational const & r, unsigned p) {
    return r.expt(p);
}

inline rational abs(rational const & r) {
  rational result(r);
  rational::m().abs(result.m_val);
  return result;
}

inline rational gcd(rational const & r1, rational const & r2) {
  rational result;
  rational::m().gcd(r1.m_val, r2.m_val, result.m_val);
  return result;
}

inline rational gcd(rational const & r1, rational const & r2, rational & a, rational & b) {
  rational result;
  rational::m().gcd(r1.m_val, r2.m_val, a.m_val, b.m_val, result.m_val);
  return result;
}

inline void swap(rational& r1, rational& r2) noexcept {
    r1.swap(r2);
}