emons branch

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

src/test/emonomials.cpp

@@ -0,0 +1,36 @@
+#include "util/lp/emonomials.h"
+
+static void display_mons(nla::emonomials const& em) {
+    std::cout << em << "\n";
+    for (auto const& m1 : em) {
+        std::cout << "factors: " << m1 << ": \n";
+        for (auto const& m2 : em.get_factors_of(m1)) {
+            std::cout << m2 << "\n";
+        }
+    }
+}
+
+void tst_emonomials() {
+    nla::var_eqs ve;
+    nla::emonomials em(ve);
+
+    unsigned x = 1;
+    unsigned y = 2;
+    unsigned z = 3;
+    unsigned u = 4;
+    unsigned v1 = 5;
+    unsigned v2 = 6;
+    unsigned v3 = 7;
+    em.add(v1, x, x);
+    em.add(v2, x, z);
+    em.add(v3, x, y, x);
+    display_mons(em);
+    em.push();
+    ve.merge_plus(x, y, nla::eq_justification({}));
+    display_mons(em);
+    em.push();
+    ve.merge_plus(x, z, nla::eq_justification({}));
+    display_mons(em);
+    em.pop(1);
+    display_mons(em);
+}

src/util/lp/CMakeLists.txt

@@ -6,6 +6,7 @@ z3_add_component(lp
     core_solver_pretty_printer.cpp
     dense_matrix.cpp
     eta_matrix.cpp
+    emonomials.cpp
     factorization.cpp
     factorization_factory_imp.cpp
     gomory.cpp

src/util/lp/emonomials.cpp

@@ -0,0 +1,286 @@
+/*++
+  Copyright (c) 2019 Microsoft Corporation
+
+  Module Name:
+
+     emonomials.cpp
+
+  Abstract:
+
+     table that associate monomials to congruence class representatives modulo a union find structure.
+
+  Author:
+    Nikolaj Bjorner (nbjorner)
+    Lev Nachmanson (levnach)
+
+  Revision History:
+
+    to replace rooted_mons.h and rooted_mon, rooted_mon_tabled
+
+  --*/
+
+#include "util/lp/emonomials.h"
+#include "util/lp/nla_defs.h"
+
+namespace nla {
+
+    void emonomials::inc_canonical() {
+        ++m_canonical;
+        if (m_canonical == 0) {
+            for (auto& svt : m_canonized) {
+                svt.m_canonical = 0;
+            }
+            ++m_canonical;
+        }
+    }
+
+    void emonomials::inc_visited() const {
+        ++m_visited;
+        if (m_visited == 0) {
+            for (auto& svt : m_canonized) {
+                svt.m_visited = 0;
+            }
+            ++m_visited;
+        }
+    }
+
+    void emonomials::push() { 
+        m_lim.push_back(m_monomials.size());
+        m_region.push_scope();
+        m_ve.push();
+    }
+
+    void emonomials::pop(unsigned n) {
+        m_ve.pop(n);
+        unsigned old_sz = m_lim[m_lim.size() - n];
+        for (unsigned i = m_monomials.size(); i-- > old_sz; ) {
+            monomial const& m = m_monomials[i];
+            m_var2index[m.var()] = UINT_MAX;
+            lpvar last_var = UINT_MAX;
+            for (lpvar v : m.vars()) {
+                if (v != last_var) {
+                    remove_var2monomials(v, i);
+                    last_var = v;
+                }
+            }
+        }
+        m_monomials.shrink(old_sz);
+        m_canonized.shrink(old_sz);
+        m_region.pop_scope(n);
+        m_lim.shrink(m_lim.size() - n);
+    }
+
+    void emonomials::remove_var2monomials(lpvar v, unsigned mIndex) {
+        cell*& cur_head = m_use_lists[v].m_head;
+        cell*& cur_tail = m_use_lists[v].m_tail;
+        cell* old_head = cur_head->m_next;
+        if (old_head == cur_head) {
+            cur_head = nullptr;
+            cur_tail = nullptr;
+        }
+        else {
+            cur_head = old_head;
+            cur_tail->m_next = old_head;
+        }
+    }
+
+    void emonomials::insert_var2monomials(lpvar v, unsigned mIndex) {
+        m_use_lists.reserve(v + 1);
+        cell*& cur_head = m_use_lists[v].m_head;
+        cell*& cur_tail = m_use_lists[v].m_tail;
+        cell* new_head = new (m_region) cell(mIndex, cur_head);
+        cur_head = new_head;
+        if (!cur_tail) cur_tail = new_head;
+        cur_tail->m_next = new_head;
+    }
+
+    void emonomials::merge_var2monomials(lpvar root, lpvar other) {
+        if (root == other) return;
+        m_use_lists.reserve(std::max(root, other) + 1);
+        cell*& root_head = m_use_lists[root].m_head;
+        cell*& root_tail = m_use_lists[root].m_tail;
+        cell* other_head = m_use_lists[other].m_head;
+        cell* other_tail = m_use_lists[other].m_tail;
+        if (root_head == nullptr) {
+            root_head = other_head;
+            root_tail = other_tail;
+        }
+        else if (other_head) {
+            // other_head -> other_tail -> root_head --> root_tail -> other_head.
+            root_tail->m_next = other_head;
+            other_tail->m_next = root_head;
+            root_head = other_head;
+        }
+        else {
+            // other_head = other_tail = nullptr
+        }
+    }
+
+    void emonomials::unmerge_var2monomials(lpvar root, lpvar other) {
+        if (root == other) return;
+        cell*& root_head = m_use_lists[root].m_head;
+        cell*& root_tail = m_use_lists[root].m_tail;
+        cell* other_head = m_use_lists[other].m_head;
+        cell* other_tail = m_use_lists[other].m_tail;
+        if (other_head == nullptr) {
+            // no-op
+        }
+        else if (root_tail == other_tail) {
+            root_head = nullptr;
+            root_tail = nullptr;
+        }
+        else {
+            root_head = other_tail->m_next;
+            root_tail->m_next = root_head;
+            other_tail->m_next = other_head;
+        }
+    }
+
+    emonomials::cell* emonomials::head(lpvar v) const {
+        v = m_ve.find(v).var();
+        m_use_lists.reserve(v + 1);
+        return m_use_lists[v].m_head;
+    }
+
+    void emonomials::set_visited(monomial const& m) const {
+        m_canonized[m_var2index[m.var()]].m_visited = m_visited;
+    }
+
+    bool emonomials::is_visited(monomial const& m) const {
+        return m_visited == m_canonized[m_var2index[m.var()]].m_visited;
+    }
+
+
+    /**
+       \brief insert a new monomial.
+
+       Assume that the variables are canonical, that is, not equal in current
+       context so another variable. To support equal in current context we
+       could track sign information in monomials, or ensure that insert_var2monomials 
+       works on the equivalence class roots.
+     */
+    void emonomials::add(lpvar v, unsigned sz, lpvar const* vs) {
+        unsigned idx = m_monomials.size();
+        m_monomials.push_back(monomial(v, sz, vs));
+        m_canonized.push_back(signed_vars_ts());
+        lpvar last_var = UINT_MAX;
+        for (unsigned i = 0; i < sz; ++i) {
+            lpvar w = vs[i];
+            SASSERT(m_ve.is_root(w));
+            if (w != last_var) {
+                insert_var2monomials(w, idx);
+                last_var = w; 
+            }
+        }
+        SASSERT(m_ve.is_root(v));
+        m_var2index.setx(v, idx, UINT_MAX);
+    }
+
+    signed_vars const& emonomials::canonize(monomial const& mon) const {
+        unsigned index = m_var2index[mon.var()];
+        if (m_canonized[index].m_canonical != m_canonical) {        
+            m_vars.reset();
+            bool sign = false;
+            for (lpvar v : mon) {
+                signed_var sv = m_ve.find(v);
+                m_vars.push_back(sv.var());
+                sign ^= sv.sign();
+            }
+            m_canonized[index].set_vars(m_vars.size(), m_vars.c_ptr());
+            m_canonized[index].set_sign(sign);
+            m_canonized[index].set_var(mon.var());
+            m_canonized[index].m_canonical = m_canonical;
+        }
+        return m_canonized[index];
+    }
+
+    bool emonomials::canonize_divides(monomial const& m1, monomial const& m2) const {
+        if (m1.size() > m2.size()) return false;
+        signed_vars const& s1 = canonize(m1);
+        signed_vars const& s2 = canonize(m2);
+        unsigned sz1 = s1.size(), sz2 = s2.size();
+        unsigned i = 0, j = 0;
+        while (true) {
+            if (i == sz1) {
+                return true;
+            }
+            else if (j == sz2) {
+                return false;
+            }
+            else if (s1[i] == s2[j]) {
+                ++i; ++j;
+            }
+            else if (s1[i] < s2[j]) {
+                return false;
+            }
+            else {
+                ++j;
+            }
+        }
+    }
+
+    void emonomials::explain_canonized(monomial const& m, lp::explanation& exp) {
+        for (lpvar v : m) {
+            signed_var w = m_ve.find(v);
+            m_ve.explain(signed_var(v, false), w, exp);
+        }
+    }
+
+    // yes, assume that monomials are non-empty.
+    emonomials::pf_iterator::pf_iterator(emonomials const& m, monomial const& mon, bool at_end):
+        m(m), m_mon(mon), m_it(iterator(m, m.head(mon[0]), at_end)), m_end(iterator(m, m.head(mon[0]), true)) {
+        fast_forward();
+    }
+
+    void emonomials::pf_iterator::fast_forward() {
+        for (; m_it != m_end; ++m_it) {
+            if (m_mon.var() != (*m_it).var() && m.canonize_divides(m_mon, *m_it) && !m.is_visited(*m_it)) {
+                m.set_visited(*m_it);
+                break;
+            }
+        }
+    }
+
+    void emonomials::merge_eh(signed_var r2, signed_var r1, signed_var v2, signed_var v1) {
+        if (!r2.sign() && m_ve.find(~r2) != m_ve.find(r1)) {
+            merge_var2monomials(r2.var(), r1.var());
+        }        
+        inc_canonical();
+    }
+
+    void emonomials::after_merge_eh(signed_var r2, signed_var r1, signed_var v2, signed_var v1) {
+        // skip
+    }
+
+    void emonomials::unmerge_eh(signed_var r2, signed_var r1) {
+        if (!r2.sign() && m_ve.find(~r2) != m_ve.find(r1)) {
+            unmerge_var2monomials(r2.var(), r1.var());            
+        }        
+        inc_canonical();
+    }
+
+   std::ostream& emonomials::display(std::ostream& out) const {
+       out << "monomials\n";
+       unsigned idx = 0;
+       for (auto const& m : m_monomials) {
+           out << (idx++) << ": " << m << "\n";
+       }
+       out << "use lists\n";
+       idx = 0;
+       for (auto const& ht : m_use_lists) {
+           cell* c = ht.m_head;
+           if (c) {
+               out << "v" << idx << ": ";
+               do {
+                   out << c->m_index << " ";
+                   c = c->m_next;
+               }
+               while (c != ht.m_head);
+               out << "\n";
+           }
+           ++idx;
+       }
+       return out;
+   }
+
+}

src/util/lp/emonomials.h

@@ -0,0 +1,246 @@
+/*++
+  Copyright (c) 2019 Microsoft Corporation
+
+  Module Name:
+
+     emonomials.h
+
+  Abstract:
+
+     table that associate monomials to congruence class representatives modulo a union find structure.
+
+  Author:
+    Nikolaj Bjorner (nbjorner)
+    Lev Nachmanson (levnach)
+
+  Revision History:
+
+    to replace rooted_mons.h and rooted_mon, rooted_mon_tabled
+
+  --*/
+
+#pragma once
+#include "util/lp/lp_utils.h"
+#include "util/lp/var_eqs.h"
+#include "util/lp/monomial.h"
+#include "util/region.h"
+
+namespace nla {
+
+    /**
+       \brief class used to summarize the coefficients to a monomial after
+       canonization with respect to current equalities.
+    */
+    class signed_vars {
+        lpvar           m_var; // variable representing original monomial
+        svector<lpvar>  m_vars;
+        bool            m_sign;
+    public:
+        signed_vars() : m_sign(false) {}
+        lpvar var() const { return m_var; }
+        svector<lpvar> const& vars() const { return m_vars; }
+        unsigned size() const { return m_vars.size(); }
+        lpvar operator[](unsigned i) const { return m_vars[i]; }
+        bool sign() const { return m_sign; }
+        void set_sign(bool s) { m_sign = s; }
+        void set_var(lpvar v) { m_var = v; }
+        void set_vars(unsigned n, lpvar const* vars) {
+            m_vars.reset();
+            m_vars.append(n, vars);
+            std::sort(m_vars.begin(), m_vars.end());
+        }
+        std::ostream& display(std::ostream& out) const {
+            out << "v" << var() << " := ";
+            if (sign()) out << "- ";
+            for (lpvar v : vars()) out << "v" << v << " ";
+            return out;
+        }
+    };
+
+    inline std::ostream& operator<<(std::ostream& out, signed_vars const& m) { return m.display(out); }
+
+
+    class emonomials : public var_eqs_merge_handler {
+        /**
+           \brief private fields used by emonomials for maintaining state of canonized monomials.
+        */
+        class signed_vars_ts : public signed_vars {
+        public:
+            signed_vars_ts(): m_canonical(0), m_visited(0) {}
+            unsigned m_canonical;
+            unsigned m_visited;
+        };
+
+        /**
+           \brief singly-lined cyclic list of monomial indices where variable occurs.
+           Each variable points to the head and tail of the cyclic list.
+           Initially, head and tail are nullptr.
+           New elements are inserted in the beginning of the list.
+           Two lists are merged when equivalence class representatives are merged, 
+           and the merge is undone when the representative variables are unmerged.
+         */
+        struct cell {
+            cell(unsigned mIndex, cell* c): m_next(c), m_index(mIndex) {}
+            cell*        m_next;
+            unsigned     m_index;
+        };
+        struct head_tail {
+            head_tail(): m_head(nullptr), m_tail(nullptr) {}
+            cell* m_head;
+            cell* m_tail;
+        };
+
+        var_eqs&                m_ve;
+        vector<monomial>        m_monomials;     // set of monomials
+        unsigned_vector         m_lim;           // backtracking point
+        unsigned_vector         m_var2index;     // var_mIndex -> mIndex
+        unsigned                m_canonical;     // timestamp of last merge
+        mutable unsigned        m_visited;       // timestamp of visited monomials during pf_iterator
+        region                  m_region;        // region for allocating linked lists
+        mutable vector<signed_vars_ts>  m_canonized;     // canonized versions of signed variables
+        mutable svector<lpvar>       m_vars;          // temporary vector of variables
+        mutable svector<head_tail>   m_use_lists;     // use list of monomials where variables occur.
+
+        void inc_canonical();
+        void inc_visited() const;
+
+        void remove_var2monomials(lpvar v, unsigned mIndex);
+        void insert_var2monomials(lpvar v, unsigned mIndex);
+        void merge_var2monomials(lpvar root, lpvar other);
+        void unmerge_var2monomials(lpvar root, lpvar other);
+
+        cell* head(lpvar v) const;
+        void set_visited(monomial const& m) const;
+        bool is_visited(monomial const& m) const;
+        
+    public:
+        /**
+           \brief emonomials builds on top of var_eqs.
+           push and pop on emonomials calls push/pop on var_eqs, so no 
+           other calls to push/pop to the var_eqs should take place. 
+         */
+        emonomials(var_eqs& ve): m_ve(ve), m_canonical(1), m_visited(0) { m_ve.set_merge_handler(this); }
+
+        /**
+           \brief push/pop scopes. 
+           The life-time of a merge is local within a scope.
+        */
+        void push();
+
+        void pop(unsigned n);
+
+        /**
+           \brief create a monomial from an equality v := vs
+        */
+        void add(lpvar v, unsigned sz, lpvar const* vs);
+        void add(lpvar v, svector<lpvar> const& vs) { add(v, vs.size(), vs.c_ptr()); }
+        void add(lpvar v, lpvar x, lpvar y) { lpvar vs[2] = { x, y }; add(v, 2, vs); }
+        void add(lpvar v, lpvar x, lpvar y, lpvar z) { lpvar vs[3] = { x, y, z }; add(v, 3, vs); }
+
+        /**
+           \brief retrieve monomial corresponding to variable v from definition v := vs
+        */
+        monomial const* var2monomial(lpvar v) const { unsigned idx = m_var2index.get(v, UINT_MAX); return idx == UINT_MAX ? nullptr : &m_monomials[idx]; }
+        
+        /**
+           \brief obtain a canonized signed monomial
+           corresponding to current equivalence class.
+        */
+        signed_vars const& canonize(monomial const& m) const;
+
+        /**
+           \brief determine if m1 divides m2 over the canonization obtained from merged variables.
+         */
+        bool canonize_divides(monomial const& m1, monomial const& m2) const;
+
+        /**
+           \brief produce explanation for monomial canonization.
+        */
+        void explain_canonized(monomial const& m, lp::explanation& exp);
+
+        /**
+           \brief iterator over monomials that are declared.
+        */
+        vector<monomial>::const_iterator begin() const { return m_monomials.begin(); }
+        vector<monomial>::const_iterator end() const { return m_monomials.end(); }
+
+        /**
+           \brief iterators over monomials where an equivalent variable is used
+        */
+        class iterator {
+            emonomials const& m;
+            cell*       m_cell;
+            bool        m_touched;            
+        public:
+            iterator(emonomials const& m, cell* c, bool at_end): m(m), m_cell(c), m_touched(at_end || c == nullptr) {}
+            monomial const& operator*() { return m.m_monomials[m_cell->m_index]; }
+            iterator& operator++() { m_touched = true; m_cell = m_cell->m_next; return *this; }
+            iterator operator++(int) { iterator tmp = *this; ++*this; return tmp; }
+            bool operator==(iterator const& other) const { return m_cell == other.m_cell && m_touched == other.m_touched; }
+            bool operator!=(iterator const& other) const { return m_cell != other.m_cell || m_touched != other.m_touched; }
+        };
+        
+        class use_list {
+            emonomials const& m;
+            lpvar     m_var;
+            cell*     head() { return m.head(m_var); } 
+        public:
+            use_list(emonomials const& m, lpvar v): m(m), m_var(v) {}
+            iterator begin() { return iterator(m, head(), false); }
+            iterator end() { return iterator(m, head(), true); }
+        };
+
+        use_list get_use_list(lpvar v) const { return use_list(*this, v); }
+
+
+        /**
+           \brief retrieve monomials m' where m is a proper factor of modulo current equalities.
+        */
+        class pf_iterator {
+            emonomials const& m;
+            monomial const&   m_mon; // canonized monomial
+            iterator          m_it;  // iterator over the first variable occurs list, ++ filters out elements that are not factors.
+            iterator          m_end;
+
+            void fast_forward();
+        public:
+            pf_iterator(emonomials const& m, monomial const& mon, bool at_end);
+            monomial const& operator*() { return *m_it; }
+            pf_iterator& operator++() { ++m_it; fast_forward(); return *this; }
+            pf_iterator operator++(int) { pf_iterator tmp = *this; ++*this; return tmp; }
+            bool operator==(pf_iterator const& other) const { return m_it == other.m_it; }
+            bool operator!=(pf_iterator const& other) const { return m_it != other.m_it; }
+        };
+
+        class factors_of {
+            emonomials const& m;
+            monomial const& mon;
+        public:
+            factors_of(emonomials const& m, monomial const& mon): m(m), mon(mon) {}
+            pf_iterator begin() { return pf_iterator(m, mon, false); }
+            pf_iterator end() { return pf_iterator(m, mon, true); }
+        };
+
+        factors_of get_factors_of(monomial const& m) const { inc_visited(); return factors_of(*this, m); }
+        
+        /**
+           \brief display state of emonomials
+        */
+        std::ostream& display(std::ostream& out) const;
+
+        /**
+           \brief
+           these are merge event handlers to interect the union-find handlers.
+           r2 becomes the new root. r2 is the root of v2, r1 is the old root of v1
+        */
+        void merge_eh(signed_var r2, signed_var r1, signed_var v2, signed_var v1) override;
+
+        void after_merge_eh(signed_var r2, signed_var r1, signed_var v2, signed_var v1) override;
+
+        void unmerge_eh(signed_var r2, signed_var r1) override;        
+
+    };
+
+    inline std::ostream& operator<<(std::ostream& out, emonomials const& m) { return m.display(out); }
+
+}

src/util/lp/incremental_vector.h

@@ -0,0 +1,76 @@
+/*++
+Copyright (c) 2017 Microsoft Corporation
+
+Module Name:
+
+    <name>
+
+Abstract:
+
+    <abstract>
+
+Author:
+
+    Lev Nachmanson (levnach)
+
+Revision History:
+
+
+--*/
+
+#pragma once
+#include <unordered_map>
+#include <set>
+#include <stack>
+#include "util/vector.h"
+namespace lp {
+template <typename B> class incremental_vector {
+    vector<unsigned> m_stack_of_vector_sizes;
+    vector<B> m_vector;
+public:
+    const B & operator[](unsigned a) const {
+        return m_vector[a];
+    }
+
+    unsigned size() const {
+        return m_vector.size();
+    }
+    
+    void push() {
+        m_stack_of_vector_sizes.push_back(m_vector.size());
+    }
+
+    void pop() {
+        pop(1);
+    }
+
+    template <typename T>  
+    void pop_tail(vector<T> & v, unsigned k) {
+        lp_assert(v.size() >= k);
+        v.shrink(v.size() - k);
+    }
+
+    template <typename T>  
+    void resize(vector<T> & v, unsigned new_size) {
+        v.resize(new_size);
+    }
+
+    void pop(unsigned k) {
+        lp_assert(m_stack_of_vector_sizes.size() >= k);
+        lp_assert(k > 0);
+        m_vector.shrink(peek_size(k));
+        unsigned new_st_size = m_stack_of_vector_sizes.size() - k;
+        m_stack_of_vector_sizes.shrink(k);
+    }
+
+    void push_back(const B & b) {
+        m_vector.push_back(b);
+    }
+
+    unsigned peek_size(unsigned k) const {
+        lp_assert(k > 0 && k <= m_stack_of_vector_sizes.size());
+        return m_stack_of_vector_sizes[m_stack_of_vector_sizes.size() - k];
+    }
+};
+}
+

src/util/lp/monomial.h

@@ -35,10 +35,22 @@ namespace nla {
         unsigned operator[](unsigned idx) const { return m_vs[idx]; }
         svector<lp::var_index>::const_iterator begin() const { return m_vs.begin(); }
         svector<lp::var_index>::const_iterator end() const { return m_vs.end(); }
-        const svector<lp::var_index> vars() const { return m_vs; }
+        const svector<lp::var_index>& vars() const { return m_vs; }
         bool empty() const { return m_vs.empty(); }
+
+        std::ostream& display(std::ostream& out) const {
+            out << "v" << var() << " := ";
+            for (auto v : *this) {
+                out << "v" << v << " ";
+            }
+            return out;
+        }
     };
 
+    inline std::ostream& operator<<(std::ostream& out, monomial const& m) {
+        return m.display(out);
+    }
+
     typedef std::unordered_map<lp::var_index, rational> variable_map_type;
     
     bool check_assignment(monomial const& m, variable_map_type & vars);

src/util/lp/nla_core.cpp

@@ -27,7 +27,8 @@ core::core(lp::lar_solver& s) :
     m_tangents(this),
     m_basics(this),
     m_order(this),
-    m_monotone(this) {
+    m_monotone(this),
+    m_emons(m_evars) {
 }
     
 bool core::compare_holds(const rational& ls, llc cmp, const rational& rs) const {
@@ -1608,25 +1609,24 @@ void core::init_search() {
     register_monomials_in_tables();
 }
 
-#if 0
+void core::init_to_refine() {
-void init_to_refine() {
     m_to_refine.clear();
     for (auto const & m : m_emons) 
         if (!check_monomial(m)) 
             m_to_refine.push_back(m.var());
     
     TRACE("nla_solver", 
-          tout << m_to_refine.size() << " mons to refine:\n");
+          tout << m_to_refine.size() << " mons to refine:\n";
-    for (unsigned v : m_to_refine) tout << m_emons.var2monomial(v) << "\n";
+          for (unsigned v : m_to_refine) tout << m_emons.var2monomial(v) << "\n";);
 }
-
+        
-std::unordered_set<lpvar> collect_vars(  const lemma& l) {
+std::unordered_set<lpvar> core::collect_vars(const lemma& l) const {
     std::unordered_set<lpvar> vars;
     auto insert_j = [&](lpvar j) { 
-                        vars.insert(j);
+        vars.insert(j);
-                        auto const* m = m_emons.var2monomial(j);
+        auto const* m = m_emons.var2monomial(j);
-                        if (m) for (lpvar k : *m) vars.insert(k);
+        if (m) for (lpvar k : *m) vars.insert(k);
-                    };
+    };
     
     for (const auto& i : current_lemma().ineqs()) {
         for (const auto& p : i.term()) {                
@@ -1641,21 +1641,6 @@ std::unordered_set<lpvar> collect_vars(  const lemma& l) {
     }
     return vars;
 }
-#endif
-
-void core::init_to_refine() {
-    m_to_refine.clear();
-    for (unsigned i = 0; i < m_monomials.size(); i++) {
-        if (!check_monomial(m_monomials[i]))
-            m_to_refine.push_back(i);
-    }
-    TRACE("nla_solver",
-          tout << m_to_refine.size() << " mons to refine:\n";
-          for (unsigned i: m_to_refine) {
-              print_monomial_with_vars(m_monomials[i], tout);
-          }
-          );
-}
 
 bool core:: divide(const rooted_mon& bc, const factor& c, factor & b) const {
     svector<lpvar> c_vars = sorted_vars(c);
@@ -1706,34 +1691,6 @@ void core::negate_factor_relation(const rational& a_sign, const factor& a, const
     mk_ineq(a_fs*a_sign, var(a), - b_fs*b_sign, var(b), cmp);
 }
 
-std::unordered_set<lpvar> core::collect_vars(const lemma& l) const {
-    std::unordered_set<lpvar> vars;
-    for (const auto& i : current_lemma().ineqs()) {
-        for (const auto& p : i.term()) {
-            lpvar j = p.var();
-            vars.insert(j);
-            auto it = m_var_to_its_monomial.find(j);
-            if (it != m_var_to_its_monomial.end()) {
-                for (lpvar k : m_monomials[it->second])
-                    vars.insert(k);
-            }
-        }
-    }
-    for (const auto& p : current_expl()) {
-        const auto& c = m_lar_solver.get_constraint(p.second);
-        for (const auto& r : c.coeffs()) {
-            lpvar j = r.second;
-            vars.insert(j);
-            auto it = m_var_to_its_monomial.find(j);
-            if (it != m_var_to_its_monomial.end()) {
-                for (lpvar k : m_monomials[it->second])
-                    vars.insert(k);
-            }
-        }
-    }
-    return vars;
-}
-
 std::ostream& core::print_lemma(std::ostream& out) const {
     print_specific_lemma(current_lemma(), out);
     return out;

src/util/lp/nla_core.h

@@ -26,6 +26,7 @@
 #include "util/lp/nla_basics_lemmas.h"
 #include "util/lp/nla_order_lemmas.h"
 #include "util/lp/nla_monotone_lemmas.h"
+#include "util/lp/emonomials.h"
 namespace nla {
 
 template <typename A, typename B>
@@ -79,7 +80,6 @@ struct core {
     var_eqs                                                          m_evars;
     vector<monomial>                                                 m_monomials;
 
-    rooted_mon_table                                                 m_rm_table;
     // this field is used for the push/pop operations
     unsigned_vector                                                  m_monomials_counts;
     lp::lar_solver&                                                  m_lar_solver;
@@ -94,6 +94,7 @@ struct core {
     basics                                                           m_basics;
     order                                                           m_order;
     monotone                                                        m_monotone;
+    emonomials                                                       m_emons;
     // methods
     core(lp::lar_solver& s);
     
@@ -115,7 +116,7 @@ struct core {
 
     rational vvr(const rooted_mon& rm) const { return vvr(m_monomials[rm.orig_index()].var()) * rm.orig_sign(); }
 
-    rational vvr(const factor& f) const {  return f.is_var()? vvr(f.index()) : vvr(m_rm_table.rms()[f.index()]); }
+   rational vvr(const factor& f) const {  return f.is_var()? vvr(f.index()) : vvr(m_rm_table.rms()[f.index()]); }
 
     lpvar var(const factor& f) const { return f.is_var()? f.index() : var(m_rm_table.rms()[f.index()]); }
 

src/util/lp/nla_defs.h

@@ -0,0 +1,103 @@
+/*++
+  Copyright (c) 2017 Microsoft Corporation
+
+  Module Name:
+
+  <name>
+
+  Abstract:
+
+  <abstract>
+
+  Author:
+  Nikolaj Bjorner (nbjorner)
+  Lev Nachmanson (levnach)
+
+  Revision History:
+
+
+  --*/
+#pragma once
+#include "util/lp/lp_types.h"
+#include "util/lp/column_info.h"
+#include "util/lp/explanation.h"
+
+namespace nla {
+typedef lp::constraint_index lpci;
+typedef lp::lconstraint_kind llc;
+typedef lp::constraint_index     lpci;
+typedef lp::explanation          expl_set;
+typedef lp::var_index            lpvar;
+
+struct from_index_dummy{};
+class signed_var {
+    unsigned m_sv;
+public:
+    // constructor, sign = true means minus
+    signed_var(lpvar v, bool sign): m_sv((v << 1) + (sign ? 1 : 0)) {}
+    // constructor
+    signed_var(unsigned sv, from_index_dummy): m_sv(sv) {}
+    bool sign() const { return 0 != (m_sv & 0x1); }
+    lpvar var() const { return m_sv >> 1; }
+    unsigned index() const { return m_sv; }        
+    void neg() { m_sv = m_sv ^ 1; }
+    friend signed_var operator~(signed_var const& sv) {
+        return signed_var(sv.var(), !sv.sign());
+    }
+    bool operator==(signed_var const& other) const {
+        return m_sv == other.m_sv;
+    }
+    bool operator!=(signed_var const& other) const {
+        return m_sv != other.m_sv;
+    }
+    rational rsign() const { return sign() ? rational::minus_one() : rational::one(); }
+
+    std::ostream& display(std::ostream& out) const {
+        return out << (sign()?"-":"") << var();
+    }
+};
+
+inline std::ostream& operator<<(std::ostream& out, signed_var const& sv) { return sv.display(out); }
+
+/*
+ *  represents definition m_v = coeff* v1*v2*...*vn, 
+ *  where m_vs = [v1, v2, .., vn]
+ */
+class monomial_coeff  {
+    svector<lp::var_index> m_vs;
+    rational m_coeff;
+public:
+    monomial_coeff(const svector<lp::var_index>& vs, rational const& coeff): m_vs(vs), m_coeff(coeff) {}
+    rational const& coeff() const { return m_coeff; }
+    const svector<lp::var_index> & vars() const { return m_vs; } 
+};
+template <typename T> bool has_zero(const T& product) {
+    for (const rational & t : product) {
+        if (t.is_zero())
+            return true;
+    }
+    return false;
+}
+template <typename T>
+bool uniform_le(const T& a, const T& b,  unsigned & strict_i) {
+    SASSERT(a.size() == b.size());
+    strict_i = -1;
+    bool z_b = false;
+        
+    for (unsigned i = 0; i < a.size(); i++) {
+        if (a[i] > b[i]){
+            return false;
+        }
+            
+        SASSERT(!a[i].is_neg());
+        if (a[i] < b[i]){
+            strict_i = i;
+        } else if (b[i].is_zero()) {
+            z_b = true;
+        }
+    }
+    if (z_b) {strict_i = -1;}
+    return true;
+}
+
+}

src/util/lp/var_eqs.cpp

@@ -23,164 +23,164 @@
 namespace nla {
     
 
-    var_eqs::var_eqs(): m_uf(m_ufctx) {}
+var_eqs::var_eqs(): m_merge_handler(nullptr), m_uf(*this), m_stack(*this) {}
     
-    void var_eqs::push() {
+void var_eqs::push() {
-        m_trail_lim.push_back(m_trail.size());
+    m_trail_lim.push_back(m_trail.size());
-        m_ufctx.get_trail_stack().push_scope();
+    get_trail_stack().push_scope();
-    }
+}
 
-    void var_eqs::pop(unsigned n) {
+void var_eqs::pop(unsigned n) {
-        unsigned old_sz = m_trail_lim[m_trail_lim.size() - n];
+    unsigned old_sz = m_trail_lim[m_trail_lim.size() - n];
-        for (unsigned i = m_trail.size(); i-- > old_sz; ) {
+    for (unsigned i = m_trail.size(); i-- > old_sz; ) {
-            auto const& sv = m_trail[i];
+        auto const& sv = m_trail[i];
-            m_eqs[sv.first.index()].pop_back();
+        m_eqs[sv.first.index()].pop_back();
-            m_eqs[sv.second.index()].pop_back();
+        m_eqs[sv.second.index()].pop_back();
-            m_eqs[(~sv.first).index()].pop_back();
+        m_eqs[(~sv.first).index()].pop_back();
-            m_eqs[(~sv.second).index()].pop_back();
+        m_eqs[(~sv.second).index()].pop_back();
+    }
+    m_trail_lim.shrink(m_trail_lim.size() - n);
+    m_trail.shrink(old_sz);
+    get_trail_stack().pop_scope(n);
+}
+
+void var_eqs::merge(signed_var v1, signed_var v2, eq_justification const& j) {
+    unsigned max_i = std::max(v1.index(), v2.index()) + 2;
+    m_eqs.reserve(max_i);
+    while (m_uf.get_num_vars() <= max_i) m_uf.mk_var();
+    m_trail.push_back(std::make_pair(v1, v2));
+    m_uf.merge(v1.index(), v2.index());
+    m_uf.merge((~v1).index(), (~v2).index());
+    m_eqs[v1.index()].push_back(justified_var(v2, j));
+    m_eqs[v2.index()].push_back(justified_var(v1, j));
+    m_eqs[(~v1).index()].push_back(justified_var(~v2, j));
+    m_eqs[(~v2).index()].push_back(justified_var(~v1, j));
+}
+
+signed_var var_eqs::find(signed_var v) const {
+    if (v.index() >= m_uf.get_num_vars()) {
+        return v;
+    }
+    unsigned idx = m_uf.find(v.index());
+    return signed_var(idx, from_index_dummy());
+}
+
+void var_eqs::explain_dfs(signed_var v1, signed_var v2, lp::explanation& e) const {
+    SASSERT(find(v1) == find(v2));
+    if (v1 == v2) {
+        return;
+    }
+    m_todo.push_back(var_frame(v1, 0));
+    m_jtrail.reset();
+    m_marked.reserve(m_eqs.size(), false);
+    SASSERT(m_marked_trail.empty());
+    m_marked[v1.index()] = true;
+    m_marked_trail.push_back(v1.index());
+    while (true) {
+        SASSERT(!m_todo.empty());
+        var_frame& f = m_todo.back();
+        signed_var v = f.m_var;
+        if (v == v2) {
+            break;
         }
-        m_trail_lim.shrink(m_trail_lim.size() - n);
+        auto const& next = m_eqs[v.index()];
-        m_trail.shrink(old_sz);
+        bool seen_all = true;
-        m_ufctx.get_trail_stack().pop_scope(n);
+        unsigned sz = next.size();
-    }
+        for (unsigned i = f.m_index; seen_all && i < sz; ++i) {
-
+            justified_var const& jv = next[i];
-    void var_eqs::merge(signed_var v1, signed_var v2, eq_justification const& j) {
+            signed_var v3 = jv.m_var;
-        unsigned max_i = std::max(v1.index(), v2.index()) + 2;
+            if (!m_marked[v3.index()]) {
-        m_eqs.reserve(max_i);
+                seen_all = false;
-        while (m_uf.get_num_vars() <= max_i) m_uf.mk_var();
+                f.m_index = i + 1;
-        m_trail.push_back(std::make_pair(v1, v2));
+                m_todo.push_back(var_frame(v3, 0));
-        m_uf.merge(v1.index(), v2.index());
+                m_jtrail.push_back(jv.m_j);
-        m_uf.merge((~v1).index(), (~v2).index());
+                m_marked_trail.push_back(v3.index());
-        m_eqs[v1.index()].push_back(justified_var(v2, j));
+                m_marked[v3.index()] = true;
-        m_eqs[v2.index()].push_back(justified_var(v1, j));
-        m_eqs[(~v1).index()].push_back(justified_var(~v2, j));
-        m_eqs[(~v2).index()].push_back(justified_var(~v1, j));
-    }
-
-    signed_var var_eqs::find(signed_var v) const {
-        if (v.index() >= m_uf.get_num_vars()) {
-            return v;
-        }
-        unsigned idx = m_uf.find(v.index());
-        return signed_var(idx, from_index());
-    }
-
-    void var_eqs::explain_dfs(signed_var v1, signed_var v2, lp::explanation& e) const {
-        SASSERT(find(v1) == find(v2));
-        if (v1 == v2) {
-            return;
-        }
-        m_todo.push_back(var_frame(v1, 0));
-        m_jtrail.reset();
-        m_marked.reserve(m_eqs.size(), false);
-        SASSERT(m_marked_trail.empty());
-        m_marked[v1.index()] = true;
-        m_marked_trail.push_back(v1.index());
-        while (true) {
-            SASSERT(!m_todo.empty());
-            var_frame& f = m_todo.back();
-            signed_var v = f.m_var;
-            if (v == v2) {
-                break;
-            }
-            auto const& next = m_eqs[v.index()];
-            bool seen_all = true;
-            unsigned sz = next.size();
-            for (unsigned i = f.m_index; seen_all && i < sz; ++i) {
-                justified_var const& jv = next[i];
-                signed_var v3 = jv.m_var;
-                if (!m_marked[v3.index()]) {
-                    seen_all = false;
-                    f.m_index = i + 1;
-                    m_todo.push_back(var_frame(v3, 0));
-                    m_jtrail.push_back(jv.m_j);
-                    m_marked_trail.push_back(v3.index());
-                    m_marked[v3.index()] = true;
-                }
-            }
-            if (seen_all) {
-                m_todo.pop_back();
-                m_jtrail.pop_back();
             }
         }
+        if (seen_all) {
+            m_todo.pop_back();
+            m_jtrail.pop_back();
+        }
+    }
         
-        for (eq_justification const& j : m_jtrail) {
+    for (eq_justification const& j : m_jtrail) {
-            j.explain(e);
+        j.explain(e);
-        }
-        m_stats.m_num_explains += m_jtrail.size();
-        m_stats.m_num_explain_calls++;
-        m_todo.reset();
-        m_jtrail.reset();
-        for (unsigned idx : m_marked_trail) {
-            m_marked[idx] = false;
-        }
-        m_marked_trail.reset();
-
-        // IF_VERBOSE(2, verbose_stream() << (double)m_stats.m_num_explains / m_stats.m_num_explain_calls << "\n");
     }
+    m_stats.m_num_explains += m_jtrail.size();
+    m_stats.m_num_explain_calls++;
+    m_todo.reset();
+    m_jtrail.reset();
+    for (unsigned idx : m_marked_trail) {
+        m_marked[idx] = false;
+    }
+    m_marked_trail.reset();
 
-    void var_eqs::explain_bfs(signed_var v1, signed_var v2, lp::explanation& e) const {
+    // IF_VERBOSE(2, verbose_stream() << (double)m_stats.m_num_explains / m_stats.m_num_explain_calls << "\n");
-        SASSERT(find(v1) == find(v2));
+}
-        if (v1 == v2) {
+
-            return;
+void var_eqs::explain_bfs(signed_var v1, signed_var v2, lp::explanation& e) const {
+    SASSERT(find(v1) == find(v2));
+    if (v1 == v2) {
+        return;
+    }
+    m_todo.push_back(var_frame(v1, 0));
+    m_jtrail.push_back(eq_justification({}));
+    m_marked.reserve(m_eqs.size(), false);
+    SASSERT(m_marked_trail.empty());
+    m_marked[v1.index()] = true;
+    m_marked_trail.push_back(v1.index());
+    unsigned head = 0;
+    for (; ; ++head) {
+        var_frame& f = m_todo[head];
+        signed_var v = f.m_var;
+        if (v == v2) {
+            break;
         }
-        m_todo.push_back(var_frame(v1, 0));
+        auto const& next = m_eqs[v.index()];
-        m_jtrail.push_back(eq_justification({}));
+        unsigned sz = next.size();
-        m_marked.reserve(m_eqs.size(), false);
+        for (unsigned i = sz; i-- > 0; ) {
-        SASSERT(m_marked_trail.empty());
+            justified_var const& jv = next[i];
-        m_marked[v1.index()] = true;
+            signed_var v3 = jv.m_var;
-        m_marked_trail.push_back(v1.index());
+            if (!m_marked[v3.index()]) {
-        unsigned head = 0;
+                m_todo.push_back(var_frame(v3, head));
-        for (; ; ++head) {
+                m_jtrail.push_back(jv.m_j);
-            var_frame& f = m_todo[head];
+                m_marked_trail.push_back(v3.index());
-            signed_var v = f.m_var;
+                m_marked[v3.index()] = true;
-            if (v == v2) {
-                break;
-            }
-            auto const& next = m_eqs[v.index()];
-            unsigned sz = next.size();
-            for (unsigned i = sz; i-- > 0; ) {
-                justified_var const& jv = next[i];
-                signed_var v3 = jv.m_var;
-                if (!m_marked[v3.index()]) {
-                    m_todo.push_back(var_frame(v3, head));
-                    m_jtrail.push_back(jv.m_j);
-                    m_marked_trail.push_back(v3.index());
-                    m_marked[v3.index()] = true;
-                }
             }
         }
+    }
         
-        while (head != 0) {
+    while (head != 0) {
-            m_jtrail[head].explain(e);
+        m_jtrail[head].explain(e);
-            head = m_todo[head].m_index;
+        head = m_todo[head].m_index;
-            ++m_stats.m_num_explains;
+        ++m_stats.m_num_explains;
-        }
+    }
-        ++m_stats.m_num_explain_calls;
+    ++m_stats.m_num_explain_calls;
         
-        m_todo.reset();
+    m_todo.reset();
-        m_jtrail.reset();
+    m_jtrail.reset();
-        for (unsigned idx : m_marked_trail) {
+    for (unsigned idx : m_marked_trail) {
-            m_marked[idx] = false;
+        m_marked[idx] = false;
-        }
-        m_marked_trail.reset();
-
-        // IF_VERBOSE(2, verbose_stream() << (double)m_stats.m_num_explains / m_stats.m_num_explain_calls << "\n");
     }
+    m_marked_trail.reset();
+
+    // IF_VERBOSE(2, verbose_stream() << (double)m_stats.m_num_explains / m_stats.m_num_explain_calls << "\n");
+}
 
 
-    std::ostream& var_eqs::display(std::ostream& out) const {
+std::ostream& var_eqs::display(std::ostream& out) const {
-        m_uf.display(out);
+    m_uf.display(out);
-        unsigned idx = 0;
+    unsigned idx = 0;
-        for (auto const& edges : m_eqs) {
+    for (auto const& edges : m_eqs) {
-            out << signed_var(idx, from_index()) << ": ";
+        out << signed_var(idx, from_index_dummy()) << ": ";
-            for (auto const& jv : edges) {
+        for (auto const& jv : edges) {
-                out << jv.m_var << " ";
+            out << jv.m_var << " ";
-            }
-            out << "\n";
-            ++idx;
         }
-        return out;
+        out << "\n";
+        ++idx;
     }
+    return out;
+}
 
 }

src/util/lp/var_eqs.h

@@ -24,41 +24,6 @@
 #include "util/rational.h"
 #include "util/lp/explanation.h"
 namespace nla {
-    
-typedef lp::var_index lpvar;
-typedef lp::constraint_index lpci;
-struct from_index{};
-
-class signed_var {
-    unsigned m_sv;
-public:
-    // constructor, sign = true means minus
-    signed_var(lpvar v, bool sign): m_sv((v << 1) + (sign ? 1 : 0)) {}
-    // constructor
-    signed_var(unsigned sv, from_index): m_sv(sv) {}
-    bool sign() const { return 0 != (m_sv & 0x1); }
-    lpvar var() const { return m_sv >> 1; }
-    unsigned index() const { return m_sv; }        
-    void neg() { m_sv = m_sv ^ 1; }
-    friend signed_var operator~(signed_var const& sv) {
-        return signed_var(sv.var(), !sv.sign());
-    }
-    bool operator==(signed_var const& other) const {
-        return m_sv == other.m_sv;
-    }
-    bool operator!=(signed_var const& other) const {
-        return m_sv != other.m_sv;
-    }
-    rational rsign() const { return sign() ? rational::minus_one() : rational::one(); }
-
-    std::ostream& display(std::ostream& out) const {
-        return out << (sign()?"-":"") << var();
-    }
-};
-
-    inline std::ostream& operator<<(std::ostream& out, signed_var const& sv) {
-        return sv.display(out);
-    }
 
 class eq_justification {
     lpci m_cs[4];
@@ -80,6 +45,16 @@ public:
     }
 };
 
+ class var_eqs_merge_handler {
+ public:
+     virtual void merge_eh(signed_var r2, signed_var r1, signed_var v2, signed_var v1) = 0;
+
+     virtual void after_merge_eh(signed_var r2, signed_var r1, signed_var v2, signed_var v1) = 0;
+
+     virtual void unmerge_eh(signed_var r2, signed_var r1) = 0;;        
+
+ };
+
 class var_eqs {
     struct justified_var {
         signed_var       m_var;
@@ -101,12 +76,14 @@ class var_eqs {
 
     typedef std::pair<signed_var, signed_var> signed_var_pair;
 
-    union_find_default_ctx  m_ufctx;
+    var_eqs_merge_handler*  m_merge_handler;    
-    union_find<>            m_uf;
+    union_find<var_eqs>         m_uf;
     svector<signed_var_pair>    m_trail;
     unsigned_vector         m_trail_lim;
     vector<justified_vars>  m_eqs;    // signed_var-index -> justified_var corresponding to edges in a graph used to extract justifications.
 
+    typedef trail_stack<var_eqs>    trail_stack_t;       
+    trail_stack_t                   m_stack;
     mutable svector<var_frame>      m_todo;
     mutable svector<bool>           m_marked;
     mutable unsigned_vector         m_marked_trail;
@@ -197,8 +174,32 @@ public:
 
     eq_class equiv_class(lpvar v) { return equiv_class(signed_var(v, false)); }
 
+
     std::ostream& display(std::ostream& out) const;
-    
+
+    // union find event handlers
+    void set_merge_handler(var_eqs_merge_handler* mh) { m_merge_handler = mh; }
+
+    trail_stack_t& get_trail_stack() { return m_stack; }
+
+    void unmerge_eh(unsigned i, unsigned j) { 
+        if (m_merge_handler) {
+            m_merge_handler->unmerge_eh(signed_var(i, from_index_dummy()), signed_var(j, from_index_dummy())); 
+        }
+    }
+    void merge_eh(unsigned r2, unsigned r1, unsigned v2, unsigned v1) {
+        if (m_merge_handler) {
+            m_merge_handler->merge_eh(signed_var(r2, from_index_dummy()), signed_var(r1, from_index_dummy()), signed_var(v2, from_index_dummy()), signed_var(v1, from_index_dummy())); 
+        }
+    }
+
+    void after_merge_eh(unsigned r2, unsigned r1, unsigned v2, unsigned v1) {
+        if (m_merge_handler) {
+            m_merge_handler->after_merge_eh(signed_var(r2, from_index_dummy()), signed_var(r1, from_index_dummy()), signed_var(v2, from_index_dummy()), signed_var(v1, from_index_dummy())); 
+        }
+    }
+
+
 }; 
 
     inline std::ostream& operator<<(var_eqs const& ve, std::ostream& out) { return ve.display(out); }