diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp
index 538e0d5e6..da1204685 100644
--- a/src/smt/theory_lra.cpp
+++ b/src/smt/theory_lra.cpp
@@ -2092,7 +2092,7 @@ public:
         for (auto const& mon : poly) {
             SASSERT(!mon.empty());
             if (mon.size() == 1) {
-                term.add_coeff_var(mon.get_coeff(), mon[0]);
+                term.add_var(mon[0]);
             }
             else {
                 // create the expression corresponding to the product.
@@ -2107,7 +2107,7 @@ public:
                 app_ref t(a.mk_mul(mul.size(), mul.c_ptr()), m);
                 VERIFY(internalize_term(t));
                 theory_var w = ctx().get_enode(t)->get_th_var(get_id());
-                term.add_coeff_var(mon.get_coeff(), lp().external_to_local(w));
+                term.add_var(lp().external_to_local(w));
             }
         }
         return term;
diff --git a/src/util/lp/CMakeLists.txt b/src/util/lp/CMakeLists.txt
index 356994e09..a1d942be6 100644
--- a/src/util/lp/CMakeLists.txt
+++ b/src/util/lp/CMakeLists.txt
@@ -32,6 +32,7 @@ z3_add_component(lp
     square_dense_submatrix.cpp
     square_sparse_matrix.cpp
     static_matrix.cpp
+    var_eqs.cpp
   COMPONENT_DEPENDENCIES
     util
     polynomial
diff --git a/src/util/lp/explanation.h b/src/util/lp/explanation.h
index c7fb22120..9388e1d53 100644
--- a/src/util/lp/explanation.h
+++ b/src/util/lp/explanation.h
@@ -18,6 +18,8 @@ Revision History:
 
 --*/
 #pragma once
+#include <unordered_set>
+#include "util/lp/lp_utils.h"
 namespace lp {
 class explanation {
     vector<std::pair<mpq, constraint_index>> m_explanation;
diff --git a/src/util/lp/lp_settings.h b/src/util/lp/lp_settings.h
index 14477af2a..fb8c103c3 100644
--- a/src/util/lp/lp_settings.h
+++ b/src/util/lp/lp_settings.h
@@ -26,11 +26,9 @@ Revision History:
 #include <iomanip>
 #include "util/lp/lp_utils.h"
 #include "util/stopwatch.h"
+#include "util/lp/lp_types.h"
 
 namespace lp {
-typedef unsigned var_index;
-typedef unsigned constraint_index;
-typedef unsigned row_index;
 
 enum class column_type  {
     free_column = 0,
diff --git a/src/util/lp/lp_types.h b/src/util/lp/lp_types.h
new file mode 100644
index 000000000..ced3fb4c5
--- /dev/null
+++ b/src/util/lp/lp_types.h
@@ -0,0 +1,25 @@
+/*++
+Copyright (c) 2017 Microsoft Corporation
+
+Module Name:
+
+    <name>
+
+Abstract:
+
+    <abstract>
+
+Author:
+
+    Lev Nachmanson (levnach)
+
+Revision History:
+
+
+--*/
+#pragma once
+namespace lp {
+typedef unsigned var_index;
+typedef unsigned constraint_index;
+typedef unsigned row_index;
+}
diff --git a/src/util/lp/monomial.h b/src/util/lp/monomial.h
index 408878dae..d5c0bf24f 100644
--- a/src/util/lp/monomial.h
+++ b/src/util/lp/monomial.h
@@ -18,7 +18,6 @@ namespace nla {
         // fields
         lp::var_index          m_v;
         svector<lp::var_index> m_vs;
-        rational               m_coeff;
     public:
         // constructors
         monomial(lp::var_index v, unsigned sz, lp::var_index const* vs):
@@ -38,7 +37,6 @@ namespace nla {
         svector<lp::var_index>::const_iterator end() const { return m_vs.end(); }
         const svector<lp::var_index> vars() const { return m_vs; }
         bool empty() const { return m_vs.empty(); }
-        const rational& get_coeff() const { return m_coeff; } 
     };
 
     typedef std::unordered_map<lp::var_index, rational> variable_map_type;
diff --git a/src/util/lp/nla_solver.cpp b/src/util/lp/nla_solver.cpp
index cf265fb46..40d674231 100644
--- a/src/util/lp/nla_solver.cpp
+++ b/src/util/lp/nla_solver.cpp
@@ -21,7 +21,7 @@
 #include <map>
 #include "util/lp/monomial.h"
 #include "util/lp/lp_utils.h"
-#include "util/lp/vars_equivalence.h"
+#include "util/lp/var_eqs.h"
 #include "util/lp/factorization.h"
 #include "util/lp/rooted_mons.h"
 namespace nla {
@@ -34,6 +34,7 @@ bool try_insert(const A& elem, B& collection) {
     collection.insert(elem);
     return true;
 }
+typedef lp::constraint_index lpci;
 
 typedef lp::lconstraint_kind llc;
 struct point {
@@ -85,7 +86,7 @@ struct solver::imp {
     };
 
     //fields    
-    vars_equivalence                                                 m_vars_equivalence;
+    var_eqs                                                          m_vars_equivalence;
     vector<monomial>                                                 m_monomials;
 
     rooted_mon_table                                                 m_rm_table;
@@ -113,7 +114,7 @@ struct solver::imp {
     
     imp(lp::lar_solver& s, reslimit& lim, params_ref const& p)
         :
-        m_vars_equivalence([this](unsigned h){return vvr(h);}),
+        m_vars_equivalence(),
         m_lar_solver(s)
         // m_limit(lim),
         // m_params(p)
@@ -193,21 +194,19 @@ struct solver::imp {
     }
 
     // the value of the factor is equal to the value of the variable multiplied
-    // by the flip_sign
-    rational flip_sign(const factor& f) const {
+    // by the canonize_sign
+    rational canonize_sign(const factor& f) const {
         return f.is_var()?
-            flip_sign_of_var(f.index()) : m_rm_table.rms()[f.index()].orig_sign();
+            canonize_sign_of_var(f.index()) : m_rm_table.rms()[f.index()].orig_sign();
     }
 
-    rational flip_sign_of_var(lpvar j) const {
-        rational sign(1);
-        m_vars_equivalence.map_to_root(j, sign);
-        return sign;
+    rational canonize_sign_of_var(lpvar j) const {
+        return m_vars_equivalence.find_sign(j);        
     }
     
     // the value of the rooted monomias is equal to the value of the variable multiplied
-    // by the flip_sign
-    rational flip_sign(const rooted_mon& m) const {
+    // by the canonize_sign
+    rational canonize_sign(const rooted_mon& m) const {
         return m.orig().sign();
     }
     
@@ -269,8 +268,9 @@ struct solver::imp {
         return product_value(m) == m_lar_solver.get_column_value_rational(m.var());
     }
     
-    void explain(const monomial& m, lp::explanation& exp) const {
-        m_vars_equivalence.explain(m, exp);
+    void explain(const monomial& m, lp::explanation& exp) const {       
+        for (lpvar j : m)
+            explain(j, exp);
     }
 
     void explain(const rooted_mon& rm, lp::explanation& exp) const {
@@ -280,9 +280,9 @@ struct solver::imp {
 
     void explain(const factor& f, lp::explanation& exp) const {
         if (f.type() == factor_type::VAR) {
-            m_vars_equivalence.explain(f.index(), exp);
+            explain(f.index(), exp);
         } else {
-            m_vars_equivalence.explain(m_monomials[m_rm_table.rms()[f.index()].orig_index()], exp);
+            explain(m_monomials[m_rm_table.rms()[f.index()].orig_index()], exp);
         }
     }
 
@@ -618,7 +618,7 @@ struct solver::imp {
         svector<lpvar> ret;
         sign = 1;
         for (lpvar v : vars) {
-            unsigned root = m_vars_equivalence.map_to_root(v, sign);
+            unsigned root = m_vars_equivalence.find(v, sign);
             SASSERT(m_vars_equivalence.is_root(root));
             TRACE("nla_solver_eq",
                   print_var(v,tout);
@@ -841,6 +841,7 @@ struct solver::imp {
         return  rat_sign(vvr(m)) != rat_sign(m);
     }
 
+    /*
     unsigned_vector eq_vars(lpvar j) const {
         TRACE("nla_solver_eq", tout << "j = "; print_var(j, tout); tout << "eqs = ";
               for(auto jj : m_vars_equivalence.eq_vars(j)) {
@@ -848,7 +849,7 @@ struct solver::imp {
               });
         return m_vars_equivalence.eq_vars(j);
     }
-
+    */
     // Monomials m and n vars have the same values, up to "sign"
     // Generate a lemma if values of m.var() and n.var() are not the same up to sign
     bool basic_sign_lemma_on_two_monomials(const monomial& m, const monomial& n, const rational& sign) {
@@ -1090,7 +1091,7 @@ struct solver::imp {
             if (!(var_has_positive_lower_bound(j) || var_has_negative_upper_bound(j))) {
                 return 0;
             }
-            m_vars_equivalence.map_to_root(j, sign);
+            sign *= m_vars_equivalence.find_sign(j);
         }
         return rat_sign(sign);
     }
@@ -1338,10 +1339,10 @@ struct solver::imp {
 
     bool vars_are_equiv(lpvar a, lpvar b) const {
         SASSERT(abs(vvr(a)) == abs(vvr(b)));
-        return m_vars_equivalence.vars_are_equiv(a, b) ||
-            (var_is_fixed(a) && var_is_fixed(b));
+        return m_vars_equivalence.vars_are_equiv(a, b);
     }
 
+    
     void explain_equiv_vars(lpvar a, lpvar b) {
         SASSERT(abs(vvr(a)) == abs(vvr(b)));
         if (m_vars_equivalence.vars_are_equiv(a, b)) {
@@ -1455,7 +1456,7 @@ struct solver::imp {
         for (auto j : f){
             lpvar var_j = var(j);
             if (not_one == var_j) continue;
-            mk_ineq(var_j, llc::NE, j.is_var()? vvr(j) : flip_sign(j) * vvr(j));   
+            mk_ineq(var_j, llc::NE, j.is_var()? vvr(j) : canonize_sign(j) * vvr(j));   
         }
 
         if (not_one == static_cast<lpvar>(-1)) {
@@ -1549,7 +1550,7 @@ struct solver::imp {
         for (auto j : f){
             lpvar var_j = var(j);
             if (not_one == var_j) continue;
-            mk_ineq(var_j, llc::NE, j.is_var()? vvr(j) : flip_sign(j) * vvr(j));   
+            mk_ineq(var_j, llc::NE, j.is_var()? vvr(j) : canonize_sign(j) * vvr(j));   
         }
 
         if (not_one == static_cast<lpvar>(-1)) {
@@ -1837,8 +1838,40 @@ struct solver::imp {
                 add_equivalence_maybe(s.terms()[i], s.get_column_upper_bound_witness(j), s.get_column_lower_bound_witness(j));
             }
         }
+        collect_equivs_of_fixed_vars();
     }
 
+    void collect_equivs_of_fixed_vars() {
+        std::unordered_map<rational, svector<lpvar> > abs_map;
+        for (lpvar j = 0; j < m_lar_solver.number_of_vars(); j++) {
+            if (!var_is_fixed(j))
+                continue;
+            rational v = abs(vvr(j));
+            auto it = abs_map.find(v);
+            if (it == abs_map.end()) {
+                abs_map[v] = svector<lpvar>();
+                abs_map[v].push_back(j);
+            } else {
+                it->second.push_back(j);
+            }
+        }
+        for (auto p : abs_map) {
+            svector<lpvar>& v = p.second;
+            lpvar head = v[0];
+            auto c0 = m_lar_solver.get_column_upper_bound_witness(head);
+            auto c1 = m_lar_solver.get_column_lower_bound_witness(head);
+            for (unsigned k = 1; k < v.size(); k++) {
+                auto c2 = m_lar_solver.get_column_upper_bound_witness(v[k]);
+                auto c3 = m_lar_solver.get_column_lower_bound_witness(v[k]);
+                if (vvr(head) == vvr(v[k]))
+                    m_vars_equivalence.merge_plus(head, v[k], eq_justification({c0, c1, c2, c3}));
+                else 
+                    m_vars_equivalence.merge_minus(head, v[k], eq_justification({c0, c1, c2, c3}));
+            }
+        }
+    }
+    
+    
     void add_equivalence_maybe(const lp::lar_term *t, lpci c0, lpci c1) {
         if (t->size() != 2)
             return;
@@ -1860,16 +1893,18 @@ struct solver::imp {
                 j = p.var();
         }
         SASSERT(j != static_cast<unsigned>(-1));
-        rational sign((seen_minus && seen_plus)? 1 : -1);
-        m_vars_equivalence.add_equiv(i, j, sign, c0, c1);
+        bool sign = (seen_minus && seen_plus)? false : true;
+        if (sign)
+            m_vars_equivalence.merge_minus(i, j, eq_justification({c0, c1}));
+        else 
+            m_vars_equivalence.merge_plus(i, j, eq_justification({c0, c1}));
     }
 
     // x is equivalent to y if x = +- y
     void init_vars_equivalence() {
-        SASSERT(m_vars_equivalence.empty());
+        /*       SASSERT(m_vars_equivalence.empty());*/
         collect_equivs();
-        m_vars_equivalence.create_tree();
-        TRACE("nla_solver_details", tout << "number of equivs = " << m_vars_equivalence.size(););
+        /*        TRACE("nla_solver_details", tout << "number of equivs = " << m_vars_equivalence.size(););*/
         
         SASSERT((settings().random_next() % 100) || tables_are_ok());
     }
@@ -1976,7 +2011,6 @@ struct solver::imp {
 
     void clear() {
         m_var_to_its_monomial.clear();
-        m_vars_equivalence.clear();
         m_rm_table.clear();
         m_monomials_containing_var.clear();
         m_lemma_vec->clear();
@@ -2046,14 +2080,14 @@ struct solver::imp {
     }
     
     void negate_factor_relation(const rational& a_sign, const factor& a, const rational& b_sign, const factor& b) {
-        rational a_fs = flip_sign(a);
-        rational b_fs = flip_sign(b);
+        rational a_fs = canonize_sign(a);
+        rational b_fs = canonize_sign(b);
         llc cmp = a_sign*vvr(a) < b_sign*vvr(b)? llc::GE : llc::LE;
         mk_ineq(a_fs*a_sign, var(a), - b_fs*b_sign, var(b), cmp);
     }
 
     void negate_var_factor_relation(const rational& a_sign, lpvar a, const rational& b_sign, const factor& b) {
-        rational b_fs = flip_sign(b);
+        rational b_fs = canonize_sign(b);
         llc cmp = a_sign*vvr(a) < b_sign*vvr(b)? llc::GE : llc::LE;
         mk_ineq(a_sign, a, - b_fs*b_sign, var(b), cmp);
     }
@@ -2068,8 +2102,8 @@ struct solver::imp {
                      const factor& b,
                      llc ab_cmp) {
          add_empty_lemma();
-        mk_ineq(rational(c_sign) * flip_sign(c), var(c), llc::LE);
-        mk_ineq(flip_sign(ac), var(ac), -flip_sign(bc), var(bc), ab_cmp);
+        mk_ineq(rational(c_sign) * canonize_sign(c), var(c), llc::LE);
+        mk_ineq(canonize_sign(ac), var(ac), -canonize_sign(bc), var(bc), ab_cmp);
         explain(ac, current_expl());
         explain(a, current_expl());
         explain(bc, current_expl());
@@ -2091,7 +2125,7 @@ struct solver::imp {
         mk_ineq(c_sign, c, llc::LE);
         explain(c, current_expl()); // this explains c == +- d
         negate_var_factor_relation(c_sign, a, d_sign, b);
-        mk_ineq(ac.var(), -flip_sign(bd), var(bd), ab_cmp);
+        mk_ineq(ac.var(), -canonize_sign(bd), var(bd), ab_cmp);
         explain(bd, current_expl());
         explain(b, current_expl());
         explain(d, current_expl());
@@ -2210,7 +2244,7 @@ struct solver::imp {
         SASSERT(abs(vvr(i)) == abs(vvr(c)));
         auto it = m_var_to_its_monomial.find(i);
         if (it == m_var_to_its_monomial.end()) {
-            i = m_vars_equivalence.map_to_root(i);
+            i = m_vars_equivalence.find(i);
             if (try_insert(i, found_vars)) {
                 r.push_back(factor(i, factor_type::VAR));
             }
@@ -2252,7 +2286,7 @@ struct solver::imp {
         TRACE("nla_solver", tout << "orig_sign = " << rm.orig_sign() << "\n";);
         rational sign = rm.orig_sign();
         for(factor f: ab)
-            sign *= flip_sign(f);
+            sign *= canonize_sign(f);
         const rational & fv = vvr(ab[0]) * vvr(ab[1]);
         const rational mv = sign * vvr(m);
         TRACE("nla_solver",
@@ -2348,7 +2382,7 @@ struct solver::imp {
     void order_lemma_on_factor_binomial_explore(const monomial& m, unsigned k) {
         SASSERT(m.size() == 2);
         lpvar c = m[k];
-        lpvar d = m_vars_equivalence.map_to_root(c);
+        lpvar d = m_vars_equivalence.find(c);
         auto it = m_rm_table.var_map().find(d);
         SASSERT(it != m_rm_table.var_map().end());
         for (unsigned bd_i : it->second) {
@@ -2359,7 +2393,7 @@ struct solver::imp {
     }
 
     void order_lemma_on_factor_binomial_rm(const monomial& ac, unsigned k, const rooted_mon& rm_bd) {
-        factor d(m_vars_equivalence.map_to_root(ac[k]), factor_type::VAR);
+        factor d(m_vars_equivalence.find(ac[k]), factor_type::VAR);
         factor b;
         if (!divide(rm_bd, d, b))
             return;
@@ -2373,7 +2407,7 @@ struct solver::imp {
         int p = (k + 1) % 2;
         lpvar a = ac[p];
         lpvar c = ac[k];
-        SASSERT(m_vars_equivalence.map_to_root(c) == d);
+        SASSERT(m_vars_equivalence.find(c) == d);
         rational acv = vvr(ac);
         rational av = vvr(a);
         rational c_sign = rrat_sign(vvr(c));
diff --git a/src/util/lp/rooted_mons.h b/src/util/lp/rooted_mons.h
index 10b75ac96..31cc208e0 100644
--- a/src/util/lp/rooted_mons.h
+++ b/src/util/lp/rooted_mons.h
@@ -54,6 +54,12 @@ struct rooted_mon {
     }
 };
 
+struct hash_svector {
+    size_t operator()(const unsigned_vector & v) const {
+        return svector_hash<unsigned_hash>()(v);
+    }
+};
+
 
 struct rooted_mon_table {
     std::unordered_map<svector<lpvar>,
diff --git a/src/util/lp/var_eqs.cpp b/src/util/lp/var_eqs.cpp
new file mode 100644
index 000000000..bd120cf17
--- /dev/null
+++ b/src/util/lp/var_eqs.cpp
@@ -0,0 +1,114 @@
+/*++
+  Copyright (c) 2017 Microsoft Corporation
+
+  Module Name:
+
+  <name>
+
+  Abstract:
+
+  <abstract>
+
+  Author:
+  Nikolaj Bjorner (nbjorner)
+  Lev Nachmanson (levnach)
+
+  Revision History:
+
+
+  --*/
+
+#include "util/lp/var_eqs.h"
+
+
+namespace nla {
+    
+
+    var_eqs::var_eqs(): m_uf(m_ufctx) {}
+    
+    void var_eqs::push() {
+        m_trail_lim.push_back(m_trail.size());
+        m_ufctx.get_trail_stack().push_scope();
+    }
+
+    void var_eqs::pop(unsigned n) {
+        unsigned old_sz = m_trail_lim[m_trail_lim.size() - n];
+        for (unsigned i = m_trail.size(); i-- > old_sz; ) {
+            auto const& sv = m_trail[i];
+            m_eqs[sv.first.index()].pop_back();
+            m_eqs[sv.second.index()].pop_back();
+            m_eqs[(~sv.first).index()].pop_back();
+            m_eqs[(~sv.second).index()].pop_back();
+        }
+        m_trail_lim.shrink(n);
+        m_trail.shrink(old_sz);
+        m_ufctx.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()) + 1;
+        m_eqs.reserve(max_i);
+        while (m_uf.get_num_vars() <= max_i) m_uf.mk_var();
+        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);
+    }
+
+    void var_eqs::explain(signed_var v1, signed_var v2, lp::explanation& e) const {
+    SASSERT(find(v1) == find(v2));
+        unsigned_vector ret;
+        if (v1 == v2) {
+            return;
+        }
+        m_todo.push_back(dfs_frame(v1, 0));
+        m_dfs_trail.reset();
+        m_marked.reserve(m_eqs.size(), false);
+        SASSERT(m_marked_trail.empty());
+        while (true) {
+            SASSERT(!m_todo.empty());
+            dfs_frame& f = m_todo.back();
+            signed_var v = f.m_var;
+            if (v == v2) {
+                break;
+            }
+            auto const& next = m_eqs[v.index()];
+            unsigned sz = next.size();
+            bool seen_all = true;
+            for (unsigned i = f.m_index; !seen_all && i < sz; ++i) {
+                justified_var const& jv = next[i];
+                if (!m_marked[jv.m_var.index()]) {
+                    seen_all = false;
+                    f.m_index = i + 1;
+                    m_todo.push_back(dfs_frame(jv.m_var, 0));
+                    m_dfs_trail.push_back(jv.m_j);
+                    m_marked[jv.m_var.index()] = true;
+                }
+            }
+            if (seen_all) {
+                m_todo.pop_back();
+                m_dfs_trail.pop_back();
+            }
+        }
+        
+        for (eq_justification const& j : m_dfs_trail) {
+            j.explain(e);
+        }
+        m_dfs_trail.reset();
+        for (unsigned idx : m_marked_trail) {
+            m_marked[idx] = false;
+        }
+        m_marked_trail.reset();
+    }
+
+}
diff --git a/src/util/lp/var_eqs.h b/src/util/lp/var_eqs.h
new file mode 100644
index 000000000..aff7ec9c6
--- /dev/null
+++ b/src/util/lp/var_eqs.h
@@ -0,0 +1,182 @@
+/*++
+  Copyright (c) 2017 Microsoft Corporation
+
+  Module Name:
+
+  <name>
+
+  Abstract:
+
+  <abstract>
+
+  Author:
+  Nikolaj Bjorner (nbjorner)
+  Lev Nachmanson (levnach)
+
+  Revision History:
+
+
+  --*/
+
+#pragma once
+#include "util/union_find.h"
+#include "util/lp/lp_types.h"
+#include "util/rational.h"
+#include "util/lp/explanation.h"
+namespace nla {
+    
+typedef lp::var_index lpvar;
+typedef lp::constraint_index lpcindex;
+
+    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
+        explicit signed_var(unsigned sv): 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;
+        }
+    };
+
+
+    class eq_justification {
+        svector<lpcindex> m_cs;
+    public:
+        eq_justification(std::initializer_list<lpcindex> cs) {
+            for (lpcindex c: cs)
+                m_cs.push_back(c);
+        }
+        void explain(lp::explanation& e) const {
+            for (lpcindex c : m_cs)
+                e.add(c);
+        }
+    };
+
+    class var_eqs {
+        struct justified_var {
+            signed_var       m_var;
+            eq_justification m_j;
+            justified_var(signed_var v, eq_justification const& j): m_var(v), m_j(j) {}
+        };
+        typedef svector<justified_var> justified_vars;
+
+        struct dfs_frame {
+            signed_var m_var;
+            unsigned   m_index;
+            dfs_frame(signed_var v, unsigned i): m_var(v), m_index(i) {}
+        };
+        typedef std::pair<signed_var, signed_var> signed_var_pair;
+
+        union_find_default_ctx  m_ufctx;
+        union_find<>            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.
+
+        mutable svector<dfs_frame>      m_todo;
+        mutable svector<bool>           m_marked;
+        mutable unsigned_vector         m_marked_trail;
+        mutable svector<eq_justification> m_dfs_trail;
+        
+    public:
+        var_eqs();
+
+        /**
+           \brief push a scope
+        */
+        void push();
+
+        /**
+           \brief pop n scopes
+         */
+        void pop(unsigned n);
+
+        /**
+           \brief merge equivalence classes for v1, v2 with justification j
+        */
+        void merge(signed_var v1, signed_var v2, eq_justification const& j);
+        void merge_plus(lpvar v1, lpvar v2, eq_justification const& j)  { merge(signed_var(v1, false), signed_var(v2, false), j); }
+        void merge_minus(lpvar v1, lpvar v2, eq_justification const& j) { merge(signed_var(v1, false), signed_var(v2, true),  j); }
+
+        /**
+           \brief find equivalence class representative for v
+         */
+        signed_var find(signed_var v) const;
+        inline lpvar find(lpvar j, rational& sign) const {
+            signed_var sv = find(signed_var(j, false));
+            sign = sv.sign()? rational(-1) : rational(1);
+            return sv.var();            
+        }
+        
+        inline rational find_sign(lpvar j) const {
+            signed_var sv = find(signed_var(j, false));        
+            return sv.sign()? rational(-1) : rational(1);
+        }
+
+        inline lpvar find(lpvar j) const {
+            signed_var sv = find(signed_var(j, false));        
+            return sv.var();
+        }
+
+        inline bool is_root(lpvar j) const {
+            signed_var sv = find(signed_var(j, false));
+            return sv.var() == j;            
+        }
+        bool vars_are_equiv(lpvar j, lpvar k) const {
+            signed_var sj = find(signed_var(j, false));
+            signed_var sk = find(signed_var(k, false));
+            return sj.var() == sk.var();
+        }
+        /**
+           \brief Returns eq_justifications for
+           \pre find(v1) == find(v2)
+        */
+        void explain(signed_var v1, signed_var v2, lp::explanation& e) const;
+        inline 
+        void explain(lpvar v1, lpvar v2, lp::explanation & e) const {
+            return explain(signed_var(v1), signed_var(v2), e);
+        }
+
+        inline void explain(lpvar j, lp::explanation& e) const {
+            signed_var s(j, false);
+            return explain(find(s), s, e); 
+        }
+
+        class iterator {
+            var_eqs& m_ve;        // context.
+            unsigned m_idx;       // index into a signed variable, same as union-find index
+            bool     m_touched;   // toggle between initial and final state
+        public:
+            iterator(var_eqs& ve, unsigned idx, bool t) : m_ve(ve), m_idx(idx), m_touched(t) {}
+            signed_var operator*() const { return signed_var(m_idx); }
+            iterator& operator++() { m_idx = m_ve.m_uf.next(m_idx); m_touched = true; return *this; }
+            bool operator==(iterator const& other) const { return m_idx == other.m_idx && m_touched == other.m_touched; }
+            bool operator!=(iterator const& other) const { return m_idx != other.m_idx || m_touched != other.m_touched; }
+        };
+
+        class eq_class {
+            var_eqs& m_ve;
+            signed_var m_v;
+        public:
+            eq_class(var_eqs& ve, signed_var v) : m_ve(ve), m_v(v) {}
+            iterator begin() { return iterator(m_ve, m_v.index(), false); }
+            iterator end() { return iterator(m_ve, m_v.index(), true); }
+        };
+
+        eq_class equiv_class(signed_var v) { return eq_class(*this, v); }
+
+        eq_class equiv_class(lpvar v) { return equiv_class(signed_var(v, false)); }
+    }; 
+}