diff --git a/src/math/dd/dd_pdd.cpp b/src/math/dd/dd_pdd.cpp
index 146b8948b..01cca363c 100644
--- a/src/math/dd/dd_pdd.cpp
+++ b/src/math/dd/dd_pdd.cpp
@@ -25,24 +25,24 @@ namespace dd {
 
     pdd_manager::pdd_manager(unsigned num_vars) {
         m_spare_entry = nullptr;
-        m_max_num_pdd_nodes = 1 << 24; // up to 16M nodes
+        m_max_num_nodes = 1 << 24; // up to 16M nodes
         m_mark_level = 0;
-        alloc_free_nodes(1024 + num_vars);
         m_disable_gc = false;
         m_is_new_node = false;
         m_mod2_semantics = false;
 
         // add dummy nodes for operations, and 0, 1 pdds.
-        const_info info;
-        init_value(info, rational::zero()); // becomes pdd_zero
-        init_value(info, rational::one());  // becomes pdd_one
-        m_nodes[0].m_refcount = max_rc;
-        m_nodes[1].m_refcount = max_rc;
-        for (unsigned i = 2; i <= pdd_no_op + 2; ++i) {
-            m_nodes.push_back(pdd_node(0,0,0));
-            m_nodes.back().m_refcount = max_rc;
-            m_nodes.back().m_index = m_nodes.size()-1;
+        for (unsigned i = 0; i < pdd_no_op; ++i) {
+            m_nodes.push_back(node());
+            m_nodes[i].m_refcount = max_rc;
+            m_nodes[i].m_index = i;
         }
+        init_value(rational::zero(), 0);
+        init_value(rational::one(), 1);
+        SASSERT(is_val(0));
+        SASSERT(is_val(1));
+
+        alloc_free_nodes(1024 + num_vars);
         
         // add variables
         for (unsigned i = 0; i < num_vars; ++i) {
@@ -70,6 +70,8 @@ namespace dd {
     pdd pdd_manager::zero() { return pdd(zero_pdd, this); }
     pdd pdd_manager::one() { return pdd(one_pdd, this); }
     
+    pdd pdd_manager::mk_or(pdd const& p, pdd const& q) { return p*q + p + q; }
+
     pdd_manager::PDD pdd_manager::apply(PDD arg1, PDD arg2, pdd_op op) {
         bool first = true;
         SASSERT(well_formed());
@@ -284,11 +286,16 @@ namespace dd {
     // q = lt(a)/lt(b), return a - b*q
     pdd_manager::PDD pdd_manager::reduce_on_match(PDD a, PDD b) {
         SASSERT(level(a) == level(b) && !is_val(a) && !is_val(b));
-        while (lm_divides(b, a)) {
-            PDD q = lt_quotient(b, a);
-            PDD r = apply(q, b, pdd_mul_op);
-            a = apply(a, r, pdd_add_op);
+        push(a);
+        while (lm_divides(b, a)) {           
+            push(lt_quotient(b, a));
+            push(apply_rec(read(1), b, pdd_mul_op));
+            push(apply_rec(a, read(1), pdd_add_op));
+            a = read(1);
+            pop(4);
+            push(a);
         }
+        pop(1);
         return a;
     }
 
@@ -315,13 +322,16 @@ namespace dd {
         if (is_val(p)) {
             if (is_val(q)) {
                 SASSERT(!val(p).is_zero());
-                return imk_val(-val(q)/val(p));
-            }     
-        }   
+                return imk_val(-val(q) / val(p));
+            }
+        }
         else if (level(p) == level(q)) {
             return lt_quotient(hi(p), hi(q));
         }
-        return apply(m_var2pdd[var(q)], lt_quotient(p, hi(q)), pdd_mul_op);
+        push(lt_quotient(p, hi(q)));
+        PDD r = apply_rec(m_var2pdd[var(q)], read(1), pdd_mul_op);
+        pop(1);
+        return r;                   
     }
 
     //
@@ -335,12 +345,10 @@ namespace dd {
 
     pdd pdd_manager::spoly(pdd const& a, pdd const& b, unsigned_vector const& p, unsigned_vector const& q, rational const& pc, rational const& qc) { 
         pdd r1 = mk_val(qc);  
-        for (unsigned i = q.size(); i-- > 0; ) r1 = mul(mk_var(q[i]), r1);
-        r1 = mul(a, r1);
+        for (unsigned i = q.size(); i-- > 0; ) r1 *= mk_var(q[i]);
         pdd r2 = mk_val(-pc);
-        for (unsigned i = p.size(); i-- > 0; ) r2 = mul(mk_var(p[i]), r2);
-        r2 = mul(b, r2);
-        return add(r1, r2);
+        for (unsigned i = p.size(); i-- > 0; ) r2 *= mk_var(p[i]);
+        return (r1*a) + (r2*b);
     }
 
     bool pdd_manager::common_factors(pdd const& a, pdd const& b, unsigned_vector& p, unsigned_vector& q, rational& pc, rational& qc) { 
@@ -502,22 +510,32 @@ namespace dd {
         }
 
         m_freeze_value = r;
-        pdd_node n(vi);
+        node n(vi);
         info.m_value_index = vi;        
         info.m_node_index = insert_node(n);
         m_mpq_table.insert(r, info);
     }
 
+    void pdd_manager::init_value(rational const& v, unsigned node_index) {
+        const_info info;
+        m_nodes[node_index].m_hi = 0;
+        m_nodes[node_index].m_lo = node_index;
+        info.m_value_index = m_values.size();
+        info.m_node_index = node_index;
+        m_mpq_table.insert(v, info);
+        m_values.push_back(v);
+    }
+
     pdd_manager::PDD pdd_manager::make_node(unsigned lvl, PDD l, PDD h) {
         m_is_new_node = false;
         if (is_zero(h)) return l;
         SASSERT(is_val(l) || level(l) < lvl);
         SASSERT(is_val(h) || level(h) <= lvl);
-        pdd_node n(lvl, l, h);
+        node n(lvl, l, h);
         return insert_node(n);
     }
 
-    pdd_manager::PDD pdd_manager::insert_node(pdd_node const& n) {
+    pdd_manager::PDD pdd_manager::insert_node(node const& n) {
         node_table::entry* e = m_node_table.insert_if_not_there2(n);
         if (e->get_data().m_index != 0) {
             unsigned result = e->get_data().m_index;
@@ -528,12 +546,11 @@ namespace dd {
         bool do_gc = m_free_nodes.empty();
         if (do_gc && !m_disable_gc) {
             gc();
-            SASSERT(n.m_hi == 0 || (!m_free_nodes.contains(n.m_hi) && !m_free_nodes.contains(n.m_lo)));            
             e = m_node_table.insert_if_not_there2(n);
             e->get_data().m_refcount = 0;      
         }
         if (do_gc) {
-            if (m_nodes.size() > m_max_num_pdd_nodes) {
+            if (m_nodes.size() > m_max_num_nodes) {
                 throw mem_out();
             }
             alloc_free_nodes(m_nodes.size()/2);
@@ -686,22 +703,28 @@ namespace dd {
     void pdd_manager::alloc_free_nodes(unsigned n) {
         for (unsigned i = 0; i < n; ++i) {
             m_free_nodes.push_back(m_nodes.size());
-            m_nodes.push_back(pdd_node());
+            m_nodes.push_back(node());
             m_nodes.back().m_index = m_nodes.size() - 1;
         }        
+        std::sort(m_free_nodes.begin(), m_free_nodes.end());
         m_free_nodes.reverse();
     }
 
-    void pdd_manager::gc() {
-        m_free_nodes.reset();
-        SASSERT(well_formed());
-        IF_VERBOSE(13, verbose_stream() << "(pdd :gc " << m_nodes.size() << ")\n";);
+    bool pdd_manager::is_reachable(PDD p) {
         svector<bool> reachable(m_nodes.size(), false);
+        compute_reachable(reachable);
+        return reachable[p];
+    }
+
+    void pdd_manager::compute_reachable(svector<bool>& reachable) {
         for (unsigned i = m_pdd_stack.size(); i-- > 0; ) {
             reachable[m_pdd_stack[i]] = true;
             m_todo.push_back(m_pdd_stack[i]);
         }
-        for (unsigned i = m_nodes.size(); i-- > 2; ) {
+        for (unsigned i = pdd_no_op; i-- > 0; ) {
+            reachable[i] = true;
+        }
+        for (unsigned i = m_nodes.size(); i-- > pdd_no_op; ) {
             if (m_nodes[i].m_refcount > 0) {
                 reachable[i] = true;
                 m_todo.push_back(i);
@@ -711,7 +734,9 @@ namespace dd {
             PDD p = m_todo.back();
             m_todo.pop_back();
             SASSERT(reachable[p]);
-            if (is_val(p)) continue;
+            if (is_val(p)) { 
+                continue;
+            }
             if (!reachable[lo(p)]) {
                 reachable[lo(p)] = true;
                 m_todo.push_back(lo(p));
@@ -721,7 +746,15 @@ namespace dd {
                 m_todo.push_back(hi(p));
             }
         }
-        for (unsigned i = m_nodes.size(); i-- > 2; ) {
+    }
+
+    void pdd_manager::gc() {
+        m_free_nodes.reset();
+        SASSERT(well_formed());
+        IF_VERBOSE(13, verbose_stream() << "(pdd :gc " << m_nodes.size() << ")\n";);
+        svector<bool> reachable(m_nodes.size(), false);
+        compute_reachable(reachable);
+        for (unsigned i = m_nodes.size(); i-- > pdd_no_op; ) {
             if (!reachable[i]) {
                 if (is_val(i)) {
                     if (m_freeze_value == val(i)) continue;
@@ -835,7 +868,7 @@ namespace dd {
                 return false;
             }
         }
-        for (pdd_node const& n : m_nodes) {
+        for (node const& n : m_nodes) {
             if (!well_formed(n)) {
                 IF_VERBOSE(0, display(verbose_stream() << n.m_index << " lo " << n.m_lo << " hi " << n.m_hi << "\n"););
                 UNREACHABLE();
@@ -845,7 +878,7 @@ namespace dd {
         return ok;
     }
 
-    bool pdd_manager::well_formed(pdd_node const& n) {
+    bool pdd_manager::well_formed(node const& n) {
         PDD lo = n.m_lo;
         PDD hi = n.m_hi;        
         if (n.is_internal() || hi == 0) return true;
@@ -856,7 +889,7 @@ namespace dd {
 
     std::ostream& pdd_manager::display(std::ostream& out) {
         for (unsigned i = 0; i < m_nodes.size(); ++i) {
-            pdd_node const& n = m_nodes[i];
+            node const& n = m_nodes[i];
             if (i != 0 && n.is_internal()) {
                 continue;
             }
@@ -870,7 +903,14 @@ namespace dd {
         return out;
     }
 
-    pdd& pdd::operator=(pdd const& other) { unsigned r1 = root; root = other.root; m->inc_ref(root); m->dec_ref(r1); return *this; }
+    pdd& pdd::operator=(pdd const& other) { 
+        unsigned r1 = root; 
+        root = other.root; 
+        m.inc_ref(root); 
+        m.dec_ref(r1); 
+        return *this; 
+    }
+
     std::ostream& operator<<(std::ostream& out, pdd const& b) { return b.display(out); }
 
 }
diff --git a/src/math/dd/dd_pdd.h b/src/math/dd/dd_pdd.h
index 4cb4b0cde..c9d41bfb4 100644
--- a/src/math/dd/dd_pdd.h
+++ b/src/math/dd/dd_pdd.h
@@ -58,15 +58,15 @@ namespace dd {
             pdd_no_op = 6
         };
 
-        struct pdd_node {
-            pdd_node(unsigned level, PDD lo, PDD hi):
+        struct node {
+            node(unsigned level, PDD lo, PDD hi):
                 m_refcount(0),
                 m_level(level),
                 m_lo(lo),
                 m_hi(hi),
                 m_index(0)
             {}
-            pdd_node(unsigned value):
+            node(unsigned value):
                 m_refcount(0),
                 m_level(0),
                 m_lo(value),
@@ -74,28 +74,29 @@ namespace dd {
                 m_index(0)
             {}
 
-            pdd_node(): m_refcount(0), m_level(0), m_lo(0), m_hi(0), m_index(0) {}
+            node(): m_refcount(0), m_level(0), m_lo(0), m_hi(0), m_index(0) {}
             unsigned m_refcount : 10;
             unsigned m_level : 22;
             PDD      m_lo;
             PDD      m_hi;
             unsigned m_index;
             unsigned hash() const { return mk_mix(m_level, m_lo, m_hi); } 
-            bool is_internal() const { return m_lo == 0 && m_hi == 0; }
+            bool is_val() const { return m_hi == 0 && (m_lo != 0 || m_index == 0); }
+            bool is_internal() const { return m_hi == 0 && m_lo == 0 && m_index != 0; }
             void set_internal() { m_lo = 0; m_hi = 0; }
         };
 
         struct hash_node {
-            unsigned operator()(pdd_node const& n) const { return n.hash(); }
+            unsigned operator()(node const& n) const { return n.hash(); }
         };
 
         struct eq_node {
-            bool operator()(pdd_node const& a, pdd_node const& b) const {
+            bool operator()(node const& a, node const& b) const {
                 return a.m_lo == b.m_lo && a.m_hi == b.m_hi && a.m_level == b.m_level;
             }
         };
         
-        typedef hashtable<pdd_node, hash_node, eq_node> node_table;
+        typedef hashtable<node, hash_node, eq_node> node_table;
 
         struct const_info {
             unsigned m_value_index;
@@ -131,7 +132,7 @@ namespace dd {
 
         typedef ptr_hashtable<op_entry, hash_entry, eq_entry> op_table;
 
-        svector<pdd_node>          m_nodes;
+        svector<node>          m_nodes;
         vector<rational>           m_values;
         op_table                   m_op_cache;
         node_table                 m_node_table;
@@ -149,14 +150,14 @@ namespace dd {
         mutable svector<double>    m_tree_size;
         bool                       m_disable_gc;
         bool                       m_is_new_node;
-        unsigned                   m_max_num_pdd_nodes;
+        unsigned                   m_max_num_nodes;
         bool                       m_mod2_semantics;
         unsigned_vector            m_free_vars;
         unsigned_vector            m_free_values;
         rational                   m_freeze_value;
 
         PDD make_node(unsigned level, PDD l, PDD r);
-        PDD insert_node(pdd_node const& n);
+        PDD insert_node(node const& n);
         bool is_new_node() const { return m_is_new_node; }
 
         PDD apply(PDD arg1, PDD arg2, pdd_op op);
@@ -169,6 +170,7 @@ namespace dd {
 
         PDD imk_val(rational const& r);       
         void init_value(const_info& info, rational const& r);
+        void init_value(rational const& v, unsigned r);
 
         void push(PDD b);
         void pop(unsigned num_scopes);
@@ -185,25 +187,28 @@ namespace dd {
 
         static const unsigned max_rc = (1 << 10) - 1;
 
-        inline bool is_zero(PDD b) const { return b == zero_pdd; } 
-        inline bool is_one(PDD b) const { return b == one_pdd; } 
-        inline bool is_val(PDD b) const { return hi(b) == 0; }
-        inline unsigned level(PDD b) const { return m_nodes[b].m_level; }
-        inline unsigned var(PDD b) const { return m_level2var[level(b)]; }
-        inline PDD lo(PDD b) const { return m_nodes[b].m_lo; }
-        inline PDD hi(PDD b) const { return m_nodes[b].m_hi; }
-        inline rational const& val(PDD b) const { SASSERT(is_val(b)); return m_values[lo(b)]; }
-        inline void inc_ref(PDD b) { if (m_nodes[b].m_refcount != max_rc) m_nodes[b].m_refcount++; SASSERT(!m_free_nodes.contains(b)); }
-        inline void dec_ref(PDD b) { if (m_nodes[b].m_refcount != max_rc) m_nodes[b].m_refcount--; SASSERT(!m_free_nodes.contains(b)); }
+        inline bool is_zero(PDD p) const { return p == zero_pdd; } 
+        inline bool is_one(PDD p) const { return p == one_pdd; } 
+        inline bool is_val(PDD p) const { return m_nodes[p].is_val(); }
+        inline bool is_internal(PDD p) const { return m_nodes[p].is_internal(); }
+        inline unsigned level(PDD p) const { return m_nodes[p].m_level; }
+        inline unsigned var(PDD p) const { return m_level2var[level(p)]; }
+        inline PDD lo(PDD p) const { return m_nodes[p].m_lo; }
+        inline PDD hi(PDD p) const { return m_nodes[p].m_hi; }
+        inline rational const& val(PDD p) const { SASSERT(is_val(p)); return m_values[lo(p)]; }
+        inline void inc_ref(PDD p) { if (m_nodes[p].m_refcount != max_rc) m_nodes[p].m_refcount++; SASSERT(!m_free_nodes.contains(p)); }
+        inline void dec_ref(PDD p) { if (m_nodes[p].m_refcount != max_rc) m_nodes[p].m_refcount--; SASSERT(!m_free_nodes.contains(p)); }
         inline PDD level2pdd(unsigned l) const { return m_var2pdd[m_level2var[l]]; }
 
-        unsigned dag_size(pdd const& b);
+        unsigned dag_size(pdd const& p);
         unsigned degree(pdd const& p);
 
+        bool is_reachable(PDD p);
+        void compute_reachable(svector<bool>& reachable);
         void try_gc();
         void reserve_var(unsigned v);
         bool well_formed();
-        bool well_formed(pdd_node const& n);
+        bool well_formed(node const& n);
 
         unsigned_vector m_p, m_q;
         rational m_pc, m_qc;
@@ -227,7 +232,7 @@ namespace dd {
         ~pdd_manager();
 
         void set_mod2_semantics() { m_mod2_semantics = true; }
-        void set_max_num_nodes(unsigned n) { m_max_num_pdd_nodes = n; }
+        void set_max_num_nodes(unsigned n) { m_max_num_nodes = n; }
         void set_level2var(unsigned_vector const& level2var);  
         unsigned_vector const& get_level2var() const { return m_level2var; }
 
@@ -241,6 +246,7 @@ namespace dd {
         pdd sub(pdd const& a, pdd const& b);
         pdd mul(pdd const& a, pdd const& b);
         pdd mul(rational const& c, pdd const& b);
+        pdd mk_or(pdd const& p, pdd const& q);
         pdd reduce(pdd const& a, pdd const& b);
 
         bool is_linear(PDD p);
@@ -264,38 +270,40 @@ namespace dd {
     class pdd {
         friend class pdd_manager;
         unsigned     root;
-        pdd_manager* m;
-        pdd(unsigned root, pdd_manager* m): root(root), m(m) { m->inc_ref(root); }
+        pdd_manager& m;
+        pdd(unsigned root, pdd_manager& m): root(root), m(m) { m.inc_ref(root); }
+        pdd(unsigned root, pdd_manager* _m): root(root), m(*_m) { m.inc_ref(root); }
     public:
-        pdd(pdd_manager& pm): root(0), m(&pm) { SASSERT(is_zero()); }
-        pdd(pdd const& other): root(other.root), m(other.m) { m->inc_ref(root); }
+        pdd(pdd_manager& pm): root(0), m(pm) { SASSERT(is_zero()); }
+        pdd(pdd const& other): root(other.root), m(other.m) { m.inc_ref(root); }
         pdd(pdd && other): root(0), m(other.m) { std::swap(root, other.root); }
         pdd& operator=(pdd const& other);
-        ~pdd() { m->dec_ref(root); }
-        pdd lo() const { return pdd(m->lo(root), m); }
-        pdd hi() const { return pdd(m->hi(root), m); }
-        unsigned var() const { return m->var(root); }
-        rational const& val() const { SASSERT(is_val()); return m->val(root); }
-        bool is_val() const { return m->is_val(root); }
-        bool is_zero() const { return m->is_zero(root); }
-        bool is_linear() const { return m->is_linear(root); }
+        ~pdd() { m.dec_ref(root); }
+        pdd lo() const { return pdd(m.lo(root), m); }
+        pdd hi() const { return pdd(m.hi(root), m); }
+        unsigned var() const { return m.var(root); }
+        rational const& val() const { SASSERT(is_val()); return m.val(root); }
+        bool is_val() const { return m.is_val(root); }
+        bool is_zero() const { return m.is_zero(root); }
+        bool is_linear() const { return m.is_linear(root); }
 
-        pdd operator+(pdd const& other) const { return m->add(*this, other); }
-        pdd operator-(pdd const& other) const { return m->sub(*this, other); }
-        pdd operator*(pdd const& other) const { return m->mul(*this, other); }
-        pdd operator*(rational const& other) const { return m->mul(other, *this); }
-        pdd operator+(rational const& other) const { return m->add(other, *this); }
-        pdd reduce(pdd const& other) const { return m->reduce(*this, other); }
-        bool different_leading_term(pdd const& other) const { return m->different_leading_term(*this, other); }
+        pdd operator+(pdd const& other) const { return m.add(*this, other); }
+        pdd operator-(pdd const& other) const { return m.sub(*this, other); }
+        pdd operator*(pdd const& other) const { return m.mul(*this, other); }
+        pdd operator*(rational const& other) const { return m.mul(other, *this); }
+        pdd operator+(rational const& other) const { return m.add(other, *this); }
+        pdd operator|(pdd const& other) const { return m.mk_or(*this, other); }
+        pdd reduce(pdd const& other) const { return m.reduce(*this, other); }
+        bool different_leading_term(pdd const& other) const { return m.different_leading_term(*this, other); }
 
-        std::ostream& display(std::ostream& out) const { return m->display(out, *this); }
+        std::ostream& display(std::ostream& out) const { return m.display(out, *this); }
         bool operator==(pdd const& other) const { return root == other.root; }
         bool operator!=(pdd const& other) const { return root != other.root; }
-        bool operator<(pdd const& b) const { return m->lt(*this, b); }
+        bool operator<(pdd const& b) const { return m.lt(*this, b); }
 
-        unsigned dag_size() const { return m->dag_size(*this); }
-        double tree_size() const { return m->tree_size(*this); }
-        unsigned degree() const { return m->degree(*this); }
+        unsigned dag_size() const { return m.dag_size(*this); }
+        double tree_size() const { return m.tree_size(*this); }
+        unsigned degree() const { return m.degree(*this); }
     };
 
     inline pdd operator*(rational const& r, pdd const& b) { return b * r; }
@@ -307,6 +315,8 @@ namespace dd {
     inline pdd operator+(pdd const& b, int x) { return b + rational(x); }
 
     inline pdd operator-(pdd const& b, int x) { return b + (-rational(x)); }
+    inline pdd& operator*=(pdd & p, pdd const& q) { p = p * q; return p; }
+    inline pdd& operator|=(pdd & p, pdd const& q) { p = p | q; return p; }
 
 
     std::ostream& operator<<(std::ostream& out, pdd const& b);
diff --git a/src/math/grobner/pdd_grobner.cpp b/src/math/grobner/pdd_grobner.cpp
index 5327e2cd6..18dfcc335 100644
--- a/src/math/grobner/pdd_grobner.cpp
+++ b/src/math/grobner/pdd_grobner.cpp
@@ -143,7 +143,10 @@ namespace dd {
     }
 
     bool grobner::basic_step() {
-        equation* e = pick_next();
+        return basic_step(pick_next());
+    }
+
+    bool grobner::basic_step(equation* e) {
         if (!e) return false;
         scoped_detach sd(*this, e);
         equation& eq = *e;
@@ -164,13 +167,89 @@ namespace dd {
                 eq = curr;
             }
         }
-        if (eq) pop_equation(eq->idx(), m_to_simplify);
+        if (eq) pop_equation(eq, m_to_simplify);
         return eq;
     }
 
+    grobner::equation* grobner::pick_linear() {
+        equation* eq = nullptr;
+        for (auto* curr : m_to_simplify) {
+            if (!eq || curr->poly().is_linear() && is_simpler(*curr, *eq)) {
+                eq = curr;
+            }
+        }
+        if (eq) pop_equation(eq, m_to_simplify);
+        return eq;
+    }
+
+    void grobner::simplify_linear() {
+        try {
+            while (simplify_linear_step()) {
+                DEBUG_CODE(invariant(););
+            }
+        }
+        catch (pdd_manager::mem_out) {
+            // done reduce
+        }
+    }
+
+    struct grobner::compare_top_var {
+        bool operator()(equation* a, equation* b) const {
+            return a->poly().var() < b->poly().var();
+        }
+    };
+
+    bool grobner::simplify_linear_step() {
+        equation_vector linear;
+        for (equation* e : m_to_simplify) {
+            if (e->poly().is_linear()) {
+                linear.push_back(e);
+            }
+        }
+        if (linear.empty()) return false;
+        std::cout << "linear eqs: " << linear.size() << "\n";
+        vector<equation_vector> use_list;
+        for (equation * e : m_to_simplify) {
+            add_to_use(e, use_list);
+        }
+        for (equation * e : m_processed) {
+            add_to_use(e, use_list);
+        }
+        compare_top_var ctv;
+        std::stable_sort(linear.begin(), linear.end(), ctv);
+        for (equation* src : linear) {
+            equation_vector const& uses = use_list[src->poly().var()];
+            bool changed_leading_term;
+            for (equation* dst : uses) {
+                if (src == dst || !simplify_source_target(*src, *dst, changed_leading_term)) {
+                    continue;                    
+                }
+                if (dst->is_processed() && changed_leading_term) {
+                    dst->set_processed(false);
+                    pop_equation(dst, m_processed);
+                    push_equation(src, m_to_simplify);
+                }
+            }            
+        } 
+        for (equation* src : linear) {
+            pop_equation(src, m_to_simplify);
+            push_equation(src, m_processed);
+            src->set_processed(true);
+        }
+        return true;
+    }
+
+    void grobner::add_to_use(equation* e, vector<equation_vector>& use_list) {
+        unsigned_vector const& fv = m.free_vars(e->poly());
+        for (unsigned v : fv) {
+            use_list.reserve(v + 1);
+            use_list[v].push_back(e);
+        }
+    }
+
     void grobner::superpose(equation const & eq) {
         for (equation* target : m_processed) {
-            retire(target);             
+            superpose(eq, *target);
         }
     }
 
@@ -343,7 +422,7 @@ namespace dd {
                 }
             }
             if (eq) {
-                pop_equation(eq->idx(), m_to_simplify);
+                pop_equation(eq, m_to_simplify);
                 m_watch[eq->poly().var()].erase(eq);
                 return eq;
             }
@@ -390,11 +469,12 @@ namespace dd {
     }
 
     void grobner::del_equation(equation& eq) {
-        pop_equation(eq.idx(), eq.is_processed() ? m_processed : m_to_simplify);
+        pop_equation(eq, eq.is_processed() ? m_processed : m_to_simplify);
         retire(&eq);
     }
 
-    void grobner::pop_equation(unsigned idx, equation_vector& v) {
+    void grobner::pop_equation(equation& eq, equation_vector& v) {
+        unsigned idx = eq.idx();
         if (idx != v.size() - 1) {
             equation* eq2 = v.back();
             eq2->set_index(idx);
diff --git a/src/math/grobner/pdd_grobner.h b/src/math/grobner/pdd_grobner.h
index 4aaa1e6b8..8fd0dc37f 100644
--- a/src/math/grobner/pdd_grobner.h
+++ b/src/math/grobner/pdd_grobner.h
@@ -99,6 +99,7 @@ public:
     void add(pdd const& p) { add(p, nullptr); }
     void add(pdd const& p, u_dependency * dep);
 
+    void simplify_linear();
     void saturate();
 
     equation_vector const& equations();
@@ -111,7 +112,9 @@ public:
 private:
     bool step();
     bool basic_step();
+    bool basic_step(equation* e);
     equation* pick_next();
+    equation* pick_linear();
     bool canceled();
     bool done();
     void superpose(equation const& eq1, equation const& eq2);
@@ -128,7 +131,6 @@ private:
     bool is_too_complex(const equation& eq) const { return is_too_complex(eq.poly()); }
     bool is_too_complex(const pdd& p) const { return p.tree_size() > m_config.m_expr_size_limit;  }
 
-
     // tuned implementation
     vector<equation_vector> m_watch;           // watch list mapping variables to vector of equations where they occur (generally a subset)
     unsigned                m_levelp1;         // index into level+1
@@ -142,11 +144,16 @@ private:
     void add_to_watch(equation& eq);
 
     void del_equation(equation& eq);    
-    void retire(equation* eq) {
-        dealloc(eq);
-    }
-    void pop_equation(unsigned idx, equation_vector& v);
+    void retire(equation* eq) { dealloc(eq); }
+    void pop_equation(equation& eq, equation_vector& v);
+    void pop_equation(equation* eq, equation_vector& v) { pop_equation(*eq, v); }
     void push_equation(equation& eq, equation_vector& v);
+    void push_equation(equation* eq, equation_vector& v) { push_equation(*eq, v); }
+
+    struct compare_top_var;
+    bool simplify_linear_step();
+    void add_to_use(equation* e, vector<equation_vector>& use_list);
+
 
     void invariant() const;
     struct scoped_detach {
diff --git a/src/test/pdd_grobner.cpp b/src/test/pdd_grobner.cpp
index c1172311a..783821efc 100644
--- a/src/test/pdd_grobner.cpp
+++ b/src/test/pdd_grobner.cpp
@@ -1,6 +1,13 @@
 #include "util/rlimit.h"
 #include "math/grobner/pdd_grobner.h"
 
+#include "ast/bv_decl_plugin.h"
+#include "ast/ast_pp.h"
+#include "ast/reg_decl_plugins.h"
+#include "tactic/goal.h"
+#include "tactic/tactic.h"
+#include "tactic/bv/bit_blaster_tactic.h"
+
 namespace dd {
     void print_eqs(ptr_vector<grobner::equation> const& eqs) {
         std::cout << "eqs\n";
@@ -65,8 +72,97 @@ namespace dd {
 
         // 
     }
+
+    void add_def(unsigned_vector const& id2var, app* e, ast_manager& m, pdd_manager& p, grobner& g) {
+        expr* a, *b;
+        pdd v1 = p.mk_var(id2var[e->get_id()]);
+        if (m.is_and(e)) {
+            pdd r = p.one();
+            for (expr* arg : *e) r *= p.mk_var(id2var[arg->get_id()]);
+            g.add(v1 + r + 1);
+        }
+        else if (m.is_or(e)) {
+            pdd r = p.zero();
+            for (expr* arg : *e) r |= p.mk_var(id2var[arg->get_id()]);
+            g.add(v1 + r + 1);
+        }
+        else if (m.is_not(e, a)) {
+            pdd v2 = p.mk_var(id2var[a->get_id()]);
+            g.add(v1 + v2 + 1);
+        }
+        else if (m.is_eq(e, a, b)) {
+            pdd v2 = p.mk_var(id2var[a->get_id()]);
+            pdd v3 = p.mk_var(id2var[b->get_id()]);
+            g.add(v1 + v2 + v3 + 1);
+        }
+        else if (is_uninterp_const(e)) {
+            // pass
+        }
+        else {
+            UNREACHABLE();
+        }
+    }
+
+    void collect_id2var(unsigned_vector& id2var, expr_ref_vector const& fmls) {
+        svector<std::pair<unsigned, unsigned>> ds;
+        unsigned maxid = 0;
+        for (expr* e : subterms(fmls)) {
+            ds.push_back(std::make_pair(to_app(e)->get_depth(), e->get_id()));
+            maxid = std::max(maxid, e->get_id());
+        }
+        std::sort(ds.begin(), ds.end());
+        unsigned v = 1;
+        id2var.resize(maxid + 1);
+        for (auto p : ds) {
+            id2var[p.second] = v++;
+        }
+    }
+
+    void test_simplify(expr_ref_vector& fmls) {
+        ast_manager& m = fmls.get_manager();
+        unsigned_vector id2var;
+
+        collect_id2var(id2var, fmls);
+        pdd_manager p(id2var.size());
+        p.set_mod2_semantics();
+        grobner g(m.limit(), p);
+
+        for (expr* e : subterms(fmls)) {
+            add_def(id2var, to_app(e), m, p, g);
+        }
+        for (expr* e : fmls) {
+            g.add(p.mk_var(id2var[e->get_id()]) + 1);
+        }
+        g.simplify_linear();
+        g.display(std::cout);
+    }
+
+    void test2() {
+        ast_manager m;
+        reg_decl_plugins(m);
+        bv_util bv(m);
+        expr_ref x(m.mk_const("x", bv.mk_sort(4)), m);
+        expr_ref y(m.mk_const("y", bv.mk_sort(4)), m);
+        expr_ref xy(bv.mk_bv_mul(x, y), m);
+        expr_ref yx(bv.mk_bv_mul(y, x), m);
+        expr_ref eq(m.mk_not(m.mk_eq(xy,yx)), m);
+        goal_ref g = alloc(goal, m);
+        g->assert_expr(eq);
+        goal_ref_buffer res;
+        tactic_ref tac = mk_bit_blaster_tactic(m);
+        (*tac)(g, res);
+        g = res[0];
+
+        expr_ref_vector fmls(m);
+        for (unsigned i = 0; i < g->size(); ++i) {
+            fmls.push_back(g->form(i));        
+        }
+        test_simplify(fmls);
+        
+    }
 }
 
 void tst_pdd_grobner() {
     dd::test1();
+    dd::test2();
 }