add toggle to use polynomial translation

Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>

src/math/lp/nla_stellensatz.cpp

@@ -224,6 +224,7 @@ namespace nla {
         if (m_active.contains(ci))
             return;
         m_active.insert(ci);
+        m_tabu.reserve(ci + 1);
         m_tabu[ci] = tabu;
     }
 
@@ -544,7 +545,7 @@ namespace nla {
 
     lbool stellensatz::factor(lp::constraint_index ci) {
         auto const &[p, k, j] = m_constraints[ci];
-        auto [vars, q] = p.var_factors();
+        auto [vars, q] = p.var_factors(); // p = vars * q
 
         auto add_new = [&](lp::constraint_index new_ci) {
             SASSERT(!constraint_is_true(new_ci));
@@ -567,7 +568,7 @@ namespace nla {
         };
 
         TRACE(arith, tout << "factor (" << ci << ") " << p << " q: " << q << " vars: " << vars << "\n");
-        if (vars.empty()) {
+        if (false && vars.empty()) {
             for (auto v : p.free_vars()) {
                 if (value(v) == 0)
                     return subst(v, pddm.mk_val(0));
@@ -577,10 +578,13 @@ namespace nla {
                     return subst(v, pddm.mk_val(rational(-1)));
             }
             return l_undef;
-        }
+        }        
+        if (vars.empty())
+            return l_undef;
         for (auto v : vars) {
             if (value(v) == 0)
-                return subst(v, pddm.mk_val(0));
+                return l_undef;
+//                return subst(v, pddm.mk_val(0));
         }
         vector<dd::pdd> muls;
         muls.push_back(q);
@@ -661,7 +665,7 @@ namespace nla {
             explain_constraint(new_lemma, m.ci, ex);
         } 
         else if (std::holds_alternative<assumption_justification>(b)) {
-            auto [t, coeff] = to_term_offset(p);
+            auto [t, coeff] = to_term_offset(p);            
             new_lemma |= ineq(t, negate(k), -coeff);
         }
         else if (std::holds_alternative<gcd_justification>(b)) {

src/math/lp/nra_solver.cpp

@@ -57,6 +57,99 @@ struct solver::imp {
         m_lp2nl.reset();
     }
 
+
+    struct eq {
+        bool operator()(unsigned_vector const &a, unsigned_vector const &b) const {
+            return a == b;
+        }
+    };
+
+    map<unsigned_vector, unsigned, svector_hash<unsigned_hash>, eq> m_vars2mon;
+    // Create polynomial definition for variable v used in setup_assignment_solver.
+    // Side-effects: updates m_vars2mon when v is a monic variable.
+    void mk_definition(unsigned v, polynomial_ref_vector &definitions) {
+        auto &pm = m_nlsat->pm();
+        polynomial::polynomial_ref p(pm);
+        if (m_nla_core.emons().is_monic_var(v)) {
+            auto const &m = m_nla_core.emons()[v];
+            auto vars = m.vars();
+            std::sort(vars.begin(), vars.end());
+            m_vars2mon.insert(vars, v);
+            for (auto v2 : vars) {
+                auto pv = definitions.get(v2);
+                if (!p)
+                    p = pv;
+                else
+                    p = pm.mul(p, pv);
+            }
+        }
+        else if (lra.column_has_term(v)) {
+            for (auto const &[w, coeff] : lra.get_term(v)) {
+                auto pw = definitions.get(w);
+                if (!p)
+                    p = pm.mul(coeff, pw);
+                else
+                    p = pm.add(p, pm.mul(coeff, pw));
+            }
+        }
+        else {
+            p = pm.mk_polynomial(v);  // nlsat var index equals v (verified above when created)
+        }
+        definitions.push_back(p);
+    }
+
+    void setup_solver_poly() {
+        m_vars2mon.reset();
+        m_coi.init();
+        auto &pm = m_nlsat->pm();
+        polynomial_ref_vector definitions(pm);
+        for (unsigned v = 0; v < lra.number_of_vars(); ++v) {
+            auto j = m_nlsat->mk_var(lra.var_is_int(v));
+            VERIFY(j == v);
+            m_lp2nl.insert(v, j);  // we don't really need this. It is going to be the identify map.
+            scoped_anum a(am());
+            am().set(a, m_nla_core.val(v).to_mpq());
+            m_values->push_back(a);
+            mk_definition(v, definitions);
+        }
+
+        // we rely on that all information encoded into the tableau is present as a constraint.
+        for (auto ci : m_coi.constraints()) {
+            auto &c = lra.constraints()[ci];
+            auto &pm = m_nlsat->pm();
+            auto k = c.kind();
+            auto rhs = c.rhs();
+            auto lhs = c.coeffs();
+            rational den = denominator(rhs);
+            for (auto [coeff, v] : lhs)
+                den = lcm(den, denominator(coeff));
+            polynomial::polynomial_ref p(pm);
+            p = pm.mk_const(-den * rhs);
+
+            for (auto [coeff, v] : lhs) {
+                polynomial_ref poly(pm);
+                poly = pm.mul(den * coeff, definitions.get(v));
+                p = p + poly;
+            }
+            auto lit = add_constraint(p, ci, k);
+        }
+    }
+
+    void setup_solver_terms() {
+        m_coi.init();
+        // add linear inequalities from lra_solver
+        for (auto ci : m_coi.constraints())
+            add_constraint(ci);
+
+        // add polynomial definitions.
+        for (auto const &m : m_coi.mons())
+            add_monic_eq(m_nla_core.emons()[m]);
+
+        // add term definitions.
+        for (unsigned i : m_coi.terms())
+            add_term(i);
+    }
+
     /**
        \brief one-shot nlsat check.
        A one shot checker is the least functionality that can 
@@ -73,22 +166,17 @@ struct solver::imp {
         reset();
         vector<nlsat::assumption, false> core;
 
-        m_coi.init();
-        // add linear inequalities from lra_solver
-        for (auto ci : m_coi.constraints())
-            add_constraint(ci);
-        
-        // add polynomial definitions.
-        for (auto const& m : m_coi.mons())
-            add_monic_eq(m_nla_core.emons()[m]);
 
-        // add term definitions.
+        
-        for (unsigned i : m_coi.terms())
+        smt_params_helper p(m_params);
-            add_term(i);
+
+        if (p.arith_nl_nra_poly())
+            setup_solver_poly();
+        else 
+            setup_solver_terms();
 
         TRACE(nra, m_nlsat->display(tout));
 
-        smt_params_helper p(m_params);
         if (p.arith_nl_log()) {
             static unsigned id = 0;
             std::stringstream strm;
@@ -248,6 +336,23 @@ struct solver::imp {
         m_nlsat->mk_clause(1, &lit, a);
     }
 
+    nlsat::literal add_constraint(polynomial::polynomial *p, unsigned idx, lp::lconstraint_kind k) {
+        polynomial::polynomial *ps[1] = {p};
+        bool is_even[1] = {false};
+        nlsat::literal lit;
+        nlsat::assumption a = this + idx;
+        switch (k) {
+        case lp::lconstraint_kind::LE: lit = ~m_nlsat->mk_ineq_literal(nlsat::atom::kind::GT, 1, ps, is_even); break;
+        case lp::lconstraint_kind::GE: lit = ~m_nlsat->mk_ineq_literal(nlsat::atom::kind::LT, 1, ps, is_even); break;
+        case lp::lconstraint_kind::LT: lit = m_nlsat->mk_ineq_literal(nlsat::atom::kind::LT, 1, ps, is_even); break;
+        case lp::lconstraint_kind::GT: lit = m_nlsat->mk_ineq_literal(nlsat::atom::kind::GT, 1, ps, is_even); break;
+        case lp::lconstraint_kind::EQ: lit = m_nlsat->mk_ineq_literal(nlsat::atom::kind::EQ, 1, ps, is_even); break;
+        default: UNREACHABLE();  // unreachable
+        }
+        m_nlsat->mk_clause(1, &lit, a);
+        return lit;
+    }
+
     bool check_monic(mon_eq const& m) {
         scoped_anum val1(am()), val2(am());
         am().set(val1, value(m.var()));

src/params/smt_params_helper.pyg

@@ -91,6 +91,7 @@ def_module_params(module_name='smt',
                           ('arith.nl.cross_nested', BOOL, True, 'enable cross-nested consistency checking'),
                           ('arith.nl.stellensatz', BOOL, False, 'enable stellensatz saturation'),
                           ('arith.nl.linearize', BOOL, True, 'enable NL linearization'),
+                          ('arith.nl.nra.poly', BOOL, False, 'use polynomial translation'),
                           ('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'),