Fstar opt2 (#10261)

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
Copilot-Session: 84b07e32-6458-4ea9-bf14-1cecfb7f1a99

src/math/lp/horner.cpp

@@ -110,6 +110,13 @@ bool horner::horner_lemmas() {
     // so the LP maximization only runs when horner is actually scheduled.
     // optimize_nl_bounds() checks arith.nl.optimize_bounds internally.
     c().optimize_nl_bounds();
+    // optimize_nl_bounds re-calibrated m_to_refine against the model it produced.
+    // If nothing remains to refine, every monomial is consistent under a feasible
+    // LP model: the nonlinear goal is satisfied. Declare it and stop.
+    if (c().to_refine().empty()) {
+        c().set_nla_satisfied();
+        return false;
+    }
     c().lp_settings().stats().m_horner_calls++;
     const auto& matrix = c().lra.A_r();
     // choose only rows that depend on m_to_refine variables

src/math/lp/monomial_bounds.cpp

@@ -677,5 +677,382 @@ namespace nla {
         return new_bound;
     }
 
+    // ================================================================
+    // max_min: incremental LP bound optimization.
+    //
+    // A direct adaptation of smt::theory_arith::max_min (see
+    // src/smt/theory_arith_aux.h).  We maximize (or minimize) a single
+    // column 'v' over the current LP tableau by a bounded-effort primal
+    // simplex walk: repeatedly pick a non-basic variable that improves the
+    // objective, ratio-test its column to find the tightest blocking basic
+    // variable, and pivot.  The tableau is left at a feasible vertex; the
+    // implied bound is then read off 'v's tableau row and rounded to respect
+    // the integrality of integer columns.
+    //
+    // Integrality is maintained during the walk (the 'maintain_integrality ==
+    // true' configuration of theory_arith): every move of a column is a multiple
+    // of the integrality quantum 'min_gain', so integer columns keep integral
+    // values throughout.  The final implied bound is additionally floored/ceiled
+    // for integer 'v'.
+    // ================================================================
+
+    static lp::impq mm_abs(lp::impq const& v) {
+        return v.is_neg() ? -v : v;
+    }
+
+    // Round 'val' down to the nearest multiple of the (integral) 'divisor'.
+    // Mirrors theory_arith::normalize_gain.  'divisor == -1' means "no quantum".
+    static void mm_round_down(lp::impq& val, rational const& divisor) {
+        if (divisor.is_one())
+            val = lp::impq(lp::floor(val));
+        else if (!divisor.is_minus_one())
+            val = lp::impq(lp::floor(val / divisor) * divisor);
+    }
+
+    lpvar monomial_bounds::mm_basic_in_row(unsigned row) const {
+        return c().lra.get_base_column_in_row(row);
+    }
+
+    // A gain is safe when the column is unbounded in the chosen direction, or
+    // the required integral quantum still fits within the maximal feasible move.
+    // Mirrors theory_arith::safe_gain.
+    bool monomial_bounds::mm_safe_gain(mm_gain const& g) const {
+        return g.unbounded || lp::impq(g.min_gain) <= g.max_gain;
+    }
+
+    // Initialize the gain for moving 'x' in direction 'inc' (increase when inc,
+    // decrease otherwise) from its own bound.  For integer columns the quantum
+    // 'min_gain' starts at 1.  Mirrors theory_arith::init_gains.
+    monomial_bounds::mm_gain monomial_bounds::mm_init_gains(lpvar x, bool inc) const {
+        auto& s = c().lra;
+        mm_gain g;
+        if (inc && s.column_has_upper_bound(x)) {
+            g.unbounded = false;
+            g.max_gain = s.column_upper_bound(x) - s.get_column_value(x);
+        }
+        else if (!inc && s.column_has_lower_bound(x)) {
+            g.unbounded = false;
+            g.max_gain = s.get_column_value(x) - s.column_lower_bound(x);
+        }
+        if (s.column_is_int(x))
+            g.min_gain = rational::one();
+        return g;
+    }
+
+    // Tighten 'g' by the room that basic variable 'x_i' (with coefficient 'a_ij'
+    // on the moving column) has before hitting a bound.  When 'x_i' is an integer
+    // column, the quantum 'min_gain' is raised to the lcm of the denominators of
+    // the involved coefficients and both gains are rounded down to that quantum,
+    // so the move keeps 'x_i' integral.  Returns true when 'max_gain' was
+    // strengthened.  Mirrors theory_arith::update_gains.
+    bool monomial_bounds::mm_update_gains(bool inc, lpvar x_i, rational const& a_ij, mm_gain& g) const {
+        auto& s = c().lra;
+        SASSERT(!a_ij.is_zero());
+        if (!mm_safe_gain(g))
+            return false;
+
+        bool decrement_x_i = (inc && a_ij.is_pos()) || (!inc && a_ij.is_neg());
+        bool bounded_i = false;
+        lp::impq max_inc;
+        if (decrement_x_i && s.column_has_lower_bound(x_i)) {
+            max_inc = mm_abs((s.get_column_value(x_i) - s.column_lower_bound(x_i)) / a_ij);
+            bounded_i = true;
+        }
+        else if (!decrement_x_i && s.column_has_upper_bound(x_i)) {
+            max_inc = mm_abs((s.column_upper_bound(x_i) - s.get_column_value(x_i)) / a_ij);
+            bounded_i = true;
+        }
+
+        bool xi_int = s.column_is_int(x_i);
+        rational den_aij(1);
+        if (xi_int)
+            den_aij = denominator(a_ij);
+        SASSERT(den_aij.is_pos() && den_aij.is_int());
+
+        // Moving 'x_i' by k requires moving the entering column by k/a_ij; to keep
+        // an integer 'x_i' integral the entering column must step in multiples of
+        // denominator(a_ij).  Accumulate that into the quantum and re-round.
+        if (xi_int && !den_aij.is_one()) {
+            if (g.min_gain.is_neg())
+                g.min_gain = den_aij;
+            else
+                g.min_gain = lcm(g.min_gain, den_aij);
+            if (!g.unbounded)
+                mm_round_down(g.max_gain, g.min_gain);
+        }
+        if (xi_int && !g.unbounded && !g.max_gain.is_int()) {
+            g.max_gain = lp::impq(lp::floor(g.max_gain));
+            mm_round_down(g.max_gain, g.min_gain);
+        }
+
+        if (bounded_i) {
+            if (xi_int) {
+                max_inc = lp::impq(lp::floor(max_inc));
+                mm_round_down(max_inc, g.min_gain);
+            }
+            if (g.unbounded) {
+                g.unbounded = false;
+                g.max_gain = max_inc;
+                return true;
+            }
+            if (g.max_gain > max_inc) {
+                g.max_gain = max_inc;
+                return true;
+            }
+        }
+        return false;
+    }
+
+    // Ratio test: for entering column 'x_j' moving in direction 'inc', find the
+    // basic variable 'x_i' that first blocks the move and the maximal gain.
+    // Returns false (unsafe) when the integrality quantum cannot be satisfied, so
+    // the caller treats 'x_j' as unusable.  Mirrors theory_arith::pick_var_to_leave.
+    bool monomial_bounds::mm_pick_var_to_leave(lpvar x_j, bool inc, rational& a_ij, mm_gain& g, lpvar& x_i) const {
+        auto& s = c().lra;
+        x_i = null_lpvar;
+        g = mm_init_gains(x_j, inc);
+        // an integer entering column must sit at an integral value to move in
+        // integral steps.
+        if (s.column_is_int(x_j) && !s.get_column_value(x_j).is_int())
+            return false;
+        for (auto const& cell : s.A_r().m_columns[x_j]) {
+            lpvar si = mm_basic_in_row(cell.var());
+            rational const& coeff_ij = s.A_r().get_val(cell);
+            if (mm_update_gains(inc, si, coeff_ij, g) ||
+                (x_i == null_lpvar && !g.unbounded)) {
+                x_i = si;
+                a_ij = coeff_ij;
+            }
+        }
+        return mm_safe_gain(g);
+    }
+
+    // Apply 'delta' to non-basic column 'j', propagating to dependent basic
+    // columns (theory_arith::update_value).
+    void monomial_bounds::mm_update_value(lpvar j, lp::impq const& delta) {
+        if (delta.is_zero())
+            return;
+        auto& s = c().lra;
+        lp::impq new_val = s.get_column_value(j) + delta;
+        s.set_value_for_nbasic_column_report(j, new_val, [](unsigned) {});
+    }
+
+    // Move (now non-basic) 'x_i' maximally towards its bound in direction 'inc'
+    // without violating other columns' bounds, in integral steps when 'x_i' is an
+    // integer column (theory_arith::move_to_bound).
+    bool monomial_bounds::mm_move_to_bound(lpvar x_i, bool inc, unsigned& best_efforts) {
+        auto& s = c().lra;
+        if (s.column_is_int(x_i) && !s.get_column_value(x_i).is_int()) {
+            ++best_efforts;
+            return false;
+        }
+        mm_gain g = mm_init_gains(x_i, inc);
+        for (auto const& cell : s.A_r().m_columns[x_i]) {
+            lpvar si = mm_basic_in_row(cell.var());
+            rational const& coeff = s.A_r().get_val(cell);
+            mm_update_gains(inc, si, coeff, g);
+        }
+        bool result = false;
+        if (mm_safe_gain(g) && !g.unbounded) {
+            lp::impq step = g.max_gain;
+            if (!inc)
+                step = -step;
+            mm_update_value(x_i, step);
+            result = !g.max_gain.is_zero();
+        }
+        if (!result)
+            ++best_efforts;
+        return result;
+    }
+
+    // Primal-simplex walk maximizing/minimizing 'v' (theory_arith::max_min).
+    void monomial_bounds::mm_optimize(lpvar v, bool maximize) {
+        auto& s = c().lra;
+        unsigned best_efforts = 0;
+        unsigned const max_efforts = 20;
+        unsigned rounds = 0;
+        unsigned const max_rounds = 200;
+
+        while (best_efforts < max_efforts && rounds < max_rounds && !c().lp_settings().get_cancel_flag()) {
+            ++rounds;
+            lpvar x_j = null_lpvar, x_i = null_lpvar;
+            rational a_ij(0);
+            mm_gain best;          // gain of the selected move
+            bool inc = false;
+            bool has_bound = false;
+
+            // Consider a candidate entering variable 'cand' whose coefficient in
+            // the objective (v expressed over the non-basic columns) is
+            // 'obj_coeff'.  Returns true to stop scanning (unbounded direction).
+            auto consider = [&](lpvar cand, rational const& obj_coeff) -> bool {
+                bool curr_inc = obj_coeff.is_pos() ? maximize : !maximize;
+                if ((curr_inc && s.column_has_upper_bound(cand)) ||
+                    (!curr_inc && s.column_has_lower_bound(cand)))
+                    has_bound = true;
+                // cannot move a variable already at the relevant bound
+                if (curr_inc && s.column_has_upper_bound(cand) &&
+                    s.get_column_value(cand) == s.column_upper_bound(cand))
+                    return false;
+                if (!curr_inc && s.column_has_lower_bound(cand) &&
+                    s.get_column_value(cand) == s.column_lower_bound(cand))
+                    return false;
+                rational curr_a(0);
+                mm_gain cur;
+                lpvar curr_xi = null_lpvar;
+                bool safe = mm_pick_var_to_leave(cand, curr_inc, curr_a, cur, curr_xi);
+                if (!safe) {
+                    // the integrality quantum cannot be met on this column
+                    has_bound = true;
+                    ++best_efforts;
+                    return false;
+                }
+                if (curr_xi == null_lpvar) {
+                    // limited only by its own bound (or fully unbounded)
+                    x_j = cand; x_i = null_lpvar; inc = curr_inc; best = cur; a_ij = curr_a;
+                    return true;
+                }
+                if (cur.max_gain > best.max_gain) {
+                    x_i = curr_xi; x_j = cand; a_ij = curr_a; best = cur; inc = curr_inc;
+                }
+                else if (cur.max_gain.is_zero() && (x_i == null_lpvar || curr_xi < x_i)) {
+                    x_i = curr_xi; x_j = cand; a_ij = curr_a; best = cur; inc = curr_inc;
+                }
+                return false;
+            };
+
+            if (!s.is_base(v)) {
+                consider(v, rational::one());
+            }
+            else {
+                unsigned ri = s.r_heading()[v];
+                rational a_v(0);
+                for (auto const& e : s.A_r().m_rows[ri])
+                    if (e.var() == v) { a_v = e.coeff(); break; }
+                for (auto const& e : s.A_r().m_rows[ri]) {
+                    if (e.var() == v)
+                        continue;
+                    // v = -(1/a_v) * sum a_e x_e, so d(v)/d(x_e) has the sign of
+                    // -a_e/a_v; only the sign steers the search direction.
+                    rational objc = -e.coeff();
+                    if (a_v.is_neg())
+                        objc.neg();
+                    if (consider(e.var(), objc))
+                        break;
+                }
+            }
+
+            if (!has_bound && x_i == null_lpvar && x_j == null_lpvar)
+                return; // objective is unbounded in the chosen direction
+            if (x_j == null_lpvar)
+                return; // optimized: no improving move remains
+
+            // a non-unit integral quantum means the exact optimum may not be
+            // reachable in integral steps: count it as best-effort progress.
+            if (best.min_gain.is_pos() && !best.min_gain.is_one())
+                ++best_efforts;
+
+            if (x_i == null_lpvar) {
+                // move x_j directly to its own bound
+                if (inc && s.column_has_upper_bound(x_j)) {
+                    if (best.max_gain.is_zero())
+                        return;
+                    mm_update_value(x_j, best.max_gain);
+                    continue;
+                }
+                if (!inc && s.column_has_lower_bound(x_j)) {
+                    if (best.max_gain.is_zero())
+                        return;
+                    mm_update_value(x_j, -best.max_gain);
+                    continue;
+                }
+                return; // unbounded
+            }
+
+            // x_j can move exactly across to its opposite bound without pivoting
+            if (s.column_has_lower_bound(x_j) && s.column_has_upper_bound(x_j) &&
+                s.column_lower_bound(x_j) != s.column_upper_bound(x_j) &&
+                (s.column_upper_bound(x_j) - s.column_lower_bound(x_j) == best.max_gain)) {
+                lp::impq step = best.max_gain;
+                if (!inc)
+                    step = -step;
+                mm_update_value(x_j, step);
+                continue;
+            }
+
+            // pivot x_j into the basis (x_i leaves); the degenerate pivot keeps
+            // the current point, then move x_i to its bound to raise v.
+            s.pivot(x_j, x_i);
+            bool inc_xi = inc ? a_ij.is_neg() : a_ij.is_pos();
+            mm_move_to_bound(x_i, inc_xi, best_efforts);
+        }
+    }
+
+    // Read the implied bound on 'v' off its final tableau row and round it to
+    // respect the integrality of integer columns (theory_arith::mk_bound_from_row
+    // + normalize_bound).  Returns the joined explanation, or nullptr if no bound
+    // is implied (e.g. a required bound on a row variable is missing).
+    u_dependency* monomial_bounds::mm_bound_from_row(lpvar v, bool maximize, rational& bound) {
+        auto& s = c().lra;
+        if (!s.is_base(v))
+            return nullptr;
+        unsigned ri = s.r_heading()[v];
+        auto const& row = s.A_r().m_rows[ri];
+        rational a_v(0);
+        for (auto const& e : row)
+            if (e.var() == v) { a_v = e.coeff(); break; }
+        if (a_v.is_zero())
+            return nullptr;
+        lp::impq acc(0);
+        u_dependency* dep = nullptr;
+        for (auto const& e : row) {
+            if (e.var() == v)
+                continue;
+            lpvar k = e.var();
+            rational ck = -e.coeff() / a_v; // v = sum ck * x_k
+            if (ck.is_zero())
+                continue;
+            bool use_upper = maximize ? ck.is_pos() : ck.is_neg();
+            if (use_upper) {
+                if (!s.column_has_upper_bound(k))
+                    return nullptr;
+                acc += s.column_upper_bound(k) * ck;
+                dep = s.join_deps(dep, s.get_column_upper_bound_witness(k));
+            }
+            else {
+                if (!s.column_has_lower_bound(k))
+                    return nullptr;
+                acc += s.column_lower_bound(k) * ck;
+                dep = s.join_deps(dep, s.get_column_lower_bound_witness(k));
+            }
+        }
+        if (s.column_is_int(v))
+            bound = maximize ? lp::floor(acc) : lp::ceil(acc);
+        else
+            bound = acc.x;
+        return dep;
+    }
+
+    u_dependency* monomial_bounds::improve_bound(lpvar j, bool is_lower, rational& bound) {
+        auto& s = c().lra;
+        if (!s.is_feasible())
+            return nullptr;
+        bool maximize = !is_lower;
+        mm_optimize(j, maximize);
+        rational b(0);
+        u_dependency* dep = mm_bound_from_row(j, maximize, b);
+        if (!dep)
+            return nullptr;
+        if (is_lower) {
+            if (s.column_has_lower_bound(j) && b <= s.column_lower_bound(j).x)
+                return nullptr;
+        }
+        else {
+            if (s.column_has_upper_bound(j) && b >= s.column_upper_bound(j).x)
+                return nullptr;
+        }
+        bound = b;
+        return dep;
+    }
+
 }
 

src/math/lp/monomial_bounds.h

@@ -51,10 +51,48 @@ namespace nla {
         bool propagate_changed_bound(monic & m);
         bool is_linear(monic const& m, lpvar& w, lpvar & fixed_to_zero);
         rational fixed_var_product(monic const& m, lpvar w);
+
+        // ----------------------------------------------------------------
+        // max_min: incremental LP bound optimization.
+        //
+        // Adapted from smt::theory_arith::max_min (theory_arith_aux.h): a
+        // bounded-variable primal-simplex walk that maximizes/minimizes a
+        // single column over the current LP tableau by pivoting, leaving the
+        // tableau at an optimal feasible vertex.  Replaces the expensive
+        // lar_solver::find_improved_bound (which rebuilds the objective and
+        // re-solves from scratch on every call) and rounds the resulting
+        // bound to respect the integrality of integer columns.
+        // ----------------------------------------------------------------
+        //
+        // Gain of a candidate move: how far a non-basic column may travel in the
+        // chosen direction.  'min_gain' is the integrality quantum -- the
+        // smallest integral step that keeps every dependent integer column at an
+        // integral value (-1 means "no quantum", i.e. a rational column);
+        // 'max_gain' is the largest feasible step (valid only when !unbounded).
+        //
+        struct mm_gain {
+            bool      unbounded = true;              // no bound blocks the move
+            rational  min_gain  = rational::minus_one();
+            lp::impq  max_gain;
+        };
+        lpvar mm_basic_in_row(unsigned row) const;
+        bool mm_safe_gain(mm_gain const& g) const;
+        mm_gain mm_init_gains(lpvar x, bool inc) const;
+        bool mm_update_gains(bool inc, lpvar x_i, rational const& a_ij, mm_gain& g) const;
+        bool mm_pick_var_to_leave(lpvar x_j, bool inc, rational& a_ij, mm_gain& g, lpvar& x_i) const;
+        bool mm_move_to_bound(lpvar x_i, bool inc, unsigned& best_efforts);
+        void mm_update_value(lpvar j, lp::impq const& delta);
+        void mm_optimize(lpvar v, bool maximize);
+        u_dependency* mm_bound_from_row(lpvar v, bool maximize, rational& bound);
     public:
         monomial_bounds(core* core);
         void generate_lemmas();
         bool tighten_lp_bounds();
         bool propagate_changed_bounds();
+
+        // Maximize (is_lower == false) or minimize (is_lower == true) column j
+        // over the LP tableau and, if the resulting bound improves j's current
+        // bound, return its explanation and set 'bound'; otherwise return nullptr.
+        u_dependency* improve_bound(lpvar j, bool is_lower, rational& bound);
     }; 
 }

src/math/lp/nla_core.cpp

@@ -1299,13 +1299,14 @@ lbool core::check(unsigned level) {
     if (m_to_refine.empty())
         return l_true;    
     init_search();
+    m_nla_satisfied = false;
 
     lbool ret = l_undef;
     bool run_grobner = need_run_grobner();
     bool run_horner = need_run_horner();
     bool run_bounds = params().arith_nl_branching();
 
-    auto no_effect = [&]() { return ret == l_undef && !done() && m_lemmas.empty() && m_literals.empty() && !m_check_feasible; };
+    auto no_effect = [&]() { return ret == l_undef && !done() && !m_nla_satisfied && m_lemmas.empty() && m_literals.empty() && !m_check_feasible; };
     
     if (no_effect())
         m_monomial_bounds.generate_lemmas();
@@ -1329,6 +1330,9 @@ lbool core::check(unsigned level) {
             return l_undef;
         if (!m_lemmas.empty() || !m_literals.empty() || m_check_feasible)
             return l_false;
+        // bound optimization proved all monomials consistent: goal satisfied.
+        if (m_nla_satisfied)
+            return l_true;
     }
 
     if (no_effect() && params().arith_nl_nra_check_assignment() && m_check_assignment_fail_cnt < params().arith_nl_nra_check_assignment_max_fail()) {
@@ -1563,8 +1567,11 @@ bool core::optimize_nl_bounds() {
 
     if (!lra.is_feasible())
         return false;
-    if (lra.find_feasible_solution() == lp::lp_status::INFEASIBLE)
+    if (lra.find_feasible_solution() == lp::lp_status::INFEASIBLE) {
+        // find_feasible_solution moved the model; keep m_to_refine in sync.
+        init_to_refine();
         return false;
+    }
 
     // Gather the candidate columns: every non-fixed leaf variable that
     // participates in a monomial (mirrors solver=2's max_min_nl_vars).
@@ -1605,7 +1612,7 @@ bool core::optimize_nl_bounds() {
             break;
         for (bool is_lower : { true, false }) {
             rational bound;
-            u_dependency* dep = lra.find_improved_bound(j, is_lower, bound);
+            u_dependency* dep = m_monomial_bounds.improve_bound(j, is_lower, bound);
             if (!dep)
                 continue;
             auto kind = is_lower ? lp::lconstraint_kind::GE : lp::lconstraint_kind::LE;
@@ -1613,12 +1620,21 @@ bool core::optimize_nl_bounds() {
         }
     }
 
-    if (improvements.empty())
+    if (improvements.empty()) {
+        // The exploratory simplex walk in improve_bound/mm_optimize mutated the
+        // LP model even though no bound was tightened.  Restore a clean feasible
+        // model and re-calibrate m_to_refine so downstream lemma passes (grobner,
+        // basic_lemma) never see a stale monomial that is now consistent.
+        lra.find_feasible_solution();
+        init_to_refine();
         return false;
+    }
 
     for (auto const& ib : improvements)
         lra.update_column_type_and_bound(ib.j, ib.kind, ib.bound, ib.dep);
     lra.find_feasible_solution();
+    // The model changed: re-calibrate m_to_refine against the new assignment.
+    init_to_refine();
     return true;
 }
 

src/math/lp/nla_core.h

@@ -90,6 +90,9 @@ class core {
     monomial_bounds          m_monomial_bounds;
     unsigned                 m_conflicts;
     bool                     m_check_feasible = false;
+    // set when bound optimization re-calibrates m_to_refine to empty: every
+    // monomial is consistent under the optimized model, so the goal is satisfied.
+    bool                     m_nla_satisfied = false;
     horner                   m_horner;
     grobner                  m_grobner;
     emonics                  m_emons;
@@ -466,6 +469,9 @@ public:
         return m_to_refine;
     }
 
+    void set_nla_satisfied() { m_nla_satisfied = true; }
+    bool nla_satisfied() const { return m_nla_satisfied; }
+
 };  // end of core
 
 struct pp_mon {

src/math/lp/nla_grobner.cpp

@@ -719,6 +719,10 @@ namespace nla {
                 TRACE(grobner, m_solver.display(tout << "horner conflict ", e) << "\n");
                 return true;
             }
+            if (c().params().arith_nl_grobner_perfect_squares() && add_perfect_square_conflict(e)) {
+                TRACE(grobner, m_solver.display(tout << "perfect-square conflict ", e) << "\n");
+                return true;
+            }
             return false;
         }
         evali.get_interval<dd::w_dep::with_deps>(e.poly(), i_wd);  
@@ -736,6 +740,156 @@ namespace nla {
         }
     }
 
+    // Helpers for perfect-square conflict detection (port of theory_arith_nl's
+    // is_inconsistent2). Variable lists are assumed sorted ascending.
+
+    // Return true and set root=sqrt(coeff) if (coeff, vars) is a perfect square
+    // a^2 * M^2 (coeff a perfect square and every variable of even multiplicity).
+    static bool ps_root(rational const& coeff, unsigned_vector const& vars, rational& root) {
+        if (vars.size() % 2 == 1)
+            return false;
+        if (!coeff.is_perfect_square(root))
+            return false;
+        unsigned i = 0, n = vars.size();
+        while (i < n) {
+            unsigned v = vars[i];
+            unsigned power = 1;
+            ++i;
+            for (; i < n && vars[i] == v; ++i)
+                ++power;
+            if (power % 2 == 1)
+                return false;
+        }
+        return true;
+    }
+
+    // Return true if v12/c12 is the cross term (-2ab)*M1*M2, given that
+    // (v1,a^2) and (v2,b^2) are perfect squares a^2*M1^2 and b^2*M2^2.
+    static bool ps_cross(rational const& a, unsigned_vector const& v1,
+                         rational const& b, unsigned_vector const& v2,
+                         rational const& c12, unsigned_vector const& v12) {
+        if (!c12.is_neg())
+            return false;
+        rational c(-2);
+        c *= a;
+        c *= b;
+        if (c12 != c)
+            return false;
+        unsigned n1 = v1.size(), n2 = v2.size(), n12 = v12.size();
+        if (n1 + n2 != n12 * 2)
+            return false;
+        unsigned i1 = 0, i2 = 0, i12 = 0;
+        while (true) {
+            bool e1 = i1 < n1, e2 = i2 < n2, e12 = i12 < n12;
+            if (!e1 && !e2 && !e12)
+                return true;
+            if (!e12)
+                return false;
+            unsigned x12 = v12[i12];
+            if (e1 && v1[i1] == x12) { i1 += 2; ++i12; }
+            else if (e2 && v2[i2] == x12) { i2 += 2; ++i12; }
+            else return false;
+        }
+    }
+
+    // r <- coeff * prod(var^power) over the (sorted) variable list.
+    template <dep_intervals::with_deps_t wd>
+    void grobner::monomial_interval(rational const& coeff, unsigned_vector const& vars, scoped_dep_interval& r) {
+        auto& di = c().m_intervals.get_dep_intervals();
+        di.set_value(r, coeff);
+        unsigned i = 0, n = vars.size();
+        while (i < n) {
+            lpvar v = vars[i];
+            unsigned power = 1;
+            ++i;
+            for (; i < n && vars[i] == v; ++i)
+                ++power;
+            scoped_dep_interval vi(di), vp(di), tmp(di);
+            c().m_intervals.set_var_interval<wd>(v, vi);
+            di.power<wd>(vi, power, vp);
+            di.mul<wd>(r, vp, tmp);
+            di.set<wd>(r, tmp);
+        }
+    }
+
+    // Perfect-square conflict detection: within a grobner equation e = 0, find a
+    // triple of monomials a^2*M1^2, b^2*M2^2, -2ab*M1*M2 forming (a*M1 - b*M2)^2.
+    // Replace such a triple by [0, oo) when that improves the lower bound, then
+    // check whether the resulting interval of e separates from zero.
+    bool grobner::add_perfect_square_conflict(const dd::solver::equation& e) {
+        vector<rational> coeffs;
+        vector<unsigned_vector> varss;
+        for (auto const& [coeff, vars] : e.poly()) {
+            unsigned_vector vs(vars);
+            std::sort(vs.begin(), vs.end());
+            coeffs.push_back(coeff);
+            varss.push_back(vs);
+        }
+        unsigned num = coeffs.size();
+        if (num < 3)
+            return false;
+        auto& di = c().m_intervals.get_dep_intervals();
+        svector<bool> deleted;
+        deleted.resize(num, false);
+        bool found = false;
+        for (unsigned i = 0; i < num; ++i) {
+            if (deleted[i])
+                continue;
+            rational a;
+            if (!ps_root(coeffs[i], varss[i], a))
+                continue;
+            for (unsigned j = i + 1; j < num && !deleted[i]; ++j) {
+                if (deleted[j])
+                    continue;
+                rational b;
+                if (!ps_root(coeffs[j], varss[j], b))
+                    continue;
+                for (unsigned k = i + 1; k < num; ++k) {
+                    if (k == j || deleted[k])
+                        continue;
+                    if (!ps_cross(a, varss[i], b, varss[j], coeffs[k], varss[k]))
+                        continue;
+                    // does [0,oo) improve on the summed interval of the triple?
+                    scoped_dep_interval Ii(di), Ij(di), Ik(di), s1(di), s2(di);
+                    monomial_interval<dd::w_dep::without_deps>(coeffs[i], varss[i], Ii);
+                    monomial_interval<dd::w_dep::without_deps>(coeffs[j], varss[j], Ij);
+                    monomial_interval<dd::w_dep::without_deps>(coeffs[k], varss[k], Ik);
+                    di.add<dd::w_dep::without_deps>(Ii, Ij, s1);
+                    di.add<dd::w_dep::without_deps>(s1, Ik, s2);
+                    if (di.lower_is_inf(s2) || rational(di.lower(s2)).is_neg()) {
+                        deleted[i] = deleted[j] = deleted[k] = true;
+                        found = true;
+                        break;
+                    }
+                }
+            }
+        }
+        if (!found)
+            return false;
+        // Rebuild the interval of e with the perfect squares replaced by [0, oo)
+        // (an unconditional lower bound, hence no dependency), tracking deps of
+        // the remaining monomials, and check whether it separates from zero.
+        scoped_dep_interval sum(di);
+        di.reset(sum);
+        di.set_lower_is_inf(sum, false);
+        di.set_lower(sum, rational::zero());
+        di.set_lower_is_open(sum, false);
+        di.set_upper_is_inf(sum, true);
+        for (unsigned i = 0; i < num; ++i) {
+            if (deleted[i])
+                continue;
+            scoped_dep_interval Ii(di), tmp(di);
+            monomial_interval<dd::w_dep::with_deps>(coeffs[i], varss[i], Ii);
+            di.add<dd::w_dep::with_deps>(sum, Ii, tmp);
+            di.set<dd::w_dep::with_deps>(sum, tmp);
+        }
+        std::function<void (const lp::explanation&)> f = [this](const lp::explanation& expl) {
+            lemma_builder lemma(m_core, "pdd-perfect-square");
+            lemma &= expl;
+        };
+        return di.check_interval_for_conflict_on_zero(sum, e.dep(), f);
+    }
+
     bool grobner::propagate_linear_equations() {
         unsigned changed = 0;
         m_mon2var.clear();

src/math/lp/nla_grobner.h

@@ -75,6 +75,12 @@ namespace nla {
         bool add_nla_conflict(dd::solver::equation const& eq);
         void check_missing_propagation(dd::solver::equation const& eq);
 
+        // perfect-square conflict detection (arith.nl.grobner_perfect_squares),
+        // a port of theory_arith_nl's is_inconsistent2.
+        bool add_perfect_square_conflict(dd::solver::equation const& eq);
+        template <dep_intervals::with_deps_t wd>
+        void monomial_interval(rational const& coeff, unsigned_vector const& vars, scoped_dep_interval& r);
+
         bool equation_is_true(dd::solver::equation const& eq);
 
         // adaptive growth (gated by arith.nl.grobner_adaptive)

src/params/smt_params_helper.pyg

@@ -93,6 +93,7 @@ def_module_params(module_name='smt',
                           ('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.grobner_expand_terms', BOOL, True, 'expand terms before computing grobner basis'),
+                          ('arith.nl.grobner_perfect_squares', BOOL, True, 'expand perfect squares with Grobner'),
                           ('arith.nl.monomial_sandwich', BOOL, False, 'derive bound on a monomial factor by pairing two LP rows that share the other factor'),
                           ('arith.nl.monomial_sandwich.max_fanout', UINT, 0, 'skip monomial sandwich when the conclusion factor appears in more than this many monomials (0 = no limit)'),
                           ('arith.nl.monomial_binomial_sign', BOOL, False, 'derive bound on a binomial-monomial factor anchored on the current LP value of the monomial; replaces order_lemma_on_binomial_sign with a deterministic factor bound conditioned on a one-sided snapshot of the monomial value'),