diff --git a/src/nlsat/nlsat_explain.cpp b/src/nlsat/nlsat_explain.cpp
index 618394ac1..c7d072232 100644
--- a/src/nlsat/nlsat_explain.cpp
+++ b/src/nlsat/nlsat_explain.cpp
@@ -45,6 +45,8 @@ namespace nlsat {
         bool                    m_minimize_cores;
         bool                    m_factor;
         bool                    m_signed_project;
+        bool m_cell_sample;
+
 
         struct todo_set {
             polynomial::cache  &    m_cache;
@@ -130,7 +132,7 @@ namespace nlsat {
         evaluator &             m_evaluator;
 
         imp(solver & s, assignment const & x2v, polynomial::cache & u, atom_vector const & atoms, atom_vector const & x2eq,
-            evaluator & ev):
+            evaluator & ev, bool is_sample):
             m_solver(s),
             m_assignment(x2v),
             m_atoms(atoms),
@@ -148,6 +150,9 @@ namespace nlsat {
             m_core1(s),
             m_core2(s),
             m_result(nullptr),
+
+            m_cell_sample(is_sample),
+
             m_evaluator(ev) {
             m_simplify_cores   = false;
             m_full_dimensional = false;
@@ -646,6 +651,51 @@ namespace nlsat {
                 m_todo.insert(p);
             }
         }
+
+        void add_sample_coeff(polynomial_ref_vector &ps, var x){
+            polynomial_ref p(m_pm);
+            polynomial_ref lc(m_pm);
+            unsigned sz = ps.size();
+            for (unsigned i = 0; i < sz; i++){
+                p = ps.get(i);
+                unsigned k = degree(p, x);
+                SASSERT(k > 0);
+                TRACE("nlsat_explain", tout << "add_lc, x: "; display_var(tout, x); tout << "\nk: " << k << "\n"; display(tout, p); tout << "\n";);
+                for(; k > 0; k--){
+                    lc = m_pm.coeff(p, x, k);
+                    add_factors(lc);
+                    if (m_pm.nonzero_const_coeff(p, x, k)){
+                        TRACE("nlsat_explain", tout << "constant coefficient, skipping...\n";);
+                        break;
+                    }
+                }
+                SASSERT(sign(lc) != 0);
+                SASSERT(!is_const(lc));
+            }
+        }
+
+        void psc_resultant_sample(polynomial_ref_vector &ps, var x, polynomial_ref_vector & samples){
+            polynomial_ref p(m_pm);
+            polynomial_ref q(m_pm);
+
+
+            // polynomial_ref_vector samples(m_pm);
+            // samples.reset();
+            // sample_poly(ps, x, samples);
+            
+            SASSERT(samples.size() <= 2);
+
+            for (unsigned i = 0; i < ps.size(); i++){
+                p = ps.get(i);
+                for (unsigned j = 0; j < samples.size(); j++){
+                    q = samples.get(j);
+                    if (!m_pm.eq(p, q)) {
+                        psc(p, q, x);
+                    }
+                }
+            }
+        }
+        
         
         /**
            \brief Add leading coefficients of the polynomials in ps.
@@ -973,6 +1023,116 @@ namespace nlsat {
             add_simple_assumption(k, p, lsign);
         }
 
+        void cac_add_cell_lits(polynomial_ref_vector & ps, var y, polynomial_ref_vector & res) {
+            res.reset();
+            SASSERT(m_assignment.is_assigned(y));
+            bool lower_inf = true;
+            bool upper_inf = true;
+            scoped_anum_vector & roots = m_roots_tmp;
+            scoped_anum lower(m_am);
+            scoped_anum upper(m_am);
+            anum const & y_val = m_assignment.value(y);
+            TRACE("nlsat_explain", tout << "adding literals for "; display_var(tout, y); tout << " -> ";
+                  m_am.display_decimal(tout, y_val); tout << "\n";);
+            polynomial_ref p_lower(m_pm);
+            unsigned i_lower = UINT_MAX;
+            polynomial_ref p_upper(m_pm);
+            unsigned i_upper = UINT_MAX;
+            polynomial_ref p(m_pm);
+            unsigned sz = ps.size();
+            for (unsigned k = 0; k < sz; k++) {
+                p = ps.get(k);
+                if (max_var(p) != y)
+                    continue;
+                roots.reset();
+                // Variable y is assigned in m_assignment. We must temporarily unassign it.
+                // Otherwise, the isolate_roots procedure will assume p is a constant polynomial.
+                m_am.isolate_roots(p, undef_var_assignment(m_assignment, y), roots);
+                unsigned num_roots = roots.size();
+//#linxi begin add_cell_lits faster
+                bool all_lt = true;
+                for (unsigned i = 0; i < num_roots; i++) {
+                    int s = m_am.compare(y_val, roots[i]);
+                    TRACE("nlsat_explain", 
+                        m_am.display_decimal(tout << "comparing root: ", roots[i]); tout << "\n";
+                        m_am.display_decimal(tout << "with y_val:", y_val); 
+                        tout << "\nsign " << s << "\n";
+                        tout << "poly: " << p << "\n";
+                        );
+                    if (s == 0) {
+                        // y_val == roots[i]
+                        // add literal
+                        // ! (y = root_i(p))
+                        add_root_literal(atom::ROOT_EQ, y, i+1, p);
+                        res.push_back(p);
+                        return;
+                    }
+                    else if (s < 0) {
+                        // y_val < roots[i]
+                        if (i > 0) {
+                            // y_val > roots[j]
+                            int j = i - 1;
+                            if (lower_inf || m_am.lt(lower, roots[j])) {
+                                lower_inf = false;
+                                m_am.set(lower, roots[j]);
+                                p_lower = p;
+                                i_lower = j + 1;
+                            }
+                        }
+                        if (upper_inf || m_am.lt(roots[i], upper)) {
+                            upper_inf = false;
+                            m_am.set(upper, roots[i]);
+                            p_upper = p;
+                            i_upper = i + 1;
+                        }
+                        all_lt = false;
+                        break;
+                    }
+                    // else if (s < 0) {
+                    //     // y_val < roots[i]
+                        
+                    //     // check if roots[i] is a better upper bound
+                    //     if (upper_inf || m_am.lt(roots[i], upper)) {
+                    //         upper_inf = false;
+                    //         m_am.set(upper, roots[i]);
+                    //         p_upper = p;
+                    //         i_upper = i+1;
+                    //     }
+                    // }
+                    // else if (s > 0) {
+                    //     // roots[i] < y_val
+
+                    //     // check if roots[i] is a better lower bound
+                    //     if (lower_inf || m_am.lt(lower, roots[i])) {
+                    //         lower_inf = false;
+                    //         m_am.set(lower, roots[i]);
+                    //         p_lower = p;
+                    //         i_lower = i+1;
+                    //     }
+                    // }
+                }
+                if (all_lt && num_roots > 0) {
+                    int j = num_roots - 1;
+                    if (lower_inf || m_am.lt(lower, roots[j])) {
+                        lower_inf = false;
+                        m_am.set(lower, roots[j]);
+                        p_lower = p;
+                        i_lower = j + 1;
+                    }
+                }
+//#linxi end add_cell_lits faster
+            }
+
+            if (!lower_inf) {
+                res.push_back(p_lower);
+                add_root_literal(m_full_dimensional ? atom::ROOT_GE : atom::ROOT_GT, y, i_lower, p_lower);
+            }
+            if (!upper_inf) {
+                res.push_back(p_upper);
+                add_root_literal(m_full_dimensional ? atom::ROOT_LE : atom::ROOT_LT, y, i_upper, p_upper);
+            }
+        }
+
         /**
            Add one or two literals that specify in which cell of variable y the current interpretation is.
            One literal is added for the cases: 
@@ -1016,6 +1176,8 @@ namespace nlsat {
                 // Otherwise, the isolate_roots procedure will assume p is a constant polynomial.
                 m_am.isolate_roots(p, undef_var_assignment(m_assignment, y), roots);
                 unsigned num_roots = roots.size();
+//#linxi begin add_cell_lits faster
+                bool all_lt = true;
                 for (unsigned i = 0; i < num_roots; i++) {
                     int s = m_am.compare(y_val, roots[i]);
                     TRACE("nlsat_explain", 
@@ -1033,27 +1195,58 @@ namespace nlsat {
                     }
                     else if (s < 0) {
                         // y_val < roots[i]
-                        
-                        // check if roots[i] is a better upper bound
+                        if (i > 0) {
+                            // y_val > roots[j]
+                            int j = i - 1;
+                            if (lower_inf || m_am.lt(lower, roots[j])) {
+                                lower_inf = false;
+                                m_am.set(lower, roots[j]);
+                                p_lower = p;
+                                i_lower = j + 1;
+                            }
+                        }
                         if (upper_inf || m_am.lt(roots[i], upper)) {
                             upper_inf = false;
                             m_am.set(upper, roots[i]);
                             p_upper = p;
-                            i_upper = i+1;
+                            i_upper = i + 1;
                         }
+                        all_lt = false;
+                        break;
                     }
-                    else if (s > 0) {
-                        // roots[i] < y_val
+                    // else if (s < 0) {
+                    //     // y_val < roots[i]
+                        
+                    //     // check if roots[i] is a better upper bound
+                    //     if (upper_inf || m_am.lt(roots[i], upper)) {
+                    //         upper_inf = false;
+                    //         m_am.set(upper, roots[i]);
+                    //         p_upper = p;
+                    //         i_upper = i+1;
+                    //     }
+                    // }
+                    // else if (s > 0) {
+                    //     // roots[i] < y_val
 
-                        // check if roots[i] is a better lower bound
-                        if (lower_inf || m_am.lt(lower, roots[i])) {
-                            lower_inf = false;
-                            m_am.set(lower, roots[i]);
-                            p_lower = p;
-                            i_lower = i+1;
-                        }
+                    //     // check if roots[i] is a better lower bound
+                    //     if (lower_inf || m_am.lt(lower, roots[i])) {
+                    //         lower_inf = false;
+                    //         m_am.set(lower, roots[i]);
+                    //         p_lower = p;
+                    //         i_lower = i+1;
+                    //     }
+                    // }
+                }
+                if (all_lt && num_roots > 0) {
+                    int j = num_roots - 1;
+                    if (lower_inf || m_am.lt(lower, roots[j])) {
+                        lower_inf = false;
+                        m_am.set(lower, roots[j]);
+                        p_lower = p;
+                        i_lower = j + 1;
                     }
                 }
+//#linxi end add_cell_lits faster
             }
 
             if (!lower_inf) {
@@ -1064,6 +1257,7 @@ namespace nlsat {
             }
         }
 
+
         /**
            \brief Return true if all polynomials in ps are univariate in x.
         */
@@ -1082,7 +1276,8 @@ namespace nlsat {
         /**
            \brief Apply model-based projection operation defined in our paper.
         */
-        void project(polynomial_ref_vector & ps, var max_x) {
+
+        void project_original(polynomial_ref_vector & ps, var max_x) {
             if (ps.empty())
                 return;
             m_todo.reset();
@@ -1094,12 +1289,12 @@ namespace nlsat {
             if (x < max_x)
                 add_cell_lits(ps, x);
             while (true) {
-                TRACE("nlsat_explain", tout << "project loop, processing var "; display_var(tout, x); tout << "\npolynomials\n";
-                      display(tout, ps); tout << "\n";);
                 if (all_univ(ps, x) && m_todo.empty()) {
                     m_todo.reset();
                     break;
                 }
+                TRACE("nlsat_explain", tout << "project loop, processing var "; display_var(tout, x); tout << "\npolynomials\n";
+                      display(tout, ps); tout << "\n";);
                 add_lc(ps, x);
                 psc_discriminant(ps, x);
                 psc_resultant(ps, x);
@@ -1109,6 +1304,72 @@ namespace nlsat {
                 add_cell_lits(ps, x);
             }
         }
+        void project_cdcac(polynomial_ref_vector & ps, var max_x) {
+            // whz
+            bool first = true;
+            
+
+            if (ps.empty())
+                return;
+            m_todo.reset();
+            for (poly* p : ps) {
+                m_todo.insert(p);
+            }
+            var x = m_todo.remove_max_polys(ps);
+            // Remark: after vanishing coefficients are eliminated, ps may not contain max_x anymore
+            
+            polynomial_ref_vector samples(m_pm);
+
+            
+            if (x < max_x){
+                cac_add_cell_lits(ps, x, samples);
+            }
+
+            while (true) {
+                if (all_univ(ps, x) && m_todo.empty()) {
+                    m_todo.reset();
+                    break;
+                }
+                TRACE("nlsat_explain", tout << "project loop, processing var "; display_var(tout, x); tout << "\npolynomials\n";
+                      display(tout, ps); tout << "\n";);
+
+
+                /**
+                 * Sample Projection
+                 * Reference:
+                 * Haokun Li and Bican Xia. 
+                 * "Solving Satisfiability of Polynomial Formulas By Sample - Cell Projection"
+                 * https://arxiv.org/abs/2003.00409 
+                 */
+
+                if (first) {
+                    add_lc(ps, x);
+                    psc_discriminant(ps, x);
+                    psc_resultant(ps, x);
+                    first = false;
+                }
+
+                else {
+                    add_lc(ps, x);
+                    // add_sample_coeff(ps, x);
+                    psc_discriminant(ps, x);
+                    psc_resultant_sample(ps, x, samples);
+                }
+                
+                if (m_todo.empty())
+                    break;
+                x = m_todo.remove_max_polys(ps);
+                cac_add_cell_lits(ps, x, samples);
+            }
+        }
+        void project(polynomial_ref_vector & ps, var max_x) {
+            if (m_cell_sample) {
+                project_cdcac(ps, max_x);
+            }
+            else {
+                project_original(ps, max_x);
+            }
+        }
 
         bool check_already_added() const {
             for (bool b : m_already_added_literal) {
@@ -1474,7 +1735,7 @@ namespace nlsat {
             var max_x = max_var(m_ps);
             TRACE("nlsat_explain", tout << "polynomials in the conflict:\n"; display(tout, m_ps); tout << "\n";);
             elim_vanishing(m_ps);
-            TRACE("nlsat_explain", tout << "elim vanishing x" << max_x << "\n"; display(tout, m_ps); tout << "\n";);
+            TRACE("nlsat_explain", tout << "elim vanishing\n"; display(tout, m_ps); tout << "\n";);
             project(m_ps, max_x);
             TRACE("nlsat_explain", tout << "after projection\n"; display(tout, m_ps); tout << "\n";);
         }
@@ -1936,8 +2197,8 @@ namespace nlsat {
     };
 
     explain::explain(solver & s, assignment const & x2v, polynomial::cache & u, 
-                     atom_vector const& atoms, atom_vector const& x2eq, evaluator & ev) {
-        m_imp = alloc(imp, s, x2v, u, atoms, x2eq, ev);
+                     atom_vector const& atoms, atom_vector const& x2eq, evaluator & ev, bool is_sample) {
+        m_imp = alloc(imp, s, x2v, u, atoms, x2eq, ev, is_sample);
     }
 
     explain::~explain() {
diff --git a/src/nlsat/nlsat_explain.h b/src/nlsat/nlsat_explain.h
index 20322a680..a19fe3982 100644
--- a/src/nlsat/nlsat_explain.h
+++ b/src/nlsat/nlsat_explain.h
@@ -34,7 +34,8 @@ namespace nlsat {
         imp * m_imp;
     public:
         explain(solver & s, assignment const & x2v, polynomial::cache & u, 
-                atom_vector const& atoms, atom_vector const& x2eq, evaluator & ev);
+                atom_vector const& atoms, atom_vector const& x2eq, evaluator & ev, bool is_sample);
+
         ~explain();
 
         void reset();
diff --git a/src/nlsat/nlsat_params.pyg b/src/nlsat/nlsat_params.pyg
index cac33fe87..dec3fac94 100644
--- a/src/nlsat/nlsat_params.pyg
+++ b/src/nlsat/nlsat_params.pyg
@@ -14,6 +14,7 @@ def_module_params('nlsat',
                           ('shuffle_vars', BOOL, False, "use a random variable order."),
                           ('inline_vars', BOOL, False, "inline variables that can be isolated from equations (not supported in incremental mode)"),
                           ('seed', UINT, 0, "random seed."),
-                          ('factor', BOOL, True, "factor polynomials produced during conflict resolution.")     
+                          ('factor', BOOL, True, "factor polynomials produced during conflict resolution."),
+                          ('cell_sample', BOOL, True, "cell sample projection"),                          
                           ))         
                 
diff --git a/src/nlsat/nlsat_solver.cpp b/src/nlsat/nlsat_solver.cpp
index bb169a05e..3d38d7a5e 100644
--- a/src/nlsat/nlsat_solver.cpp
+++ b/src/nlsat/nlsat_solver.cpp
@@ -253,7 +253,7 @@ namespace nlsat {
             m_simplify(s, m_atoms, m_clauses, m_learned, m_pm),
             m_display_var(m_perm),
             m_display_assumption(nullptr),
-            m_explain(s, m_assignment, m_cache, m_atoms, m_var2eq, m_evaluator),
+            m_explain(s, m_assignment, m_cache, m_atoms, m_var2eq, m_evaluator, nlsat_params(c.m_params).cell_sample()),
             m_scope_lvl(0),
             m_lemma(s),
             m_lazy_clause(s),