From caf314e569bcfd291df592daa900dc5d38a58c5e Mon Sep 17 00:00:00 2001
From: Nikolaj Bjorner <nbjorner@microsoft.com>
Date: Tue, 20 Jun 2023 17:05:48 -0700
Subject: [PATCH] experiment with better gcd test for non-linear arihmetic

Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>
---
 src/math/dd/dd_pdd.cpp            |   9 ++
 src/math/dd/dd_pdd.h              |   1 +
 src/math/dd/pdd_interval.h        |  42 ++++++++
 src/math/interval/dep_intervals.h |  28 +++++
 src/math/lp/int_solver.cpp        |   2 +
 src/math/lp/lar_solver.h          |   8 +-
 src/math/lp/lp_settings.h         |   2 +
 src/math/lp/nla_core.h            |   2 +-
 src/math/lp/nla_grobner.cpp       | 165 ++++++++++++++++++++++++++----
 src/math/lp/nla_grobner.h         |   9 +-
 src/math/lp/nla_intervals.cpp     |  32 +++++-
 src/math/lp/nla_intervals.h       |   2 +
 12 files changed, 274 insertions(+), 28 deletions(-)

diff --git a/src/math/dd/dd_pdd.cpp b/src/math/dd/dd_pdd.cpp
index 291af3b5c..e08c42b33 100644
--- a/src/math/dd/dd_pdd.cpp
+++ b/src/math/dd/dd_pdd.cpp
@@ -1959,5 +1959,14 @@ namespace dd {
         return out;
     }
 
+    pdd& operator+=(pdd & p, pdd_monomial const& m) {
+        pdd q(p.manager());
+        q = m.coeff;
+        for (auto v : m.vars)
+            q *= p.manager().mk_var(v);
+        return p += q;
+    }
+
+
 
 }
diff --git a/src/math/dd/dd_pdd.h b/src/math/dd/dd_pdd.h
index 6dee7977f..1fd6f068f 100644
--- a/src/math/dd/dd_pdd.h
+++ b/src/math/dd/dd_pdd.h
@@ -536,6 +536,7 @@ namespace dd {
     };
 
     std::ostream& operator<<(std::ostream& out, pdd_monomial const& m);
+    pdd& operator+=(pdd& p, pdd_monomial const& m);
 
     class pdd_iterator {
         friend class pdd;
diff --git a/src/math/dd/pdd_interval.h b/src/math/dd/pdd_interval.h
index c5b4c9e1c..f3e2aafa7 100644
--- a/src/math/dd/pdd_interval.h
+++ b/src/math/dd/pdd_interval.h
@@ -29,6 +29,7 @@ typedef dep_intervals::with_deps_t w_dep;
 class pdd_interval {
     dep_intervals& m_dep_intervals;
     std::function<void (unsigned, bool, scoped_dep_interval&)> m_var2interval;
+    std::function<void (unsigned, bool, vector<scoped_dep_interval>&)> m_var2intervals;
 
     // retrieve intervals after distributing multiplication over addition.
     template <w_dep wd>
@@ -65,6 +66,9 @@ public:
     std::function<void (unsigned, bool, scoped_dep_interval&)>& var2interval() { return m_var2interval; } // setter
     const std::function<void (unsigned, bool, scoped_dep_interval&)>& var2interval() const { return m_var2interval; } // getter
 
+    std::function<void (unsigned, bool, vector<scoped_dep_interval>&)>& var2intervals() { return m_var2intervals; } // setter
+    const std::function<void (unsigned, bool, vector<scoped_dep_interval>&)>& var2intervals() const { return m_var2intervals; } // getter
+
     template <w_dep wd>
     void get_interval(pdd const& p, scoped_dep_interval& ret) {
         if (p.is_val()) {
@@ -91,5 +95,43 @@ public:
         m_dep_intervals.set_interval_for_scalar(i, rational::one());        
         get_interval_distributed<wd>(p, i, ret);
     }
+
+
+    //
+    // produce an explanation for a range using weaker bounds.
+    //
+    // lo_interval := interval(lo)
+    // hi_bound    := bound - lo_interval
+    // hi_interval := explain(var*hi, hi_bound);
+    // lo_bound    := bound - hi_interval
+    // lo_interval := explain(lo, lo_bound);
+    // return lo_interval + hi_interval
+    // 
+    void explain(pdd const& p, interval const& bound, scoped_dep_interval& ret) {
+        if (p.is_val()) {
+            m_dep_intervals.set_interval_for_scalar(ret, p.val());
+            return;
+        }
+        scoped_dep_interval a(m()), t(m()), hi(m()), lo(m());
+        scoped_dep_interval lo_interval(m()), lo_bound(m());
+        scoped_dep_interval hi_interval(m()), hi_bound(m());
+        if (!p.hi().is_val()) {
+            m_var2interval(p.var(), true, a);
+            get_interval<w_dep::with_deps>(p.hi(), hi);
+            m_dep_intervals.mul<dep_intervals::with_deps>(hi, a, hi_interval);
+            m_dep_intervals.sub(bound, hi_interval, lo_bound);
+            explain(p.lo(), lo_bound, lo_interval);
+            m_dep_intervals.add<dep_intervals::with_deps>(lo_interval, hi_interval, ret);
+        }
+        else {
+            get_interval<w_dep::without_deps>(p.lo(), lo_interval);
+            m_dep_intervals.sub(bound, lo_interval, hi_bound);
+            m_dep_intervals.div(hi_bound, p.hi().val(), hi_bound);
+            vectro<scoped_dep_interval> as;
+            m_var2intervals(p.var(), true, as);
+            // use hi_bound to adjust for variable bound.
+        }
+        
+    }
 };
 }
diff --git a/src/math/interval/dep_intervals.h b/src/math/interval/dep_intervals.h
index d641a294d..e5523e09f 100644
--- a/src/math/interval/dep_intervals.h
+++ b/src/math/interval/dep_intervals.h
@@ -195,6 +195,8 @@ public:
     void add(const rational& r, interval& a) const;
     void mul(const interval& a, const interval& b, interval& c) { m_imanager.mul(a, b, c); }
     void add(const interval& a, const interval& b, interval& c) { m_imanager.add(a, b, c); }
+    void sub(const interval& a, const interval& b, interval& c) { m_imanager.sub(a, b, c); }
+    void div(const interval& a, const mpq& b, interval& c) { m_imanager.div(a, b, c); }
 
     template <enum with_deps_t wd>
     void set(interval& a, const interval& b) const {
@@ -271,6 +273,30 @@ public:
         }
     }
 
+    template <enum with_deps_t wd>
+    void sub(const interval& a, const interval& b, interval& c) {
+        if (wd == with_deps) {
+            interval_deps_combine_rule comb_rule;
+            sub(a, b, c, comb_rule);
+            combine_deps(a, b, comb_rule, c);
+        }
+        else {
+            sub(a, b, c);
+        }
+    }
+    
+    template <enum with_deps_t wd>
+    void div(const interval& a, const mpq& b, interval& c) {
+        if (wd == with_deps) {
+            interval_deps_combine_rule comb_rule;
+            div(a, b, c, comb_rule);
+            combine_deps(a, comb_rule, c);
+        }
+        else {
+            div(a, b, c);
+        }
+    }
+
     template <enum with_deps_t wd>
     void copy_upper_bound(const interval& a, interval& i) const {
         SASSERT(a.m_upper_inf == false);
@@ -413,6 +439,8 @@ public:
 private:
     void mul(const interval& a, const interval& b, interval& c, interval_deps_combine_rule& deps) { m_imanager.mul(a, b, c, deps); }
     void add(const interval& a, const interval& b, interval& c, interval_deps_combine_rule& deps) { m_imanager.add(a, b, c, deps); }
+    void sub(const interval& a, const interval& b, interval& c, interval_deps_combine_rule& deps) { m_imanager.sub(a, b, c, deps); }
+    void div(const interval& a, const mpq& b, interval& c, interval_deps_combine_rule& deps) { m_imanager.div(a, b, c, deps); }
    
     void combine_deps(interval const& a, interval const& b, interval_deps_combine_rule const& deps, interval& i) const {
         m_config.add_deps(a, b, deps, i);
diff --git a/src/math/lp/int_solver.cpp b/src/math/lp/int_solver.cpp
index b4c6bff99..7516b35ea 100644
--- a/src/math/lp/int_solver.cpp
+++ b/src/math/lp/int_solver.cpp
@@ -52,6 +52,8 @@ namespace lp {
         unsigned non_int_before, non_int_after;
 
         non_int_before = count_non_int();
+
+
         unsigned num = lra.A_r().column_count();
         for (unsigned j : lra.r_basis()) 
             if (!lra.get_value(j).is_int())
diff --git a/src/math/lp/lar_solver.h b/src/math/lp/lar_solver.h
index 77fcbfc14..b3b7e09aa 100644
--- a/src/math/lp/lar_solver.h
+++ b/src/math/lp/lar_solver.h
@@ -368,7 +368,8 @@ public:
 
     bool is_fixed_at_bound(column_index const& j, vector<std::tuple<explanation, column_index, bool, mpq>>& bounds);
     bool has_fixed_at_bound(vector<std::tuple<explanation, column_index, bool, mpq>>& bounds);
-    
+
+
     bool is_fixed(column_index const& j) const { return column_is_fixed(j); }    
     inline column_index to_column_index(unsigned v) const { return column_index(external_to_column_index(v)); }
     bool external_is_used(unsigned) const;
@@ -388,6 +389,11 @@ public:
             m_term_register.local_to_external(idx) : m_var_register.local_to_external(idx);
     }
     bool column_corresponds_to_term(unsigned) const;
+    const lar_term & column_to_term(unsigned j) const {
+        SASSERT(column_corresponds_to_term(j));
+        return get_term(column2tv(to_column_index(j)));
+    }
+
     inline unsigned row_count() const { return A_r().row_count(); }
     bool var_is_registered(var_index vj) const;
     void clear_inf_set() {
diff --git a/src/math/lp/lp_settings.h b/src/math/lp/lp_settings.h
index c213333e0..07df58bcf 100644
--- a/src/math/lp/lp_settings.h
+++ b/src/math/lp/lp_settings.h
@@ -120,6 +120,7 @@ struct statistics {
     unsigned m_cross_nested_forms;
     unsigned m_grobner_calls;
     unsigned m_grobner_conflicts;
+    unsigned m_grobner_propagations;
     unsigned m_offset_eqs;
     unsigned m_fixed_eqs;
     statistics() { reset(); }
@@ -142,6 +143,7 @@ struct statistics {
         st.update("arith-horner-cross-nested-forms", m_cross_nested_forms);
         st.update("arith-grobner-calls", m_grobner_calls);
         st.update("arith-grobner-conflicts", m_grobner_conflicts);
+        st.update("arith-grobner-propagations", m_grobner_propagations);
         st.update("arith-offset-eqs", m_offset_eqs);
         st.update("arith-fixed-eqs", m_fixed_eqs);
 
diff --git a/src/math/lp/nla_core.h b/src/math/lp/nla_core.h
index 938bcbe83..e627f9d0a 100644
--- a/src/math/lp/nla_core.h
+++ b/src/math/lp/nla_core.h
@@ -174,7 +174,7 @@ public:
     }
 
     bool need_run_grobner() const { 
-        return m_nla_settings.run_grobner && lp_settings().stats().m_nla_calls % m_nla_settings.grobner_frequency == 0; 
+        return m_nla_settings.run_grobner; // && lp_settings().stats().m_nla_calls % m_nla_settings.grobner_frequency == 0; 
     }
 
     void set_active_vars_weights(nex_creator&);
diff --git a/src/math/lp/nla_grobner.cpp b/src/math/lp/nla_grobner.cpp
index 974c48d14..0c3c98f1d 100644
--- a/src/math/lp/nla_grobner.cpp
+++ b/src/math/lp/nla_grobner.cpp
@@ -24,8 +24,8 @@ namespace nla {
         common(c),
         m_pdd_manager(m_core.m_lar_solver.number_of_vars()),
         m_solver(m_core.m_reslim, m_pdd_manager),
-        m_lar_solver(m_core.m_lar_solver)
-
+        m_lar_solver(m_core.m_lar_solver),
+        m_quota(m_core.m_nla_settings.grobner_quota*2)
     {}
 
     lp::lp_settings& grobner::lp_settings() {
@@ -33,19 +33,26 @@ namespace nla {
     }
 
     void grobner::operator()() {
-        unsigned& quota = c().m_nla_settings.grobner_quota;
-        if (quota == 1)
+        if (m_quota <= c().m_nla_settings.grobner_quota) {
+            ++m_quota;
             return;
+        }
 
         lp_settings().stats().m_grobner_calls++;
         find_nl_cluster();        
         configure();
         m_solver.saturate();
 
+        m_quota += 2;
+
         if (is_conflicting())
             return;
 
         try {
+
+            if (propagate_gcd())
+                return;
+
             if (propagate_bounds())
                 return;
 
@@ -59,10 +66,12 @@ namespace nla {
             
         }
 
-        if (quota > 1)
-            quota--;
+        m_quota -= 3;
 
-        IF_VERBOSE(2, verbose_stream() << "grobner miss, quota " << quota << "\n");
+        TRACE("grobner", tout << "saturated\n");
+
+
+        IF_VERBOSE(2, verbose_stream() << "grobner miss, quota " << m_quota << "\n");
         IF_VERBOSE(4, diagnose_pdd_miss(verbose_stream()));
 
 #if 0
@@ -89,28 +98,37 @@ namespace nla {
         return conflicts > 0;
     }
 
-    bool grobner::propagate_bounds() {
+    bool grobner::propagate(std::function<bool(dd::solver::equation const&)>& fn) {
         unsigned changed = 0;
-        for (auto eq : m_solver.equations()) 
-            if (propagate_bounds(*eq) && ++changed >= m_solver.number_of_conflicts_to_report())
-                return true;
+        for (auto eq : m_solver.equations()) {
+            if (fn(*eq))
+                ++changed;
+            if (changed >= m_solver.number_of_conflicts_to_report())
+                break;            
+        }
+        lp_settings().stats().m_grobner_propagations += changed;
         return changed > 0;
     }
 
+
+    bool grobner::propagate_bounds() {
+        std::function<bool(dd::solver::equation const&)> p = [&](auto const& eq) { return propagate_bounds(eq); };
+        return propagate(p);
+    }
+
     bool grobner::propagate_eqs() {
-        unsigned changed = 0;
-        for (auto eq : m_solver.equations())
-            if (propagate_fixed(*eq) && ++changed >= m_solver.number_of_conflicts_to_report())
-                return true;
-        return changed > 0;
+        std::function<bool(dd::solver::equation const&)> p = [&](auto const& eq) { return propagate_fixed(eq); };
+        return propagate(p);
     }
 
     bool grobner::propagate_factorization() {
-        unsigned changed = 0;
-        for (auto eq : m_solver.equations())
-            if (propagate_factorization(*eq) && ++changed >= m_solver.number_of_conflicts_to_report())
-                return true;
-        return changed > 0;
+        std::function<bool(dd::solver::equation const&)> p = [&](auto const& eq) { return propagate_factorization(eq); };
+        return propagate(p);
+    }
+
+    bool grobner::propagate_gcd() {
+        std::function<bool(dd::solver::equation const&)> p = [&](auto const& eq) { return propagate_gcd(eq); };
+        return propagate(p);
     }
 
     /**
@@ -189,6 +207,108 @@ namespace nla {
         return true;
     }
 
+    //
+    // v1: 
+    // p + mx * q = 0
+    // such that -mx < lb(p) <= ub(p) < mx
+    // => propagate p = 0
+    //
+    // other versions TBD
+    // p + mx * q = 0
+    // such that ub(p) - lb(p) <= 2mx - 1
+    // => propagate p = k
+    // 
+    // p + mx * q = 0
+    // such that ub(p) does not divide mx
+    // => propagte p <= floor(ub(p)/mx)
+    //
+    // also the gcd tests from int_gcd_test
+    // 
+    bool grobner::propagate_gcd(const dd::solver::equation& eq) {
+        dd::pdd const& p = eq.poly();
+        rational denom(1), mx(0);
+
+        auto is_fixed = [&](dd::pdd_monomial const& m) {
+            for (auto v : m.vars)
+                if (!c().var_is_fixed(v))
+                    return false;
+            return true;
+        };
+        
+        if (p.is_val())
+            return false;
+        for (auto const& m : p) {
+            for (auto v : m.vars)
+                if (!m_lar_solver.column_is_int(v))
+                    return false;            
+            denom = lcm(denom, denominator(m.coeff));
+            if (is_fixed(m))
+                continue;
+            mx = std::max(mx, abs(m.coeff));
+        }
+        dd::pdd lo(p.manager());
+        lo = 0;
+        for (auto const& m : p) 
+            if (abs(m.coeff) < mx || is_fixed(m)) 
+                lo += m;
+        
+        if (lo.is_zero())
+            return false;
+        
+        auto& di = c().m_intervals.get_dep_intervals();
+        dd::pdd_interval eval(di);
+        eval.var2interval() = [this](lpvar j, bool deps, scoped_dep_interval& a) {
+            if (deps) c().m_intervals.set_var_interval<dd::w_dep::with_deps>(j, a);
+            else c().m_intervals.set_var_interval<dd::w_dep::without_deps>(j, a);
+        };
+        scoped_dep_interval i(di), i_wd(di);
+        eval.get_interval<dd::w_dep::without_deps>(lo, i);    
+
+        if (di.lower_is_inf(i))
+            return false;
+        if (di.upper_is_inf(i))
+            return false;
+        if (di.lower(i) <= -mx)
+            return false;
+        if (di.upper(i) >= mx)
+            return false;
+
+        // TODO - relax bound dependencies to weakest that admit within interval -mx, mx.
+        eval.get_interval<dd::w_dep::with_deps>(lo, i_wd);
+
+        lp::lar_term lo_t;
+        rational k(0);
+        for (auto const& m : lo) {
+            unsigned j;
+            auto coeff = m.coeff * denom;
+            if (m.vars.empty())
+                k -= coeff;
+            else if (m.vars.size() == 1) {
+                lo_t.add_monomial(coeff, m.vars[0]);
+            }
+            else if (c().find_canonical_monic_of_vars(m.vars, j)) {
+                verbose_stream() << "canonical monic\n";
+                lo_t.add_monomial(coeff, j);
+            }
+            else
+                return false;
+        }
+
+        c().m_intervals.display(verbose_stream(), i); verbose_stream() << "\n";
+        c().print_ineq(ineq(lo_t, lp::EQ, k), verbose_stream()) << "\n";
+
+        new_lemma lemma(c(), "pdd-gcd");
+        add_dependencies(lemma, eq);
+        svector<lp::constraint_index> deps;
+        di.get_lower_dep(i_wd, deps);
+        di.get_upper_dep(i_wd, deps);
+        lp::explanation e(deps);
+        lemma &= e;
+        lemma |= ineq(lo_t, lp::EQ, k);        
+
+        verbose_stream() << lemma << "\n";
+        return true;
+    }
 
     void grobner::add_dependencies(new_lemma& lemma, const dd::solver::equation& eq) {
         lp::explanation ex;
@@ -346,7 +466,8 @@ namespace nla {
                 continue;
             m_rows.insert(row);
             unsigned k = m_lar_solver.get_base_column_in_row(row);
-            if (m_lar_solver.column_is_free(k) && k != j)
+            CTRACE("grobner", m_lar_solver.column_is_free(k) && k != j, tout << "skip " << k << "\n");
+            if (false && m_lar_solver.column_is_free(k) && k != j)
                 continue;
             CTRACE("grobner", matrix.m_rows[row].size() > c().m_nla_settings.grobner_row_length_limit,
                    tout << "ignore the row " << row << " with the size " << matrix.m_rows[row].size() << "\n";); 
diff --git a/src/math/lp/nla_grobner.h b/src/math/lp/nla_grobner.h
index 902ad3a46..a1c877fb4 100644
--- a/src/math/lp/nla_grobner.h
+++ b/src/math/lp/nla_grobner.h
@@ -23,6 +23,7 @@ namespace nla {
         dd::solver               m_solver;
         lp::lar_solver&          m_lar_solver;
         lp::u_set                m_rows;
+        unsigned                 m_quota = 0;
 
         lp::lp_settings& lp_settings();
 
@@ -30,6 +31,9 @@ namespace nla {
         bool is_conflicting();
         bool is_conflicting(const dd::solver::equation& eq);
 
+
+        bool propagate(std::function<bool(dd::solver::equation const&)>& fn);
+        
         bool propagate_bounds();
         bool propagate_bounds(const dd::solver::equation& eq);
 
@@ -38,7 +42,10 @@ namespace nla {
 
         bool propagate_factorization();
         bool propagate_factorization(const dd::solver::equation& eq);
-                
+
+        bool propagate_gcd();
+        bool propagate_gcd(const dd::solver::equation& eq);
+        
         void add_dependencies(new_lemma& lemma, const dd::solver::equation& eq);
 
         // setup
diff --git a/src/math/lp/nla_intervals.cpp b/src/math/lp/nla_intervals.cpp
index 4ffbcb7e3..df5862f62 100644
--- a/src/math/lp/nla_intervals.cpp
+++ b/src/math/lp/nla_intervals.cpp
@@ -274,8 +274,32 @@ void intervals::set_var_interval(lpvar v, interval& b) {
         m_dep_intervals.set_upper_is_inf(b, true);
         if (wd == e_with_deps::with_deps) b.m_upper_dep = nullptr;
     }
+
+    if (ls().column_corresponds_to_term(v)) {
+        auto const& lt = ls().column_to_term(v);
+        scoped_dep_interval ti(m_dep_intervals), r(m_dep_intervals);
+        if (interval_from_lar_term<wd>(lt, ti)) {
+            m_dep_intervals.intersect<wd>(b, ti, r);
+            m_dep_intervals.set<wd>(b, r);
+        }
+    }
 }
 
+template <dep_intervals::with_deps_t wd>
+bool intervals::interval_from_lar_term(lp::lar_term const& t, scoped_dep_interval& i) {
+    m_dep_intervals.set_value(i, rational::zero());
+    scoped_dep_interval vi(m_dep_intervals);
+    for (auto const& arg : t) {
+        set_var_interval<wd>(arg.column().index(), vi);
+        m_dep_intervals.mul<wd>(arg.coeff(), vi, vi);
+        m_dep_intervals.add<wd>(vi, i, i);
+        if (m_dep_intervals.is_inf(i))
+            return false;
+    }
+    return true;
+}
+
+
 template <e_with_deps wd>
 bool intervals::interval_from_term(const nex& e, scoped_dep_interval& i) {
     rational a, b;
@@ -360,7 +384,7 @@ bool intervals::conflict_u_l(const interval& a, const interval& b) const {
 template <e_with_deps wd, typename T>
 bool intervals::interval_of_sum(const nex_sum& e, scoped_dep_interval& a, const std::function<void (const T&)>& f) {
     TRACE("nla_intervals_details", tout << "e=" << e << "\n";);
-    if(! interval_of_sum_no_term<wd>(e, a, f)) {
+    if (!interval_of_sum_no_term<wd>(e, a, f)) {
         return false;
     }
     TRACE("nla_intervals_details", tout << "a = "; display(tout, a););
@@ -379,12 +403,14 @@ bool intervals::interval_of_sum(const nex_sum& e, scoped_dep_interval& a, const
                     if (conflict_u_l(a, i_from_term)) {
                         get_dep_intervals().linearize(a.get().m_upper_dep, expl);
                         get_dep_intervals().linearize(r.get().m_lower_dep, expl);                        
-                    } else {
+                    }
+                    else {
                         get_dep_intervals().linearize(r.get().m_upper_dep, expl);
                         get_dep_intervals().linearize(a.get().m_lower_dep, expl);                        
                     }                        
                     f(expl);
-                } else {
+                }
+                else {
                     // need to recalculate the interval with dependencies
                     scoped_dep_interval sa(get_dep_intervals());
                     interval_of_sum<e_with_deps::with_deps>(e, sa, f);
diff --git a/src/math/lp/nla_intervals.h b/src/math/lp/nla_intervals.h
index 0545e9933..d39f1c8d8 100644
--- a/src/math/lp/nla_intervals.h
+++ b/src/math/lp/nla_intervals.h
@@ -55,6 +55,8 @@ public:
     template <dep_intervals::with_deps_t wd>
     bool interval_from_term(const nex& e, scoped_dep_interval& i); 
 
+    template <dep_intervals::with_deps_t wd>
+    bool interval_from_lar_term(lp::lar_term const& t, scoped_dep_interval& i); 
 
     template <dep_intervals::with_deps_t wd, typename T>
     bool interval_of_sum_no_term(const nex_sum& e, scoped_dep_interval&, const std::function<void (const T&)>& f );