diff --git a/src/math/lp/gomory.cpp b/src/math/lp/gomory.cpp
index a21279714..20d58bcc4 100644
--- a/src/math/lp/gomory.cpp
+++ b/src/math/lp/gomory.cpp
@@ -68,7 +68,7 @@ class gomory::imp {
             m_k.addmul(new_a, upper_bound(j).x);
             m_ex->push_justification(column_upper_bound_constraint(j));
         }
-        m_t.add_coeff_var(new_a, j);
+        m_t.add_monomial(new_a, j);
         m_lcm_den = lcm(m_lcm_den, denominator(new_a));
         TRACE("gomory_cut_detail", tout << "new_a = " << new_a << ", k = " << m_k << ", lcm_den = " << m_lcm_den << "\n";);
     }
@@ -99,7 +99,7 @@ class gomory::imp {
             m_ex->push_justification(column_upper_bound_constraint(j));
         }
         TRACE("gomory_cut_detail_real", tout << a << "*v" << j << " k: " << m_k << "\n";);
-        m_t.add_coeff_var(new_a, j);
+        m_t.add_monomial(new_a, j);
     }
 
     lia_move report_conflict_from_gomory_cut() {
@@ -124,12 +124,12 @@ class gomory::imp {
                 if (!m_k.is_int())
                     m_k = ceil(m_k);
                 // switch size
-                m_t.add_coeff_var(- mpq(1), v);
+                m_t.add_monomial(- mpq(1), v);
                 m_k.neg();
             } else {
                 if (!m_k.is_int())
                     m_k = floor(m_k);
-                m_t.add_coeff_var(mpq(1), v);
+                m_t.add_monomial(mpq(1), v);
             }
         } else {
             m_lcm_den = lcm(m_lcm_den, denominator(m_k));
@@ -145,7 +145,7 @@ class gomory::imp {
             }
             // negate everything to return -pol <= -m_k
             for (const auto & pi: pol)
-                m_t.add_coeff_var(-pi.first, pi.second);
+                m_t.add_monomial(-pi.first, pi.second);
             m_k.neg();
         }
         TRACE("gomory_cut_detail", tout << "k = " << m_k << std::endl;);
diff --git a/src/math/lp/hnf_cutter.h b/src/math/lp/hnf_cutter.h
index 65c545613..c155f3336 100644
--- a/src/math/lp/hnf_cutter.h
+++ b/src/math/lp/hnf_cutter.h
@@ -165,7 +165,7 @@ public:
     void fill_term(const vector<mpq> & row, lar_term& t) {
         for (unsigned j = 0; j < row.size(); j++) {
             if (!is_zero(row[j]))
-                t.add_coeff_var(row[j], m_var_register.local_to_external(j));
+                t.add_monomial(row[j], m_var_register.local_to_external(j));
         }
     }
 #ifdef Z3DEBUG
diff --git a/src/math/lp/horner.cpp b/src/math/lp/horner.cpp
index 2d5531f18..fa265fca1 100644
--- a/src/math/lp/horner.cpp
+++ b/src/math/lp/horner.cpp
@@ -39,7 +39,7 @@ bool horner::row_has_monomial_to_refine(const T& row) const {
 // Returns true if the row has at least two monomials sharing a variable
 template <typename T>
 bool horner::row_is_interesting(const T& row) const {
-    TRACE("nla_solver_details", m_core->print_term(row, tout););
+    TRACE("nla_solver_details", m_core->print_row(row, tout););
     if (row.size() > m_core->m_nla_settings.horner_row_length_limit()) {
         TRACE("nla_solver_details", tout << "disregard\n";);
         return false;
diff --git a/src/math/lp/int_solver.cpp b/src/math/lp/int_solver.cpp
index e17f84a9c..774648b99 100644
--- a/src/math/lp/int_solver.cpp
+++ b/src/math/lp/int_solver.cpp
@@ -1000,7 +1000,7 @@ lia_move int_solver::create_branch_on_column(int j) {
     TRACE("check_main_int", tout << "branching" << std::endl;);
     lp_assert(m_t.is_empty());
     lp_assert(j != -1);
-    m_t.add_coeff_var(mpq(1), m_lar_solver->adjust_column_index_to_term_index(j));
+    m_t.add_monomial(mpq(1), m_lar_solver->adjust_column_index_to_term_index(j));
     if (is_free(j)) {
         m_upper = true;
         m_k = mpq(0);
diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp
index bc8a886ca..b8640d7af 100644
--- a/src/math/lp/lar_solver.cpp
+++ b/src/math/lp/lar_solver.cpp
@@ -603,7 +603,7 @@ lar_term lar_solver::get_term_to_maximize(unsigned j_or_term) const {
     }
     if (j_or_term < m_mpq_lar_core_solver.m_r_x.size()) {
         lar_term r;
-        r.add_coeff_var(one_of_type<mpq>(), j_or_term);
+        r.add_monomial(one_of_type<mpq>(), j_or_term);
         return r;
     } 
     return lar_term(); // return an empty term
diff --git a/src/math/lp/lar_solver.h b/src/math/lp/lar_solver.h
index 9b01a110c..4ca0fc2d2 100644
--- a/src/math/lp/lar_solver.h
+++ b/src/math/lp/lar_solver.h
@@ -56,9 +56,9 @@ class lar_solver : public column_namer {
             using std::string;
             size_t seed = 0;
             int i = 0;
-            for (const auto& p : t.coeffs()) {
-                hash_combine(seed, p.first);
-                hash_combine(seed, p.second);
+            for (const auto p : t) {
+                hash_combine(seed, p.var());
+                hash_combine(seed, p.coeff());
                 if (i++ > 10)
                     break;
             }
@@ -69,7 +69,7 @@ class lar_solver : public column_namer {
     struct term_comparer {
         bool operator()(const lar_term &a, const lar_term& b) const
         {
-            return a.coeffs() == b.coeffs();            
+            return a == b;            
         }
     };
 
diff --git a/src/math/lp/lar_term.h b/src/math/lp/lar_term.h
index 2330a647c..6bc8ba8b0 100644
--- a/src/math/lp/lar_term.h
+++ b/src/math/lp/lar_term.h
@@ -19,7 +19,7 @@
   --*/
 #pragma once
 #include "math/lp/indexed_vector.h"
-#include <map>
+#include "util/map.h"
 
 namespace lp {
 class lar_term {
@@ -70,7 +70,7 @@ public:
     vector<std::pair<mpq, lpvar>> coeffs_as_vector() const {
         vector<std::pair<mpq, lpvar>> ret;
         for (const auto & p :  m_coeffs) {
-            ret.push_back(std::make_pair(p.second, p.first));
+            ret.push_back(std::make_pair(p.m_value, p.m_key));
         }
         return ret;
     }
@@ -95,34 +95,34 @@ public:
         m_coeffs.insert(k, b);
     }
     
-    bool contains(lpvar j) const {
-        return m_coeffs.find(j) != m_coeffs.end();
+    bool contains(unsigned j) const {
+        return m_coeffs.contains(j);
     }
 
     void negate() {
         for (auto & t : m_coeffs)
-            t.second.neg();
+            t.m_value.neg();
     }
 
     template <typename T>
     T apply(const vector<T>& x) const {
         T ret(0);
         for (const auto & t : m_coeffs) {
-            ret += t.second * x[t.first];
+            ret += t.m_value * x[t.m_key];
         }
         return ret;
     }
    
     void clear() {
-        m_coeffs.clear();
+        m_coeffs.reset();
     }
 
     struct ival {
-        lpvar m_var;
+        unsigned m_var;
         const mpq & m_coeff;
-        ival(lpvar var, const mpq & val) : m_var(var), m_coeff(val) {
+        ival(unsigned var, const mpq & val) : m_var(var), m_coeff(val) {
         }
-        lpvar var() const { return m_var;}
+        unsigned var() const { return m_var;}
         const mpq & coeff() const { return m_coeff; }
     };
     
@@ -138,19 +138,20 @@ public:
         typedef std::forward_iterator_tag iterator_category;
 
         reference operator*() const {
-            return ival(m_it->first, m_it->second);
+            return ival(m_it->m_key, m_it->m_value);
         }
         
         self_type operator++() {  self_type i = *this; m_it++; return i;  }
         self_type operator++(int) { m_it++; return *this; }
 
-        const_iterator(std::map<lpvar, mpq>::const_iterator it) : m_it(it) {}
+        const_iterator(u_map<mpq>::iterator it) : m_it(it) {}
         bool operator==(const self_type &other) const {
             return m_it == other.m_it;
         }
         bool operator!=(const self_type &other) const { return !(*this == other); }
     };
 
+    
     bool is_normalized() const {
         lpvar min_var = -1;
         mpq c;
@@ -169,23 +170,22 @@ public:
         return false;        
     }
 
-    // a is the coefficient by which we diveded the term to normalize it
+    // a is the coefficient by which we divided the term to normalize it
     lar_term get_normalized_by_min_var(mpq& a) const {
-        a = m_coeffs.begin()->second;
+        a = m_coeffs.begin()->m_value;
         if (a.is_one()) {
             return *this;
         }
         lar_term r;
         auto it = m_coeffs.begin();
-        r.add_var(it->first);
+        r.add_var(it->m_key);
         it++;
         for(;it != m_coeffs.end(); it++) {
-            r.add_coeff_var(it->second / a, it->first);
+            r.add_monomial(it->m_value / a, it->m_key);
         }
         return r;        
     }
-    const_iterator begin() const { return const_iterator(m_coeffs.begin());}
-    const_iterator end() const { return const_iterator(m_coeffs.end()); }
-
+    const_iterator begin() const { return m_coeffs.begin();}
+    const_iterator end() const { return m_coeffs.end(); }
 };
 }
diff --git a/src/math/lp/nla_common.cpp b/src/math/lp/nla_common.cpp
index 72c2a974e..7feb0183e 100644
--- a/src/math/lp/nla_common.cpp
+++ b/src/math/lp/nla_common.cpp
@@ -169,7 +169,7 @@ template <typename T> u_dependency* common::create_sum_from_row(const T& row,
                                                                 u_dependency_manager* dep_manager
                                                                 ) {
 
-    TRACE("nla_horner", tout << "row="; m_core->print_term(row, tout) << "\n";);
+    TRACE("nla_horner", tout << "row="; m_core->print_row(row, tout) << "\n";);
     u_dependency * dep = nullptr;
     SASSERT(row.size() > 1);
     sum.reset();
diff --git a/src/math/lp/nla_core.cpp b/src/math/lp/nla_core.cpp
index 058acdfee..93c2caf83 100644
--- a/src/math/lp/nla_core.cpp
+++ b/src/math/lp/nla_core.cpp
@@ -69,7 +69,7 @@ lp::lar_term core::subs_terms_to_columns(const lp::lar_term& t) const {
         lpvar j = p.var();
         if (m_lar_solver.is_term(j))
             j = m_lar_solver.map_term_index_to_column_index(j);
-        r.add_coeff_var(p.coeff(), j);
+        r.add_monomial(p.coeff(), j);
     }
     return r;
 } 
@@ -453,8 +453,8 @@ void core::mk_ineq_no_expl_check(lp::lar_term& t, llc cmp, const rational& rs) {
 
 void core::mk_ineq(const rational& a, lpvar j, const rational& b, lpvar k, llc cmp, const rational& rs) {
     lp::lar_term t;
-    t.add_coeff_var(a, j);
-    t.add_coeff_var(b, k);
+    t.add_monomial(a, j);
+    t.add_monomial(b, k);
     mk_ineq(t, cmp, rs);
 }
 
@@ -484,7 +484,7 @@ void core:: mk_ineq(lpvar j, llc cmp, const rational& rs) {
 
 void core:: mk_ineq(const rational& a, lpvar j, llc cmp, const rational& rs) {
     lp::lar_term t;        
-    t.add_coeff_var(a, j);
+    t.add_monomial(a, j);
     mk_ineq(t, cmp, rs);
 }
 
@@ -1157,7 +1157,7 @@ void core::add_abs_bound(lpvar v, llc cmp) {
 void core::add_abs_bound(lpvar v, llc cmp, rational const& bound) {
     SASSERT(!val(v).is_zero());
     lp::lar_term t;  // t = abs(v)
-    t.add_coeff_var(rrat_sign(val(v)), v);
+    t.add_monomial(rrat_sign(val(v)), v);
 
     switch (cmp) {
     case llc::GT:
@@ -1638,7 +1638,7 @@ void core::display_matrix_of_m_rows(std::ostream & out) const {
     out << m_rows.size() << " rows" <<"\n";
     out << "the matrix\n";          
     for (const auto & r : matrix.m_rows) {
-        print_term(r, out) << std::endl;
+        print_row(r, out) << std::endl;
     }
 }
 
diff --git a/src/math/lp/nla_core.h b/src/math/lp/nla_core.h
index 435d37914..633026ee2 100644
--- a/src/math/lp/nla_core.h
+++ b/src/math/lp/nla_core.h
@@ -391,7 +391,16 @@ public:
     lbool  test_check(vector<lemma>& l);
     lpvar map_to_root(lpvar) const;
     std::ostream& print_terms(std::ostream&) const;
-    std::ostream& print_term( const lp::lar_term&, std::ostream&) const;
+    std::ostream& print_term(const lp::lar_term&, std::ostream&) const;
+    template <typename T>
+    std::ostream& print_row(const T & row , std::ostream& out) const {
+        vector<std::pair<rational, lpvar>> v;
+        for (auto p : row) {
+            v.push_back(std::make_pair(p.coeff(), p.var()));
+        }
+        return lp::print_linear_combination_customized(v, [this](lpvar j) { return var_str(j); },
+        out);        
+    }
     void run_pdd_grobner();
     void find_nl_cluster();
     void prepare_rows_and_active_vars();
diff --git a/src/math/lp/nla_intervals.cpp b/src/math/lp/nla_intervals.cpp
index e2e0eae1e..f524137e6 100644
--- a/src/math/lp/nla_intervals.cpp
+++ b/src/math/lp/nla_intervals.cpp
@@ -138,13 +138,13 @@ lp::lar_term intervals::expression_to_normalized_term(const nex_sum* e, rational
 
     if (a.is_one()) {
         for (auto& p : v) {
-            t.add_coeff_var(p.first, p.second);
+            t.add_monomial(p.first, p.second);
         }
     } else {
         for (unsigned k = 0; k < v.size(); k++) {
             auto& p = v[k];
             if (k != a_index) 
-                t.add_coeff_var(p.first/a, p.second);
+                t.add_monomial(p.first/a, p.second);
             else
                 t.add_var(p.second);
         }
diff --git a/src/math/lp/nla_tangent_lemmas.cpp b/src/math/lp/nla_tangent_lemmas.cpp
index d654303cd..9606b609f 100644
--- a/src/math/lp/nla_tangent_lemmas.cpp
+++ b/src/math/lp/nla_tangent_lemmas.cpp
@@ -63,8 +63,8 @@ void tangents::generate_tang_plane(const rational & a, const rational& b, const
 #endif
 
     lp::lar_term t;
-    t.add_coeff_var(-a, jy);
-    t.add_coeff_var(-b, jx);
+    t.add_monomial(-a, jy);
+    t.add_monomial(-b, jx);
     t.add_var(j);
     c().mk_ineq(t, below? llc::GT : llc::LT, - a*b);
 }
diff --git a/src/test/lp/gomory_test.h b/src/test/lp/gomory_test.h
index b26cfab0b..fdf62d992 100644
--- a/src/test/lp/gomory_test.h
+++ b/src/test/lp/gomory_test.h
@@ -82,7 +82,7 @@ struct gomory_test {
             expl.push_justification(column_upper_bound_constraint(x_j), new_a);
         }
         TRACE("gomory_cut_detail_real", tout << a << "*v" << x_j << " k: " << k << "\n";);
-        pol.add_coeff_var(new_a, x_j);
+        pol.add_monomial(new_a, x_j);
     }
     
     void int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & k, lar_term & t, explanation& expl, mpq & lcm_den, const mpq& f_0, const mpq& one_minus_f_0) {
@@ -122,7 +122,7 @@ struct gomory_test {
             expl.push_justification(column_upper_bound_constraint(x_j), new_a);
         }
         TRACE("gomory_cut_detail", tout << "new_a: " << new_a << " k: " << k << "\n";);
-        t.add_coeff_var(new_a, x_j);
+        t.add_monomial(new_a, x_j);
         lcm_den = lcm(lcm_den, denominator(new_a));
     }
 
@@ -145,12 +145,12 @@ struct gomory_test {
                 if (!k.is_int())
                     k = ceil(k);
                 // switch size
-                t.add_coeff_var(- mpq(1), v);
+                t.add_monomial(- mpq(1), v);
                 k.neg();
             } else {
                 if (!k.is_int())
                     k = floor(k);
-                t.add_coeff_var(mpq(1), v);
+                t.add_monomial(mpq(1), v);
             }
         } else {
             TRACE("gomory_cut_detail", tout << "pol.size() > 1" << std::endl;);
@@ -177,7 +177,7 @@ struct gomory_test {
             
             // negate everything to return -pol <= -k
             for (const auto & pi: pol)
-                t.add_coeff_var(-pi.first, pi.second);
+                t.add_monomial(-pi.first, pi.second);
             k.neg();
         }
         TRACE("gomory_cut_detail", tout << "k = " << k << std::endl;);
diff --git a/src/test/lp/nla_solver_test.cpp b/src/test/lp/nla_solver_test.cpp
index e2f2f2e27..d8ab6d978 100644
--- a/src/test/lp/nla_solver_test.cpp
+++ b/src/test/lp/nla_solver_test.cpp
@@ -118,10 +118,10 @@ void test_basic_lemma_for_mon_neutral_from_factors_to_monomial_0() {
     i0.m_term.add_var(lp_ac);
     ineq i1(llc::EQ, lp::lar_term(), rational(0));
     i1.m_term.add_var(lp_bde);
-    i1.m_term.add_coeff_var(-rational(1), lp_abcde);
+    i1.m_term.add_monomial(-rational(1), lp_abcde);
     ineq i2(llc::EQ, lp::lar_term(), rational(0));
     i2.m_term.add_var(lp_abcde);
-    i2.m_term.add_coeff_var(-rational(1), lp_bde);
+    i2.m_term.add_monomial(-rational(1), lp_bde);
     bool found0 = false;
     bool found1 = false;
     bool found2 = false;
@@ -428,7 +428,7 @@ void test_horner() {
 
     lp::lar_term t;
     t.add_var(lp_c);
-    t.add_coeff_var(rational(-1), lp_b);
+    t.add_monomial(rational(-1), lp_b);
     lpvar lp_c_min_b = s.add_term(t.coeffs_as_vector(), c_min_b);
     
     reslimit l;
@@ -573,7 +573,7 @@ void test_order_lemma_params(bool var_equiv, int sign) {
     if (var_equiv) { // make k equivalent to j
         lp::lar_term t;
         t.add_var(lp_k);
-        t.add_coeff_var(-rational(1), lp_j);
+        t.add_monomial(-rational(1), lp_j);
         lpvar kj = s.add_term(t.coeffs_as_vector(), -1);
         s.add_var_bound(kj, llc::LE, rational(0));
         s.add_var_bound(kj, llc::GE, rational(0));
@@ -630,15 +630,15 @@ void test_order_lemma_params(bool var_equiv, int sign) {
 
     SASSERT(nla.get_core()->test_check(lemma) == l_false);
     // lp::lar_term t;
-    // t.add_coeff_var(lp_bde);
-    // t.add_coeff_var(lp_acd);
+    // t.add_monomial(lp_bde);
+    // t.add_monomial(lp_acd);
     // ineq q(llc::EQ, t, rational(0));
    
     nla.get_core()->print_lemma(std::cout);
     // SASSERT(q == lemma.back());
     // ineq i0(llc::EQ, lp::lar_term(), rational(0));
-    // i0.m_term.add_coeff_var(lp_bde);
-    // i0.m_term.add_coeff_var(rational(1), lp_acd);
+    // i0.m_term.add_monomial(lp_bde);
+    // i0.m_term.add_monomial(rational(1), lp_acd);
     // bool found = false;
     // for (const auto& k : lemma){
     //     if (k == i0) {
diff --git a/src/test/pdd.cpp b/src/test/pdd.cpp
index 7152096bd..53abf991e 100644
--- a/src/test/pdd.cpp
+++ b/src/test/pdd.cpp
@@ -1,90 +1,106 @@
 #include "math/dd/dd_pdd.h"
 
 namespace dd {
-    static void test1() {
-        pdd_manager m(3);
-        pdd v0 = m.mk_var(0);
-        pdd v1 = m.mk_var(1);
-        pdd v2 = m.mk_var(2);
-        std::cout << v0 << "\n";
-        std::cout << v1 << "\n";
-        std::cout << v2 << "\n";
-        pdd c1 = v0 * v1 * v2;
-        pdd c2 = v2 * v0 * v1;
-        std::cout << c1 << "\n";
-        SASSERT(c1 == c2);
+static void test1() {
+    pdd_manager m(3);
+    pdd v0 = m.mk_var(0);
+    pdd v1 = m.mk_var(1);
+    pdd v2 = m.mk_var(2);
+    std::cout << v0 << "\n";
+    std::cout << v1 << "\n";
+    std::cout << v2 << "\n";
+    pdd c1 = v0 * v1 * v2;
+    pdd c2 = v2 * v0 * v1;
+    std::cout << c1 << "\n";
+    SASSERT(c1 == c2);
 
-        c1 = v0 + v1 + v2;
-        c2 = v2 + v1 + v0;
-        std::cout << c1 << "\n";
-        SASSERT(c1 == c2);
+    c1 = v0 + v1 + v2;
+    c2 = v2 + v1 + v0;
+    std::cout << c1 << "\n";
+    SASSERT(c1 == c2);
 
-        c1 = (v0+v1) * v2;
-        c2 = (v0*v2) + (v1*v2);
-        std::cout << c1 << "\n";
-        SASSERT(c1 == c2);
-        c1 = (c1 + 3) + 1;
-        c2 = (c2 + 1) + 3;
-        std::cout << c1 << "\n";
-        SASSERT(c1 == c2);
-        c1 = v0 - v1;
-        c2 = v1 - v0;
-        std::cout << c1 << " " << c2 << "\n";
+    c1 = (v0+v1) * v2;
+    c2 = (v0*v2) + (v1*v2);
+    std::cout << c1 << "\n";
+    SASSERT(c1 == c2);
+    c1 = (c1 + 3) + 1;
+    c2 = (c2 + 1) + 3;
+    std::cout << c1 << "\n";
+    SASSERT(c1 == c2);
+    c1 = v0 - v1;
+    c2 = v1 - v0;
+    std::cout << c1 << " " << c2 << "\n";
 
-        c1 = v1*v2;
-        c2 = (v0*v2) + (v2*v2);
-        pdd c3 = m.zero();
-        VERIFY(m.try_spoly(c1, c2, c3));
-        std::cout << c1 << " " << c2 << " spoly: " << c3 << "\n";
+    c1 = v1*v2;
+    c2 = (v0*v2) + (v2*v2);
+    pdd c3 = m.zero();
+    VERIFY(m.try_spoly(c1, c2, c3));
+    std::cout << c1 << " " << c2 << " spoly: " << c3 << "\n";
 
-        c1 = v1*v2;
-        c2 = (v0*v2) + (v1*v1);
-        VERIFY(m.try_spoly(c1, c2, c3));
-        std::cout << c1 << " " << c2 << " spoly: " << c3 << "\n";
+    c1 = v1*v2;
+    c2 = (v0*v2) + (v1*v1);
+    VERIFY(m.try_spoly(c1, c2, c3));
+    std::cout << c1 << " " << c2 << " spoly: " << c3 << "\n";
 
-        c1 = (v0*v1) - (v0*v0);
-        c2 = (v0*v1*(v2 + v0)) + v2;
-        c3 = c2.reduce(c1);
-        std::cout << c1 << " " << c2 << " reduce: " << c3 << "\n";
-    }
-
-    static void test2() {
-        std::cout << "\ntest2\n";
-        // a(b^2)cd + abc + bcd + bc + cd + 3 reduce by  bc
-        pdd_manager m(4);
-        pdd a = m.mk_var(0);
-        pdd b = m.mk_var(1);
-        pdd c = m.mk_var(2);
-        pdd d = m.mk_var(3);
-        pdd e = (a * b * b * c * d) + (2*a*b*c) + (b*c*d) + (b*c) + (c*d) + 3;
-        std::cout << e << "\n";
-        pdd f = b * c;
-        pdd r_ef = m.reduce(e, f);
-        m.display(std::cout);
-        std::cout << "result of reduce " << e << " by " << f << " is " << r_ef <<  "\n";
-        pdd r_fe = m.reduce(f, e);
-        std::cout << "result of reduce " << f << " by " << e << " is " << r_fe << "\n" ;
-        VERIFY(r_fe == f);
+    c1 = (v0*v1) - (v0*v0);
+    c2 = (v0*v1*(v2 + v0)) + v2;
+    c3 = c2.reduce(c1);
+    std::cout << c1 << " " << c2 << " reduce: " << c3 << "\n";
+}
+
+static void test2() {
+    std::cout << "\ntest2\n";
+    // a(b^2)cd + abc + bcd + bc + cd + 3 reduce by  bc
+    pdd_manager m(4);
+    pdd a = m.mk_var(0);
+    pdd b = m.mk_var(1);
+    pdd c = m.mk_var(2);
+    pdd d = m.mk_var(3);
+    pdd e = (a * b * b * c * d) + (2*a*b*c) + (b*c*d) + (b*c) + (c*d) + 3;
+    std::cout << e << "\n";
+    pdd f = b * c;
+    pdd r_ef = m.reduce(e, f);
+    m.display(std::cout);
+    std::cout << "result of reduce " << e << " by " << f << " is " << r_ef <<  "\n";
+    pdd r_fe = m.reduce(f, e);
+    std::cout << "result of reduce " << f << " by " << e << " is " << r_fe << "\n" ;
+    VERIFY(r_fe == f);
+}
+
+static void test3() {
+    std::cout << "\ntest3\n";
+    pdd_manager m(4);
+    pdd a = m.mk_var(0);
+    pdd b = m.mk_var(1);
+    pdd c = m.mk_var(2);
+    pdd d = m.mk_var(3);
+        
+    pdd e = a + c;
+    for (unsigned i = 0; i < 5; i++) {
+        e = e * e;
     }
+    e = e * b;
+    std::cout << e << "\n";
+}
 
 static void test_reset() {
     std::cout << "\ntest reset\n";
     pdd_manager m(4);
-        pdd a = m.mk_var(0);
-        pdd b = m.mk_var(1);
-        pdd c = m.mk_var(2);
-        pdd d = m.mk_var(3);
-        std::cout << (a + b)*(c + d) << "\n";
+    pdd a = m.mk_var(0);
+    pdd b = m.mk_var(1);
+    pdd c = m.mk_var(2);
+    pdd d = m.mk_var(3);
+    std::cout << (a + b)*(c + d) << "\n";
 
-        unsigned_vector l2v;
-        for (unsigned i = 0; i < 4; ++i) 
-            l2v.push_back(3 - i);
-        m.reset(l2v);
-        a = m.mk_var(0);
-        b = m.mk_var(1);
-        c = m.mk_var(2);
-        d = m.mk_var(3);
-        std::cout << (a + b)*(c + d) << "\n";
+    unsigned_vector l2v;
+    for (unsigned i = 0; i < 4; ++i) 
+        l2v.push_back(3 - i);
+    m.reset(l2v);
+    a = m.mk_var(0);
+    b = m.mk_var(1);
+    c = m.mk_var(2);
+    d = m.mk_var(3);
+    std::cout << (a + b)*(c + d) << "\n";
 }
 
 static void test5() {
@@ -110,25 +126,6 @@ static void test5() {
     SASSERT(e == f);        
 }
 
-    static void test_reset() {
-        std::cout << "\ntest reset\n";
-        pdd_manager m(4);
-        pdd a = m.mk_var(0);
-        pdd b = m.mk_var(1);
-        pdd c = m.mk_var(2);
-        pdd d = m.mk_var(3);
-        std::cout << (a + b)*(c + d) << "\n";
-
-        unsigned_vector l2v;
-        for (unsigned i = 0; i < 4; ++i) 
-            l2v.push_back(3 - i);
-        m.reset(l2v);
-        a = m.mk_var(0);
-        b = m.mk_var(1);
-        c = m.mk_var(2);
-        d = m.mk_var(3);
-        std::cout << (a + b)*(c + d) << "\n";
-    }
 
     void test_iterator() {
         std::cout << "test iterator\n";
@@ -150,6 +147,7 @@ void tst_pdd() {
     dd::test1();
     dd::test2();
     dd::test3();
+    dd::test5();
     dd::test_reset();
     dd::test_iterator();
 }