diff --git a/src/util/lp/emonomials.cpp b/src/util/lp/emonomials.cpp
index d7875aaca..c5621ce19 100644
--- a/src/util/lp/emonomials.cpp
+++ b/src/util/lp/emonomials.cpp
@@ -233,7 +233,7 @@ namespace nla {
         if (m_cg_table.find(v, w)) {
             SASSERT(w != v);
             unsigned idxr = m_var2index[w];
-            unsigned idxl = m_canonized[idxr].m_next;
+            unsigned idxl = m_canonized[idxr].m_prev;
             m_canonized[idx].m_next  = idxr;
             m_canonized[idx].m_prev  = idxl;
             m_canonized[idxr].m_prev = idx;
@@ -324,13 +324,22 @@ namespace nla {
 
     // 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)) {
+        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();
+    }
+
+    emonomials::pf_iterator::pf_iterator(emonomials const& m, lpvar v, bool at_end):
+        m(m), m_mon(nullptr), m_it(iterator(m, m.head(v), at_end)), m_end(iterator(m, m.head(v), 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)) {
+            if (m_mon && m_mon->var() != (*m_it).var() && m.canonize_divides(*m_mon, *m_it) && !m.is_visited(*m_it)) {
+                m.set_visited(*m_it);
+                break;
+            }
+            if (!m_mon && !m.is_visited(*m_it)) {
                 m.set_visited(*m_it);
                 break;
             }
diff --git a/src/util/lp/emonomials.h b/src/util/lp/emonomials.h
index 847d6d2a3..e34e72c96 100644
--- a/src/util/lp/emonomials.h
+++ b/src/util/lp/emonomials.h
@@ -271,13 +271,15 @@ namespace nla {
         */
         class pf_iterator {
             emonomials const& m;
-            monomial const&   m_mon; // canonized monomial
+            monomial const*   m_mon; // monomial
+            lpvar             m_var;
             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);
+            pf_iterator(emonomials const& m, lpvar v, 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; }
@@ -287,15 +289,17 @@ namespace nla {
 
         class factors_of {
             emonomials const& m;
-            monomial const& mon;
+            monomial const* mon;
+            lpvar           m_var;
         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(emonomials const& m, monomial const& mon): m(m), mon(&mon), m_var(UINT_MAX) {}
+            factors_of(emonomials const& m, lpvar v): m(m), mon(nullptr), m_var(v) {}
+            pf_iterator begin() { if (mon) return pf_iterator(m, *mon, false); return pf_iterator(m, m_var, false); }
+            pf_iterator end() { if (mon) return pf_iterator(m, *mon, true); return pf_iterator(m, m_var, true); }
         };
 
         factors_of get_factors_of(monomial const& m) const { inc_visited(); return factors_of(*this, m); }
-        factors_of get_factors_of(lpvar v) const { return get_factors_of(var2monomial(v)); }
+        factors_of get_factors_of(lpvar v) const { inc_visited(); return factors_of(*this, v); }
        
         signed_vars const* find_canonical(svector<lpvar> const& vars) const;
 
diff --git a/src/util/lp/factorization.cpp b/src/util/lp/factorization.cpp
index da30c6cf3..48be88383 100644
--- a/src/util/lp/factorization.cpp
+++ b/src/util/lp/factorization.cpp
@@ -24,27 +24,23 @@ bool const_iterator_mon::get_factors(factor& k, factor& j, rational& sign) const
     std::sort(j_vars.begin(), j_vars.end());
 
     if (k_vars.size() == 1) {
-        k.var() = k_vars[0];
-        k.type() = factor_type::VAR;
+        k.set(k_vars[0], factor_type::VAR);
     } else {
         unsigned i;
         if (!m_ff->find_rm_monomial_of_vars(k_vars, i)) {
             return false;
         }
-        k.var() = i;
-        k.type() = factor_type::RM;
+        k.set(i, factor_type::RM);
     }
 
     if (j_vars.size() == 1) {
-        j.var() = j_vars[0];
-        j.type() = factor_type::VAR;
+        j.set(j_vars[0], factor_type::VAR);
     } else {
         unsigned i;
         if (!m_ff->find_rm_monomial_of_vars(j_vars, i)) {
             return false;
         }
-        j.var() = i;
-        j.type() = factor_type::RM;
+        j.set(i, factor_type::RM);
     }
     return true;
 }
@@ -90,22 +86,10 @@ bool const_iterator_mon::operator==(const const_iterator_mon::self_type &other)
         m_full_factorization_returned == other.m_full_factorization_returned &&
         m_mask == other.m_mask;
 }
+
 bool const_iterator_mon::operator!=(const const_iterator_mon::self_type &other) const { return !(*this == other); }
             
 factorization const_iterator_mon::create_binary_factorization(factor j, factor k) const {
-    // todo : the current explanation is an overkill
-    // std::function<void (expl_set&)> explain = [&](expl_set& exp){
-    //                                               const imp & impl = m_ff->m_impf;
-    //                                               unsigned mon_index = 0;
-    //                                               if (impl.m_var_to_its_monomial.find(k, mon_index)) {
-    //                                                   impl.add_explanation_of_reducing_to_rooted_monomial(impl.m_monomials[mon_index], exp);
-    //                                               }
-    //                                               if (impl.m_var_to_its_monomial.find(j, mon_index)) {
-    //                                                   impl.add_explanation_of_reducing_to_rooted_monomial(impl.m_monomials[mon_index], exp);
-    //                                               }
-
-    //                                               impl.add_explanation_of_reducing_to_rooted_monomial(m_ff->m_mon, exp);
-    //                                           };
     factorization f(nullptr);
     f.push_back(j);
     f.push_back(k);  
diff --git a/src/util/lp/factorization.h b/src/util/lp/factorization.h
index 56808824e..ff6908113 100644
--- a/src/util/lp/factorization.h
+++ b/src/util/lp/factorization.h
@@ -30,16 +30,14 @@ typedef unsigned lpvar;
 enum class factor_type { VAR, RM }; // RM stands for rooted monomial 
 
 class factor {
-    unsigned     m_var;
+    lpvar        m_var;
     factor_type  m_type;
 public:
-    factor() {}
-    explicit factor(unsigned j) : factor(j, factor_type::VAR) {}
-    factor(unsigned i, factor_type t) : m_var(i), m_type(t) {}
+    factor(): m_var(UINT_MAX), m_type(factor_type::VAR) {}
+    explicit factor(lpvar v, factor_type t) : m_var(v), m_type(t) {}
     unsigned var() const { return m_var; }
-    unsigned& var() { return m_var; }
     factor_type type() const { return m_type; }
-    factor_type& type() { return m_type; }
+    void set(lpvar v, factor_type t) { m_var = v; m_type = t; }
     bool is_var() const { return m_type == factor_type::VAR; }
     bool operator==(factor const& other) const {
         return m_var == other.var() && m_type == other.type();
@@ -51,12 +49,12 @@ public:
 
 
 class factorization {
-    vector<factor>        m_vars;
+    svector<factor>       m_vars;
     const monomial*       m_mon;
 public:
     factorization(const monomial* m): m_mon(m) {
         if (m != nullptr) {
-            for(lpvar j : *m)
+            for (lpvar j : *m)
                 m_vars.push_back(factor(j, factor_type::VAR));
         }
     }
@@ -66,7 +64,7 @@ public:
     size_t size() const { return m_vars.size(); }
     const factor* begin() const { return m_vars.begin(); }
     const factor* end() const { return m_vars.end(); }
-    void push_back(factor v) {
+    void push_back(factor const& v) {
         SASSERT(!is_mon());
         m_vars.push_back(v);
     }
@@ -122,7 +120,7 @@ struct factorization_factory {
         // repeating a pair twice, that is why m_mask is shorter by one then m_vars
         
         return
-            m_vars.size() != 2?
+            m_vars.size() != 2 ?
             svector<bool>(m_vars.size() - 1, false)
             :
             svector<bool>(1, true); // init mask as in the end() since the full iteration will do the job
diff --git a/src/util/lp/nla_core.cpp b/src/util/lp/nla_core.cpp
index 717ab2a48..198260903 100644
--- a/src/util/lp/nla_core.cpp
+++ b/src/util/lp/nla_core.cpp
@@ -182,7 +182,8 @@ std::ostream & core::print_factor(const factor& f, std::ostream& out) const {
 std::ostream & core::print_factor_with_vars(const factor& f, std::ostream& out) const {
     if (f.is_var()) {
         print_var(f.var(), out);
-    } else {
+    } 
+    else {
         out << " RM = " << m_emons.canonical[f.var()];
         out << "\n orig mon = " << m_emons[f.var()];
     }
@@ -191,8 +192,7 @@ std::ostream & core::print_factor_with_vars(const factor& f, std::ostream& out)
 
 std::ostream& core::print_monomial(const monomial& m, std::ostream& out) const {
     out << "( [" << m.var() << "] = " << m_lar_solver.get_variable_name(m.var()) << " = " << vvr(m.var()) << " = ";
-    print_product(m.vars(), out);
-    out << ")\n";
+    print_product(m.vars(), out) << ")\n";
     return out;
 }
 
@@ -207,13 +207,9 @@ std::ostream& core::print_tangent_domain(const point &a, const point &b, std::os
 }
 
 std::ostream& core::print_bfc(const bfc& m, std::ostream& out) const {
-    out << "( x = "; print_factor(m.m_x, out); out <<  ", y = "; print_factor(m.m_y, out); out <<  ")";
-    return out;
+    return out << "( x = "; print_factor(m.m_x, out); out <<  ", y = "; print_factor(m.m_y, out); out <<  ")";
 }
 
-std::ostream& core::print_monomial(unsigned i, std::ostream& out) const {
-    return print_monomial(m_emons[i], out);
-}
 
 std::ostream& core::print_monomial_with_vars(lpvar v, std::ostream& out) const {
     return print_monomial_with_vars(m_emons[v], out);
@@ -221,8 +217,7 @@ std::ostream& core::print_monomial_with_vars(lpvar v, std::ostream& out) const {
 
 template <typename T>
 std::ostream& core::print_product_with_vars(const T& m, std::ostream& out) const {
-    print_product(m, out);
-    out << '\n';
+    print_product(m, out) << "\n";
     for (unsigned k = 0; k < m.size(); k++) {
         print_var(m[k], out);
     }
@@ -772,8 +767,9 @@ std::ostream & core::print_ineqs(const lemma& l, std::ostream & out) const {
     
 std::ostream & core::print_factorization(const factorization& f, std::ostream& out) const {
     if (f.is_mon()){
-        print_monomial(*f.mon(), out << "is_mon ");
-    } else {
+        out << "is_mon " << pp_mon(*this, *f.mon());
+    } 
+    else {
         for (unsigned k = 0; k < f.size(); k++ ) {
             print_factor(f[k], out);
             if (k < f.size() - 1)
@@ -830,7 +826,7 @@ void core::trace_print_monomial_and_factorization(const signed_vars& rm, const f
     print_product(rm.vars(), out);
     out << "\n";
         
-    print_monomial(rm.var(), out << "mon:  ") << "\n";
+    out << "mon:  " << pp_mon(*this, rm.var()) << "\n";
     out << "value: " << vvr(rm) << "\n";
     print_factorization(f, out << "fact: ") << "\n";
 }
@@ -1483,7 +1479,7 @@ bool core::divide(const signed_vars& bc, const factor& c, factor & b) const {
     TRACE("nla_solver_div", tout << "b_vars = "; print_product(b_vars, tout););
     SASSERT(b_vars.size() > 0);
     if (b_vars.size() == 1) {
-        b = factor(b_vars[0]);
+        b = factor(b_vars[0], factor_type::VAR);
         return true;
     }
     signed_vars const* sv = m_emons.find_canonical(b_vars);
@@ -1670,8 +1666,8 @@ bool core::find_bfc_to_refine(bfc& bf, lpvar &j, rational& sign, const signed_va
             const monomial & m = m_emons[rm.var()];
             j = m.var();
             rm_found = nullptr;
-            bf.m_x = factor(m[0]);
-            bf.m_y = factor(m[1]);
+            bf.m_x = factor(m[0], factor_type::VAR);
+            bf.m_y = factor(m[1], factor_type::VAR);
             return true;
         }
                 
diff --git a/src/util/lp/nla_core.h b/src/util/lp/nla_core.h
index 0843b6347..3b72f85b9 100644
--- a/src/util/lp/nla_core.h
+++ b/src/util/lp/nla_core.h
@@ -151,7 +151,6 @@ public:
     void explain_var_separated_from_zero(lpvar j);
     void explain_fixed_var(lpvar j);
 
-
     std::ostream & print_ineq(const ineq & in, std::ostream & out) const;
     std::ostream & print_var(lpvar j, std::ostream & out) const;
     std::ostream & print_monomials(std::ostream & out) const;    
@@ -163,7 +162,6 @@ public:
     std::ostream & print_factor_with_vars(const factor& f, std::ostream& out) const;
     std::ostream& print_monomial(const monomial& m, std::ostream& out) const;
     std::ostream& print_bfc(const bfc& m, std::ostream& out) const;
-    std::ostream& print_monomial(unsigned i, std::ostream& out) const;
     std::ostream& print_monomial_with_vars(unsigned i, std::ostream& out) const;
     template <typename T>
     std::ostream& print_product_with_vars(const T& m, std::ostream& out) const;
@@ -343,5 +341,15 @@ public:
     
     lbool  test_check(vector<lemma>& l);
 };  // end of core
+
+    struct pp_mon {
+        core const& c;
+        monomial const& m;
+        pp_mon(core const& c, monomial const& m): c(c), m(m) {}
+        pp_mon(core const& c, lpvar v): c(c), m(c.m_emons[v]) {}
+    };
+
+    inline std::ostream& operator<<(std::ostream& out, pp_mon const& p) { return p.c.print_monomial(p.m, out); }
+
 } // end of namespace nla
 
diff --git a/src/util/lp/nla_order_lemmas.cpp b/src/util/lp/nla_order_lemmas.cpp
index 72aafa773..3c6fce009 100644
--- a/src/util/lp/nla_order_lemmas.cpp
+++ b/src/util/lp/nla_order_lemmas.cpp
@@ -276,7 +276,7 @@ void order::order_lemma_on_factor_binomial_rm(const monomial& ac, unsigned k, co
 }
 
 void order::order_lemma_on_binomial_ac_bd(const monomial& ac, unsigned k, const signed_vars& bd, const factor& b, lpvar d) {
-    TRACE("nla_solver",  print_monomial(ac, tout << "ac=");
+    TRACE("nla_solver",  tout << "ac=" << pp_mon(c(), ac);
           tout << "\nrm=" << bd;
           print_factor(b, tout << ", b="); print_var(d, tout << ", d=") << "\n";);
     int p = (k + 1) % 2;
@@ -312,7 +312,7 @@ void order::order_lemma_on_factor_binomial_explore(const monomial& m, unsigned k
 }
 
 void order::order_lemma_on_binomial(const monomial& ac) {
-    TRACE("nla_solver", print_monomial(ac, tout););
+    TRACE("nla_solver", tout << pp_mon(c(), ac););
     SASSERT(!check_monomial(ac) && ac.size() == 2);
     const rational & mult_val = vvr(ac[0]) * vvr(ac[1]);
     const rational acv = vvr(ac);