add option for gcd-test to grobner

src/math/dd/dd_pdd.cpp

@@ -1208,6 +1208,35 @@ namespace dd {
         return is_binary(p.root); 
     }
 
+    /*
+    * Retrieve variables that occur as powers.
+    */
+    void pdd_manager::get_powers(pdd const& p, svector<std::pair<unsigned, unsigned>>& powers) {
+        powers.reset();
+        init_mark();
+        m_todo.push_back(p.root);
+        while (!m_todo.empty()) {
+            PDD r = m_todo.back();
+            m_todo.pop_back();
+            if (is_val(r) || is_marked(r)) 
+                continue;
+            set_mark(r);
+            unsigned v = var(r);
+            unsigned d = 1;
+            while (!is_val(r) && var(hi(r)) == v) {
+                d++;
+                m_todo.push_back(lo(r));
+                r = hi(r);
+            }
+            if (d > 1)
+                powers.push_back({ v, d });
+            if (!is_val(r)) {
+                m_todo.push_back(hi(r));
+                m_todo.push_back(lo(r));
+            }
+        }
+    }
+
     /**
        Determine if p is a monomial.
     */
@@ -1986,6 +2015,42 @@ namespace dd {
     pdd_iterator pdd::begin() const { return pdd_iterator(*this, true); }
     pdd_iterator pdd::end() const { return pdd_iterator(*this, false); }
 
+    pdd_coeff_iterator pdd::pdd_coeffients::begin() const { return m_pdd.begin_coeff(); }
+    pdd_coeff_iterator pdd::pdd_coeffients::end() const { return m_pdd.end_coeff(); }
+
+    void pdd_coeff_iterator::next() {
+        auto& m = m_pdd.manager();
+        while (!m_nodes.empty()) {
+            auto p = m_nodes.back();
+            m_nodes.pop_back();
+            SASSERT(!m.is_val(p));
+            unsigned n = m.lo(p);
+            if (m.is_val(n) && m.val(n).is_zero())
+                continue;                
+            while (!m.is_val(n)) {
+                m_nodes.push_back(n);
+                n = m.hi(n);
+            }
+            m_value = { m.val(n), m_nodes.empty()};
+            return;
+        }
+        at_end = true;
+    }
+
+    void pdd_coeff_iterator::first() {
+        unsigned n = m_pdd.root;
+        auto& m = m_pdd.manager();
+        while (!m.is_val(n)) {
+            m_nodes.push_back(n);
+            n = m.hi(n);
+        }
+        at_end = false;
+        m_value = { m.val(n), m_nodes.empty() };
+    }
+
+    pdd_coeff_iterator pdd::begin_coeff() const { return pdd_coeff_iterator(*this, true); }
+    pdd_coeff_iterator pdd::end_coeff() const { return pdd_coeff_iterator(*this, false); }
+
     std::ostream& operator<<(std::ostream& out, pdd_monomial const& m) {
         if (!m.coeff.is_one()) {
             out << m.coeff;

src/math/dd/dd_pdd.h

@@ -46,6 +46,7 @@ namespace dd {
     class pdd_manager;
     class pdd_iterator;
     class pdd_linear_iterator;
+    class pdd_coeff_iterator;
 
     class pdd_manager {
     public:
@@ -55,6 +56,7 @@ namespace dd {
         friend pdd;
         friend pdd_iterator;
         friend pdd_linear_iterator;
+        friend pdd_coeff_iterator;
 
         typedef unsigned PDD;
         typedef vector<std::pair<rational,unsigned_vector>> monomials_t;
@@ -366,6 +368,8 @@ namespace dd {
         bool is_binary(PDD p);
         bool is_binary(pdd const& p);
 
+        void get_powers(pdd const& p, svector<std::pair<unsigned, unsigned>>& powers);
+
         bool is_monomial(PDD p);
 
         bool is_univariate(PDD p);
@@ -406,6 +410,7 @@ namespace dd {
         friend class pdd_manager;
         friend class pdd_iterator;
         friend class pdd_linear_iterator;
+        friend class pdd_coeff_iterator;
         unsigned     root;
         pdd_manager* m;
         pdd(unsigned root, pdd_manager& pm): root(root), m(&pm) { m->inc_ref(root); }
@@ -440,6 +445,7 @@ namespace dd {
         bool is_unary() const { return !is_val() && lo().is_zero() && hi().is_val(); }
         bool is_offset() const { return !is_val() && lo().is_val() && hi().is_one(); }
         bool is_binary() const { return m->is_binary(root); }
+        void get_powers(svector<std::pair<unsigned, unsigned>>& powers) const { m->get_powers(*this, powers); }
         bool is_monomial() const { return m->is_monomial(root); }
         bool is_univariate() const { return m->is_univariate(root); }
         bool is_univariate_in(unsigned v) const { return m->is_univariate_in(root, v); }
@@ -508,6 +514,20 @@ namespace dd {
         pdd_iterator begin() const;
         pdd_iterator end() const;
 
+
+        pdd_coeff_iterator begin_coeff() const;
+        pdd_coeff_iterator end_coeff() const;
+
+        class pdd_coeffients {
+            friend class pdd;
+            pdd const& m_pdd;
+            pdd_coeffients(pdd const& p) : m_pdd(p) {}
+        public:
+            pdd_coeff_iterator begin() const;
+            pdd_coeff_iterator end() const;
+        };
+        pdd_coeffients coefficients() const { return pdd_coeffients(*this); }
+
         class pdd_linear_monomials {
             friend class pdd;
             pdd const& m_pdd;
@@ -593,6 +613,31 @@ namespace dd {
         bool operator!=(pdd_linear_iterator const& other) const { return m_next != other.m_next; }
     };
 
+    // iterator class that visits each coefficient with information
+    // if it is a constant or is used as a coefficient to a monomial.
+
+    struct coeff_value {
+        rational coeff;
+        bool is_constant; // true if it is a constant term.
+    };
+
+    class pdd_coeff_iterator {
+        friend class pdd;
+        pdd m_pdd;
+        unsigned_vector m_nodes;
+        coeff_value m_value;
+        bool at_end = true;
+        pdd_coeff_iterator(pdd const& p, bool at_start): m_pdd(p) { if (at_start) first(); }
+        void first();
+        void next();
+    public:
+        coeff_value operator*() const { return m_value; }
+        coeff_value const* operator->() const { return &m_value; }
+        pdd_coeff_iterator& operator++() { next(); return *this; }
+        pdd_coeff_iterator operator++(int) { auto tmp = *this; next(); return tmp; }
+        bool operator!=(pdd_coeff_iterator const& other) const { return at_end != other.at_end; }
+    };
+
     class val_pp {
         pdd_manager const& m;
         rational const& val;

src/math/lp/nla_grobner.cpp

@@ -31,6 +31,7 @@ namespace nla {
     void grobner::updt_params(params_ref const& p) {
         smt_params_helper ph(p);
         m_config.m_propagate_quotients = ph.arith_nl_grobner_propagate_quotients();
+        m_config.m_gcd_test = ph.arith_nl_grobner_gcd_test();
     }
 
     lp::lp_settings& grobner::lp_settings() {
@@ -81,6 +82,9 @@ namespace nla {
 
             if (propagate_quotients())
                 return;
+
+            if (propagate_gcd_test())
+                return;
             
         }
         catch (...) {
@@ -99,14 +103,6 @@ namespace nla {
         IF_VERBOSE(5, diagnose_pdd_miss(verbose_stream()));
     }
 
-    dd::solver::equation_vector const& grobner::core_equations(bool all_eqs) {
-        flet<bool> _add_all(m_add_all_eqs, all_eqs);
-        find_nl_cluster();        
-        if (!configure()) 
-            throw dd::pdd_manager::mem_out();
-        return m_solver.equations();
-    }
-
     bool grobner::is_conflicting() {
         for (auto eq : m_solver.equations()) {
             if (is_conflicting(*eq)) {
@@ -228,6 +224,76 @@ namespace nla {
         return {t, offset};
     }
 
+    bool grobner::propagate_gcd_test() {
+        if (!m_config.m_gcd_test)
+            return false;
+        unsigned changed = 0;
+        for (auto eq : m_solver.equations())
+            if (propagate_gcd_test(*eq) && ++changed >= m_solver.number_of_conflicts_to_report())
+                return true;
+        return changed > 0;
+    }
+
+    // x^2 - y = 0 => y mod 4 = { 0, 1 }
+    // x^k - y = 0 => y mod k^2 = {0, 1, .., (k-1)^k }
+    bool grobner::propagate_gcd_test(dd::solver::equation const& eq) {
+        dd::pdd p = eq.poly();
+        if (p.is_linear())
+            return false;
+        if (p.is_val())
+            return false;
+        auto v = p.var();
+        if (!c().var_is_int(v))
+            return false;
+        for (auto v : p.free_vars())
+            if (!c().var_is_int(v))
+                return false;
+        dd::pdd_eval eval;
+        eval.var2val() = [&](unsigned j) { return val(j); };
+        if (eval(p) == 0)
+            return false;
+        p.get_powers(m_powers);
+        if (m_powers.empty())
+            return false;
+        rational d(1);
+        for (auto const& value : p.coefficients())
+            d = lcm(d, denominator(value.coeff));
+        if (d != 1)
+            p *= d;
+        auto& m = p.manager();
+        auto fails_gcd_test = [&](dd::pdd const& q) {
+            rational g(0);
+            rational k(0);
+            for (auto const& value : q.coefficients()) {
+                if (value.is_constant) {
+                    k += value.coeff;
+                    continue;
+                }
+                g = gcd(g, value.coeff);
+                if (g == 1)
+                    return false;                
+            }
+            return k != 0 && !divides(g, k);
+        };
+        for (auto [x, k] : m_powers) {
+            SASSERT(k > 1);
+            bool gcd_fail = true;
+            dd::pdd kx = m.mk_var(x) * m.mk_val(k); 
+            for (unsigned r = 0; gcd_fail && r < k; r++) {
+                dd::pdd kx_plus_r = kx + m.mk_val(r);
+                auto q = p.subst_pdd(x, kx_plus_r);
+                if (!fails_gcd_test(q))
+                    gcd_fail = false;
+            }
+            if (gcd_fail) {
+                lemma_builder lemma(c(), "pdd-power");
+                add_dependencies(lemma, eq);
+                return true;
+            }
+        }
+        return false;
+    }
+
     bool grobner::propagate_quotients() {
         if (!m_config.m_propagate_quotients)
             return false;
@@ -249,8 +315,6 @@ namespace nla {
     // z < 0 & x < 0 => -x <= -z
     bool grobner::propagate_quotients(dd::solver::equation const& eq) {
         dd::pdd const& p = eq.poly();
-        dd::pdd_eval eval;
-        eval.var2val() = [&](unsigned j) { return val(j); };
         if (p.is_linear())
             return false;
         if (p.is_val())
@@ -261,6 +325,8 @@ namespace nla {
         for (auto v : p.free_vars())
             if (!c().var_is_int(v))
                 return false;
+        dd::pdd_eval eval;
+        eval.var2val() = [&](unsigned j) { return val(j); };
         if (eval(p) == 0)
             return false;
         tracked_uint_set nl_vars;
@@ -577,8 +643,12 @@ namespace nla {
                 for (const factor& fc: fcn) 
                     q.push_back(var(fc));        
         }
-
-        if (c().var_is_fixed(j))
+        else if (c().lra.column_has_term(j)) {
+            const lp::lar_term& t = c().lra.get_term(j);
+            for (auto k : t)
+                q.push_back(k.j());
+        }
+        else if (c().var_is_fixed(j))
             return;
         const auto& matrix = lra.A_r();
         for (auto & s : matrix.m_columns[j]) {
@@ -587,12 +657,18 @@ namespace nla {
                 continue;
             m_rows.insert(row);
             unsigned k = lra.get_base_column_in_row(row);
-            // grobner bassis does not know about integer constraints
-            if (lra.column_is_free(k) && !m_add_all_eqs && k != j)
-                continue;
             // a free column over the reals can be assigned
-            if (lra.column_is_free(k) && k != j && !lra.var_is_int(k)) 
-                continue;
+            if (lra.column_is_free(k) && k != j) {
+                if (!lra.var_is_int(k))
+                    continue;
+                // free integer columns are ignored unless m_add_all_eqs is set or we are doing gcd test.
+                if (!m_add_all_eqs && !m_config.m_gcd_test)
+                    continue;
+                // a free integer column with integer coefficients can be assigned.
+                if (!m_add_all_eqs && all_of(c().lra.get_row(row), [&](auto& ri) { return ri.coeff().is_int();}))
+                    continue;
+            }
+
             CTRACE(grobner, matrix.m_rows[row].size() > c().params().arith_nl_grobner_row_length_limit(),
                    tout << "ignore the row " << row << " with the size " << matrix.m_rows[row].size() << "\n";);
             // limits overhead of grobner equations, unless this is for extracting a complete COI of the non-satisfied subset.
@@ -618,6 +694,7 @@ namespace nla {
         while (!vars.empty()) {
             j = vars.back();
             vars.pop_back();
+
             if (c().params().arith_nl_grobner_subs_fixed() > 0 && c().var_is_fixed_to_zero(j)) {
                 r = m_pdd_manager.mk_val(val_of_fixed_var_with_deps(j, zero_dep));
                 dep = zero_dep;

src/math/lp/nla_grobner.h

@@ -21,6 +21,7 @@ namespace nla {
     class grobner : common {
         struct config {
             bool m_propagate_quotients = false;
+            bool m_gcd_test = false;
         };
         dd::pdd_manager          m_pdd_manager;
         dd::solver               m_solver;
@@ -51,6 +52,10 @@ namespace nla {
         bool propagate_quotients();
         bool propagate_quotients(dd::solver::equation const& eq);
 
+        svector<std::pair<unsigned, unsigned>> m_powers;
+        bool propagate_gcd_test();
+        bool propagate_gcd_test(dd::solver::equation const& eq);
+
         std::pair<lp::lar_term, rational> linear_to_term(dd::pdd q);
         
         void add_dependencies(lemma_builder& lemma, dd::solver::equation const& eq);
@@ -74,7 +79,7 @@ namespace nla {
         bool is_solved(dd::pdd const& p, unsigned& v, dd::pdd& r);
         void add_eq(dd::pdd& p, u_dependency* dep);        
         const rational& val_of_fixed_var_with_deps(lpvar j, u_dependency*& dep);
-        dd::pdd pdd_expr(const rational& c, lpvar j, u_dependency*& dep);                
+        dd::pdd pdd_expr(const rational& c, lpvar j, u_dependency*& dep);  
 
         void display_matrix_of_m_rows(std::ostream& out) const;
         std::ostream& diagnose_pdd_miss(std::ostream& out);
@@ -83,6 +88,6 @@ namespace nla {
         grobner(core *core);        
         void operator()();
         void updt_params(params_ref const& p);
-        dd::solver::equation_vector const& core_equations(bool all_eqs);
+        // dd::solver::equation_vector const& core_equations(bool all_eqs);
     }; 
 }

src/math/lp/nla_pp.cpp

@@ -19,13 +19,13 @@ std::ostream& core::print_product(const T& m, std::ostream& out) const {
     bool first = true;
     for (lpvar v : m) {
         if (!first)
-            out << "*";
+            out << " * ";
         else
             first = false;
         if (lp_settings().print_external_var_name())
-            out << "(" << lra.get_variable_name(v) << "=" << val(v) << ")";
+            out << lra.get_variable_name(v);
         else
-            out << "(j" << v << " = " << val(v) << ")";
+            out << "j" << v;
     }
     return out;
 }
@@ -69,13 +69,13 @@ std::ostream& core::print_factor_with_vars(const factor& f, std::ostream& out) c
 
 std::ostream& core::print_monic(const monic& m, std::ostream& out) const {
     if (lp_settings().print_external_var_name())
-        out << "([" << m.var() << "] = " << lra.get_variable_name(m.var()) << " = " << val(m.var()) << " = ";
+        out << "[" << m.var() << "] := " << lra.get_variable_name(m.var()) << " := ";
     else
-        out << "(j" << m.var() << " = " << val(m.var()) << " = ";
-    print_product(m.vars(), out) << ")\n";
+        out << "j" << m.var() << " := ";
+    print_product(m.vars(), out) << "\n";
     return out;
 }
-
+ 
 std::ostream& core::print_bfc(const factorization& m, std::ostream& out) const {
     SASSERT(m.size() == 2);
     out << "( x = " << pp(m[0]) << "* y = " << pp(m[1]) << ")";
@@ -85,6 +85,7 @@ std::ostream& core::print_bfc(const factorization& m, std::ostream& out) const {
 std::ostream& core::print_monic_with_vars(lpvar v, std::ostream& out) const {
     return print_monic_with_vars(m_emons[v], out);
 }
+
 template <typename T>
 std::ostream& core::print_product_with_vars(const T& m, std::ostream& out) const {
     print_product(m, out) << "\n";
@@ -141,9 +142,8 @@ std::ostream& core::print_var(lpvar j, std::ostream& out) const {
 }
 
 std::ostream& core::print_monics(std::ostream& out) const {
-    for (auto& m : m_emons) {
-        print_monic_with_vars(m, out);
-    }
+    for (auto& m : m_emons) 
+        print_monic_with_vars(m, out);    
     return out;
 }
 
@@ -312,8 +312,7 @@ std::ostream& core::display_row(std::ostream& out, lp::row_strip<lp::mpq> const&
 }
 
 std::ostream& core::display(std::ostream& out) {
-    print_monics(out);
-    for (unsigned i = 0; i < lra.row_count(); ++i)
-        display_row(out, lra.get_row(i));
+    for (auto& m : m_emons)
+        print_monic(m, out);
     return out;
 }

src/params/smt_params_helper.pyg

@@ -79,6 +79,7 @@ def_module_params(module_name='smt',
                           ('arith.nl.grobner_max_simplified', UINT, 10000, 'grobner\'s maximum number of simplifications'),
                           ('arith.nl.grobner_cnfl_to_report', UINT, 1, 'grobner\'s maximum number of conflicts to report'),
                           ('arith.nl.grobner_propagate_quotients', BOOL, False, 'detect conflicts x*y + z = 0 where x doesn\'t divide z'),
+                          ('arith.nl.grobner_gcd_test', BOOL, False, 'detect gcd conflicts for polynomial powers x^k - y = 0'),
                           ('arith.nl.gr_q', UINT, 10, 'grobner\'s quota'),
                           ('arith.nl.grobner_subs_fixed', UINT, 1, '0 - no subs, 1 - substitute, 2 - substitute fixed zeros only'),   
 	                  ('arith.nl.delay', UINT, 10, 'number of calls to final check before invoking bounded nlsat check'),