mul-saturation wip

fixup conflict explanations in mul_saturation, add parameter to enable it, add statistics

src/math/lp/lp_settings.cpp

@@ -34,6 +34,7 @@ void lp::lp_settings::updt_params(params_ref const& _p) {
     report_frequency = p.arith_rep_freq();
     m_simplex_strategy = static_cast<lp::simplex_strategy_enum>(p.arith_simplex_strategy());
     m_nlsat_delay = p.arith_nl_delay();
+    m_enable_stellensatz = p.arith_nl_stellensatz();
     auto eps = p.arith_epsilon();
     m_epsilon = rational(std::max(1, (int)(100000*eps)), 100000); 
     m_dio = lp_p.dio();

src/math/lp/lp_settings.h

@@ -137,6 +137,7 @@ struct statistics {
     unsigned m_bounds_tightening_conflicts = 0;
     unsigned m_bounds_tightenings = 0;
     unsigned m_nla_throttled_lemmas = 0;
+    unsigned m_nla_stellensatz = 0;
 
     ::statistics m_st = {};
 
@@ -176,6 +177,7 @@ struct statistics {
         st.update("arith-bounds-tightening-conflicts", m_bounds_tightening_conflicts);
         st.update("arith-bounds-tightenings", m_bounds_tightenings);
         st.update("arith-nla-throttled-lemmas", m_nla_throttled_lemmas);
+        st.update("arith-nla-stellensatz", m_nla_stellensatz);
         st.copy(m_st);
     }
 };
@@ -260,6 +262,7 @@ private:
     bool             m_dio_run_gcd = true;
 public:
     bool m_enable_relevancy = false;
+    bool m_enable_stellensatz = true;
     unsigned dio_calls_period() const { return m_dio_calls_period; }
     unsigned & dio_calls_period() { return m_dio_calls_period; }
     bool print_external_var_name() const { return m_print_external_var_name; }

src/math/lp/lp_types.h

@@ -40,6 +40,7 @@ namespace lp {
         EQ = 0,
         NE = 3
     };
+
     typedef unsigned lpvar;
     const lpvar null_lpvar = UINT_MAX;
     const constraint_index null_ci = UINT_MAX;

src/math/lp/nla_core.cpp

@@ -1333,7 +1333,7 @@ lbool core::check() {
             return l_false;
     }
 
-    if (false && no_effect())
+    if (no_effect() && lp_settings().m_enable_stellensatz)
         ret = m_mul_saturate.saturate();
     
     if (no_effect() && should_run_bounded_nlsat()) 

src/math/lp/nla_mul_saturate.cpp

@@ -1,19 +1,30 @@
 /*++
   Copyright (c) 2025 Microsoft Corporation
 
-  given a monic m = x * y * z ... with evaluation val(m) != val(x) * val(y) * val(z) ...
+  given a monic m = x * y * z ... used in a constraint that is false under the current evaluation of x,y,z
 
   saturate constraints with respect to m
   in other words, if a constraint contains x*y + p >= 0, 
-  then include the constraint z >= 0 => x*y*z + z*p >= 0
+  then include the constraint x*y*z + z*p >= 0
   assuming current value of z is non-negative.
-  Check if the system with new constraints is LP feasible.
+  Check if the system with new constraints is LP (and MIP) feasible.
-  If it is not, then produce a lemma that explains the infeasibility.
+  If it is not, then produce a lemma that explains infeasibility.
 
   Strategy 1: The lemma is in terms of the original constraints and bounds.
   Strategy 2: Attempt to eliminate new monomials from the lemma by relying on Farkas multipliers.
               If it succeeds to eliminate new monomials we have a lemma that is a linear 
               combination of existing variables. 
+              The idea is that a conflict over the new system is a set of multipliers lambda, such
+              that lambda * A is separated from b for the constraints Axy >= b.
+              The coefficients in lambda are non-negative.
+              They correspond to variables x, y where x are variables from the input constraints
+              and y are new variables introduced for new monomials.
+              We can test if lambda allows eliminating y by taking a subset of lambda indices where
+              A has rows containing y_i for some fresh y_i. Replace those rows r_j by the partial sum of
+              of rows multiplied by lambda_j. The sum r_j1 * lambda_j1 + .. + r_jk * lambda_jk does not
+              contain y_i. Repeat the same process with other variables y_i'. If a sum is
+              generated without any y, it is a linear consequence of the new constraints but not
+              necessarily derivable with the old constraints.
   Strategy 3: The lemma uses the new constraints.
 
   --*/
@@ -38,7 +49,8 @@ namespace nla {
     }
 
     void mul_saturate::init_solver() {
-        local_solver = alloc(lp::lar_solver);
+        lra_solver = alloc(lp::lar_solver);
+        int_solver = alloc(lp::int_solver, *lra_solver);
         m_vars2mon.reset();
         m_mon2vars.reset();
         m_values.reset();
@@ -50,6 +62,7 @@ namespace nla {
         auto const& lra = c().lra_solver();
         auto sz = lra.number_of_vars();
         for (unsigned v = 0; v < sz; ++v) {
+            // Declare v into lra_solver
             unsigned w;
             if (m_coi.mons().contains(v))
                 init_monomial(v);
@@ -59,24 +72,28 @@ namespace nla {
                 auto const& t = lra.get_term(v);
                 // Assumption: variables in coefficients are always declared before term variable.
                 SASSERT(all_of(t, [&](auto p) { return p.j() < v; }));
-                w = local_solver->add_term(t.coeffs_as_vector(), v);
+                w = lra_solver->add_term(t.coeffs_as_vector(), v);
             } 
             else
-                w = local_solver->add_var(v, lra.var_is_int(v));
+                w = lra_solver->add_var(v, lra.var_is_int(v));
 
+            // assert bounds on v in the new solver.
             VERIFY(w == v);
             if (lra.column_has_lower_bound(v)) {
-                auto lo_dep = lra.get_column_lower_bound_witness(v);
                 auto lo_bound = lra.get_lower_bound(v);
                 auto k = lo_bound.y > 0 ? lp::lconstraint_kind::GT : lp::lconstraint_kind::GE;
-                auto ci = local_solver->add_var_bound(v, k, lo_bound.x);
+                auto rhs = lo_bound.x;
+                auto dep = lra.get_column_lower_bound_witness(v);
+                auto ci = lra_solver->add_var_bound(v, k, rhs);
+                m_ci2dep.setx(ci, dep, nullptr);
             }
             if (lra.column_has_upper_bound(v)) {
-                auto hi_dep = lra.get_column_upper_bound_witness(v);
                 auto hi_bound = lra.get_upper_bound(v);
                 auto k = hi_bound.y < 0 ? lp::lconstraint_kind::LT : lp::lconstraint_kind::LE;
-                auto ci = local_solver->add_var_bound(v, k, hi_bound.x);
+                auto rhs = hi_bound.x;
-                m_ci2dep.setx(ci, hi_dep, nullptr);
+                auto dep = lra.get_column_upper_bound_witness(v);
+                auto ci = lra_solver->add_var_bound(v, k, rhs);
+                m_ci2dep.setx(ci, dep, nullptr);
             }
         }
     }
@@ -96,81 +113,106 @@ namespace nla {
     void mul_saturate::add_lemma(lp::explanation const& ex1) {
         auto& lra = c().lra_solver();
         lp::explanation ex2;
-        for (auto p : ex1) {
-            lp::constraint_index src_ci;
-            if (!m_new_mul_constraints.find(p.ci(), src_ci))
-                src_ci = p.ci();
-            auto dep = m_ci2dep.get(src_ci, nullptr);
-            local_solver->push_explanation(dep, ex2);
-        }
-        for (auto [v, sign, dep] : m_var_signs) {            
-            if (!dep) {
-                verbose_stream() << "TODO explain non-implied bound\n";
-                continue;
-            }
-            local_solver->push_explanation(dep, ex2);
-        }
         lemma_builder new_lemma(c(), "stellensatz");
+        for (auto p : ex1) {
+            auto dep = m_ci2dep.get(p.ci(), nullptr);
+            lra_solver->push_explanation(dep, ex2);
+            if (!m_new_mul_constraints.contains(p.ci()))
+                continue;
+            
+            auto const& bounds = m_new_mul_constraints[p.ci()];
+            for (auto const& b : bounds) {
+                if (std::holds_alternative<u_dependency*>(b)) {
+                    auto dep = *std::get_if<u_dependency*>(&b);
+                    lra_solver->push_explanation(dep, ex2);
+                } 
+                else {
+                    auto const &[v, k, rhs] = *std::get_if<bound>(&b);    
+                    new_lemma |= ineq(v, negate(k), rhs);
+                }
+            }
+        }
+
         new_lemma &= ex2;
-        IF_VERBOSE(1, verbose_stream() << "unsat \n" << new_lemma << "\n");
+        IF_VERBOSE(5, verbose_stream() << "unsat \n" << new_lemma << "\n");
         TRACE(arith, tout << "unsat\n" << new_lemma << "\n");
+        c().lra_solver().settings().stats().m_nla_stellensatz++;
     }
 
+
     lbool mul_saturate::solve(lp::explanation& ex) {
-        for (auto const& [new_ci, old_ci] : m_new_mul_constraints)
+        lbool r = solve_lra(ex);
-            local_solver->activate(new_ci);
+        if (r != l_true)
-        auto st = local_solver->solve();
+            return r;
-        lbool r = l_undef;
+        r = solve_lia(ex);
+        if (r != l_true)
+            return r;
+        // if (r == l_true) check if solution satisfies constraints
+        // variables outside the slice have values from the outer solver.
+        return l_undef;
+    }
+
+    lbool mul_saturate::solve_lra(lp::explanation& ex) {         
+        auto st = lra_solver->solve();
         if (st == lp::lp_status::INFEASIBLE) {
-            local_solver->get_infeasibility_explanation(ex);
+            lra_solver->get_infeasibility_explanation(ex);
-            r = l_false;
+            return l_false;
         } 
-        if (st == lp::lp_status::OPTIMAL || st == lp::lp_status::FEASIBLE) {
+        else if (lra_solver->is_feasible())
-            // TODO: check model just in case it got lucky.
+            return l_true;
-            IF_VERBOSE(1, verbose_stream() << "saturation is LP feasible, maybe it is a model of the NLA problem\n");
+        else
+            return l_undef;
+    }
+
+    lbool mul_saturate::solve_lia(lp::explanation& ex) {
+        switch (int_solver->check(&ex)) {
+        case lp::lia_move::sat:
+            return l_true;
+        case lp::lia_move::conflict:
+            return l_false;
+        default: // TODO: an option is to perform (bounded) search here to get an LIA verdict.
+            return l_undef;
         }
-        IF_VERBOSE(0, display(verbose_stream()));
+        return l_undef;
-        return r;
     }
 
     // record new monomials that are created and recursively down-saturate with respect to these.
     // this is a simplistic pass
     void mul_saturate::add_multiply_constraints() {
         m_new_mul_constraints.reset();
-        m_seen_vars.reset();
-        m_var_signs.reset();
         m_to_refine.reset();
         vector<svector<lp::constraint_index>> var2cs; 
 
-        for (auto ci : local_solver->constraints().indices()) {
+        // current approach: only resolve against var2cs, which is initialized
-            auto const& con = local_solver->constraints()[ci];
+        // with monomials in the input.
+
+        for (auto ci : lra_solver->constraints().indices()) {
+            auto const& con = lra_solver->constraints()[ci];
             for (auto [coeff, v] : con.coeffs()) {
                 if (v >= var2cs.size())
                     var2cs.resize(v + 1);
                 var2cs[v].push_back(ci);
             }
-
             // insert monomials to be refined
             insert_monomials_from_constraint(ci);            
         }
 
         for (unsigned it = 0; it < m_to_refine.size(); ++it) {
             auto j = m_to_refine[it];
-            verbose_stream() << "refining " << j << " := " << m_mon2vars[j] << "\n";
             auto vars = m_mon2vars[j];
             for (auto v : vars) {
                 if (v >= var2cs.size())
                     continue;
                 auto cs = var2cs[v];
                 for (auto ci : cs) {
-                    for (auto [coeff, u] : local_solver->constraints()[ci].coeffs()) {
+                    for (auto [coeff, u] : lra_solver->constraints()[ci].coeffs()) {
                         if (u == v)
                             add_multiply_constraint(ci, j, v);
                     }                    
                 }
             }
         }
-        IF_VERBOSE(0,
+        IF_VERBOSE(5,
                    c().lra_solver().display(verbose_stream() << "original\n");
                    c().display(verbose_stream());
                    display(verbose_stream() << "saturated\n"));
@@ -178,52 +220,60 @@ namespace nla {
 
     // multiply by remaining vars
     void mul_saturate::add_multiply_constraint(lp::constraint_index old_ci, lp::lpvar mi, lpvar x) {
-        lp::lar_base_constraint const& con = local_solver->constraints()[old_ci];
+        lp::lar_base_constraint const& con = lra_solver->constraints()[old_ci];
         auto &lra = c().lra_solver();
         auto const& lhs = con.coeffs();
         auto const& rhs = con.rhs();
         auto k = con.kind();
         if (k == lp::lconstraint_kind::NE || k == lp::lconstraint_kind::EQ)
             return;  // not supported
-        auto sign = false;
+        auto sign = false, strict = true;
         svector<lpvar> vars;
         bool first = true;
-        for (auto v : c().emon(mi).vars()) {
+        for (auto v : m_mon2vars[mi]) {
             if (v != x || !first)
                 vars.push_back(v);
             else
                 first = false;
         }
-        // compute sign of vars
+        SASSERT(!first); // v was a member and was removed
+
+        vector<bound_justification> bounds;
+        // compute bounds constraints and sign of vars
+
+        auto add_bound = [&](lpvar v, lp::lconstraint_kind k, int r) { 
+            bound b(v, k, rational(r));
+            bounds.push_back(b);
+        };
         for (auto v : vars) {
-            if (m_seen_vars.contains(v))
-                continue;
-            m_seen_vars.insert(v);
             // retrieve bounds of v
             // if v has positive lower bound add as positive
             // if v has negative upper bound add as negative
-            // otherwise, soft-fail (for now unsound)
+            // otherwise look at the current value of v and add bounds assumption based on current sign.
-            // proper signs of variables from old tableau should be extracted using lra_solver()
-            // instead of local_solver.
-            // TODO is to also add case where lower or upper bound is zero and then the sign
-            // of the inequality is relaxed if it is strict.
             if (lra.number_of_vars() > v && lra.column_has_lower_bound(v) && lra.get_lower_bound(v).is_pos()) {
-                m_var_signs.push_back({v, false, lra.get_column_lower_bound_witness(v)});
+                bounds.push_back(lra.get_column_lower_bound_witness(v));
             } 
             else if (lra.number_of_vars() > v && lra.column_has_upper_bound(v) && lra.get_upper_bound(v).is_neg()) {
-                m_var_signs.push_back({v, true, lra.get_column_upper_bound_witness(v)});
+                bounds.push_back(lra.get_column_upper_bound_witness(v));                
                 sign = !sign;
             } 
             else if (m_values[v].is_neg()) {
-                m_var_signs.push_back({v, true, nullptr});
+                if (lra.var_is_int(v))
+                    add_bound(v, lp::lconstraint_kind::LE, -1);
+                else
+                    add_bound(v, lp::lconstraint_kind::LT, 0);
                 sign = !sign;
             } 
             else if (m_values[v].is_pos()) {
-                m_var_signs.push_back({v, false, nullptr});
+                if (lra.var_is_int(v)) 
+                    add_bound(v, lp::lconstraint_kind::GE, 1);                
+                else 
+                    add_bound(v, lp::lconstraint_kind::GT, 0);                
             } 
             else {
-                IF_VERBOSE(0, verbose_stream() << "not separated from 0\n");
+                SASSERT(m_values[v] == 0);
-                return;
+                strict = false;
+                add_bound(v, lp::lconstraint_kind::GE, 0);                
             }
         }
         lp::lar_term new_lhs;
@@ -246,34 +296,35 @@ namespace nla {
             new_rhs = 0;
         }
 
-        if (sign) {
+        if (sign) 
+            k = swap_side(k);
+        
+        if (!strict) {
             switch (k) {
-            case lp::lconstraint_kind::LE:
+            case lp::lconstraint_kind::GT:
                 k = lp::lconstraint_kind::GE;
                 break;
             case lp::lconstraint_kind::LT:
-                k = lp::lconstraint_kind::GT;
-                break;
-            case lp::lconstraint_kind::GE:
                 k = lp::lconstraint_kind::LE;
                 break;
-            case lp::lconstraint_kind::GT:
-                k = lp::lconstraint_kind::LT;
-                break;
-            default:
-                break;
             }
         }
-        display_constraint(verbose_stream(), old_ci) << " -> ";
+
-        display_constraint(verbose_stream(), new_lhs.coeffs_as_vector(), k, new_rhs) << "\n";
+        auto ti = lra_solver->add_term(new_lhs.coeffs_as_vector(), lra_solver->number_of_vars());
-        auto ti = local_solver->add_term(new_lhs.coeffs_as_vector(), local_solver->number_of_vars());
+        auto new_ci = lra_solver->add_var_bound(ti, k, new_rhs);
-        auto new_ci = local_solver->add_var_bound(ti, k, new_rhs);
+        IF_VERBOSE(4, display_constraint(verbose_stream(), old_ci) << " -> ";
+                   display_constraint(verbose_stream(), new_lhs.coeffs_as_vector(), k, new_rhs) << "\n");
         insert_monomials_from_constraint(new_ci);
         m_values.push_back(term_value);
         SASSERT(m_values.size() - 1 == ti);
-        m_new_mul_constraints.insert(new_ci, old_ci);
+        m_new_mul_constraints.insert(new_ci, bounds);
+        m_ci2dep.setx(new_ci, m_ci2dep.get(old_ci, nullptr), nullptr);
     }
 
+    // Insert a (new) monomial representing product of vars.
+    // if the length of vars is 1 then the new monomial is vars[0]
+    // if the same monomial was previous defined, return the previously defined monomial.
+    // otherwise create a fresh variable representing the monomial.
     lpvar mul_saturate::add_monomial(svector<lpvar> const& vars) {
         lpvar v;
         if (vars.size() == 1)
@@ -300,23 +351,25 @@ namespace nla {
     }
 
     lpvar mul_saturate::add_var(bool is_int) {
-        auto v = local_solver->number_of_vars();
+        auto v = lra_solver->number_of_vars();
-        auto w = local_solver->add_var(v, is_int);
+        auto w = lra_solver->add_var(v, is_int);
-        VERIFY(v == w);
+        SASSERT(v == w);
         return w;
     }
 
+    // if a constraint is false, then insert _all_ monomials from that constraint.
+    // other approach: use a lex ordering on monomials and insert the lex leading monomial.   
     void mul_saturate::insert_monomials_from_constraint(lp::constraint_index ci) {
         if (constraint_is_true(ci))
             return;
-        auto const& con = local_solver->constraints()[ci];
+        auto const& con = lra_solver->constraints()[ci];
         for (auto [coeff, v] : con.coeffs()) 
             if (c().is_monic_var(v))
                 m_to_refine.insert(v);        
     }
 
     bool mul_saturate::constraint_is_true(lp::constraint_index ci) {
-        auto const& con = local_solver->constraints()[ci];
+        auto const& con = lra_solver->constraints()[ci];
         auto lhs = -con.rhs();
         for (auto const& [coeff, v] : con.coeffs())
             lhs += coeff * m_values[v];
@@ -340,7 +393,7 @@ namespace nla {
     }
 
     std::ostream& mul_saturate::display(std::ostream& out) const {
-        local_solver->display(out);
+        lra_solver->display(out);
         for (auto const& [vars, v] : m_vars2mon) {
             out << "j" << v << " := ";
             display_product(out, vars);
@@ -362,7 +415,7 @@ namespace nla {
     }
 
     std::ostream& mul_saturate::display_constraint(std::ostream& out, lp::constraint_index ci) const {
-        auto const& con = local_solver->constraints()[ci];
+        auto const& con = lra_solver->constraints()[ci];
         return display_constraint(out, con.coeffs(), con.kind(), con.rhs());
     }
 

src/math/lp/nla_mul_saturate.h

@@ -5,6 +5,8 @@
 #pragma once
 
 #include "math/lp/nla_coi.h"
+#include "math/lp/int_solver.h"
+#include <variant>
 
 namespace nla {
 
@@ -12,18 +14,19 @@ namespace nla {
     class lar_solver;    
     class mul_saturate : common {
 
-        struct var_sign {
+
+        struct bound {
             lpvar v = lp::null_lpvar;
-            bool is_neg = false;
+            lp::lconstraint_kind k;
-            u_dependency* dep = nullptr;
+            rational rhs;           
         };
+        using bound_justification = std::variant<u_dependency*, bound>;
+
         coi m_coi;
-        // source of multiplication constraint
+        u_map<vector<bound_justification>> m_new_mul_constraints;
-        u_map<lp::constraint_index> m_new_mul_constraints;
-        svector<var_sign> m_var_signs;
-        tracked_uint_set m_seen_vars;
         indexed_uint_set m_to_refine;
-        scoped_ptr<lp::lar_solver> local_solver;
+        scoped_ptr<lp::lar_solver> lra_solver;
+        scoped_ptr<lp::int_solver> int_solver;
         ptr_vector<u_dependency> m_ci2dep;
         vector<rational> m_values;
         struct eq {
@@ -51,6 +54,8 @@ namespace nla {
 
         // solving
         lbool solve(lp::explanation& ex);
+        lbool solve_lra(lp::explanation &ex);
+        lbool solve_lia(lp::explanation &ex);
 
         // lemmas
         void add_lemma(lp::explanation const& ex);

src/params/smt_params_helper.pyg

@@ -88,6 +88,7 @@ def_module_params(module_name='smt',
                           ('arith.nl.propagate_linear_monomials', BOOL, True, 'propagate linear monomials'),
                           ('arith.nl.optimize_bounds', BOOL, True, 'enable bounds optimization'),
                           ('arith.nl.cross_nested', BOOL, True, 'enable cross-nested consistency checking'),
+                          ('arith.nl.stellensatz', BOOL, False, 'enable stellensatz saturation'),
                           ('arith.nl.log', BOOL, False, 'Log lemmas sent to nra solver'),
                           ('arith.propagate_eqs', BOOL, True, 'propagate (cheap) equalities'),
                           ('arith.epsilon', DOUBLE, 1.0, 'initial value of epsilon used for model generation of infinitesimals'),

src/smt/theory_lra.cpp

@@ -3467,15 +3467,15 @@ public:
 
     void set_conflict_or_lemma(literal_vector const& core, bool is_conflict) {
         reset_evidence();
-        for (literal lit : core) {
+        for (literal lit : core) 
             m_core.push_back(lit);
-        }
+        
         // lp().shrink_explanation_to_minimum(m_explanation); // todo, enable when perf is fixed
         ++m_num_conflicts;
         ++m_stats.m_conflicts;
         for (auto ev : m_explanation)
             set_evidence(ev.ci(), m_core, m_eqs);
-        if (m_eqs.empty() && all_of(m_core, [&](literal l) { return ctx().get_assignment(l) == l_false; }))
+        if (all_of(m_core, [&](literal l) { return ctx().get_assignment(l) == l_false; }))
             is_conflict = true;
         TRACE(arith_conflict,
               tout << "@" << ctx().get_scope_level() << (is_conflict ? " conflict":" lemma");