diff --git a/src/nlsat/CMakeLists.txt b/src/nlsat/CMakeLists.txt
index 34dfa755c..7faae8b19 100644
--- a/src/nlsat/CMakeLists.txt
+++ b/src/nlsat/CMakeLists.txt
@@ -1,6 +1,7 @@
 z3_add_component(nlsat
   SOURCES
     nlsat_clause.cpp
+    nlsat_common.cpp
     nlsat_evaluator.cpp
     nlsat_explain.cpp
     nlsat_interval_set.cpp
diff --git a/src/nlsat/levelwise.cpp b/src/nlsat/levelwise.cpp
index fef36d738..c11f18cf3 100644
--- a/src/nlsat/levelwise.cpp
+++ b/src/nlsat/levelwise.cpp
@@ -10,6 +10,8 @@
 #include "math/polynomial/algebraic_numbers.h"
 #include "nlsat_common.h"
 #include <queue>
+#include "math/polynomial/polynomial_cache.h"
+#include "math/polynomial/polynomial.h"
 
 namespace nlsat {
     
@@ -30,17 +32,22 @@ namespace nlsat {
     struct levelwise::impl {
         // Utility: call fn for each distinct irreducible factor of poly
         template<typename Func>
-        void for_each_distinct_factor(const polynomial::polynomial_ref& poly, Func&& fn) {
-            polynomial::factors factors(m_pm);
-            factor(poly, factors);
-            for (unsigned i = 0; i < factors.distinct_factors(); ++i) 
-                fn(factors[i]);            
+        void for_each_distinct_factor(polynomial::polynomial_ref& poly, Func&& fn) {
+            polynomial::polynomial_ref_vector factors(m_pm);
+            ::nlsat::factor(poly, m_cache, factors);
+            for (unsigned i = 0; i < factors.size(); i++) {
+                polynomial_ref pr(m_pm);
+                pr = factors.get(i);
+                fn(pr);
+            }
         }
         template<typename Func>
-        void for_first_distinct_factor(const polynomial::polynomial_ref& poly, Func&& fn) {
-            polynomial::factors factors(m_pm);
-            factor(poly, factors);
-            fn(factors[0]);            
+        void for_first_distinct_factor(polynomial::polynomial_ref& poly, Func&& fn) {
+            polynomial::polynomial_ref_vector factors(m_pm);
+            ::nlsat::factor(poly, m_cache, factors);
+            polynomial_ref pr(m_pm);
+            pr = factors.get(0);
+            fn(pr);
         }
 
         // todo: consider to key polynomials in a set by using m_pm.eq
@@ -98,11 +105,11 @@ namespace nlsat {
         std::vector<root_function>   m_E; // the ordered root functions on a level
         assignment const & sample() const { return m_solver.sample();}
         assignment & sample() { return m_solver.sample(); }
-
+        polynomial::cache &          m_cache; 
        // max_x plays the role of n in algorith 1 of the levelwise paper.
 
-        impl(solver& solver, polynomial_ref_vector const& ps, var max_x, assignment const& s, pmanager& pm, anum_manager& am)
-            : m_solver(solver), m_P(ps),  m_n(max_x),  m_pm(pm), m_am(am) {
+        impl(solver& solver, polynomial_ref_vector const& ps, var max_x, assignment const& s, pmanager& pm, anum_manager& am, polynomial::cache & cache)
+            : m_solver(solver), m_P(ps),  m_n(max_x),  m_pm(pm), m_am(am), m_cache(cache) {
             TRACE(lws, tout  << "m_n:" << m_n << "\n";);
             m_I.reserve(m_n); // cannot just resize bcs of the absence of the default constructor of root_function_interval
             for (unsigned i = 0; i < m_n; ++i)
@@ -136,9 +143,18 @@ namespace nlsat {
         }
 
         bool is_irreducible(poly* p) {
-            polynomial::factors factors(m_pm);
-            factor(polynomial_ref(p, m_pm), factors);
-            return factors.total_factors() == 1;            
+            polynomial_ref_vector factors(m_pm);
+            polynomial_ref pref(p, m_pm);
+            ::nlsat::factor(pref, m_cache, factors);
+            unsigned num_factors = factors.size();
+            CTRACE(lws, num_factors != 1, ::nlsat::display(tout, m_solver, p) << std::endl;
+                   tout << "{";
+                   tout << "num_factors:" << num_factors << "\n";
+                   ::nlsat::display(tout, m_solver, factors);
+                   tout << "}\n";
+                );
+            
+            return factors.size() == 1;            
         }
         
         /*
@@ -338,6 +354,7 @@ namespace nlsat {
              }
 
              collect_E(p_non_null);
+
              std::sort(m_E.begin(), m_E.end(), [&](root_function const& a, root_function const& b){
                  return m_am.lt(a.val, b.val);
              });
@@ -411,7 +428,6 @@ namespace nlsat {
                 for_each_distinct_factor(disc, [&](polynomial::polynomial_ref f) {
                     if (coeffs_are_zeroes_on_sample(f, m_pm, sample(), m_am)) {
                         m_fail = true; // ambiguous multiplicity -- not handled yet
-                        NOT_IMPLEMENTED_YET();
                         return;
                     }
                     unsigned lvl = max_var(f);
@@ -714,6 +730,7 @@ or
             */
                 mk_prop(sample_holds, level_t(m_level - 1)); 
                 mk_prop(repr, level_t(m_level - 1));
+                mk_prop(ir_ord, level_t(m_level));
                 mk_prop(an_del, p.poly);
             }
         }
@@ -863,8 +880,8 @@ or
         }
     };
 // constructor
-    levelwise::levelwise(nlsat::solver& solver, polynomial_ref_vector const& ps, var n, assignment const& s, pmanager& pm, anum_manager& am)
-        : m_impl(new impl(solver, ps, n, s, pm, am)) {}
+    levelwise::levelwise(nlsat::solver& solver, polynomial_ref_vector const& ps, var n, assignment const& s, pmanager& pm, anum_manager& am, polynomial::cache& cache)
+        : m_impl(new impl(solver, ps, n, s, pm, am, cache)) {}
 
     levelwise::~levelwise() { delete m_impl; }
 
diff --git a/src/nlsat/levelwise.h b/src/nlsat/levelwise.h
index 0cccb31fa..b9cb53da5 100644
--- a/src/nlsat/levelwise.h
+++ b/src/nlsat/levelwise.h
@@ -1,5 +1,6 @@
 #pragma once
 #include "nlsat_types.h"
+#include "math/polynomial/polynomial_cache.h"
 
 namespace nlsat {
 
@@ -39,7 +40,7 @@ namespace nlsat {
         impl* m_impl;
     public:
     // Construct with polynomials ps, maximal variable max_x, current sample s, polynomial manager pm, and algebraic-number manager am
-        levelwise(nlsat::solver& solver, polynomial_ref_vector const& ps, var max_x, assignment const& s, pmanager& pm, anum_manager& am);
+        levelwise(nlsat::solver& solver, polynomial_ref_vector const& ps, var max_x, assignment const& s, pmanager& pm, anum_manager& am, polynomial::cache & cache);
         ~levelwise();
 
         levelwise(levelwise const&) = delete;
diff --git a/src/nlsat/nlsat_common.h b/src/nlsat/nlsat_common.h
index 44ba3fc74..a0540e4d2 100644
--- a/src/nlsat/nlsat_common.h
+++ b/src/nlsat/nlsat_common.h
@@ -13,8 +13,11 @@
   --*/
 #pragma once
 
+#include "nlsat/nlsat_assignment.h"
 #include "nlsat/nlsat_solver.h"
 #include "nlsat/nlsat_scoped_literal_vector.h"
+#include "math/polynomial/polynomial_cache.h"
+
 namespace nlsat {
 
     inline std::ostream& display(std::ostream& out, pmanager& pm, polynomial_ref const& p, display_var_proc const& proc) {
@@ -119,5 +122,9 @@ namespace nlsat {
             out << '\n';
         }
     }
+        /**
+           \brief Wrapper for factorization
+        */
+    void factor(polynomial_ref & p, polynomial::cache& cache, polynomial_ref_vector & fs);
 
 } // namespace nlsat
diff --git a/src/nlsat/nlsat_explain.cpp b/src/nlsat/nlsat_explain.cpp
index 93193fb73..1af4bbbdc 100644
--- a/src/nlsat/nlsat_explain.cpp
+++ b/src/nlsat/nlsat_explain.cpp
@@ -188,16 +188,6 @@ namespace nlsat {
         }
 
         
-        /**
-           \brief Wrapper for factorization
-        */
-        void factor(polynomial_ref & p, polynomial_ref_vector & fs) {
-            // TODO: add params, caching
-            TRACE(nlsat_factor, tout << "factor\n" << p << "\n";);
-            fs.reset();
-            m_cache.factor(p.get(), fs);
-        }
-
         /**
            \brief Wrapper for psc chain computation
         */
@@ -257,7 +247,7 @@ namespace nlsat {
             // Then, I assert p_i1 * ... * p_im  != 0
             {
                 restore_factors _restore(m_factors, m_factors_save);
-                factor(p, m_factors);
+                factor(p, m_cache, m_factors);
                 unsigned num_factors = m_factors.size();
                 m_zero_fs.reset();
                 m_is_even.reset();
@@ -576,7 +566,7 @@ namespace nlsat {
                 return;
             if (m_factor) {
                 restore_factors _restore(m_factors, m_factors_save);
-                factor(p, m_factors);
+                factor(p, m_cache, m_factors);
                 TRACE(nlsat_explain, display(tout << "adding factors of\n", m_solver, p); tout << "\n" << m_factors << "\n";);
                 polynomial_ref f(m_pm);
                 for (unsigned i = 0; i < m_factors.size(); i++) {
@@ -667,7 +657,7 @@ namespace nlsat {
         // this function also explains the value 0, if met
         bool coeffs_are_zeroes(polynomial_ref &s) {
             restore_factors _restore(m_factors, m_factors_save);
-            factor(s, m_factors);
+            factor(s, m_cache, m_factors);
             unsigned num_factors = m_factors.size();
             m_zero_fs.reset();
             m_is_even.reset();
@@ -1179,8 +1169,8 @@ namespace nlsat {
             }
         }
 
-        bool levelwise_single_cell(polynomial_ref_vector & ps, var max_x) {
-            levelwise lws(m_solver, ps, max_x, sample(), m_pm, m_am);
+        bool levelwise_single_cell(polynomial_ref_vector & ps, var max_x, polynomial::cache & cache) {
+            levelwise lws(m_solver, ps, max_x, sample(), m_pm, m_am, cache);
             auto cell = lws.single_cell();
             if (lws.failed()) {
                 return false;
@@ -1216,9 +1206,6 @@ namespace nlsat {
          */
         void project_cdcac(polynomial_ref_vector & ps, var max_x) {
             TRACE(nlsat_explain, tout << "max_x:" << max_x << std::endl;);
-            if (max_x == 0) {
-                std::cout << "*";
-            }
             if (ps.empty()) {
                 TRACE(nlsat_explain, tout << "ps.empty\n";);
                 return;
@@ -1238,7 +1225,7 @@ namespace nlsat {
             var x = m_todo.extract_max_polys(ps);
 
             TRACE(nlsat_explain, tout << "m_solver.apply_levelwise():" << m_solver.apply_levelwise() << "\n";);
-            if (m_solver.apply_levelwise() && levelwise_single_cell(ps, max_x))
+            if (m_solver.apply_levelwise() && levelwise_single_cell(ps, max_x, m_cache))
                 return;
 
             polynomial_ref_vector samples(m_pm);