diff --git a/src/util/lp/factorization.cpp b/src/util/lp/factorization.cpp
index beb442791..ec112b6de 100644
--- a/src/util/lp/factorization.cpp
+++ b/src/util/lp/factorization.cpp
@@ -27,7 +27,7 @@ bool const_iterator_mon::get_factors(factor& k, factor& j, rational& sign) const
         k.set(k_vars[0], factor_type::VAR);
     } else {
         unsigned i;
-        if (!m_ff->find_rm_monomial_of_vars(k_vars, i)) {
+        if (!m_ff->find_canonical_monomial_of_vars(k_vars, i)) {
             return false;
         }
         k.set(i, factor_type::MON);
@@ -37,7 +37,7 @@ bool const_iterator_mon::get_factors(factor& k, factor& j, rational& sign) const
         j.set(j_vars[0], factor_type::VAR);
     } else {
         unsigned i;
-        if (!m_ff->find_rm_monomial_of_vars(j_vars, i)) {
+        if (!m_ff->find_canonical_monomial_of_vars(j_vars, i)) {
             return false;
         }
         j.set(i, factor_type::MON);
@@ -92,7 +92,14 @@ bool const_iterator_mon::operator!=(const const_iterator_mon::self_type &other)
 factorization const_iterator_mon::create_binary_factorization(factor j, factor k) const {
     factorization f(nullptr);
     f.push_back(j);
-    f.push_back(k);  
+    f.push_back(k);
+    // Let m by *m_ff->m_monomial, the monomial we factorize
+    // We have canonize_sign(m)*m.vars() = m.rvars()
+    // Let s = canonize_sign(f). Then we have f[0]*f[1] = s*m.rvars()
+    // s*canonize_sign(m)*val(m).
+    // Therefore val(m) = sign*val((f[0])*val(f[1]), where sign = canonize_sign(f)*canonize_sign(m).
+    // We apply this sign to the first factor.
+    f[0].sign() ^= (m_ff->canonize_sign(f)^m_ff->canonize_sign(*m_ff->m_monomial));
     return f;
 }
 
diff --git a/src/util/lp/factorization.h b/src/util/lp/factorization.h
index c0ad8ab62..860acebf1 100644
--- a/src/util/lp/factorization.h
+++ b/src/util/lp/factorization.h
@@ -54,26 +54,30 @@ public:
 
 
 class factorization {
-    svector<factor>       m_vars;
+    svector<factor>       m_factors;
     const monomial*       m_mon;
 public:
     factorization(const monomial* m): m_mon(m) {
         if (m != nullptr) {
             for (lpvar j : m->vars())
-                m_vars.push_back(factor(j, factor_type::VAR));
+                m_factors.push_back(factor(j, factor_type::VAR));
         }
     }
-    bool is_mon() const { return m_mon != nullptr;}
-    bool is_empty() const { return m_vars.empty(); }
-    factor operator[](unsigned k) const { return m_vars[k]; }
-    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 const& v) {
-        SASSERT(!is_mon());
-        m_vars.push_back(v);
+    bool is_mon() const {
+        return m_mon != nullptr;
     }
-    const monomial* mon() const { return m_mon; }
+    bool is_empty() const { return m_factors.empty(); }
+    const factor& operator[](unsigned k) const { return m_factors[k]; }
+    factor& operator[](unsigned k) { return m_factors[k]; }
+    size_t size() const { return m_factors.size(); }
+    const factor* begin() const { return m_factors.begin(); }
+    const factor* end() const { return m_factors.end(); }
+    void push_back(factor const& v) {
+        m_factors.push_back(v);
+    }
+    const monomial& mon() const { return *m_mon; }
+    void set_mon(const monomial* m) { m_mon = m; }
+
 };
 
 struct const_iterator_mon {
@@ -112,9 +116,9 @@ struct factorization_factory {
     const svector<lpvar>&  m_vars;
     const monomial*       m_monomial;
     // returns true if found
-    virtual bool find_rm_monomial_of_vars(const svector<lpvar>& vars, unsigned& i) const = 0;
-    virtual const monomial* find_monomial_of_vars(const svector<lpvar>& vars) const = 0;
-    
+    virtual bool find_canonical_monomial_of_vars(const svector<lpvar>& vars, unsigned& i) const = 0;
+    virtual bool canonize_sign(const monomial& m) const = 0;
+    virtual bool canonize_sign(const factorization& m) const = 0;
 
     factorization_factory(const svector<lpvar>& vars, const monomial* m) :
         m_vars(vars), m_monomial(m) {
@@ -140,8 +144,7 @@ struct factorization_factory {
         auto it = const_iterator_mon(mask, this);
         it.m_full_factorization_returned = true;
         return it;
-    } 
-   
+    }   
 };
 
 }
diff --git a/src/util/lp/factorization_factory_imp.cpp b/src/util/lp/factorization_factory_imp.cpp
index 0228eaa2a..faed2ea2e 100644
--- a/src/util/lp/factorization_factory_imp.cpp
+++ b/src/util/lp/factorization_factory_imp.cpp
@@ -25,10 +25,14 @@ factorization_factory_imp::factorization_factory_imp(const monomial& rm, const c
     factorization_factory(rm.rvars(), &s.m_emons[rm.var()]),
     m_core(s), m_mon(s.m_emons[rm.var()]), m_rm(rm) { }
         
-bool factorization_factory_imp::find_rm_monomial_of_vars(const svector<lpvar>& vars, unsigned & i) const {
-    return m_core.find_rm_monomial_of_vars(vars, i);
+bool factorization_factory_imp::find_canonical_monomial_of_vars(const svector<lpvar>& vars, unsigned & i) const {
+    return m_core.find_canonical_monomial_of_vars(vars, i);
 }
-const monomial* factorization_factory_imp::find_monomial_of_vars(const svector<lpvar>& vars) const {
-    return m_core.find_monomial_of_vars(vars);
+bool factorization_factory_imp::canonize_sign(const monomial& m) const {
+    return m_core.canonize_sign(m);
 }
+bool factorization_factory_imp::canonize_sign(const factorization& f) const {
+    return m_core.canonize_sign(f);
+}
+
 }
diff --git a/src/util/lp/factorization_factory_imp.h b/src/util/lp/factorization_factory_imp.h
index ae1c8d9a0..7278d9af8 100644
--- a/src/util/lp/factorization_factory_imp.h
+++ b/src/util/lp/factorization_factory_imp.h
@@ -28,7 +28,8 @@ namespace nla {
         const monomial& m_rm;
         
         factorization_factory_imp(const monomial& rm, const core& s);
-        bool find_rm_monomial_of_vars(const svector<lpvar>& vars, unsigned & i) const;
-        const monomial* find_monomial_of_vars(const svector<lpvar>& vars) const;
-    };
+        bool find_canonical_monomial_of_vars(const svector<lpvar>& vars, unsigned & i) const;
+        virtual bool canonize_sign(const monomial& m) const;
+        virtual bool canonize_sign(const factorization& m) const;
+ };
 }
diff --git a/src/util/lp/monomial.h b/src/util/lp/monomial.h
index ead38b4f6..48e2f4a35 100644
--- a/src/util/lp/monomial.h
+++ b/src/util/lp/monomial.h
@@ -70,7 +70,7 @@ public:
 };
 
  inline std::ostream& operator<<(std::ostream& out, monomial const& m) {
-     return out << m.var() << " := " << m.vars() << " r " << m.rsign() << " * " << m.rvars();
+     return out << m.var() << " := " << m.vars() << " r ( " << sign_to_rat(m.rsign()) << " * " << m.rvars() << ")";
  }
 
 
diff --git a/src/util/lp/nla_basics_lemmas.cpp b/src/util/lp/nla_basics_lemmas.cpp
index b5bf55f01..8588d0673 100644
--- a/src/util/lp/nla_basics_lemmas.cpp
+++ b/src/util/lp/nla_basics_lemmas.cpp
@@ -446,20 +446,20 @@ void basics::generate_pl_on_mon(const monomial& m, unsigned factor_index) {
     
 // none of the factors is zero and the product is not zero
 // -> |fc[factor_index]| <= |rm|
-void basics::generate_pl(const monomial& rm, const factorization& fc, int factor_index) {
-    TRACE("nla_solver", tout << "factor_index = " << factor_index << ", rm = "
-          << pp_mon(c(), rm);
+void basics::generate_pl(const monomial& m, const factorization& fc, int factor_index) {
+    TRACE("nla_solver", tout << "factor_index = " << factor_index << ", m = "
+          << pp_mon(c(), m);
           tout << ", fc = "; c().print_factorization(fc, tout);
-          tout << "orig mon = "; c().print_monomial(c().m_emons[rm.var()], tout););
+          tout << "orig mon = "; c().print_monomial(c().emons()[m.var()], tout););
     if (fc.is_mon()) {
-        generate_pl_on_mon(*fc.mon(), factor_index);
+        generate_pl_on_mon(m, factor_index);
         return;
     }
     add_empty_lemma();
     int fi = 0;
-    rational rmv = val(rm);
-    rational sm = rational(nla::rat_sign(rmv));
-    unsigned mon_var = var(rm);
+    rational mv = val(m);
+    rational sm = rational(nla::rat_sign(mv));
+    unsigned mon_var = var(m);
     c().mk_ineq(sm, mon_var, llc::LT);
     for (factor f : fc) {
         if (fi++ != factor_index) {
@@ -468,14 +468,14 @@ void basics::generate_pl(const monomial& rm, const factorization& fc, int factor
             lpvar j = var(f);
             rational jv = val(j);
             rational sj = rational(nla::rat_sign(jv));
-            SASSERT(sm*rmv < sj*jv);
+            SASSERT(sm*mv < sj*jv);
             c().mk_ineq(sj, j, llc::LT);
             c().mk_ineq(sm, mon_var, -sj, j, llc::GE );
         }
     }
     if (!fc.is_mon()) {
         explain(fc);
-        explain(rm);
+        explain(m);
     }
     TRACE("nla_solver", c().print_lemma(tout); );
 }
@@ -691,8 +691,8 @@ bool basics::basic_lemma_for_mon_neutral_monomial_to_factor_model_based(const mo
 
 void basics::basic_lemma_for_mon_neutral_model_based(const monomial& rm, const factorization& f) {
     if (f.is_mon()) {
-        basic_lemma_for_mon_neutral_monomial_to_factor_model_based_fm(*f.mon());
-        basic_lemma_for_mon_neutral_from_factors_to_monomial_model_based_fm(*f.mon());
+        basic_lemma_for_mon_neutral_monomial_to_factor_model_based_fm(f.mon());
+        basic_lemma_for_mon_neutral_from_factors_to_monomial_model_based_fm(f.mon());
     }
     else {
         basic_lemma_for_mon_neutral_monomial_to_factor_model_based(rm, f);
@@ -701,9 +701,10 @@ void basics::basic_lemma_for_mon_neutral_model_based(const monomial& rm, const f
 }
 // use the fact
 // 1 * 1 ... * 1 * x * 1 ... * 1 = x
-bool basics::basic_lemma_for_mon_neutral_from_factors_to_monomial_model_based(const monomial& rm, const factorization& f) {
-    rational sign = sign_to_rat(rm.rsign());
-    TRACE("nla_solver_bl", tout << pp_rmon(_(), rm) <<"\nf = "; c().print_factorization(f, tout); tout << "sign = " << sign << '\n'; );
+bool basics::basic_lemma_for_mon_neutral_from_factors_to_monomial_model_based(const monomial& m, const factorization& f) {
+    rational sign = sign_to_rat(m.rsign());
+    SASSERT(m.rsign() == canonize_sign(f));
+    TRACE("nla_solver_bl", tout << pp_rmon(_(), m) <<"\nf = "; c().print_factorization(f, tout); tout << "sign = " << sign << '\n'; );
     lpvar not_one = -1;
     for (auto j : f){
         TRACE("nla_solver_bl", tout << "j = "; c().print_factor_with_vars(j, tout););
@@ -728,13 +729,13 @@ bool basics::basic_lemma_for_mon_neutral_from_factors_to_monomial_model_based(co
 
     if (not_one + 1) {
         // we found the only not_one
-        if (val(rm) == val(not_one) * sign) {
+        if (val(m) == val(not_one) * sign) {
             TRACE("nla_solver", tout << "the whole is equal to the factor" << std::endl;);
             return false;
         }
     } else {
         // we have +-ones only in the factorization
-        if (val(rm) == sign) {
+        if (val(m) == sign) {
             return false;
         }
     }
@@ -746,19 +747,20 @@ bool basics::basic_lemma_for_mon_neutral_from_factors_to_monomial_model_based(co
     for (auto j : f){
         lpvar var_j = var(j);
         if (not_one == var_j) continue;
-        c().mk_ineq(var_j, llc::NE, j.is_var()? val(j) : sign_to_rat(c().canonize_sign(j)) * val(j));   
+        TRACE("nla_solver_bl", tout << "j = "; c().print_factor_with_vars(j, tout);); 
+        c().mk_ineq(var_j, llc::NE, val(var_j));   
     }
 
     if (not_one == static_cast<lpvar>(-1)) {
-        c().mk_ineq(rm.var(), llc::EQ, sign);
+        c().mk_ineq(m.var(), llc::EQ, sign);
     } else {
-        c().mk_ineq(rm.var(), -sign, not_one, llc::EQ);
+        c().mk_ineq(m.var(), -sign, not_one, llc::EQ);
     }
-    explain(rm);
+    explain(m);
     explain(f);
     TRACE("nla_solver",
           c().print_lemma(tout);
-          tout << "rm = " << pp_rmon(c(), rm);
+          tout << "m = " << pp_rmon(c(), m);
           );
     return true;
 }
@@ -777,7 +779,7 @@ void basics::basic_lemma_for_mon_non_zero_model_based_mf(const factorization& f)
     if (zero_j == -1) { return; } 
     add_empty_lemma();
     c().mk_ineq(zero_j, llc::NE);
-    c().mk_ineq(f.mon()->var(), llc::EQ);
+    c().mk_ineq(f.mon().var(), llc::EQ);
     TRACE("nla_solver", c().print_lemma(tout););
 }
 
diff --git a/src/util/lp/nla_common.cpp b/src/util/lp/nla_common.cpp
index fe156bdb6..97e3942c6 100644
--- a/src/util/lp/nla_common.cpp
+++ b/src/util/lp/nla_common.cpp
@@ -47,6 +47,7 @@ template <typename T> bool common::canonize_sign(const T& t) const {
 template bool common::canonize_sign<monomial>(const monomial&) const;
 template bool common::canonize_sign<factor>(const factor&) const;
 template bool common::canonize_sign<lpvar>(const lpvar&) const;
+template bool common::canonize_sign<factorization>(const factorization&) const;
 
 void common::mk_ineq(lp::lar_term& t, llc cmp, const rational& rs){
     c().mk_ineq(t, cmp, rs);
diff --git a/src/util/lp/nla_core.cpp b/src/util/lp/nla_core.cpp
index 02ad4ef40..6fd94ab3a 100644
--- a/src/util/lp/nla_core.cpp
+++ b/src/util/lp/nla_core.cpp
@@ -93,7 +93,7 @@ svector<lpvar> core::sorted_rvars(const factor& f) const {
 // the value of the factor is equal to the value of the variable multiplied
 // by the canonize_sign
 bool core::canonize_sign(const factor& f) const {
-    return f.sign() ^ (f.is_var()? canonize_sign(f.var()) : m_emons[f.var()].rsign());
+    return f.sign() ^ (f.is_var()? canonize_sign(f.var()) : canonize_sign(m_emons[f.var()]));
 }
 
 bool core::canonize_sign(lpvar j) const {
@@ -113,6 +113,14 @@ bool core::canonize_sign(const monomial& m) const {
     return m.rsign();
 }
 
+bool core::canonize_sign(const factorization& f) const {
+    bool r = false;
+    for (const factor & a : f) {
+        r ^= canonize_sign(a);
+    }
+    return r;
+}
+
 void core::add(lpvar v, unsigned sz, lpvar const* vs) {
     m_emons.add(v, sz, vs);
 }
@@ -180,7 +188,7 @@ std::ostream & core::print_factor(const factor& f, std::ostream& out) const {
         out << "VAR,  ";
         print_var(f.var(), out);
     } else {
-        out << "MON, ";
+        out << "MON, v" << m_emons[f.var()] << " = ";
         print_product(m_emons[f.var()].rvars(), out);
     }
     out << "\n";
@@ -192,7 +200,7 @@ std::ostream & core::print_factor_with_vars(const factor& f, std::ostream& out)
         print_var(f.var(), out);
     } 
     else {
-        out << " RM = " << pp_rmon(*this, m_emons[f.var()]);
+        out << " MON = " << pp_rmon(*this, m_emons[f.var()]);
     }
     return out;
 }
@@ -216,11 +224,12 @@ std::ostream& core::print_tangent_domain(const point &a, const point &b, std::os
     return out;
 }
 
-std::ostream& core::print_bfc(const bfc& m, std::ostream& out) const {
+std::ostream& core::print_bfc(const factorization& m, std::ostream& out) const {
+    SASSERT(m.size() == 2);
     out << "( x = ";
-    print_factor(m.m_x, out);
-    out <<  ", y = ";
-    print_factor(m.m_y, out); out <<  ")";
+    print_factor(m[0], out);
+    out <<  "* y = ";
+    print_factor(m[1], out); out <<  ")";
     return out;
 }
 
@@ -748,7 +757,8 @@ std::ostream & core::print_var(lpvar j, std::ostream & out) const {
     }
         
     m_lar_solver.print_column_info(j, out);
-    out << "root=" << m_evars.find(j) << "\n";
+    signed_var jr = m_evars.find(j);
+    out << "root=" << (jr.sign()? "-v":"v") << jr.var() << "\n";
     return out;
 }
 
@@ -783,7 +793,7 @@ 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()){
-        out << "is_mon " << pp_mon(*this, *f.mon());
+        out << "is_mon " << pp_mon(*this, f.mon());
     } 
     else {
         for (unsigned k = 0; k < f.size(); k++ ) {
@@ -795,18 +805,12 @@ std::ostream & core::print_factorization(const factorization& f, std::ostream& o
     return out;
 }
     
-bool core:: find_rm_monomial_of_vars(const svector<lpvar>& vars, unsigned & i) const {
+bool core::find_canonical_monomial_of_vars(const svector<lpvar>& vars, unsigned & i) const {
     SASSERT(vars_are_roots(vars));
     monomial const* sv = m_emons.find_canonical(vars);
     return sv && (i = sv->var(), true);
 }
 
-const monomial* core::find_monomial_of_vars(const svector<lpvar>& vars) const {
-    monomial const* sv = m_emons.find_canonical(vars);
-    return sv ? &m_emons[sv->var()] : nullptr;
-}
-
-
 void core::explain_existing_lower_bound(lpvar j) {
     SASSERT(has_lower_bound(j));
     current_expl().add(m_lar_solver.get_column_lower_bound_witness(j));
@@ -1571,13 +1575,17 @@ void core::add_abs_bound(lpvar v, llc cmp, rational const& bound) {
 */
 
     
-bool core::find_bfc_to_refine_on_rmonomial(const monomial& rm, bfc & bf) {
-    for (auto factorization : factorization_factory_imp(rm, *this)) {
-        if (factorization.size() == 2) {
-            auto a = factorization[0];
-            auto b = factorization[1];
-            if (val(rm) != val(a) * val(b)) {
-                bf = bfc(a, b);
+bool core::find_bfc_to_refine_on_monomial(const monomial& m, factorization & bf) {
+    for (auto f : factorization_factory_imp(m, *this)) {
+        if (f.size() == 2) {
+            auto a = f[0];
+            auto b = f[1];
+            if (val(m) != val(a) * val(b)) {
+                bf = f;
+                TRACE("nla_solver", tout << "found bf";
+                      tout << ":m:" << pp_rmon(*this, m) << "\n";
+                      tout << "bf:"; print_bfc(bf, tout););
+                      
                 return true;
             }
         }
@@ -1585,30 +1593,23 @@ bool core::find_bfc_to_refine_on_rmonomial(const monomial& rm, bfc & bf) {
     return false;
 }
     
-bool core::find_bfc_to_refine(bfc& bf, lpvar &j, rational& sign, const monomial*& rm_found){
-    rm_found = nullptr;
+bool core::find_bfc_to_refine(const monomial* & m, factorization & bf){
+    m = nullptr;
+    // todo: randomise loop
     for (unsigned i: m_to_refine) {
-        const auto& rm = m_emons[i];
-        SASSERT (!check_monomial(m_emons[rm.var()]));
-        if (rm.size() == 2) {
-            sign = rational(1);
-            const monomial & m = m_emons[rm.var()];
-            j = m.var();
-            rm_found = nullptr;
-            bf.m_x = factor(m.vars()[0], factor_type::VAR);
-            bf.m_y = factor(m.vars()[1], factor_type::VAR);
+        m = &m_emons[i];
+        SASSERT (!check_monomial(*m));
+        if (m->size() == 2) {
+            bf.set_mon(m);
+            bf.push_back(factor(m->vars()[0], factor_type::VAR));
+            bf.push_back(factor(m->vars()[1], factor_type::VAR));
             return true;
         }
                 
-        rm_found = &rm;
-        if (find_bfc_to_refine_on_rmonomial(rm, bf)) {
-            j = rm.var();
-            sign = sign_to_rat(rm.rsign());
-            TRACE("nla_solver", tout << "found bf";
-                  tout << ":rm:" << pp_rmon(*this, rm) << "\n";
-                  tout << "bf:"; print_bfc(bf, tout);
-                  tout << ", product = " << val(rm) << ", but should be =" << val(bf.m_x)*val(bf.m_y);
-                  tout << ", j == "; print_var(j, tout) << "\n";);
+        if (find_bfc_to_refine_on_monomial(*m, bf)) {
+            TRACE("nla_solver",
+                  tout << "bf = "; print_factorization(bf, tout);
+                  tout << "\nval(*m) = " << val(*m) << ", should be = (val(bf[0])=" << val(bf[0]) << ")*(val(bf[1]) = " << val(bf[1]) << ") = " << val(bf[0])*val(bf[1]) << "\n";);
             return true;
         } 
     }
diff --git a/src/util/lp/nla_core.h b/src/util/lp/nla_core.h
index 0a159d97a..cf182e219 100644
--- a/src/util/lp/nla_core.h
+++ b/src/util/lp/nla_core.h
@@ -117,6 +117,7 @@ public:
     // the value of the factor is equal to the value of the variable multiplied
     // by the canonize_sign
     bool canonize_sign(const factor& f) const;
+    bool canonize_sign(const factorization& f) const;
 
     bool canonize_sign(lpvar j) const;
     
@@ -159,7 +160,7 @@ public:
     std::ostream & print_factor(const factor& f, std::ostream& out) const;
     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_bfc(const factorization& m, 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;
@@ -238,9 +239,7 @@ public:
 
     bool var_is_fixed(lpvar j) const;
         
-    bool find_rm_monomial_of_vars(const svector<lpvar>& vars, unsigned & i) const;
-
-    const monomial* find_monomial_of_vars(const svector<lpvar>& vars) const;
+    bool find_canonical_monomial_of_vars(const svector<lpvar>& vars, unsigned & i) const;
 
     bool var_has_positive_lower_bound(lpvar j) const;
 
@@ -320,9 +319,9 @@ public:
     void add_abs_bound(lpvar v, llc cmp);
     void add_abs_bound(lpvar v, llc cmp, rational const& bound);
 
-    bool  find_bfc_to_refine_on_rmonomial(const monomial& rm, bfc & bf);
+    bool  find_bfc_to_refine_on_monomial(const monomial&, factorization & bf);
     
-    bool  find_bfc_to_refine(bfc& bf, lpvar &j, rational& sign, const monomial*& rm_found  );
+    bool  find_bfc_to_refine(const monomial* & m, factorization& bf);
     void generate_simple_sign_lemma(const rational& sign, const monomial& m);
 
     void negate_relation(unsigned j, const rational& a);
@@ -349,7 +348,6 @@ struct pp_rmon {
     pp_rmon(core const& c, monomial const& m): c(c), m(m) {}
     pp_rmon(core const& c, lpvar v): c(c), m(c.m_emons[v]) {}
 };
-inline rational sign_to_rat(bool s) { return rational(s? -1 : 1); }
 inline std::ostream& operator<<(std::ostream& out, pp_mon const& p) { return p.c.print_monomial(p.m, out); }
 inline std::ostream& operator<<(std::ostream& out, pp_rmon const& p) { return p.c.print_monomial_with_vars(p.m, out); }
 
diff --git a/src/util/lp/nla_defs.h b/src/util/lp/nla_defs.h
index aafe6de11..4d3f13069 100644
--- a/src/util/lp/nla_defs.h
+++ b/src/util/lp/nla_defs.h
@@ -99,5 +99,5 @@ bool uniform_le(const T& a, const T& b,  unsigned & strict_i) {
     if (z_b) {strict_i = -1;}
     return true;
 }
-
+inline rational sign_to_rat(bool s) { return rational(s? -1 : 1); }
 }
diff --git a/src/util/lp/nla_order_lemmas.cpp b/src/util/lp/nla_order_lemmas.cpp
index f3ced5bc0..4916c3a32 100644
--- a/src/util/lp/nla_order_lemmas.cpp
+++ b/src/util/lp/nla_order_lemmas.cpp
@@ -49,7 +49,7 @@ void order::order_lemma_on_rmonomial(const monomial& m) {
         if (ac.size() != 2)
             continue;
         if (ac.is_mon())
-            order_lemma_on_binomial(*ac.mon());
+            order_lemma_on_binomial(ac.mon());
         else
             order_lemma_on_factorization(m, ac);
         if (done())
diff --git a/src/util/lp/nla_tangent_lemmas.cpp b/src/util/lp/nla_tangent_lemmas.cpp
index 479a7e1e9..941ab0619 100644
--- a/src/util/lp/nla_tangent_lemmas.cpp
+++ b/src/util/lp/nla_tangent_lemmas.cpp
@@ -29,10 +29,22 @@ tangents::tangents(core * c) : common(c) {}
 std::ostream& tangents::print_tangent_domain(const point &a, const point &b, std::ostream& out) const {
     return out << "(" << a <<  ", " << b << ")";
 }
+unsigned tangents::find_binomial_to_refine() {
+    unsigned start = c().random();
+    unsigned sz = c().m_to_refine.size();
+    for (unsigned k = 0; k < sz; k++) {
+        lpvar j = c().m_to_refine[(k + start) % sz];
+        if (c().emons()[j].size() == 2)
+            return j;
+    }
+    return -1;
+}
 
-void tangents::generate_simple_tangent_lemma(const monomial& m) {
-    if (m.size() != 2)
-        return;
+void tangents::generate_simple_tangent_lemma() {
+    lpvar j = find_binomial_to_refine();
+    if (!is_set(j)) return;
+    const monomial& m = c().emons()[j];
+    SASSERT(m.size() != 2);
     TRACE("nla_solver", tout << "m:" << pp_mon(c(), m) << std::endl;);
     c().add_empty_lemma();
     const rational v = c().product_value(m.vars());
@@ -66,65 +78,70 @@ void tangents::generate_simple_tangent_lemma(const monomial& m) {
 }
     
 void tangents::tangent_lemma() {
-    bfc bf;
-    lpvar j;
-    rational sign;
-    const monomial* rm = nullptr;
-        
-    if (c().find_bfc_to_refine(bf, j, sign, rm)) {
-        tangent_lemma_bf(bf, j, sign, rm);
-    } else {
-        TRACE("nla_solver", tout << "cannot find a bfc to refine\n"; );
-        if (rm != nullptr)
-            generate_simple_tangent_lemma(*rm);
+    if (!c().m_settings.run_tangents()) {
+        TRACE("nla_solver", tout << "not generating tangent lemmas\n";);
+        return;
+    }
+    factorization bf(nullptr);
+    const monomial* m;
+    if (c().find_bfc_to_refine(m, bf)) {
+        tangent_lemma_bf(*m, bf);
     }
 }
 
-void tangents::generate_explanations_of_tang_lemma(const monomial& rm, const bfc& bf, lp::explanation& exp) {
+void tangents::generate_explanations_of_tang_lemma(const monomial& rm, const factorization& bf, lp::explanation& exp) {
     // here we repeat the same explanation for each lemma
     c().explain(rm, exp);
-    c().explain(bf.m_x, exp);
-    c().explain(bf.m_y, exp);
+    c().explain(bf[0], exp);
+    c().explain(bf[1], exp);
 }
 
-void tangents::generate_tang_plane(const rational & a, const rational& b, const factor& x, const factor& y, bool below, lpvar j, const rational& j_sign) {
+void tangents::generate_tang_plane(const rational & a, const rational& b, const factor& x, const factor& y, bool below, lpvar j) {
     lpvar jx = var(x);
     lpvar jy = var(y);
     add_empty_lemma();
     c().negate_relation(jx, a);
     c().negate_relation(jy, b);
-    bool sbelow = j_sign.is_pos()? below: !below;
 #if Z3DEBUG
     int mult_sign = nla::rat_sign(a - val(jx))*nla::rat_sign(b - val(jy));
-    SASSERT((mult_sign == 1) == sbelow);
+    SASSERT((mult_sign == 1) == below);
     // If "mult_sign is 1"  then (a - x)(b-y) > 0 and ab - bx - ay + xy > 0
     // or -ab + bx + ay < xy or -ay - bx + xy > -ab
-    // j_sign*val(j) stands for xy. So, finally we have  -ay - bx + j_sign*j > - ab
+    // val(j) stands for xy. So, finally we have  -ay - bx + j > - ab
 #endif
 
     lp::lar_term t;
     t.add_coeff_var(-a, jy);
     t.add_coeff_var(-b, jx);
-    t.add_coeff_var( j_sign, j);
-    c().mk_ineq(t, sbelow? llc::GT : llc::LT, - a*b);
-}  
-void tangents::tangent_lemma_bf(const bfc& bf, lpvar j, const rational& sign, const monomial* rm){
+    t.add_var(j);
+    c().mk_ineq(t, below? llc::GT : llc::LT, - a*b);
+}
+
+void tangents::tangent_lemma_bf(const monomial& m, const factorization& bf){
     point a, b;
-    point xy (val(bf.m_x), val(bf.m_y));
+    point xy (val(bf[0]), val(bf[1]));
     rational correct_mult_val =  xy.x * xy.y;
-    rational v = val(j) * sign;
+    lpvar j =m.var();
+    // We have canonize_sign(m)*m.vars() = m.rvars()
+    // Let s = canonize_sign(bf). Then we have bf[1]*bf[1] = s*m.rvars()
+    // s*canonize_sign(m)*val(m).
+    // Therefore sign*val(m) = val((bf[0])*val(bf[1]), where sign = canonize_sign(bf)*canonize_sign(m)
+    
+    SASSERT(canonize_sign(bf) == canonize_sign(m));
+    rational v = val(j);    
     bool below = v < correct_mult_val;
-    TRACE("nla_solver", tout << "rm = " << rm << ", below = " << below << std::endl; );
+    TRACE("nla_solver", tout << "below = " << below << std::endl; );
     get_tang_points(a, b, below, v, xy);
-    TRACE("nla_solver", tout << "sign = " << sign << ", tang domain = "; print_tangent_domain(a, b, tout); tout << std::endl;);
+    TRACE("nla_solver", tout << "tang domain = "; print_tangent_domain(a, b, tout); tout << std::endl;);
     unsigned lemmas_size_was = c().m_lemma_vec->size();
-    generate_two_tang_lines(bf, xy, sign, j);
-    generate_tang_plane(a.x, a.y, bf.m_x, bf.m_y, below, j, sign);
-    generate_tang_plane(b.x, b.y, bf.m_x, bf.m_y, below, j, sign);
-    // if rm == nullptr there is no need to explain equivs since we work on a monomial and not on a rooted monomial
-    if (rm != nullptr) { 
+    rational sign(1);
+    generate_two_tang_lines(bf, xy, j);
+    generate_tang_plane(a.x, a.y, bf[0], bf[1], below, j);
+    generate_tang_plane(b.x, b.y, bf[0], bf[1], below, j);
+
+    if (!bf.is_mon()) { 
         lp::explanation expl;
-        generate_explanations_of_tang_lemma(*rm, bf, expl);
+        generate_explanations_of_tang_lemma(m, bf, expl);
         for (unsigned i = lemmas_size_was; i < c().m_lemma_vec->size(); i++) {
             auto &l = ((*c().m_lemma_vec)[i]);
             l.expl().add(expl);
@@ -135,14 +152,14 @@ void tangents::tangent_lemma_bf(const bfc& bf, lpvar j, const rational& sign, co
               c().print_specific_lemma((*c().m_lemma_vec)[i], tout); );
 }
     
-void tangents::generate_two_tang_lines(const bfc & bf, const point& xy, const rational& sign, lpvar j) {
+void tangents::generate_two_tang_lines(const factorization & bf, const point& xy, lpvar j) {
     add_empty_lemma();
-    c().mk_ineq(var(bf.m_x), llc::NE, xy.x);
-    c().mk_ineq(sign, j, - xy.x, var(bf.m_y), llc::EQ);
+    c().mk_ineq(var(bf[0]), llc::NE, xy.x);
+    c().mk_ineq(j, - xy.x, var(bf[1]), llc::EQ);
         
     add_empty_lemma();
-    c().mk_ineq(var(bf.m_y), llc::NE, xy.y);
-    c().mk_ineq(sign, j, - xy.y, var(bf.m_x), llc::EQ);
+    c().mk_ineq(var(bf[1]), llc::NE, xy.y);
+    c().mk_ineq(j, - xy.y, var(bf[0]), llc::EQ);
 }
 
 // Get two planes tangent to surface z = xy, one at point a,  and another at point b.
@@ -168,7 +185,7 @@ void tangents::push_tang_point(point &a, const point& xy, bool below, const rati
     while (steps--) {
         del *= rational(2);
         point na = xy + del;
-        TRACE("nla_solver", tout << "del = " << del << std::endl;);
+        TRACE("nla_solver_tp", tout << "del = " << del << std::endl;);
         if (!plane_is_correct_cut(na, xy, correct_val, val, below)) {
             TRACE("nla_solver_tp", tout << "exit";tout << std::endl;);
             return;
diff --git a/src/util/lp/nla_tangent_lemmas.h b/src/util/lp/nla_tangent_lemmas.h
index c2e608a3c..29f919edd 100644
--- a/src/util/lp/nla_tangent_lemmas.h
+++ b/src/util/lp/nla_tangent_lemmas.h
@@ -52,16 +52,16 @@ class tangents : common {
 public:
     tangents(core *core);
     void tangent_lemma();
-private:
+private:    
+    lpvar find_binomial_to_refine();
+    void generate_simple_tangent_lemma();    
 
-    void generate_simple_tangent_lemma(const monomial&);    
-
-    void generate_explanations_of_tang_lemma(const monomial& rm, const bfc& bf, lp::explanation& exp);
+    void generate_explanations_of_tang_lemma(const monomial& m, const factorization& bf, lp::explanation& exp);
     
-    void tangent_lemma_bf(const bfc& bf, lpvar j, const rational& sign, const monomial* rm);
-    void generate_tang_plane(const rational & a, const rational& b, const factor& x, const factor& y, bool below, lpvar j, const rational& j_sign);
+    void tangent_lemma_bf(const monomial& m,const factorization& bf);
+    void generate_tang_plane(const rational & a, const rational& b, const factor& x, const factor& y, bool below, lpvar j);
     
-    void generate_two_tang_lines(const bfc & bf, const point& xy, const rational& sign, lpvar j);
+    void generate_two_tang_lines(const factorization & bf, const point& xy, lpvar j);
     // Get two planes tangent to surface z = xy, one at point a,  and another at point b.
     // One can show that these planes still create a cut.
     void get_initial_tang_points(point &a, point &b, const point& xy, bool below) const;