mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2026-07-18 04:55:45 +00:00
Code Activity

Start resolving merge conflicts with master branch

Browse source
This commit is contained in:
copilot-swe-agent[bot] 2026-07-15 04:14:50 +00:00 • committed by GitHub
commit dcf81a2592
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
1020 changed files with 80693 additions and 20317 deletions
Download patch file Download diff file Expand all files Collapse all files

2
src/math/grobner/pdd_simplifier.cpp
View file

@ -494,7 +494,7 @@ namespace dd {
hash(unsigned_vector& vars):vars(vars) {}
bool operator()(mon const& m) const {
return unsigned_ptr_hash(vars.data() + m.offset, m.sz, 1);
};
}
};
struct eq {
unsigned_vector& vars;

4
src/math/lp/bound_analyzer_on_row.h
View file

@ -297,7 +297,7 @@ namespace lp {
void limit_j(unsigned bound_j, const mpq& u, bool coeff_before_j_is_pos, bool is_lower_bound, bool strict) {
auto* lar = &m_bp.lp();
const auto& row = this->m_row;
auto* row = &this->m_row;
auto explain = [row, bound_j, coeff_before_j_is_pos, is_lower_bound, strict, lar]() {
(void) strict;
TRACE(bound_analyzer, tout << "explain_bound_on_var_on_coeff, bound_j = " << bound_j << ", coeff_before_j_is_pos = " << coeff_before_j_is_pos << ", is_lower_bound = " << is_lower_bound << ", strict = " << strict << "\n";);
@ -305,7 +305,7 @@ namespace lp {
int j_sign = (coeff_before_j_is_pos ? 1 : -1) * bound_sign;
u_dependency* ret = nullptr;
for (auto const& r : row) {
for (auto const& r : *row) {
unsigned j = r.var();
if (j == bound_j)
continue;

17
src/math/lp/dioph_eq.cpp
View file

@ -238,18 +238,13 @@ namespace lp {
r.c() -= b.c();
return r;
}
#if Z3DEBUG
friend bool operator==(const term_o& a, const term_o& b) {
friend bool eq(const term_o& a, const term_o& b) {
term_o t = a.clone();
t += mpq(-1) * b;
return t.c() == mpq(0) && t.size() == 0;
}
friend bool operator!=(const term_o& a, const term_o& b) {
return ! (a == b);
}
#endif
term_o& operator+=(const term_o& t) {
for (const auto& p : t) {
add_monomial(p.coeff(), p.j());
@ -585,7 +580,7 @@ namespace lp {
const lar_term* m_t;
undo_add_term(imp& s, const lar_term* t) : m_s(s), m_t(t) {}
void undo() {
void undo() override {
m_s.undo_add_term_method(m_t);
}
};
@ -714,8 +709,8 @@ namespace lp {
while (column.size() > 1) {
auto& c = column.back();
SASSERT(c.var() != last_row_index);
m_l_matrix.pivot_row_to_row_given_cell(last_row_index, c, j);
m_changed_rows.insert(c.var());
m_l_matrix.pivot_row_to_row_given_cell(last_row_index, c, j);
}
}
@ -1541,7 +1536,7 @@ namespace lp {
term_o t1 = open_ml(t0);
t1.add_monomial(mpq(1), j);
term_o rs = fix_vars(t1);
if (ls != rs) {
if (!eq(ls, rs)) {
TRACE(dio, tout << "ls:"; print_term_o(ls, tout) << "\n";
tout << "rs:"; print_term_o(rs, tout) << "\n";);
return false;
@ -2351,7 +2346,7 @@ namespace lp {
return false;
}
bool ret = ls == fix_vars(open_ml(m_l_matrix.m_rows[ei]));
bool ret = eq(ls, fix_vars(open_ml(m_l_matrix.m_rows[ei])));
if (!ret) {
CTRACE(dio, !ret,
{

42
src/math/lp/hnf.h
View file

@ -29,6 +29,7 @@ Revision History:
#pragma once
#include "math/lp/numeric_pair.h"
#include "util/ext_gcd.h"
#include <functional>
namespace lp {
namespace hnf_calc {
@ -128,16 +129,21 @@ bool prepare_pivot_for_lower_triangle(M &m, unsigned r) {
template <typename M>
void pivot_column_non_fractional(M &m, unsigned r, bool & overflow, const mpq & big_number) {
SASSERT(!is_zero(m[r][r]));
// rows <= r are not written below, so these pivot entries are loop-invariant
const auto & mrr = m[r][r];
const mpq * denom = r > 0 ? &m[r - 1][r - 1] : nullptr;
for (unsigned j = r + 1; j < m.column_count(); ++j) {
for (unsigned i = r + 1; i < m.row_count(); ++i) {
if (
(m[i][j] = (r > 0) ? (m[r][r]*m[i][j] - m[i][r]*m[r][j]) / m[r-1][r-1] :
(m[r][r]*m[i][j] - m[i][r]*m[r][j]))
>= big_number) {
const auto & mrj = m[r][j];
for (unsigned i = r + 1; i < m.row_count(); ++i) {
auto & mij = m[i][j];
mij = mrr * mij - m[i][r] * mrj;
if (denom)
mij /= *denom;
if (mij >= big_number) {
overflow = true;
return;
}
SASSERT(is_integer(m[i][j]));
SASSERT(is_integer(mij));
}
}
}
@ -221,6 +227,8 @@ class hnf {
unsigned m_j;
mpq m_R;
mpq m_half_R;
bool m_cancelled = false;
std::function<bool()> m_cancel_flag;
mpq mod_R_balanced(const mpq & a) const {
mpq t = a % m_R;
return t > m_half_R? t - m_R : (t < - m_half_R? t + m_R : t);
@ -574,6 +582,10 @@ private:
void calculate_by_modulo() {
for (m_i = 0; m_i < m_m; m_i ++) {
if (m_cancel_flag && m_cancel_flag()) {
m_cancelled = true;
return;
}
process_row_modulo();
SASSERT(is_pos(m_W[m_i][m_i]));
m_R /= m_W[m_i][m_i];
@ -583,7 +595,7 @@ private:
}
public:
hnf(M & A, const mpq & d) :
hnf(M & A, const mpq & d, std::function<bool()> cancel_flag = nullptr) :
#ifdef Z3DEBUG
m_H(A),
m_A_orig(A),
@ -594,7 +606,8 @@ public:
m_n(A.column_count()),
m_d(d),
m_R(m_d),
m_half_R(floor(m_R / 2))
m_half_R(floor(m_R / 2)),
m_cancel_flag(cancel_flag)
{
if (m_m == 0 || m_n == 0 || is_zero(m_d))
return;
@ -605,15 +618,18 @@ public:
#endif
calculate_by_modulo();
#ifdef Z3DEBUG
CTRACE(hnf_calc, m_H != m_W,
tout << "A = "; m_A_orig.print(tout, 4); tout << std::endl;
tout << "H = "; m_H.print(tout, 4); tout << std::endl;
tout << "W = "; m_W.print(tout, 4); tout << std::endl;);
SASSERT (m_H == m_W);
if (!m_cancelled) {
CTRACE(hnf_calc, m_H != m_W,
tout << "A = "; m_A_orig.print(tout, 4); tout << std::endl;
tout << "H = "; m_H.print(tout, 4); tout << std::endl;
tout << "W = "; m_W.print(tout, 4); tout << std::endl;);
SASSERT (m_H == m_W);
}
#endif
}
const M & W() const { return m_W; }
bool is_cancelled() const { return m_cancelled; }
};

4
src/math/lp/hnf_cutter.cpp
View file

@ -199,7 +199,9 @@ branch y_i >= ceil(y0_i) is impossible.
shrink_explanation(basis_rows);
}
hnf<general_matrix> h(m_A, d);
hnf<general_matrix> h(m_A, d, [this]() { return m_settings.get_cancel_flag(); });
if (h.is_cancelled())
return lia_move::undef;
vector<mpq> b = create_b(basis_rows);
#ifdef Z3DEBUG
SASSERT(m_A * x0 == b);

273
src/math/lp/int_cube.cpp
View file

@ -16,6 +16,9 @@ Author:
Revision History:
--*/
#include <algorithm>
#include <unordered_map>
#include "math/lp/int_solver.h"
#include "math/lp/lar_solver.h"
#include "math/lp/int_cube.h"
@ -81,32 +84,264 @@ namespace lp {
SASSERT(lp_status::OPTIMAL == lra.get_status() || lp_status::FEASIBLE == lra.get_status());
}
impq int_cube::get_cube_delta_for_term(const lar_term& t) const {
if (t.size() == 2) {
bool seen_minus = false;
bool seen_plus = false;
for(lar_term::ival p : t) {
if (!lia.column_is_int(p.j()))
goto usual_delta;
const mpq & c = p.coeff();
if (c == one_of_type<mpq>()) {
seen_plus = true;
} else if (c == -one_of_type<mpq>()) {
seen_minus = true;
} else {
goto usual_delta;
// The largest cube test of Bromberger and Weidenbach:
// maximize x_e subject to Ax + a'(x_e/2) <= b, x_e >= 0, where a'_i = ||a_i||_1,
// with the 1-norm taken over the integer variables of the row.
// The solution is the center z of a largest cube contained in the polyhedron.
// If the maximal edge length is at least 1, then the rounding of z is
// an integer solution; otherwise the rounding is checked, and possibly repaired,
// against the original constraints.
lia_move int_cube::find_largest_cube() {
lia.settings().stats().m_lcube_calls++;
TRACE(cube,
for (unsigned j = 0; j < lra.number_of_vars(); ++j)
lia.display_column(tout, j);
tout << lra.constraints();
);
lra.push();
// The edge rows are ephemeral: suppress the add-term callback,
// dioph_eq's reaction to it is not undone by pop().
auto add_term_cb = lra.m_add_term_callback;
lra.m_add_term_callback = nullptr;
unsigned x_e = lra.add_var(UINT_MAX, false); // the edge length of the cube
lra.add_var_bound(x_e, lconstraint_kind::GE, mpq(0));
bool ok = add_cube_edge_rows(x_e);
lra.m_add_term_callback = add_term_cb;
if (!ok) {
lra.pop();
lra.set_status(lp_status::OPTIMAL);
return lia_move::undef;
}
lp_status st = lra.find_feasible_solution();
if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) {
TRACE(cube, tout << "cannot find a feasible solution";);
lra.pop();
lra.move_non_basic_columns_to_bounds();
// it can happen that we found an integer solution here
return !lra.r_basis_has_inf_int()? lia_move::sat: lia_move::undef;
}
impq e; // the maximal edge length
st = lra.maximize_term(x_e, e, /*fix_int_cols*/ false);
if (lia.settings().get_cancel_flag()) {
lra.pop();
return lia_move::undef;
}
if (st == lp_status::UNBOUNDED) {
// infinite lattice width: the polyhedron contains cubes of arbitrary edge length
lra.add_var_bound(x_e, lconstraint_kind::GE, mpq(1));
st = lra.find_feasible_solution();
if (st != lp_status::FEASIBLE && st != lp_status::OPTIMAL) {
lra.pop();
return lia_move::undef;
}
lra.pop();
return sat_after_rounding();
}
TRACE(cube, tout << "max edge length = " << e << "\n";);
if (e >= impq(mpq(1))) {
lra.pop();
return sat_after_rounding();
}
// the largest cube is smaller than the unit cube:
// the rounded center is only a candidate
lra.pop();
return round_and_repair();
}
bool int_cube::add_cube_edge_rows(unsigned x_e) {
// snapshot the term columns: add_edge_rows_for_term appends to lra.terms()
svector<unsigned> term_columns;
for (const lar_term* t : lra.terms())
term_columns.push_back(t->j());
for (unsigned j : term_columns)
if (!add_edge_rows_for_term(j, x_e)) {
TRACE(cube, tout << "cannot add the edge rows";);
return false;
}
return true;
}
// i is the column index having the term
bool int_cube::add_edge_rows_for_term(unsigned i, unsigned x_e) {
if (!lra.column_associated_with_row(i))
return true;
const lar_term& t = lra.get_term(i);
impq delta = get_cube_delta_for_term(t);
TRACE(cube, lra.print_term_as_indices(t, tout); tout << ", delta = " << delta << "\n";);
if (is_zero(delta))
return true;
if (!is_zero(delta.y))
// the infinitesimal delta does not scale with x_e: tighten statically,
// it is sound for any edge length
return lra.tighten_term_bounds_by_delta(i, delta);
if (lra.column_has_upper_bound(i)) {
impq u = lra.get_upper_bound(i); // copy: add_term invalidates bound references
vector<std::pair<mpq, unsigned>> coeffs = {{mpq(1), i}, {delta.x, x_e}};
unsigned s = lra.add_term(coeffs, UINT_MAX);
lra.add_var_bound(s, is_zero(u.y) ? lconstraint_kind::LE : lconstraint_kind::LT, u.x);
}
if (lra.column_has_lower_bound(i)) {
impq l = lra.get_lower_bound(i); // copy: add_term invalidates bound references
vector<std::pair<mpq, unsigned>> coeffs = {{mpq(1), i}, {-delta.x, x_e}};
unsigned s = lra.add_term(coeffs, UINT_MAX);
lra.add_var_bound(s, is_zero(l.y) ? lconstraint_kind::GE : lconstraint_kind::GT, l.x);
}
return true;
}
lia_move int_cube::sat_after_rounding() {
lra.round_to_integer_solution();
lra.set_status(lp_status::FEASIBLE);
SASSERT(lia.settings().get_cancel_flag() || lia.is_feasible());
TRACE(cube, tout << "largest cube success";);
lia.settings().stats().m_lcube_success++;
return lia_move::sat;
}
lia_move int_cube::round_and_repair() {
lra.backup_x(); // remember the cube center
vector<flip_candidate> flips;
for (unsigned j = 0; j < lra.column_count(); ++j) {
if (!lra.column_is_int(j) || lra.column_has_term(j))
continue;
const impq& v = lra.get_column_value(j);
if (v.is_int())
continue;
flips.push_back({j, floor(v), false});
}
lra.round_to_integer_solution();
for (auto& f : flips)
f.m_at_hi = lra.get_column_value(f.m_j).x > f.m_lo;
if (repair_rounded_candidate(flips)) {
lra.set_status(lp_status::FEASIBLE);
SASSERT(lia.settings().get_cancel_flag() || lia.is_feasible());
TRACE(cube, tout << "largest cube success";);
lia.settings().stats().m_lcube_success++;
return lia_move::sat;
}
// return to the cube center: an interior point of the polyhedron
lra.restore_x();
lra.set_status(lp_status::FEASIBLE);
return lia_move::undef;
}
// Checks the rounded center against the original constraints. On failure
// searches the vertices of the lattice cell around the center greedily:
// flip a coordinate between floor and ceiling to maximally decrease the
// total bound violation, within a budget.
bool int_cube::repair_rounded_candidate(vector<flip_candidate>& flips) {
vector<bounded_row> rows;
for (const lar_term* t : lra.terms()) {
unsigned j = t->j();
if (!lra.column_associated_with_row(j))
continue;
if (!lra.column_has_upper_bound(j) && !lra.column_has_lower_bound(j))
continue;
bounded_row r;
r.m_j = j;
r.m_val = t->apply(lra.r_x());
rows.push_back(r);
}
auto row_violation = [&](unsigned ri, const impq& v) {
impq w;
unsigned j = rows[ri].m_j;
if (lra.column_has_upper_bound(j) && v > lra.get_upper_bound(j))
w += v - lra.get_upper_bound(j);
if (lra.column_has_lower_bound(j) && v < lra.get_lower_bound(j))
w += lra.get_lower_bound(j) - v;
return w;
};
impq violation;
for (unsigned ri = 0; ri < rows.size(); ++ri)
violation += row_violation(ri, rows[ri].m_val);
if (is_zero(violation))
return true; // the rounded center fits as it is
if (flips.empty())
return false;
std::unordered_map<unsigned, unsigned> flip_of_var;
for (unsigned fi = 0; fi < flips.size(); ++fi)
flip_of_var[flips[fi].m_j] = fi;
// occurrences of the flip candidates in the bounded rows
vector<vector<std::pair<unsigned, mpq>>> occs(flips.size());
for (unsigned ri = 0; ri < rows.size(); ++ri) {
const lar_term& t = lra.get_term(rows[ri].m_j);
for (lar_term::ival p : t) {
auto it = flip_of_var.find(p.j());
if (it != flip_of_var.end())
occs[it->second].push_back({ri, p.coeff()});
}
}
unsigned budget = std::min(2 * flips.size(), lia.settings().lcube_flips());
bool flipped = false;
while (!is_zero(violation) && budget-- > 0) {
unsigned best_fi = UINT_MAX;
impq best_gain;
for (unsigned fi = 0; fi < flips.size(); ++fi) {
if (occs[fi].empty())
continue;
mpq step = flips[fi].m_at_hi ? mpq(-1) : mpq(1);
impq gain;
for (const auto& o : occs[fi]) {
const impq& v = rows[o.first].m_val;
gain += row_violation(o.first, v + impq(step * o.second)) - row_violation(o.first, v);
}
if (gain < best_gain) {
best_gain = gain;
best_fi = fi;
}
}
if (seen_minus && seen_plus)
return zero_of_type<impq>();
return impq(0, 1);
if (best_fi == UINT_MAX)
return false; // no flip decreases the violation
mpq step = flips[best_fi].m_at_hi ? mpq(-1) : mpq(1);
for (const auto& o : occs[best_fi])
rows[o.first].m_val += impq(step * o.second);
flips[best_fi].m_at_hi = !flips[best_fi].m_at_hi;
violation += best_gain;
flipped = true;
TRACE(cube, tout << "flipped column " << flips[best_fi].m_j << ", violation = " << violation << "\n";);
}
if (!is_zero(violation))
return false;
// apply the repaired candidate
for (const auto& f : flips)
lra.set_column_value(f.m_j, impq(f.m_at_hi ? f.m_lo + 1 : f.m_lo));
for (const lar_term* t : lra.terms()) {
unsigned j = t->j();
if (!lra.column_associated_with_row(j))
continue;
lra.set_column_value(j, t->apply(lra.r_x()));
}
if (flipped)
lia.settings().stats().m_lcube_flip_success++;
return true;
}
impq int_cube::get_cube_delta_for_term(const lar_term& t) const {
if (t.size() == 2) {
bool seen_minus = false, seen_plus = false, all_ok = true;
for (lar_term::ival p : t) {
if (!lia.column_is_int(p.j())) { all_ok = false; break; }
const mpq& c = p.coeff();
if (c == one_of_type<mpq>()) seen_plus = true;
else if (c == -one_of_type<mpq>()) seen_minus = true;
else { all_ok = false; break; }
}
if (all_ok) {
if (seen_minus && seen_plus)
return zero_of_type<impq>();
return impq(0, 1);
}
}
usual_delta:
mpq delta = zero_of_type<mpq>();
for (lar_term::ival p : t)
if (lia.column_is_int(p.j()))
delta += abs(p.coeff());
delta *= mpq(1, 2);
return impq(delta);
}

29
src/math/lp/int_cube.h
View file

@ -10,9 +10,15 @@ Abstract:
Cube finder
This routine attempts to find a feasible integer solution
by tightnening bounds and running an LRA solver on the
by tightnening bounds and running an LRA solver on the
tighter system.
find_largest_cube() implements the largest cube test of
Bromberger and Weidenbach (Fast Cube Tests for LIA Constraint
Solving, IJCAR 2016): a fresh variable x_e for the cube edge
length is introduced and maximized; the center of the largest
cube is rounded to a candidate integer solution.
Author:
Nikolaj Bjorner (nbjorner)
Lev Nachmanson (levnach)
@ -21,7 +27,10 @@ Revision History:
--*/
#pragma once
#include "util/vector.h"
#include "math/lp/lia_move.h"
#include "math/lp/numeric_pair.h"
#include "math/lp/lar_term.h"
namespace lp {
class int_solver;
@ -29,12 +38,30 @@ namespace lp {
class int_cube {
class int_solver& lia;
class lar_solver& lra;
// a fractional integer coordinate of the cube center:
// the candidate value is m_lo or m_lo + 1
struct flip_candidate {
unsigned m_j = 0;
mpq m_lo;
bool m_at_hi = false;
};
// a term column with at least one bound, tracked during the repair
struct bounded_row {
unsigned m_j = 0;
impq m_val;
};
bool tighten_term_for_cube(unsigned i);
bool tighten_terms_for_cube();
void find_feasible_solution();
impq get_cube_delta_for_term(const lar_term& t) const;
bool add_edge_rows_for_term(unsigned i, unsigned x_e);
bool add_cube_edge_rows(unsigned x_e);
lia_move sat_after_rounding();
lia_move round_and_repair();
bool repair_rounded_candidate(vector<flip_candidate>& flips);
public:
int_cube(int_solver& lia);
lia_move operator()();
lia_move find_largest_cube();
};
}

72
src/math/lp/int_solver.cpp
View file

@ -44,7 +44,11 @@ namespace lp {
dioph_eq m_dio;
int_gcd_test m_gcd;
unsigned m_initial_dio_calls_period;
unsigned m_lcube_period;
// The number of consecutive genuine dio calls that returned undef, reset on a dio
// conflict. Drives the decision to start running Gomory with dio.
unsigned m_dio_undef_in_a_row = 0;
bool column_is_int_inf(unsigned j) const {
return lra.column_is_int(j) && (!lia.value_is_int(j));
}
@ -52,7 +56,8 @@ namespace lp {
imp(int_solver& lia): lia(lia), lra(lia.lra), lrac(lia.lrac), m_hnf_cutter(lia), m_dio(lia), m_gcd(lia) {
m_hnf_cut_period = settings().hnf_cut_period();
m_initial_dio_calls_period = settings().dio_calls_period();
}
m_lcube_period = settings().m_int_find_cube_period;
}
bool has_lower(unsigned j) const {
switch (lrac.m_column_types()[j]) {
@ -176,14 +181,16 @@ namespace lp {
if (r == lia_move::conflict) {
m_dio.explain(*this->m_ex);
lia.settings().dio_calls_period() = m_initial_dio_calls_period;
lia.settings().dio_enable_gomory_cuts() = false;
m_dio_undef_in_a_row = 0;
lia.settings().stop_running_gomory_with_dio(); // dio was productive: stop running Gomory
lia.settings().set_run_gcd_test(false);
return lia_move::conflict;
}
if (r == lia_move::undef) {
lia.settings().dio_calls_period() *= 2;
if (lra.settings().dio_calls_period() >= 16) {
lia.settings().dio_enable_gomory_cuts() = true;
lia.settings().dio_calls_period() *= 2; // throttle dio scheduling on failure
++m_dio_undef_in_a_row;
if (m_dio_undef_in_a_row >= lra.settings().dio_gomory_enable_period()) {
lia.settings().start_running_gomory_with_dio(); // dio persistently unproductive: start running Gomory
lia.settings().set_run_gcd_test(true);
}
}
@ -191,26 +198,66 @@ namespace lp {
}
lp_settings& settings() { return lra.settings(); }
// Decide whether a periodic heuristic fires on this call. When
// random_hammers is enabled the gate is drawn at random with the same
// 1/period expected rate instead of a deterministic "every k-th call"
// modulus: a deterministic period can phase-lock with the search on
// some families and drown the solver in conflicts while another handler
// is starved; randomizing the gate breaks that resonance.
bool hit_period(unsigned period) {
if (period <= 1)
return true;
if (settings().random_hammers())
return settings().random_next(period) == 0;
return m_number_of_calls % period == 0;
}
bool should_find_cube() {
return m_number_of_calls % settings().m_int_find_cube_period == 0;
return hit_period(settings().m_int_find_cube_period);
}
// The largest cube test is throttled exponentially: when the polyhedron
// does not contain a large enough cube it is unlikely to contain one
// later, after more constraints are added, so each failure doubles the
// period and a success resets it.
bool should_find_lcube() {
return settings().lcube() && hit_period(m_lcube_period);
}
lia_move find_lcube() {
lia_move r = int_cube(lia).find_largest_cube();
if (r == lia_move::undef) {
if (m_lcube_period < (1u << 30))
m_lcube_period *= 2;
}
else
m_lcube_period = settings().m_int_find_cube_period;
return r;
}
bool should_gomory_cut() {
bool dio_allows_gomory = !settings().dio() || settings().dio_enable_gomory_cuts() ||
m_dio.some_terms_are_ignored();
return dio_allows_gomory && m_number_of_calls % settings().m_int_gomory_cut_period == 0;
return dio_allows_gomory && hit_period(settings().m_int_gomory_cut_period);
}
bool should_solve_dioph_eq() {
return lia.settings().dio() && (m_number_of_calls % settings().dio_calls_period() == 0);
bool ret = lia.settings().dio() && hit_period(settings().dio_calls_period());
if (!ret && lia.settings().dio_calls_period() > m_initial_dio_calls_period) {
unsigned dec = settings().dio_calls_period_decrease();
unsigned& period = lia.settings().dio_calls_period();
period = period > m_initial_dio_calls_period + dec ? period - dec : m_initial_dio_calls_period;
}
return ret;
}
// HNF
bool should_hnf_cut() {
return (!settings().dio() || settings().dio_enable_hnf_cuts())
&& settings().enable_hnf() && m_number_of_calls % settings().hnf_cut_period() == 0;
&& settings().enable_hnf() && hit_period(settings().hnf_cut_period());
}
lia_move hnf_cut() {
@ -245,11 +292,12 @@ namespace lp {
++m_number_of_calls;
if (r == lia_move::undef) r = patch_basic_columns();
if (r == lia_move::undef && should_find_cube()) r = int_cube(lia)();
if (r == lia_move::undef && should_find_cube()) r = int_cube(lia)();
if (r == lia_move::undef && should_find_lcube()) r = find_lcube();
if (r == lia_move::undef) lra.move_non_basic_columns_to_bounds();
if (r == lia_move::undef && should_hnf_cut()) r = hnf_cut();
if (r == lia_move::undef && should_gomory_cut()) r = gomory(lia).get_gomory_cuts(2);
if (r == lia_move::undef && should_solve_dioph_eq()) r = solve_dioph_eq();
if (r == lia_move::undef && should_gomory_cut()) r = gomory(lia).get_gomory_cuts(2);
if (r == lia_move::undef) r = int_branch(lia)();
if (settings().get_cancel_flag()) r = lia_move::undef;
return r;

34
src/math/lp/lar_constraints.h
View file

@ -39,7 +39,17 @@ inline std::string lconstraint_kind_string(lconstraint_kind t) {
class lar_base_constraint {
lconstraint_kind m_kind;
mpq m_right_side;
// Right-hand side as a delta-rational pair (x, y) = x + y*delta. The
// rational part x is the ordinary bound value returned by rhs(); the
// infinitesimal part y (bound_eps) is non-zero only for the delta-rational
// bounds that validate strict optimization suprema/infima (see
// opt_solver::maximize_objective). The strict kinds already carry a
// matching-sign infinitesimal in update_bound_with_* (LT -> upper bound
// (r, -1); GT -> lower bound (r, +1)); y is needed for the OPPOSITE-sign
// case that no kind yields: a lower bound (r, -1), i.e. objvar >= r - delta
// (GE gives (r, 0), GT gives (r, +1)), which is how a maximize supremum
// r - delta is asserted. y is zero for all ordinary constraints.
impq m_right_side;
bool m_active;
bool m_is_auxiliary;
unsigned m_j;
@ -53,10 +63,17 @@ public:
virtual ~lar_base_constraint() = default;
lconstraint_kind kind() const { return m_kind; }
mpq const& rhs() const { return m_right_side; }
// First (rational) component of the right-hand side pair.
mpq const& rhs() const { return m_right_side.x; }
// Whole right-hand side pair (rational value + infinitesimal, see below).
impq const& rhs_impq() const { return m_right_side; }
unsigned column() const { return m_j; }
u_dependency* dep() const { return m_dep; }
// Second (infinitesimal) component of the right-hand side pair.
mpq const& bound_eps() const { return m_right_side.y; }
void set_bound_eps(mpq const& e) { m_right_side.y = e; }
void activate() { m_active = true; }
void deactivate() { m_active = false; }
bool is_active() const { return m_active; }
@ -181,11 +198,24 @@ public:
return add(new (m_region) lar_var_constraint(j, k, mk_dep(), rhs));
}
constraint_index add_var_constraint(lpvar j, lconstraint_kind k, mpq const& rhs, mpq const& eps) {
auto* c = new (m_region) lar_var_constraint(j, k, mk_dep(), rhs);
c->set_bound_eps(eps);
return add(c);
}
constraint_index add_term_constraint(unsigned j, const lar_term* t, lconstraint_kind k, mpq const& rhs) {
auto* dep = mk_dep();
return add(new (m_region) lar_term_constraint(j, t, k, dep, rhs));
}
constraint_index add_term_constraint(unsigned j, const lar_term* t, lconstraint_kind k, mpq const& rhs, mpq const& eps) {
auto* dep = mk_dep();
auto* c = new (m_region) lar_term_constraint(j, t, k, dep, rhs);
c->set_bound_eps(eps);
return add(c);
}
// future behavior uses activation bit.
bool is_active(constraint_index ci) const { return m_constraints[ci]->is_active(); }

9
src/math/lp/lar_core_solver.h
View file

@ -81,10 +81,15 @@ public:
void backup_x() { m_backup_x = m_r_x; }
void restore_x() {
m_r_x = m_backup_x;
m_r_x.reserve(m_m());
unsigned n = m_r_A.column_count();
unsigned backup_sz = m_backup_x.size();
unsigned copy_sz = std::min(backup_sz, n);
for (unsigned j = 0; j < copy_sz; j++)
m_r_x[j] = m_backup_x[j];
}
unsigned backup_x_size() const { return m_backup_x.size(); }
vector<impq> const& r_x() const { return m_r_x; }
impq& r_x(unsigned j) { return m_r_x[j]; }
impq const& r_x(unsigned j) const { return m_r_x[j]; }

185
src/math/lp/lar_solver.cpp
View file

@ -36,7 +36,12 @@ namespace lp {
struct term_comparer {
bool operator()(const lar_term& a, const lar_term& b) const {
return a == b;
if (a.size() != b.size()) return false;
for (const auto& p : a) {
auto const* e = b.coeffs().find_core(p.j());
if (!e || e->get_data().m_value != p.coeff()) return false;
}
return true;
}
};
@ -58,6 +63,7 @@ namespace lp {
unsigned_vector m_row_bounds_to_replay;
u_dependency_manager m_dependencies;
svector<constraint_index> m_tmp_dependencies;
ptr_vector<u_dependency> m_tmp_witnesses;
u_dependency* m_crossed_bounds_deps = nullptr;
lpvar m_crossed_bounds_column = null_lpvar;
@ -467,6 +473,8 @@ namespace lp {
return ret;
}
lp_status lar_solver::solve() {
if (m_imp->m_status == lp_status::INFEASIBLE || m_imp->m_status == lp_status::CANCELLED)
return m_imp->m_status;
@ -857,7 +865,7 @@ namespace lp {
}
lp_status lar_solver::maximize_term(unsigned j,
impq& term_max) {
impq& term_max, bool fix_int_cols) {
TRACE(lar_solver, print_values(tout););
SASSERT(get_core_solver().m_r_solver.calc_current_x_is_feasible_include_non_basis());
lar_term term = get_term_to_maximize(j);
@ -872,6 +880,11 @@ namespace lp {
return lp_status::UNBOUNDED;
}
if (!fix_int_cols) {
set_status(lp_status::OPTIMAL);
return lp_status::OPTIMAL;
}
impq opt_val = term_max;
bool change = false;
@ -1113,8 +1126,28 @@ namespace lp {
void lar_solver::explain_fixed_column(unsigned j, explanation& ex) {
SASSERT(column_is_fixed(j));
auto* deps = get_bound_constraint_witnesses_for_column(j);
for (auto ci : flatten(deps))
const column& ul = m_imp->m_columns[j];
m_imp->m_tmp_dependencies.reset();
m_imp->m_dependencies.linearize(ul.lower_bound_witness(), ul.upper_bound_witness(), m_imp->m_tmp_dependencies);
for (auto ci : m_imp->m_tmp_dependencies)
ex.push_back(ci);
}
// Linearize the bound witnesses of all fixed columns in the row together, so the
// mark bits walk each dependency sub-DAG shared between columns only once.
void lar_solver::explain_fixed_in_row(unsigned row, explanation& ex) {
auto& witnesses = m_imp->m_tmp_witnesses;
witnesses.reset();
for (auto const& c : get_row(row)) {
if (!column_is_fixed(c.var()))
continue;
const column& ul = m_imp->m_columns[c.var()];
witnesses.push_back(ul.lower_bound_witness());
witnesses.push_back(ul.upper_bound_witness());
}
m_imp->m_tmp_dependencies.reset();
m_imp->m_dependencies.linearize(witnesses, m_imp->m_tmp_dependencies);
for (auto ci : m_imp->m_tmp_dependencies)
ex.push_back(ci);
}
@ -1305,10 +1338,14 @@ namespace lp {
if (m_imp->m_settings.get_cancel_flag())
return true;
std::unordered_map<lpvar, mpq> var_map;
get_model_do_not_care_about_diff_vars(var_map);
// Compute the strict-bounds delta once per model: it flattens both the
// model (var_map) and the eps component of any delta-rational bound in
// constraint_holds, so the two must use the very same value.
mpq delta = get_core_solver().find_delta_for_strict_bounds(m_imp->m_settings.m_epsilon);
get_model_do_not_care_about_diff_vars(var_map, delta);
for (auto const& c : m_imp->m_constraints.active()) {
if (!constraint_holds(c, var_map)) {
if (!constraint_holds(c, var_map, delta)) {
TRACE(lar_solver,
m_imp->m_constraints.display(tout, c) << "\n";
for (auto p : c.coeffs()) {
@ -1320,14 +1357,21 @@ namespace lp {
return true;
}
bool lar_solver::constraint_holds(const lar_base_constraint& constr, std::unordered_map<lpvar, mpq>& var_map) const {
bool lar_solver::constraint_holds(const lar_base_constraint& constr, std::unordered_map<lpvar, mpq>& var_map, const mpq& delta) const {
mpq left_side_val = get_left_side_val(constr, var_map);
// Account for a delta-rational bound rhs + eps*delta (eps != 0 only for
// the bounds that validate strict optimization optima). 'delta' is the
// same strict-bounds delta that flattened var_map, so the comparison is
// exact over the reals.
mpq rhs = constr.rhs();
if (!constr.bound_eps().is_zero())
rhs += constr.bound_eps() * delta;
switch (constr.kind()) {
case LE: return left_side_val <= constr.rhs();
case LT: return left_side_val < constr.rhs();
case GE: return left_side_val >= constr.rhs();
case GT: return left_side_val > constr.rhs();
case EQ: return left_side_val == constr.rhs();
case LE: return left_side_val <= rhs;
case LT: return left_side_val < rhs;
case GE: return left_side_val >= rhs;
case GT: return left_side_val > rhs;
case EQ: return left_side_val == rhs;
default:
UNREACHABLE();
}
@ -1550,6 +1594,10 @@ namespace lp {
void lar_solver::get_model_do_not_care_about_diff_vars(std::unordered_map<lpvar, mpq>& variable_values) const {
mpq delta = get_core_solver().find_delta_for_strict_bounds(m_imp->m_settings.m_epsilon);
get_model_do_not_care_about_diff_vars(variable_values, delta);
}
void lar_solver::get_model_do_not_care_about_diff_vars(std::unordered_map<lpvar, mpq>& variable_values, const mpq& delta) const {
for (unsigned i = 0; i < get_core_solver().r_x().size(); ++i) {
const impq& rp = get_core_solver().r_x(i);
variable_values[i] = rp.x + delta * rp.y;
@ -2109,12 +2157,12 @@ namespace lp {
void lar_solver::activate_check_on_equal(constraint_index ci, unsigned& equal_column) {
auto const& c = m_imp->m_constraints[ci];
update_column_type_and_bound_check_on_equal(c.column(), c.rhs(), ci, equal_column);
update_column_type_and_bound_check_on_equal(c.column(), c.rhs_impq(), ci, equal_column);
}
void lar_solver::activate(constraint_index ci) {
auto const& c = m_imp->m_constraints[ci];
update_column_type_and_bound(c.column(), c.rhs(), ci);
update_column_type_and_bound(c.column(), c.rhs_impq(), ci);
}
mpq lar_solver::adjust_bound_for_int(lpvar j, lconstraint_kind& k, const mpq& bound) {
@ -2158,6 +2206,24 @@ namespace lp {
return ci;
}
// Variant that attaches an infinitesimal coefficient 'eps' to the bound, so
// that activating the resulting constraint asserts the delta-rational bound
// (right_side, eps). Used to faithfully validate strict optimization optima
// (e.g. a maximize supremum r - delta is validated as a lower bound
// (r, -1)). Only supported for plain column bounds (no term column).
constraint_index lar_solver::mk_var_bound(lpvar j, lconstraint_kind kind, const mpq& right_side, const mpq& eps) {
TRACE(lar_solver, tout << "j = " << get_variable_name(j) << " " << lconstraint_kind_string(kind) << " " << right_side << " + " << eps << "*eps" << std::endl;);
mpq rs = adjust_bound_for_int(j, kind, right_side);
SASSERT(bound_is_integer_for_integer_column(j, rs));
constraint_index ci;
if (!column_has_term(j))
ci = m_imp->m_constraints.add_var_constraint(j, kind, rs, eps);
else
ci = m_imp->m_constraints.add_term_constraint(j, m_imp->m_columns[j].term(), kind, rs, eps);
SASSERT(sizes_are_correct());
return ci;
}
bool lar_solver::compare_values(lpvar j, lconstraint_kind k, const mpq& rhs) {
return compare_values(get_column_value(j), k, rhs);
}
@ -2176,7 +2242,7 @@ namespace lp {
}
void lar_solver::update_column_type_and_bound(unsigned j,
const mpq& right_side,
const impq& right_side,
constraint_index constr_index) {
TRACE(lar_solver_feas, tout << "j = " << j << " was " << (this->column_is_feasible(j)?"feas":"non-feas") << std::endl;);
m_imp->m_constraints.activate(constr_index);
@ -2240,7 +2306,10 @@ namespace lp {
ls.add_var_bound(tv, c.kind(), c.rhs());
}
void lar_solver::update_column_type_and_bound(unsigned j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
// SASSERT(validate_bound(j, kind, right_side, dep));
update_column_type_and_bound(j, kind, impq(right_side), dep);
}
void lar_solver::update_column_type_and_bound(unsigned j, lconstraint_kind kind, const impq& right_side, u_dependency* dep) {
// SASSERT(validate_bound(j, kind, right_side.x, dep));
TRACE(
lar_solver_feas,
tout << "j" << j << " " << lconstraint_kind_string(kind) << " " << right_side << std::endl;
@ -2254,11 +2323,16 @@ namespace lp {
}
});
bool was_fixed = column_is_fixed(j);
mpq rs = adjust_bound_for_int(j, kind, right_side);
// adjust_bound_for_int operates on the rational part (and may sharpen
// the kind for integer columns); the infinitesimal part y is carried
// through unchanged. y is non-zero only for the delta-rational bounds
// that validate strict optimization optima, which target real columns.
mpq rs = adjust_bound_for_int(j, kind, right_side.x);
impq bound(rs, right_side.y);
if (column_has_upper_bound(j))
update_column_type_and_bound_with_ub(j, kind, rs, dep);
update_column_type_and_bound_with_ub(j, kind, bound, dep);
else
update_column_type_and_bound_with_no_ub(j, kind, rs, dep);
update_column_type_and_bound_with_no_ub(j, kind, bound, dep);
if (!was_fixed && column_is_fixed(j) && m_fixed_var_eh)
m_fixed_var_eh(j);
@ -2287,7 +2361,7 @@ namespace lp {
}
void lar_solver::update_column_type_and_bound_check_on_equal(unsigned j,
const mpq& right_side,
const impq& right_side,
constraint_index constr_index,
unsigned& equal_to_j) {
update_column_type_and_bound(j, right_side, constr_index);
@ -2303,17 +2377,7 @@ namespace lp {
return m_imp->m_constraints.add_term_constraint(j, m_imp->m_columns[j].term(), kind, rs);
}
struct lar_solver::scoped_backup {
lar_solver& m_s;
scoped_backup(lar_solver& s) : m_s(s) {
m_s.get_core_solver().backup_x();
}
~scoped_backup() {
m_s.get_core_solver().restore_x();
}
};
void lar_solver::update_column_type_and_bound_with_ub(unsigned j, lp::lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
void lar_solver::update_column_type_and_bound_with_ub(unsigned j, lp::lconstraint_kind kind, const impq& right_side, u_dependency* dep) {
SASSERT(column_has_upper_bound(j));
if (column_has_lower_bound(j)) {
update_bound_with_ub_lb(j, kind, right_side, dep);
@ -2323,7 +2387,7 @@ namespace lp {
}
}
void lar_solver::update_column_type_and_bound_with_no_ub(unsigned j, lp::lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
void lar_solver::update_column_type_and_bound_with_no_ub(unsigned j, lp::lconstraint_kind kind, const impq& right_side, u_dependency* dep) {
SASSERT(!column_has_upper_bound(j));
if (column_has_lower_bound(j)) {
update_bound_with_no_ub_lb(j, kind, right_side, dep);
@ -2333,18 +2397,18 @@ namespace lp {
}
}
void lar_solver::update_bound_with_ub_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
void lar_solver::update_bound_with_ub_lb(lpvar j, lconstraint_kind kind, const impq& right_side, u_dependency* dep) {
SASSERT(column_has_lower_bound(j) && column_has_upper_bound(j));
SASSERT(get_core_solver().m_column_types[j] == column_type::boxed ||
get_core_solver().m_column_types[j] == column_type::fixed);
mpq y_of_bound(0);
mpq y_of_bound(right_side.y);
switch (kind) {
case LT:
y_of_bound = -1;
y_of_bound += -1;
Z3_fallthrough;
case LE: {
auto up = numeric_pair<mpq>(right_side, y_of_bound);
auto up = numeric_pair<mpq>(right_side.x, y_of_bound);
if (up < get_lower_bound(j)) {
set_crossed_bounds_column_and_deps(j, true, dep);
}
@ -2356,10 +2420,10 @@ namespace lp {
break;
}
case GT:
y_of_bound = 1;
y_of_bound += 1;
Z3_fallthrough;
case GE: {
auto low = numeric_pair<mpq>(right_side, y_of_bound);
auto low = numeric_pair<mpq>(right_side.x, y_of_bound);
if (low > get_upper_bound(j)) {
set_crossed_bounds_column_and_deps(j, false, dep);
}
@ -2372,7 +2436,7 @@ namespace lp {
break;
}
case EQ: {
auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
auto v = numeric_pair<mpq>(right_side.x, zero_of_type<mpq>());
if (v > get_upper_bound(j))
set_crossed_bounds_column_and_deps(j, false, dep);
else if (v < get_lower_bound(j))
@ -2393,17 +2457,17 @@ namespace lp {
get_core_solver().m_column_types[j] = column_type::fixed;
}
void lar_solver::update_bound_with_no_ub_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
void lar_solver::update_bound_with_no_ub_lb(lpvar j, lconstraint_kind kind, const impq& right_side, u_dependency* dep) {
SASSERT(column_has_lower_bound(j) && !column_has_upper_bound(j));
SASSERT(get_core_solver().m_column_types[j] == column_type::lower_bound);
mpq y_of_bound(0);
mpq y_of_bound(right_side.y);
switch (kind) {
case LT:
y_of_bound = -1;
y_of_bound += -1;
Z3_fallthrough;
case LE: {
auto up = numeric_pair<mpq>(right_side, y_of_bound);
auto up = numeric_pair<mpq>(right_side.x, y_of_bound);
if (up < get_lower_bound(j)) {
set_crossed_bounds_column_and_deps(j, true, dep);
}
@ -2414,9 +2478,9 @@ namespace lp {
break;
}
case GT:
y_of_bound = 1;
y_of_bound += 1;
case GE: {
auto low = numeric_pair<mpq>(right_side, y_of_bound);
auto low = numeric_pair<mpq>(right_side.x, y_of_bound);
if (low < get_lower_bound(j)) {
return;
}
@ -2424,7 +2488,7 @@ namespace lp {
break;
}
case EQ: {
auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
auto v = numeric_pair<mpq>(right_side.x, zero_of_type<mpq>());
if (v < get_lower_bound(j)) {
set_crossed_bounds_column_and_deps(j, true, dep);
}
@ -2441,28 +2505,28 @@ namespace lp {
}
}
void lar_solver::update_bound_with_ub_no_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
void lar_solver::update_bound_with_ub_no_lb(lpvar j, lconstraint_kind kind, const impq& right_side, u_dependency* dep) {
SASSERT(!column_has_lower_bound(j) && column_has_upper_bound(j));
SASSERT(get_core_solver().m_column_types[j] == column_type::upper_bound);
mpq y_of_bound(0);
mpq y_of_bound(right_side.y);
switch (kind) {
case LT:
y_of_bound = -1;
y_of_bound += -1;
Z3_fallthrough;
case LE:
{
auto up = numeric_pair<mpq>(right_side, y_of_bound);
auto up = numeric_pair<mpq>(right_side.x, y_of_bound);
if (up >= get_upper_bound(j))
return;
set_upper_bound_witness(j, dep, up);
}
break;
case GT:
y_of_bound = 1;
y_of_bound += 1;
Z3_fallthrough;
case GE:
{
auto low = numeric_pair<mpq>(right_side, y_of_bound);
auto low = numeric_pair<mpq>(right_side.x, y_of_bound);
if (low > get_upper_bound(j)) {
set_crossed_bounds_column_and_deps(j, false, dep);
}
@ -2474,7 +2538,7 @@ namespace lp {
break;
case EQ:
{
auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
auto v = numeric_pair<mpq>(right_side.x, zero_of_type<mpq>());
if (v > get_upper_bound(j)) {
set_crossed_bounds_column_and_deps(j, false, dep);
}
@ -2491,30 +2555,30 @@ namespace lp {
}
}
void lar_solver::update_bound_with_no_ub_no_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep) {
void lar_solver::update_bound_with_no_ub_no_lb(lpvar j, lconstraint_kind kind, const impq& right_side, u_dependency* dep) {
SASSERT(!column_has_lower_bound(j) && !column_has_upper_bound(j));
mpq y_of_bound(0);
mpq y_of_bound(right_side.y);
switch (kind) {
case LT:
y_of_bound = -1;
y_of_bound += -1;
Z3_fallthrough;
case LE: {
auto up = numeric_pair<mpq>(right_side, y_of_bound);
auto up = numeric_pair<mpq>(right_side.x, y_of_bound);
set_upper_bound_witness(j, dep, up);
get_core_solver().m_column_types[j] = column_type::upper_bound;
} break;
case GT:
y_of_bound = 1;
y_of_bound += 1;
Z3_fallthrough;
case GE: {
auto low = numeric_pair<mpq>(right_side, y_of_bound);
auto low = numeric_pair<mpq>(right_side.x, y_of_bound);
set_lower_bound_witness(j, dep, low);
get_core_solver().m_column_types[j] = column_type::lower_bound;
} break;
case EQ: {
auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
auto v = numeric_pair<mpq>(right_side.x, zero_of_type<mpq>());
set_upper_bound_witness(j, dep, v);
set_lower_bound_witness(j, dep, v);
get_core_solver().m_column_types[j] = column_type::fixed;
@ -3002,4 +3066,3 @@ namespace lp {
}
} // namespace lp

50
src/math/lp/lar_solver.h
View file

@ -72,7 +72,6 @@ class lar_solver : public column_namer {
void clear_columns_with_changed_bounds();
struct scoped_backup;
public:
const indexed_uint_set& columns_with_changed_bounds() const;
void insert_to_columns_with_changed_bounds(unsigned j);
@ -89,19 +88,20 @@ class lar_solver : public column_namer {
void add_bound_negation_to_solver(lar_solver& ls, lpvar j, lconstraint_kind kind, const mpq& right_side);
void add_constraint_to_validate(lar_solver& ls, constraint_index ci);
bool m_validate_blocker = false;
void update_column_type_and_bound_check_on_equal(unsigned j, const mpq& right_side, constraint_index ci, unsigned&);
void update_column_type_and_bound(unsigned j, const mpq& right_side, constraint_index ci);
void update_column_type_and_bound_check_on_equal(unsigned j, const impq& right_side, constraint_index ci, unsigned&);
void update_column_type_and_bound(unsigned j, const impq& right_side, constraint_index ci);
public:
bool validate_blocker() const { return m_validate_blocker; }
bool & validate_blocker() { return m_validate_blocker; }
void update_column_type_and_bound(unsigned j, lconstraint_kind kind, const impq& right_side, u_dependency* dep);
void update_column_type_and_bound(unsigned j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
private:
void update_column_type_and_bound_with_ub(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
void update_column_type_and_bound_with_no_ub(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
void update_bound_with_ub_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
void update_bound_with_no_ub_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
void update_bound_with_ub_no_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
void update_bound_with_no_ub_no_lb(lpvar j, lconstraint_kind kind, const mpq& right_side, u_dependency* dep);
void update_column_type_and_bound_with_ub(lpvar j, lconstraint_kind kind, const impq& right_side, u_dependency* dep);
void update_column_type_and_bound_with_no_ub(lpvar j, lconstraint_kind kind, const impq& right_side, u_dependency* dep);
void update_bound_with_ub_lb(lpvar j, lconstraint_kind kind, const impq& right_side, u_dependency* dep);
void update_bound_with_no_ub_lb(lpvar j, lconstraint_kind kind, const impq& right_side, u_dependency* dep);
void update_bound_with_ub_no_lb(lpvar j, lconstraint_kind kind, const impq& right_side, u_dependency* dep);
void update_bound_with_no_ub_no_lb(lpvar j, lconstraint_kind kind, const impq& right_side, u_dependency* dep);
void remove_non_fixed_from_fixed_var_table();
constraint_index add_var_bound_on_constraint_for_term(lpvar j, lconstraint_kind kind, const mpq& right_side);
void set_crossed_bounds_column_and_deps(unsigned j, bool lower_bound, u_dependency* dep);
@ -148,7 +148,7 @@ class lar_solver : public column_namer {
numeric_pair<mpq> get_basic_var_value_from_row(unsigned i);
bool all_constrained_variables_are_registered(const vector<std::pair<mpq, lpvar>>& left_side);
bool all_constraints_hold() const;
bool constraint_holds(const lar_base_constraint& constr, std::unordered_map<lpvar, mpq>& var_map) const;
bool constraint_holds(const lar_base_constraint& constr, std::unordered_map<lpvar, mpq>& var_map, const mpq& delta) const;
static void register_in_map(std::unordered_map<lpvar, mpq>& coeffs, const lar_base_constraint& cn, const mpq& a);
static void register_monoid_in_map(std::unordered_map<lpvar, mpq>& coeffs, const mpq& a, unsigned j);
bool the_left_sides_sum_to_zero(const vector<std::pair<mpq, unsigned>>& evidence) const;
@ -206,7 +206,9 @@ public:
set_column_value(j, v);
}
lp_status maximize_term(unsigned j_or_term, impq& term_max);
// fix_int_cols: after maximizing try to move the integer columns to integer values;
// pass false to keep the optimal (possibly fractional) vertex intact, e.g., for the largest cube test
lp_status maximize_term(unsigned j_or_term, impq& term_max, bool fix_int_cols);
core_solver_pretty_printer<lp::mpq, lp::impq> pp(std::ostream& out) const;
@ -270,6 +272,7 @@ public:
bool fixed_base_removed_correctly() const;
#endif
constraint_index mk_var_bound(lpvar j, lconstraint_kind kind, const mpq& right_side);
constraint_index mk_var_bound(lpvar j, lconstraint_kind kind, const mpq& right_side, const mpq& eps);
void activate_check_on_equal(constraint_index, lpvar&);
void activate(constraint_index);
void random_update(unsigned sz, lpvar const* vars);
@ -437,7 +440,28 @@ public:
statistics& stats();
void backup_x() { get_core_solver().backup_x(); }
void restore_x() { get_core_solver().restore_x(); }
void restore_x() {
auto& cs = get_core_solver();
unsigned backup_sz = cs.backup_x_size();
unsigned current_sz = cs.m_n();
CTRACE(lar_solver_restore, backup_sz != current_sz,
tout << "restore_x: backup_sz=" << backup_sz
<< " current_sz=" << current_sz << "\n";);
cs.restore_x();
if (backup_sz < current_sz) {
// New columns were added after backup.
// Recalculate basic variable values from non-basic ones
// to restore the Ax=0 tableau invariant, then snap
// non-basic columns to their bounds and find a feasible solution.
for (unsigned i = 0; i < A_r().row_count(); i++)
set_column_value(r_basis()[i], get_basic_var_value_from_row(i));
move_non_basic_columns_to_bounds();
}
else {
SASSERT(ax_is_correct());
SASSERT(cs.m_r_solver.calc_current_x_is_feasible_include_non_basis());
}
}
void updt_params(params_ref const& p);
column_type get_column_type(unsigned j) const { return get_core_solver().m_column_types()[j]; }
@ -456,6 +480,7 @@ public:
void get_model(std::unordered_map<lpvar, mpq>& variable_values) const;
void get_rid_of_inf_eps();
void get_model_do_not_care_about_diff_vars(std::unordered_map<lpvar, mpq>& variable_values) const;
void get_model_do_not_care_about_diff_vars(std::unordered_map<lpvar, mpq>& variable_values, const mpq& delta) const;
std::string get_variable_name(lpvar vi) const override;
void set_variable_name(lpvar vi, const std::string&);
unsigned number_of_vars() const;
@ -499,6 +524,7 @@ public:
}
void explain_fixed_column(unsigned j, explanation& ex);
void explain_fixed_in_row(unsigned row, explanation& ex);
u_dependency* join_deps(u_dependency* a, u_dependency *b) { return dep_manager().mk_join(a, b); }
const constraint_set & constraints() const;
void push();

4
src/math/lp/lar_term.h
View file

@ -129,8 +129,8 @@ public:
add_monomial(a, v1);
add_monomial(b, v2);
}
bool operator==(const lar_term & a) const { return false; } // take care not to create identical terms
bool operator!=(const lar_term & a) const { return ! (*this == a);}
bool operator==(const lar_term & a) const = delete; // take care not to create identical terms
bool operator!=(const lar_term & a) const = delete;
// some terms get used in add constraint
// it is the same as the offset in the m_constraints

14
src/math/lp/lp_api.h
View file

@ -36,15 +36,23 @@ namespace lp_api {
rational m_value;
bound_kind m_bound_kind;
lp::constraint_index m_constraints[2];
// Infinitesimal coefficient of the asserted (positive-literal) bound
// value: the bound means v (>=|<=) m_value + m_eps*delta. Non-zero
// only for the delta-rational bounds used to validate strict
// optimization optima (e.g. a lower bound (r, -1) for a maximize
// supremum). Kept so get_value reports the true delta value and no
// spurious rational fixed-variable equality is propagated to EUF.
rational m_eps;
public:
bound(Literal bv, theory_var v, lp::lpvar vi, bool is_int, rational const& val, bound_kind k, lp::constraint_index ct, lp::constraint_index cf) :
bound(Literal bv, theory_var v, lp::lpvar vi, bool is_int, rational const& val, bound_kind k, lp::constraint_index ct, lp::constraint_index cf, rational const& eps = rational::zero()) :
m_bv(bv),
m_var(v),
m_column_index(vi),
m_is_int(is_int),
m_value(val),
m_bound_kind(k) {
m_bound_kind(k),
m_eps(eps) {
m_constraints[0] = cf;
m_constraints[1] = ct;
}
@ -67,7 +75,7 @@ namespace lp_api {
inf_rational get_value(bool is_true) const {
if (is_true != get_lit().sign())
return inf_rational(m_value); // v >= value or v <= value
return inf_rational(m_value, m_eps); // v >= value (+ eps*delta) or v <= value (+ eps*delta)
if (m_is_int) {
SASSERT(m_value.is_int());
rational const& offset = (m_bound_kind == lower_t) ? rational::minus_one() : rational::one();

16
src/math/lp/lp_bound_propagator.h
View file

@ -248,22 +248,16 @@ public:
void explain_fixed_in_row(unsigned row, explanation& ex) {
TRACE(eq, tout << lp().get_row(row) << std::endl);
for (const auto& c : lp().get_row(row))
if (lp().column_is_fixed(c.var()))
lp().explain_fixed_column(c.var(), ex);
lp().explain_fixed_in_row(row, ex);
}
unsigned explain_fixed_in_row_and_get_base(unsigned row, explanation& ex) {
unsigned base = UINT_MAX;
TRACE(eq, tout << lp().get_row(row) << std::endl);
for (const auto& c : lp().get_row(row)) {
if (lp().column_is_fixed(c.var())) {
lp().explain_fixed_column(c.var(), ex);
}
else if (lp().is_base(c.var())) {
lp().explain_fixed_in_row(row, ex);
unsigned base = UINT_MAX;
for (const auto& c : lp().get_row(row))
if (!lp().column_is_fixed(c.var()) && lp().is_base(c.var()))
base = c.var();
}
}
return base;
}

8
src/math/lp/lp_core_solver_base_def.h
View file

@ -277,13 +277,11 @@ pivot_column_tableau(unsigned j, unsigned piv_row_index) {
m_A.m_rows[c.var()][c.offset()].offset() = pivot_col_cell_index;
}
while (column.size() > 1) {
auto & c = column.back();
auto& c = column.back();
SASSERT(c.var() != piv_row_index);
if(! m_A.pivot_row_to_row_given_cell(piv_row_index, c, j)) {
return false;
}
if (m_touched_rows!= nullptr)
if (m_touched_rows != nullptr)
m_touched_rows->insert(c.var());
m_A.pivot_row_to_row_given_cell(piv_row_index, c, j);
}
if (m_settings.simplex_strategy() == simplex_strategy_enum::tableau_costs)

8
src/math/lp/lp_params_helper.pyg
View file

@ -5,9 +5,15 @@ def_module_params(module_name='lp',
params=(('dio', BOOL, True, 'use Diophantine equalities'),
('dio_branching_period', UINT, 100, 'Period of calling branching on undef in Diophantine handler'),
('dio_cuts_enable_gomory', BOOL, False, 'enable Gomory cuts together with Diophantine cuts, only relevant when dioph_eq is true'),
('dio_gomory_enable_period', UINT, 16, 'number of consecutive unproductive (undef) Diophantine-handler calls after which the controller starts running Gomory cuts and the gcd test alongside dio; a dio conflict resets the count and stops them; set very large to never start them this way so Gomory follows dio_cuts_enable_gomory only'),
('dio_cuts_enable_hnf', BOOL, True, 'enable hnf cuts together with Diophantine cuts, only relevant when dioph_eq is true'),
('dio_ignore_big_nums', BOOL, True, 'Ignore the terms with big numbers in the Diophantine handler, only relevant when dioph_eq is true'),
('dio_calls_period', UINT, 1, 'Period of calling the Diophantine handler in the final_check()'),
('dio_run_gcd', BOOL, False, 'Run the GCD heuristic if dio is on, if dio is disabled the option is not used'),
('dio_calls_period_decrease', UINT, 2, 'Amount by which dio_calls_period is decreased on each final_check() call where the Diophantine handler is not triggered, until it returns to its initial value'),
('dio_run_gcd', BOOL, False, 'Run the GCD heuristic if dio is on, if dio is disabled the option is not used'),
('lcube', BOOL, True, 'use the largest cube test for integer feasibility'),
('lcube_flips', UINT, 16, 'maximal number of coordinate flips when repairing the rounded largest cube center, only relevant when lcube is true'),
('int_hammer_period', UINT, 4, 'period (in final_check calls) for the integer cut/cube heuristics (find_cube, hnf, gomory); a smaller value calls them more often'),
('random_hammers', BOOL, True, 'draw the periodic integer heuristic gates (find_cube, lcube, hnf, gomory, dio) at random with the same 1/period rate instead of a deterministic every-k-th-call modulus'),
))

2
src/math/lp/lp_primal_core_solver.h
View file

@ -267,7 +267,7 @@ namespace lp {
unsigned j, const T &m, X &theta, bool &unlimited) {
SASSERT(m > 0 && this->m_column_types[j] == column_type::upper_bound);
limit_inf_on_bound_m_pos(m, this->m_x[j], this->m_upper_bounds[j], theta, unlimited);
};
}
void get_bound_on_variable_and_update_leaving_precisely(
unsigned j, vector<unsigned> &leavings, T m, X &t,

12
src/math/lp/lp_settings.cpp
View file

@ -37,11 +37,21 @@ void lp::lp_settings::updt_params(params_ref const& _p) {
auto eps = p.arith_epsilon();
m_epsilon = rational(std::max(1, (int)(100000*eps)), 100000);
m_dio = lp_p.dio();
m_dio_enable_gomory_cuts = lp_p.dio_cuts_enable_gomory();
m_dio_cuts_enable_gomory = lp_p.dio_cuts_enable_gomory();
m_dio_gomory_enable_period = lp_p.dio_gomory_enable_period();
m_dio_enable_hnf_cuts = lp_p.dio_cuts_enable_hnf();
m_dump_bound_lemmas = p.arith_dump_bound_lemmas();
m_dio_ignore_big_nums = lp_p.dio_ignore_big_nums();
m_dio_calls_period = lp_p.dio_calls_period();
m_dio_calls_period_decrease = lp_p.dio_calls_period_decrease();
m_dio_run_gcd = lp_p.dio_run_gcd();
m_random_hammers = lp_p.random_hammers();
m_lcube = lp_p.lcube();
m_lcube_flips = lp_p.lcube_flips();
unsigned hammer_period = lp_p.int_hammer_period();
SASSERT(hammer_period != 0);
m_int_find_cube_period = hammer_period;
m_int_gomory_cut_period = hammer_period;
m_hnf_cut_period = hammer_period;
m_max_conflicts = p.max_conflicts();
}

35
src/math/lp/lp_settings.h
View file

@ -112,6 +112,9 @@ struct statistics {
unsigned m_gcd_conflicts = 0;
unsigned m_cube_calls = 0;
unsigned m_cube_success = 0;
unsigned m_lcube_calls = 0;
unsigned m_lcube_success = 0;
unsigned m_lcube_flip_success = 0;
unsigned m_patches = 0;
unsigned m_patches_success = 0;
unsigned m_hnf_cutter_calls = 0;
@ -152,6 +155,9 @@ struct statistics {
st.update("arith-gcd-conflict", m_gcd_conflicts);
st.update("arith-cube-calls", m_cube_calls);
st.update("arith-cube-success", m_cube_success);
st.update("arith-lcube-calls", m_lcube_calls);
st.update("arith-lcube-success", m_lcube_success);
st.update("arith-lcube-flip-success", m_lcube_flip_success);
st.update("arith-patches", m_patches);
st.update("arith-patches-success", m_patches_success);
st.update("arith-hnf-calls", m_hnf_cutter_calls);
@ -252,15 +258,27 @@ private:
bool m_print_external_var_name = false;
bool m_propagate_eqs = false;
bool m_dio = false;
bool m_dio_enable_gomory_cuts = false;
bool m_dio_cuts_enable_gomory = false;
bool m_run_gomory_with_dio = false;
unsigned m_dio_gomory_enable_period = 16;
bool m_dio_enable_hnf_cuts = true;
bool m_dump_bound_lemmas = false;
bool m_dio_ignore_big_nums = false;
unsigned m_dio_calls_period = 4;
unsigned m_dio_calls_period_decrease = 2;
bool m_dio_run_gcd = true;
bool m_random_hammers = true;
bool m_lcube = true;
unsigned m_lcube_flips = 16;
public:
bool lcube() const { return m_lcube; }
unsigned lcube_flips() const { return m_lcube_flips; }
unsigned dio_calls_period() const { return m_dio_calls_period; }
unsigned & dio_calls_period() { return m_dio_calls_period; }
unsigned dio_calls_period_decrease() const { return m_dio_calls_period_decrease; }
unsigned & dio_calls_period_decrease() { return m_dio_calls_period_decrease; }
bool random_hammers() const { return m_random_hammers; }
bool & random_hammers() { return m_random_hammers; }
bool print_external_var_name() const { return m_print_external_var_name; }
bool propagate_eqs() const { return m_propagate_eqs;}
unsigned hnf_cut_period() const { return m_hnf_cut_period; }
@ -268,8 +286,19 @@ public:
unsigned random_next() { return m_rand(); }
unsigned random_next(unsigned u ) { return m_rand(u); }
bool dio() { return m_dio; }
bool & dio_enable_gomory_cuts() { return m_dio_enable_gomory_cuts; }
bool dio_enable_gomory_cuts() const { return m_dio && m_dio_enable_gomory_cuts; }
// Static config: did the user request Gomory cuts up front? (lp.dio_cuts_enable_gomory)
bool dio_cuts_enable_gomory() const { return m_dio_cuts_enable_gomory; }
// dio_calls_period at which the Diophantine back-off starts running Gomory (lp.dio_gomory_enable_period)
unsigned dio_gomory_enable_period() const { return m_dio_gomory_enable_period; }
// Runtime flag owned by the Diophantine controller, kept separate from the static
// config above so toggling it never clobbers the user's parameter: once dio has
// backed off enough it starts running Gomory cuts alongside dio, and a productive
// dio conflict stops them again.
void start_running_gomory_with_dio() { m_run_gomory_with_dio = true; }
void stop_running_gomory_with_dio() { m_run_gomory_with_dio = false; }
// Effective state read by should_gomory_cut(): allowed if either the user enabled it
// statically or the dio controller started running it, guarded by dio being active.
bool dio_enable_gomory_cuts() const { return m_dio && (m_dio_cuts_enable_gomory || m_run_gomory_with_dio); }
bool dio_run_gcd() const { return m_dio && m_dio_run_gcd; }
bool dio_enable_hnf_cuts() const { return m_dio && m_dio_enable_hnf_cuts; }
bool dio_ignore_big_nums() const { return m_dio_ignore_big_nums; }

193
src/math/lp/monomial_bounds.cpp
View file

@ -312,6 +312,12 @@ namespace nla {
}
dep.mul<dep_intervals::with_deps>(product, vi, product);
}
if (do_propagate_down && c().params().arith_nl_monomial_sandwich() &&
propagate_shared_factor(m))
return true;
if (c().params().arith_nl_monomial_binomial_sign() &&
propagate_binomial_sign(m))
return true;
return do_propagate_up && propagate_value(product, m.var());
}
@ -501,11 +507,196 @@ namespace nla {
}
lpvar monomial_bounds::non_fixed_var(monic const& m) {
for (lpvar v : m)
for (lpvar v : m)
if (!c().var_is_fixed(v))
return v;
return null_lpvar;
}
/**
* Dual-row shared-factor sandwich. For a binary monomial m = u*v, find LP
* term columns whose term has shape a_m * m + a_v * v (exactly two
* variables, both factors of m). The term column's bound is a sound
* interval for (a_m * m + a_v * v). Substituting m = u*v yields
* v * (a_m * u + a_v); dividing by the interval on v (sign-determined)
* gives an interval on (a_m * u + a_v), and an affine shift gives an
* interval on u. The derived interval is fed to the existing
* propagate_value path so the lemma channel and integer rounding are
* shared with the rest of the propagation pipeline.
*/
bool monomial_bounds::propagate_shared_factor(monic const& m) {
if (m.size() != 2)
return false;
lpvar f0 = m.vars()[0], f1 = m.vars()[1];
if (f0 == f1)
return false;
unsigned const fanout_limit = c().params().arith_nl_monomial_sandwich_max_fanout();
auto try_pair = [&](lpvar u, lpvar v) -> bool {
// Skip if u participates in too many monomials: tightening such a
// factor cascades through ord-binom / monotonicity on every monic
// that contains it.
if (fanout_limit > 0) {
unsigned fanout = 0;
for (auto const& m1 : c().emons().get_use_list(u)) {
(void)m1;
if (++fanout > fanout_limit)
return false;
}
}
scoped_dep_interval vi(dep);
var2interval(v, vi);
if (!dep.separated_from_zero(vi))
return false;
auto& lra = c().lra;
unsigned const ROW_CAP = 16;
unsigned scanned = 0;
for (auto const& cell : lra.A_r().m_columns[m.var()]) {
if (++scanned > ROW_CAP)
break;
unsigned basic = lra.get_base_column_in_row(cell.var());
if (basic == m.var() || basic == v || basic == u)
continue;
if (!lra.column_has_term(basic))
continue;
auto const& term = lra.get_term(basic);
if (term.size() != 2 ||
!term.contains(m.var()) || !term.contains(v))
continue;
rational const& a_m = term.get_coeff(m.var());
rational const& a_v = term.get_coeff(v);
if (a_m.is_zero())
continue;
// Term value = a_m*m + a_v*v; bound on basic bounds the term.
// Substituting m = u*v: term = v * (a_m*u + a_v).
scoped_dep_interval bi(dep);
var2interval(basic, bi);
scoped_dep_interval inner(dep);
dep.div<dep_intervals::with_deps>(bi, vi, inner);
scoped_dep_interval shift(dep);
dep.set_value(shift, -a_v);
scoped_dep_interval scaled(dep);
dep.add<dep_intervals::with_deps>(inner, shift, scaled);
scoped_dep_interval u_int(dep);
dep.mul<dep_intervals::with_deps>(rational::one() / a_m, scaled, u_int);
TRACE(nla_solver, tout << "sandwich shared-factor basic=" << basic
<< " m=" << m.var() << " v=" << v << " u=" << u
<< " a_m=" << a_m << " a_v=" << a_v << "\n";);
if (propagate_value(u_int, u))
return true; // one lemma per call to keep the channel quiet
}
return false;
};
return try_pair(f1, f0) || try_pair(f0, f1);
}
/**
* Sign-pinned binomial bound. For a binary monomial m = u*v in m_to_refine,
* use the current LP value mv = val(m.var()) as a one-sided anchor on the
* monomial value variable, and derive a deterministic interval for u via
* sign-aware division by v.
*
* Direction is chosen by the disagreement: if val(m.var()) > val(u)*val(v)
* the LP placed the monomial above the factor product, so we condition on
* "m.var() >= mv"; otherwise on "m.var() <= mv". The resulting clause is
* structurally analogous to a propagate_value lemma plus one extra
* snapshot literal on m.var(): under the asserted bounds on v, the clause
* reduces to a 2-disjunct (snapshot literal | factor bound).
*
* Targets the case ord-binom currently handles: factors have determined
* signs, m.var() may have no LP bound at all. The clause is sound modulo
* the monomial definition (the same condition propagate_down,
* propagate_shared_factor and ord-binom rely on).
*/
bool monomial_bounds::propagate_binomial_sign(monic const& m) {
if (m.size() != 2)
return false;
lpvar f0 = m.vars()[0], f1 = m.vars()[1];
if (f0 == f1)
return false;
rational const mv = c().val(m.var());
rational const fp = c().val(f0) * c().val(f1);
if (mv == fp)
return false;
bool const below = mv > fp; // LP placed m.var() too high
llc const anchor_cmp = below ? llc::LT : llc::GT;
auto try_anchor = [&](lpvar u, lpvar v) -> bool {
// Throttle once per (m.var(), u, v, direction) tuple. Without it
// each new val(m.var()) snapshot would re-emit and the search
// would cascade across model changes the same way ord-binom does.
if (c().throttle().insert_new(
nla_throttle::MONOMIAL_BINOMIAL_SIGN,
m.var(), u, v, below))
return false;
scoped_dep_interval vi(dep);
var2interval(v, vi);
if (!dep.separated_from_zero(vi))
return false;
// Synthesize a one-sided interval for m.var() at mv. No deps;
// the snapshot literal goes into the lemma body directly.
scoped_dep_interval mi_anchor(dep);
if (below) {
dep.set_lower(mi_anchor, mv);
dep.set_lower_is_inf(mi_anchor, false);
dep.set_lower_is_open(mi_anchor, false);
dep.set_upper_is_inf(mi_anchor, true);
} else {
dep.set_upper(mi_anchor, mv);
dep.set_upper_is_inf(mi_anchor, false);
dep.set_upper_is_open(mi_anchor, false);
dep.set_lower_is_inf(mi_anchor, true);
}
scoped_dep_interval u_int(dep);
dep.div<dep_intervals::with_deps>(mi_anchor, vi, u_int);
bool emitted = false;
if (should_propagate_lower(u_int, u, 1)) {
auto const& lower = dep.lower(u_int);
if (!is_too_big(lower)) {
auto cmp = dep.lower_is_open(u_int) ? llc::GT : llc::GE;
lp::explanation ex;
dep.get_lower_dep(u_int, ex);
lemma_builder lemma(c(), "binomial sign anchor");
lemma &= ex;
lemma |= ineq(m.var(), anchor_cmp, mv);
lemma |= ineq(u, cmp, lower);
emitted = true;
}
}
if (should_propagate_upper(u_int, u, 1)) {
auto const& upper = dep.upper(u_int);
if (!is_too_big(upper)) {
auto cmp = dep.upper_is_open(u_int) ? llc::LT : llc::LE;
lp::explanation ex;
dep.get_upper_dep(u_int, ex);
lemma_builder lemma(c(), "binomial sign anchor");
lemma &= ex;
lemma |= ineq(m.var(), anchor_cmp, mv);
lemma |= ineq(u, cmp, upper);
emitted = true;
}
}
return emitted;
};
return try_anchor(f1, f0) || try_anchor(f0, f1);
}
}

2
src/math/lp/monomial_bounds.h
View file

@ -33,6 +33,8 @@ namespace nla {
u_dependency* explain_fixed(monic const& m, rational const& k);
lp::explanation get_explanation(u_dependency* dep);
bool propagate_down(monic const& m, dep_interval& mi, lpvar v, unsigned power, dep_interval& product);
bool propagate_shared_factor(monic const& m);
bool propagate_binomial_sign(monic const& m);
void analyze_monomial(monic const& m, unsigned& num_free, lpvar& free_v, unsigned& power) const;
bool is_free(lpvar v) const;
bool is_zero(lpvar v) const;

3
src/math/lp/nla_coi.cpp
View file

@ -1,4 +1,3 @@
/*++
Copyright (c) 2025 Microsoft Corporation
@ -85,4 +84,4 @@ namespace nla {
}
}
}
}
}

6
src/math/lp/nla_coi.h
View file

@ -1,4 +1,3 @@
/*++
Copyright (c) 2025 Microsoft Corporation
@ -14,6 +13,9 @@
#pragma once
#include "util/uint_set.h"
#include "util/vector.h"
namespace nla {
class core;
@ -40,4 +42,4 @@ namespace nla {
indexed_uint_set const &vars() { return m_var_set; }
};
}
}

12
src/math/lp/nla_core.cpp
View file

@ -539,6 +539,7 @@ bool core::is_octagon_term(const lp::lar_term& t, bool & sign, lpvar& i, lpvar &
bool seen_minus = false;
bool seen_plus = false;
i = null_lpvar;
j = null_lpvar;
for(lp::lar_term::ival p : t) {
const auto & c = p.coeff();
if (c == 1) {
@ -1330,7 +1331,14 @@ lbool core::check(unsigned level) {
return l_false;
}
if (no_effect() && params().arith_nl_nra_check_assignment() && m_check_assignment_fail_cnt < params().arith_nl_nra_check_assignment_max_fail()) {
scoped_limits sl(m_reslim);
sl.push_child(&m_nra_lim);
ret = m_nra.check_assignment();
if (ret != l_true)
++m_check_assignment_fail_cnt;
}
if (no_effect() && should_run_bounded_nlsat())
ret = bounded_nlsat();
@ -1581,4 +1589,4 @@ void core::refine_pseudo_linear(monic const& m) {
}
SASSERT(nlvar != null_lpvar);
lemma |= ineq(lp::lar_term(m.var(), rational(-prod), nlvar), llc::EQ, rational(0));
}
}

13
src/math/lp/nla_core.h
View file

@ -63,6 +63,7 @@ class core {
unsigned m_nlsat_delay = 0;
unsigned m_nlsat_delay_bound = 0;
unsigned m_check_assignment_fail_cnt = 0;
bool should_run_bounded_nlsat();
lbool bounded_nlsat();
@ -94,6 +95,8 @@ class core {
emonics m_emons;
svector<lpvar> m_add_buffer;
mutable indexed_uint_set m_active_var_set;
// hook installed by theory_lra for creating a multiplication definition
std::function<lpvar(unsigned, lpvar const*)> m_add_mul_def_hook;
reslimit m_nra_lim;
@ -212,9 +215,12 @@ public:
void deregister_monic_from_tables(const monic & m, unsigned i);
void add_monic(lpvar v, unsigned sz, lpvar const* vs);
void add_idivision(lpvar q, lpvar x, lpvar y) { m_divisions.add_idivision(q, x, y); }
void add_rdivision(lpvar q, lpvar x, lpvar y) { m_divisions.add_rdivision(q, x, y); }
void add_bounded_division(lpvar q, lpvar x, lpvar y) { m_divisions.add_bounded_division(q, x, y); }
void add_idivision(lpvar q, lpvar x, lpvar y, lpvar r) { m_divisions.add_idivision(q, x, y, r); }
void add_rdivision(lpvar q, lpvar x, lpvar y, lpvar r) { m_divisions.add_rdivision(q, x, y, r); }
void add_bounded_division(lpvar q, lpvar x, lpvar y, lpvar r) { m_divisions.add_bounded_division(q, x, y, r); }
void add_divisibility(lpvar r, lpvar x, lpvar y, lpvar d) { m_divisions.add_divisibility(r, x, y, d); }
void set_add_mul_def_hook(std::function<lpvar(unsigned, lpvar const*)> const& f) { m_add_mul_def_hook = f; }
lpvar add_mul_def(unsigned sz, lpvar const* vs) { SASSERT(m_add_mul_def_hook); lpvar v = m_add_mul_def_hook(sz, vs); add_monic(v, sz, vs); return v; }
void set_relevant(std::function<bool(lpvar)>& is_relevant) { m_relevant = is_relevant; }
bool is_relevant(lpvar v) const { return !m_relevant || m_relevant(v); }
@ -478,4 +484,3 @@ inline std::ostream& operator<<(std::ostream& out, pp_factorization const& f) {
inline std::ostream& operator<<(std::ostream& out, pp_var const& v) { return v.c.print_var(v.v, out); }
} // end of namespace nla

192
src/math/lp/nla_divisions.cpp
View file

@ -18,29 +18,36 @@ Description:
namespace nla {
void divisions::add_idivision(lpvar q, lpvar x, lpvar y) {
if (x == null_lpvar || y == null_lpvar || q == null_lpvar)
void divisions::add_idivision(lpvar q, lpvar x, lpvar y, lpvar r) {
if (x == null_lpvar || y == null_lpvar || q == null_lpvar || r == null_lpvar)
return;
m_idivisions.push_back({q, x, y});
m_idivisions.push_back({q, x, y, r});
m_core.trail().push(push_back_vector(m_idivisions));
}
void divisions::add_rdivision(lpvar q, lpvar x, lpvar y) {
if (x == null_lpvar || y == null_lpvar || q == null_lpvar)
void divisions::add_rdivision(lpvar q, lpvar x, lpvar y, lpvar r) {
if (x == null_lpvar || y == null_lpvar || q == null_lpvar || r == null_lpvar)
return;
m_rdivisions.push_back({ q, x, y });
m_rdivisions.push_back({ q, x, y, r });
m_core.trail().push(push_back_vector(m_rdivisions));
}
void divisions::add_bounded_division(lpvar q, lpvar x, lpvar y) {
if (x == null_lpvar || y == null_lpvar || q == null_lpvar)
void divisions::add_bounded_division(lpvar q, lpvar x, lpvar y, lpvar r) {
if (x == null_lpvar || y == null_lpvar || q == null_lpvar || r == null_lpvar)
return;
if (m_core.lra.column_has_term(x) || m_core.lra.column_has_term(y) || m_core.lra.column_has_term(q))
return;
m_bounded_divisions.push_back({ q, x, y });
m_bounded_divisions.push_back({ q, x, y, r });
m_core.trail().push(push_back_vector(m_bounded_divisions));
}
void divisions::add_divisibility(lpvar r, lpvar x, lpvar y, lpvar d) {
if (x == null_lpvar || y == null_lpvar || r == null_lpvar || d == null_lpvar)
return;
m_divisibility.push_back({ r, x, y, d });
m_core.trail().push(push_back_vector(m_divisibility));
}
typedef lp::lar_term term;
// y1 >= y2 > 0 & x1 <= x2 => x1/y1 <= x2/y2
@ -111,7 +118,7 @@ namespace nla {
return false;
};
for (auto const & [r, x, y] : m_idivisions) {
for (auto const & [r, x, y, md] : m_idivisions) {
if (!c.is_relevant(r))
continue;
auto xval = c.val(x);
@ -120,7 +127,7 @@ namespace nla {
// idiv semantics
if (!xval.is_int() || !yval.is_int() || yval == 0 || rval == div(xval, yval))
continue;
for (auto const& [q2, x2, y2] : m_idivisions) {
for (auto const& [q2, x2, y2, md2] : m_idivisions) {
if (q2 == r)
continue;
if (!c.is_relevant(q2))
@ -133,7 +140,7 @@ namespace nla {
}
}
for (auto const& [r, x, y] : m_rdivisions) {
for (auto const& [r, x, y, md] : m_rdivisions) {
if (!c.is_relevant(r))
continue;
auto xval = c.val(x);
@ -142,7 +149,7 @@ namespace nla {
// / semantics
if (yval == 0 || rval == xval / yval)
continue;
for (auto const& [q2, x2, y2] : m_rdivisions) {
for (auto const& [q2, x2, y2, md2] : m_rdivisions) {
if (q2 == r)
continue;
if (!c.is_relevant(q2))
@ -154,7 +161,10 @@ namespace nla {
return;
}
}
check_mod_mult();
check_linear_divisibility();
check_mod_congruence();
}
// if p is bounded, q a value, r = eval(p):
@ -163,11 +173,11 @@ namespace nla {
void divisions::check_bounded_divisions() {
core& c = m_core;
unsigned offset = c.random(), sz = m_bounded_divisions.size();
unsigned offset = c.random(), sz = m_bounded_divisions.size();
for (unsigned j = 0; j < sz; ++j) {
unsigned i = (offset + j) % sz;
auto [q, x, y] = m_bounded_divisions[i];
auto [q, x, y, r] = m_bounded_divisions[i];
if (!c.is_relevant(q))
continue;
auto xv = c.val(x);
@ -188,9 +198,9 @@ namespace nla {
rational lo = yv * div_v;
if (xv > hi) {
lemma_builder lemma(c, "y = yv & x <= yv * div(xv, yv) + yv - 1 => div(p, y) <= div(xv, yv)");
lemma |= ineq(y, llc::NE, yv);
lemma |= ineq(x, llc::GT, hi);
lemma |= ineq(q, llc::LE, div_v);
lemma |= ineq(y, llc::NE, yv);
lemma |= ineq(x, llc::GT, hi);
lemma |= ineq(q, llc::LE, div_v);
return;
}
if (xv < lo) {
@ -201,5 +211,147 @@ namespace nla {
return;
}
}
}
}
// mod(factor, p) = 0 => mod(factor * k, p) = 0
// For each division (q, x, y, r) where x is a monic m = f1 * f2 * ... * fk,
// if some factor fi has mod(fi, p) = 0 (fixed), then mod(x, p) = 0.
void divisions::check_mod_mult() {
core& c = m_core;
unsigned offset = c.random(), sz = m_bounded_divisions.size();
for (unsigned j = 0; j < sz; ++j) {
unsigned i = (offset + j) % sz;
auto [q, x, y, r] = m_bounded_divisions[i];
if (!c.is_relevant(q))
continue;
if (c.var_is_fixed_to_zero(r))
continue;
if (c.val(r).is_zero())
continue;
if (!c.is_monic_var(x))
continue;
auto yv = c.val(y);
if (yv <= 0 || !yv.is_int())
continue;
auto const& m = c.emons()[x];
for (lpvar f : m.vars()) {
for (auto const& [q2, x2, y2, r2] : m_bounded_divisions) {
if (x2 != f)
continue;
if (c.val(y2) != yv)
continue;
if (!c.var_is_fixed_to_zero(r2))
continue;
// mod(factor, p) = 0 => mod(product, p) = 0
lemma_builder lemma(c, "mod(factor, p) = 0 => mod(factor * k, p) = 0");
lemma |= ineq(r2, llc::NE, 0);
lemma |= ineq(r, llc::EQ, 0);
return;
}
}
}
}
// Linear divisibility closure:
// mod(a, y) = 0 & x = c * a (c an integer constant) => mod(x, y) = 0.
// The emitted clause
// (x - c*a != 0) \/ (mod(a, y) != 0) \/ (mod(x, y) = 0)
// is a tautology for every integer c (under the Euclidean semantics of mod),
// so the choice of c/a from the current model can never be unsound. We only
// emit it when all three literals are false in the current model, which makes
// the clause a real conflict/propagation and guarantees progress.
void divisions::check_linear_divisibility() {
core& c = m_core;
unsigned sz = m_divisibility.size();
for (unsigned i = 0; i < sz; ++i) {
auto const& [rx, x, y, dx] = m_divisibility[i];
if (!c.is_relevant(rx))
continue;
if (c.val(rx).is_zero()) // mod(x, y) already 0 in model: nothing to refute
continue;
auto xval = c.val(x);
if (xval.is_zero())
continue;
for (unsigned j = 0; j < sz; ++j) {
if (i == j)
continue;
auto const& [ra, a, y2, da] = m_divisibility[j];
if (y2 != y && c.val(y2) != c.val(y)) // same divisor (by column or value)
continue;
if (!c.is_relevant(ra))
continue;
if (!c.val(ra).is_zero()) // need mod(a, y) = 0 in model
continue;
auto aval = c.val(a);
if (aval.is_zero())
continue;
rational cc = xval / aval;
if (!cc.is_int() || cc.is_zero())
continue;
if (xval != cc * aval) // ensure x = c*a holds exactly in the model
continue;
lemma_builder lemma(c, "mod(a,y) = 0 & x = c*a => mod(x,y) = 0");
lemma |= ineq(term(x, -cc, a), llc::NE, 0); // x - c*a != 0
lemma |= ineq(ra, llc::NE, 0); // mod(a, y) != 0
lemma |= ineq(rx, llc::EQ, 0); // mod(x, y) = 0
return;
}
}
}
// Modular congruence over a shared (possibly symbolic) divisor.
//
// For each divisibility fact we have the Euclidean identities (asserted by
// theory_lra::mk_idiv_mod_axioms):
// x = y * div(x,y) + mod(x,y), 0 <= mod(x,y) < |y|.
// For two facts (rx = mod(x,y), dx = div(x,y)) and (rs = mod(s,y), ds = div(s,y))
// sharing divisor y, subtracting the identities gives, for every integer delta,
// div(x,y) - div(s,y) = delta => mod(x,y) - mod(s,y) = (x - s) - delta*y.
// This is a tautology (entailed by the two identities) for any fixed integer
// delta, so choosing delta from the current model can never be unsound. We emit
// the clause
// (div(x,y) - div(s,y) != delta) \/ (mod(x,y) - mod(s,y) - (x - s) + delta*y = 0)
// only when the equality literal is false in the model (delta taken as the model
// value of div(x,y) - div(s,y)), which makes the clause a real propagation and
// guarantees progress. This discharges linear congruences with a symbolic
// modulus (e.g. mod(i + s, n) = i + mod(s, n)) that the nonlinear core does not
// otherwise isolate.
void divisions::check_mod_congruence() {
core& c = m_core;
unsigned sz = m_divisibility.size();
for (unsigned i = 0; i < sz; ++i) {
auto const& [rx, x, y, dx] = m_divisibility[i];
if (!c.is_relevant(rx))
continue;
auto yval = c.val(y);
if (yval.is_zero()) // mod/div uninterpreted when the divisor is 0
continue;
for (unsigned j = i + 1; j < sz; ++j) {
auto const& [rs, s, y2, ds] = m_divisibility[j];
if (!c.is_relevant(rs))
continue;
if (y2 != y && c.val(y2) != yval) // same divisor (by column or value)
continue;
rational delta = c.val(dx) - c.val(ds);
rational lhs = c.val(rx) - c.val(rs);
rational rhs = (c.val(x) - c.val(s)) - delta * yval;
if (lhs == rhs) // residue equation already holds: nothing to propagate
continue;
lemma_builder lemma(c, "y != 0 & y = y2 & div(x,y) - div(s,y) = delta => mod(x,y) - mod(s,y) = (x - s) - delta*y");
lemma |= ineq(y, llc::EQ, 0); // y = 0 (guard: mod/div uninterpreted when divisor is 0)
if (y2 != y)
lemma |= ineq(term(y, rational(-1), y2), llc::NE, 0); // y != y2 (guard: divisors must coincide symbolically)
lemma |= ineq(term(dx, rational(-1), ds), llc::NE, delta); // div(x,y) - div(s,y) != delta
term t;
t.add_monomial(rational::one(), rx);
t.add_monomial(rational(-1), rs);
t.add_monomial(rational(-1), x);
t.add_monomial(rational::one(), s);
t.add_monomial(delta, y);
lemma |= ineq(t, llc::EQ, 0); // mod(x,y) - mod(s,y) - x + s + delta*y = 0
return;
}
}
}
}

20
src/math/lp/nla_divisions.h
View file

@ -22,16 +22,22 @@ namespace nla {
class divisions {
core& m_core;
vector<std::tuple<lpvar, lpvar, lpvar>> m_idivisions;
vector<std::tuple<lpvar, lpvar, lpvar>> m_rdivisions;
vector<std::tuple<lpvar, lpvar, lpvar>> m_bounded_divisions;
vector<std::tuple<lpvar, lpvar, lpvar, lpvar>> m_idivisions;
vector<std::tuple<lpvar, lpvar, lpvar, lpvar>> m_rdivisions;
vector<std::tuple<lpvar, lpvar, lpvar, lpvar>> m_bounded_divisions;
// divisibility facts (r, x, y, d) meaning r = mod(x, y) and d = div(x, y)
vector<std::tuple<lpvar, lpvar, lpvar, lpvar>> m_divisibility;
public:
divisions(core& c):m_core(c) {}
void add_idivision(lpvar q, lpvar x, lpvar y);
void add_rdivision(lpvar q, lpvar x, lpvar y);
void add_bounded_division(lpvar q, lpvar x, lpvar y);
void add_idivision(lpvar q, lpvar x, lpvar y, lpvar r);
void add_rdivision(lpvar q, lpvar x, lpvar y, lpvar r);
void add_bounded_division(lpvar q, lpvar x, lpvar y, lpvar r);
void add_divisibility(lpvar r, lpvar x, lpvar y, lpvar d);
void check();
void check_bounded_divisions();
void check_mod_mult();
void check_linear_divisibility();
void check_mod_congruence();
};
}

194
src/math/lp/nla_grobner.cpp
View file

@ -10,6 +10,8 @@ Author:
Nikolaj Bjorner (nbjorner)
--*/
#include <algorithm>
#include <climits>
#include "util/uint_set.h"
#include "params/smt_params_helper.hpp"
#include "math/lp/nla_core.h"
@ -77,56 +79,75 @@ namespace nla {
if (!configure())
return;
bool productive = false;
try {
if (propagate_gcd_test())
return;
productive = true;
}
catch (...) {
}
m_solver.saturate();
TRACE(grobner, m_solver.display(tout));
if (!productive) {
m_solver.saturate();
TRACE(grobner, m_solver.display(tout));
if (m_delay_base > 0)
--m_delay_base;
try {
if (m_delay_base > 0)
--m_delay_base;
if (is_conflicting())
return;
try {
productive = is_conflicting()
|| propagate_quotients()
|| propagate_gcd_test()
|| propagate_eqs()
|| propagate_factorization()
|| propagate_linear_equations();
}
catch (...) {
if (propagate_quotients())
return;
if (propagate_gcd_test())
return;
if (propagate_eqs())
return;
if (propagate_factorization())
return;
if (propagate_linear_equations())
return;
}
catch (...) {
}
}
// DEBUG_CODE(for (auto e : m_solver.equations()) check_missing_propagation(*e););
if (c().params().arith_nl_grobner_adaptive())
update_growth_boost(productive);
if (productive)
return;
// for (auto e : m_solver.equations()) check_missing_propagation(*e);
++m_delay_base;
if (m_quota > 0)
--m_quota;
--m_quota;
IF_VERBOSE(5, verbose_stream() << "grobner miss, quota " << m_quota << "\n");
IF_VERBOSE(5, diagnose_pdd_miss(verbose_stream()));
}
void grobner::update_growth_boost(bool productive) {
// Bumping is conservative: requires two consecutive productive runs
// before any boost; misses decay toward unit by 1/4 per call.
unsigned const unit = m_config.m_adaptive_unit;
unsigned const cap = m_config.m_adaptive_max;
if (productive) {
++m_hit_streak;
if (m_hit_streak >= m_config.m_adaptive_bump_after) {
unsigned next = m_growth_boost + (m_growth_boost >> 1);
m_growth_boost = std::min(next, cap);
m_hit_streak = 0;
}
}
else {
m_hit_streak = 0;
if (m_growth_boost > unit) {
unsigned excess = m_growth_boost - unit;
m_growth_boost -= (excess + 3) / 4;
if (m_growth_boost < unit)
m_growth_boost = unit;
}
}
IF_VERBOSE(5, verbose_stream() << "grobner adaptive boost " << m_growth_boost
<< "/" << unit << (productive ? " (hit)" : " (miss)") << "\n");
}
bool grobner::is_conflicting() {
for (auto eq : m_solver.equations()) {
if (is_conflicting(*eq)) {
@ -210,8 +231,6 @@ namespace nla {
if (vars.empty() || !q.is_linear())
return false;
// IF_VERBOSE(0, verbose_stream() << "factored " << q << " : " << vars << "\n");
auto [t, offset] = linear_to_term(q);
vector<ineq> ineqs;
@ -226,7 +245,6 @@ namespace nla {
add_dependencies(lemma, eq);
for (auto const& i : ineqs)
lemma |= i;
//lemma.display(verbose_stream());
return true;
}
@ -368,6 +386,70 @@ namespace nla {
nl_vars.insert(j);
}
// mod_residue: derive v's residue mod M from polynomial divisibility.
//
// Common case. Given polynomial
// p = M*v1 + v - M*v2*v3 = 0,
// every monomial except v is M-divisible, so v ≡ 0 (mod M).
// Combined with 0 ≤ v < M, this forces v = 0.
// Emit: dependencies => (v < 0) ∨ (v ≥ M) ∨ (v = 0).
//
// General case. For a linear monomial c_v*v in p with c0 the constant
// term, require c_i/c_v integer for every non-v monomial and c0/c_v
// integer (call it K). Let M = gcd(|c_i/c_v|) over non-v monomials.
// Then p/c_v gives v + M*Q + K = 0 with Q integer, so v ≡ -K (mod M).
// With target = (-K) mod M ∈ [0, M-1], emit
// dependencies => (v < 0) ∨ (v ≥ M) ∨ (v = target).
for (auto const& mv : p) {
if (mv.vars.size() != 1)
continue;
lpvar vv = mv.vars[0];
if (!c().var_is_int(vv))
continue;
rational c_v = mv.coeff;
SASSERT(c_v != 0);
rational M(0); // 0 sentinel: "no non-v non-constant monomial seen yet".
rational c0(0);
bool ok = true;
for (auto const& mi : p) {
if (mi.vars.size() == 1 && mi.vars[0] == vv)
continue; // skip the mv monomial itself
if (mi.vars.empty()) {
c0 = mi.coeff;
continue;
}
rational quot = mi.coeff / c_v;
if (!quot.is_int()) { ok = false; break; }
rational a = abs(quot);
SASSERT(a != 0);
M = M == 0 ? a : gcd(M, a);
if (M == 1) { ok = false; break; } // trivial modulus, abort
}
if (!ok || M == 0)
continue;
rational K = c0 / c_v;
if (!K.is_int())
continue;
rational target = mod(-K, M); // Euclidean: result in [0, M-1].
SASSERT(target >= 0 && target < M);
// Skip if the lemma is already satisfied by the current model:
// any of (v < 0), (v ≥ M), (v = target) trivially holding means
// emission would be redundant. Without this guard, the lemma
// re-emits every Grobner round on the same polynomial.
rational v_val = c().val(vv);
if (v_val < 0 || v_val >= M || v_val == target)
continue;
lemma_builder lemma(c(), "grobner-mod-residue");
add_dependencies(lemma, eq);
lemma |= ineq(vv, llc::LT, rational::zero());
lemma |= ineq(vv, llc::GE, M);
lemma |= ineq(vv, llc::EQ, target);
TRACE(grobner, lemma.display(tout << "mod_residue v=" << vv
<< " M=" << M << " c_v=" << c_v << " c0=" << c0
<< " target=" << target << "\n"));
return true;
}
bool found_lemma = false;
for (auto v : nl_vars) {
auto& m = p.manager();
@ -559,25 +641,27 @@ namespace nla {
}
TRACE(grobner, m_solver.display(tout));
#if 0
IF_VERBOSE(2, m_pdd_grobner.display(verbose_stream()));
dd::pdd_eval eval(m_pdd_manager);
eval.var2val() = [&](unsigned j){ return val(j); };
for (auto* e : m_pdd_grobner.equations()) {
dd::pdd p = e->poly();
rational v = eval(p);
if (p.is_linear() && !eval(p).is_zero()) {
IF_VERBOSE(0, verbose_stream() << "violated linear constraint " << p << "\n");
}
}
#endif
struct dd::solver::config cfg;
cfg.m_max_steps = m_solver.equations().size();
cfg.m_max_simplified = c().params().arith_nl_grobner_max_simplified();
cfg.m_eqs_growth = c().params().arith_nl_grobner_eqs_growth();
cfg.m_expr_size_growth = c().params().arith_nl_grobner_expr_size_growth();
cfg.m_expr_degree_growth = c().params().arith_nl_grobner_expr_degree_growth();
if (c().params().arith_nl_grobner_adaptive() && m_growth_boost != m_config.m_adaptive_unit) {
// Wider intermediate to prevent overflow when a user param is
// close to UINT_MAX; clamp before assigning back to the unsigned
// config fields.
uint64_t const unit = m_config.m_adaptive_unit;
uint64_t const boost = m_growth_boost;
auto scale = [unit, boost](unsigned x) -> unsigned {
uint64_t y = (static_cast<uint64_t>(x) * boost) / unit;
return y > UINT_MAX ? UINT_MAX : static_cast<unsigned>(y);
};
cfg.m_eqs_growth = scale(cfg.m_eqs_growth);
cfg.m_expr_size_growth = scale(cfg.m_expr_size_growth);
cfg.m_expr_degree_growth = scale(cfg.m_expr_degree_growth);
cfg.m_max_simplified = scale(cfg.m_max_simplified);
}
cfg.m_number_of_conflicts_to_report = c().params().arith_nl_grobner_cnfl_to_report();
m_solver.set(cfg);
m_solver.adjust_cfg();
@ -588,14 +672,12 @@ namespace nla {
std::ostream& grobner::diagnose_pdd_miss(std::ostream& out) {
// m_pdd_grobner.display(out);
dd::pdd_eval eval;
eval.var2val() = [&](unsigned j){ return val(j); };
for (auto* e : m_solver.equations()) {
dd::pdd p = e->poly();
rational v = eval(p);
if (!v.is_zero()) {
if (v != 0) {
out << p << " := " << v << "\n";
}
}
@ -701,7 +783,15 @@ namespace nla {
lp::lpvar j = c().lra.add_term(coeffs, UINT_MAX);
c().lra.update_column_type_and_bound(j, lp::lconstraint_kind::EQ, offset, e.dep());
c().m_check_feasible = true;
c().m_check_feasible = true;
TRACE(nla_solver,
// Print the term as installed (post subst_known_terms), not the
// pre-add_term coeffs vector. add_term normalizes/substitutes
// term-column references, so coeffs and the resulting row can
// diverge if any var is itself a term-column.
tout << "grobner-linear-eq: ";
c().lra.print_term(c().lra.get_term(j), tout);
tout << " = " << offset << "\n";);
return true;
}

13
src/math/lp/nla_grobner.h
View file

@ -12,6 +12,7 @@
#include "math/lp/nla_intervals.h"
#include "math/lp/nex.h"
#include "math/lp/cross_nested.h"
#include "util/params.h"
#include "util/uint_set.h"
#include "math/grobner/pdd_solver.h"
@ -23,7 +24,13 @@ namespace nla {
bool m_propagate_quotients = false;
bool m_gcd_test = false;
bool m_expand_terms = false;
// Adaptive growth (gated by arith.nl.grobner_adaptive). m_growth_boost
// is in fixed-point units of 1/m_adaptive_unit (m_adaptive_unit == 1.0x).
unsigned m_adaptive_unit = 16;
unsigned m_adaptive_max = 4 * 16;
unsigned m_adaptive_bump_after = 2;
};
config m_config;
dd::pdd_manager m_pdd_manager;
dd::solver m_solver;
lp::lar_solver& lra;
@ -31,8 +38,9 @@ namespace nla {
unsigned m_quota = 0;
unsigned m_delay_base = 0;
unsigned m_delay = 0;
unsigned m_growth_boost = m_config.m_adaptive_unit;
unsigned m_hit_streak = 0;
bool m_add_all_eqs = false;
config m_config;
std::unordered_map<unsigned_vector, lpvar, hash_svector> m_mon2var;
lp::lp_settings& lp_settings();
@ -69,6 +77,9 @@ namespace nla {
bool equation_is_true(dd::solver::equation const& eq);
// adaptive growth (gated by arith.nl.grobner_adaptive)
void update_growth_boost(bool productive);
// setup
bool configure();
void set_level2var();

4
src/math/lp/nla_order_lemmas.cpp
View file

@ -81,9 +81,11 @@ void order::order_lemma_on_binomial(const monic& ac) {
*/
void order::order_lemma_on_binomial_sign(const monic& xy, lpvar x, lpvar y, int sign) {
if (!c().params().arith_nl_order_binomial_sign())
return;
if (!c().var_is_int(x) && val(x).is_big())
return;
SASSERT(!_().mon_has_zero(xy.vars()));
int sy = rat_sign(val(y));

2
src/math/lp/nla_pp.cpp
View file

@ -432,4 +432,4 @@ std::ostream& core::display_constraint_smt(std::ostream& out, unsigned id, lp::l
out << (evaluation ? "true" : "false");
out << "\n";
return out;
}
}

16
src/math/lp/nla_solver.cpp
View file

@ -20,16 +20,20 @@ namespace nla {
m_core->add_monic(v, sz, vs);
}
void solver::add_idivision(lpvar q, lpvar x, lpvar y) {
m_core->add_idivision(q, x, y);
void solver::add_idivision(lpvar q, lpvar x, lpvar y, lpvar r) {
m_core->add_idivision(q, x, y, r);
}
void solver::add_rdivision(lpvar q, lpvar x, lpvar y) {
m_core->add_rdivision(q, x, y);
void solver::add_rdivision(lpvar q, lpvar x, lpvar y, lpvar r) {
m_core->add_rdivision(q, x, y, r);
}
void solver::add_bounded_division(lpvar q, lpvar x, lpvar y) {
m_core->add_bounded_division(q, x, y);
void solver::add_bounded_division(lpvar q, lpvar x, lpvar y, lpvar r) {
m_core->add_bounded_division(q, x, y, r);
}
void solver::add_divisibility(lpvar r, lpvar x, lpvar y, lpvar d) {
m_core->add_divisibility(r, x, y, d);
}
void solver::set_relevant(std::function<bool(lpvar)>& is_relevant) {

7
src/math/lp/nla_solver.h
View file

@ -28,9 +28,10 @@ namespace nla {
~solver();
const auto& monics_with_changed_bounds() const { return m_core->monics_with_changed_bounds(); }
void add_monic(lpvar v, unsigned sz, lpvar const* vs);
void add_idivision(lpvar q, lpvar x, lpvar y);
void add_rdivision(lpvar q, lpvar x, lpvar y);
void add_bounded_division(lpvar q, lpvar x, lpvar y);
void add_idivision(lpvar q, lpvar x, lpvar y, lpvar r);
void add_rdivision(lpvar q, lpvar x, lpvar y, lpvar r);
void add_bounded_division(lpvar q, lpvar x, lpvar y, lpvar r);
void add_divisibility(lpvar r, lpvar x, lpvar y, lpvar d);
void check_bounded_divisions();
void set_relevant(std::function<bool(lpvar)>& is_relevant);
void updt_params(params_ref const& p);

62
src/math/lp/nla_tangent_lemmas.cpp
View file

@ -18,6 +18,12 @@ class tangent_imp {
rational m_correct_v;
// "below" means that the incorrect value is less than the correct one, that is m_v < m_correct_v
bool m_below;
// pl is in the strict interior of the bound box (model-driven points
// get_initial_points + push_point); McCormick at the box corners
// requires non-strict inequality because the tangent meets the surface
// along the box's edges (xy = pl.y*x + pl.x*y - pl.x*pl.y at x = pl.x
// or y = pl.y).
bool m_pl_strict_interior = true;
rational m_v; // the monomial value
lpvar m_j; // the monic variable
const monic& m_m;
@ -89,7 +95,10 @@ private:
t.add_monomial(- m_y.rat_sign()*pl.x, m_jy);
t.add_monomial(- m_x.rat_sign()*pl.y, m_jx);
t.add_var(m_j);
lemma |= ineq(t, m_below? llc::GT : llc::LT, - pl.x*pl.y);
llc cmp = m_below
? (m_pl_strict_interior ? llc::GT : llc::GE)
: (m_pl_strict_interior ? llc::LT : llc::LE);
lemma |= ineq(t, cmp, - pl.x*pl.y);
explain(lemma);
}
@ -164,14 +173,61 @@ private:
return a.x * m_xy.y + a.y * m_xy.x - a.x * a.y;
}
// McCormick at box corners: choose m_a, m_b at the corners of
// [x_lo, x_hi] x [y_lo, y_hi] that bound xy from the side dictated by
// m_below. Returns false if either factor has an unbounded side, the
// box is degenerate, or the current LP value of a factor coincides with
// a chosen corner — generate_plane's negate_relation requires
// val(j) != corner_coord (SASSERT in debug; trivially-true literal in
// release). The caller falls back to the model-driven point selection in
// these cases.
bool set_box_corners() {
if (!c().has_lower_bound(m_jx) || !c().has_upper_bound(m_jx))
return false;
if (!c().has_lower_bound(m_jy) || !c().has_upper_bound(m_jy))
return false;
rational const& x_lo = c().get_lower_bound(m_jx);
rational const& x_hi = c().get_upper_bound(m_jx);
rational const& y_lo = c().get_lower_bound(m_jy);
rational const& y_hi = c().get_upper_bound(m_jy);
if (x_lo == x_hi || y_lo == y_hi)
return false;
// negate_relation requires the model value to be strictly separated
// from the corner coordinate it's compared to. If LP currently sits
// exactly at a box edge, fall back.
rational const& vx = c().val(m_jx);
rational const& vy = c().val(m_jy);
if (vx == x_lo || vx == x_hi || vy == y_lo || vy == y_hi)
return false;
if (m_below) {
// Under-approximation: tangents at (x_lo, y_lo) and (x_hi, y_hi)
// bound xy from below across the box.
m_a = point(x_lo, y_lo);
m_b = point(x_hi, y_hi);
} else {
// Over-approximation: anti-diagonal corners.
m_a = point(x_lo, y_hi);
m_b = point(x_hi, y_lo);
}
m_pl_strict_interior = false;
return true;
}
void get_points() {
if (c().params().arith_nl_tangents_box_corners() && set_box_corners()) {
// Box corners are extremes; pushing further moves out of the box
// and would invalidate the McCormick property.
TRACE(nla_solver, tout << "xy = " << m_xy << ", box-corner points: ";
print_tangent_domain(tout) << std::endl;);
return;
}
get_initial_points();
TRACE(nla_solver, tout << "xy = " << m_xy << ", correct val = " << m_correct_v;
print_tangent_domain(tout << "\ntang points:") << std::endl;);
push_point(m_a);
push_point(m_a);
push_point(m_b);
TRACE(nla_solver,
tout << "pushed a = " << m_a << std::endl
tout << "pushed a = " << m_a << std::endl
<< "pushed b = " << m_b << std::endl
<< "tang_plane(a) = " << tang_plane(m_a) << " , val = " << m_a << ", "
<< "tang_plane(b) = " << tang_plane(m_b) << " , val = " << m_b << std::endl;);

5
src/math/lp/nla_throttle.h
View file

@ -18,9 +18,10 @@ class nla_throttle {
public:
enum throttle_kind {
ORDER_LEMMA, // order lemma (9 params)
BINOMIAL_SIGN_LEMMA, // binomial sign (6 params)
BINOMIAL_SIGN_LEMMA, // binomial sign (6 params)
MONOTONE_LEMMA, // monotonicity (2 params)
TANGENT_LEMMA // tangent lemma (5 params: monic_var, x_var, y_var, below, plane_type)
TANGENT_LEMMA, // tangent lemma (5 params: monic_var, x_var, y_var, below, plane_type)
MONOMIAL_BINOMIAL_SIGN // monomial binomial sign anchor (4 params: monic_var, u, v, below)
};
private:

11
src/math/lp/nla_types.h
View file

@ -41,12 +41,11 @@ namespace nla {
ineq(const lp::lar_term& term, lp::lconstraint_kind cmp, const rational& rs) : m_cmp(cmp), m_term(term), m_rs(rs) {}
ineq(lpvar v, lp::lconstraint_kind cmp, int i): m_cmp(cmp), m_term(v), m_rs(rational(i)) {}
ineq(lpvar v, lp::lconstraint_kind cmp, rational const& r): m_cmp(cmp), m_term(v), m_rs(r) {}
bool operator==(const ineq& a) const {
return m_cmp == a.m_cmp && m_term == a.m_term && m_rs == a.m_rs;
}
const lp::lar_term& term() const { return m_term; };
lp::lconstraint_kind cmp() const { return m_cmp; };
const rational& rs() const { return m_rs; };
bool operator==(const ineq& a) const = delete;
bool operator!=(const ineq& a) const = delete;
const lp::lar_term& term() const { return m_term; }
lp::lconstraint_kind cmp() const { return m_cmp; }
const rational& rs() const { return m_rs; }
};
class lemma {

337
src/math/lp/nra_solver.cpp
View file

@ -11,6 +11,7 @@
#include "math/lp/nra_solver.h"
#include "math/lp/nla_coi.h"
#include "nlsat/nlsat_solver.h"
#include "nlsat/nlsat_assignment.h"
#include "math/polynomial/polynomial.h"
#include "math/polynomial/algebraic_numbers.h"
#include "util/map.h"
@ -35,6 +36,13 @@ struct solver::imp {
scoped_ptr<scoped_anum_vector> m_values; // values provided by LRA solver
scoped_ptr<scoped_anum> m_tmp1, m_tmp2;
nla::coi m_coi;
svector<lp::constraint_index> m_literal2constraint;
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;
nla::core& m_nla_core;
imp(lp::lar_solver& s, reslimit& lim, params_ref const& p, nla::core& nla_core):
@ -56,8 +64,10 @@ struct solver::imp {
m_lp2nl.reset();
}
// Create polynomial definition for variable v used in setup_assignment_solver.
// Side-effects: updates m_vars2mon when v is a monic variable.
// Create polynomial definition for variable v used in setup_solver_poly.
// The definition recursively expands monic and term variables into
// polynomials in leaf variables, scaled by an integer denominator
// tracked in `denominators` to keep the coefficients integral.
void mk_definition(unsigned v, polynomial_ref_vector &definitions, vector<rational>& denominators) {
auto &pm = m_nlsat->pm();
polynomial::polynomial_ref p(pm);
@ -214,20 +224,9 @@ struct solver::imp {
out.close();
}
lbool r = l_undef;
statistics& st = m_nla_core.lp_settings().stats().m_st;
try {
r = m_nlsat->check();
}
catch (z3_exception&) {
if (m_limit.is_canceled()) {
r = l_undef;
}
else {
m_nlsat->collect_statistics(st);
throw;
}
}
lbool r = m_nlsat->check();
m_nlsat->collect_statistics(st);
TRACE(nra, tout << "nra result " << r << "\n");
CTRACE(nra, false,
@ -273,7 +272,220 @@ struct solver::imp {
break;
}
return r;
}
}
void setup_assignment_solver() {
SASSERT(need_check());
reset();
m_literal2constraint.reset();
m_vars2mon.reset();
m_coi.init();
auto &pm = m_nlsat->pm();
polynomial_ref_vector definitions(pm);
vector<rational> denominators;
// Create an NLSAT polyvar for each LRA variable (identity mapping),
// seed the assignment from the current LRA model, populate
// m_vars2mon, and build the inlined polynomial definition of v.
//
// The definition expands monic and term variables into polynomials
// over leaf variables. Each definition is scaled by denominators[v]
// so that all coefficients stay integral; the scaling cancels on
// both sides of every constraint we build below (just like in
// setup_solver_poly).
//
// This "de-linearized" representation is what the linear-cell
// construction in NLSAT needs: a cell built around a constraint
// polynomial that mentions several multiplications at once can
// yield a lemma constraining all of them simultaneously, which is
// strictly stronger than the per-multiplication lemmas we would
// get from asserting `v_mon - v1*...*vk = 0` separately.
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);
scoped_anum a(am());
am().set(a, m_nla_core.val(v).to_mpq());
m_values->push_back(a);
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);
}
mk_definition(v, definitions, denominators);
}
// Substitute each variable in the LRA constraint by its definition
// and rescale to keep integer coefficients. Symbolically:
//
// v == definitions[v] / denominators[v]
//
// sum(coeff_v * v) k rhs
// == sum((coeff_v / denominators[v]) * definitions[v]) k rhs
//
// We pick den := lcm of all denominators(coeff_v / denominators[v])
// together with denominator(rhs), so that den * coeff_v / denominators[v]
// and den * rhs are all integers. The relation kind k is preserved
// because den > 0.
for (auto ci : m_coi.constraints()) {
auto &c = lra.constraints()[ci];
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 / denominators[v]));
polynomial::polynomial_ref p(pm);
p = pm.mk_const(-den * rhs);
for (auto [coeff, v] : lhs) {
polynomial_ref poly(pm);
poly = definitions.get(v);
poly = poly * constant(den * coeff / denominators[v]);
p = p + poly;
}
auto lit = add_constraint(p, ci, k);
m_literal2constraint.setx(lit.index(), ci, lp::null_ci);
}
definitions.reset();
}
void process_polynomial_check_assignment(polynomial::polynomial const* p, rational& bound, const u_map<lp::lpvar>& nl2lp, lp::lar_term& t) {
polynomial::manager& pm = m_nlsat->pm();
for (unsigned i = 0; i < pm.size(p); ++i) {
polynomial::monomial* m = pm.get_monomial(p, i);
auto& coeff = pm.coeff(p, i);
unsigned num_vars = pm.size(m);
// add mon * coeff to t;
switch (num_vars) {
case 0:
bound -= coeff;
break;
case 1: {
auto v = nl2lp[pm.get_var(m, 0)];
t.add_monomial(coeff, v);
break;
}
default: {
svector<lp::lpvar> vars;
for (unsigned j = 0; j < num_vars; ++j)
vars.push_back(nl2lp[pm.get_var(m, j)]);
std::sort(vars.begin(), vars.end());
lp::lpvar v;
if (m_vars2mon.contains(vars))
v = m_vars2mon[vars];
else
v = m_nla_core.add_mul_def(vars.size(), vars.data());
t.add_monomial(coeff, v);
break;
}
}
}
}
u_map<lp::lpvar> reverse_lp2nl() {
u_map<lp::lpvar> nl2lp;
for (auto [j, x] : m_lp2nl)
nl2lp.insert(x, j);
return nl2lp;
}
lbool check_assignment() {
setup_assignment_solver();
lbool r = l_undef;
statistics &st = m_nla_core.lp_settings().stats().m_st;
nlsat::literal_vector clause;
nlsat::assignment rvalues(m_nlsat->am());
for (auto [j, x] : m_lp2nl) {
scoped_anum a(am());
am().set(a, m_nla_core.val(j).to_mpq());
rvalues.set(x, a);
}
r = m_nlsat->check(rvalues, clause);
m_nlsat->collect_statistics(st);
switch (r) {
case l_true:
m_nla_core.set_use_nra_model(true);
lra.init_model();
for (lp::constraint_index ci : lra.constraints().indices())
if (!check_constraint(ci))
return l_undef;
for (auto const& m : m_nla_core.emons())
if (!check_monic(m))
return l_undef;
m_nla_core.set_use_nra_model(true);
break;
case l_false:
r = add_lemma(clause);
break;
default:
break;
}
return r;
}
lbool add_lemma(nlsat::literal_vector const &clause) {
u_map<lp::lpvar> nl2lp = reverse_lp2nl();
lbool result = l_false;
{
nla::lemma_builder lemma(m_nla_core, __FUNCTION__);
for (nlsat::literal l : clause) {
if (m_literal2constraint.get((~l).index(), lp::null_ci) != lp::null_ci) {
auto ci = m_literal2constraint[(~l).index()];
lp::explanation ex;
ex.push_back(ci);
lemma &= ex;
continue;
}
nlsat::atom *a = m_nlsat->bool_var2atom(l.var());
if (a->is_root_atom()) {
result = l_undef;
break;
}
SASSERT(a->is_ineq_atom());
auto &ia = *to_ineq_atom(a);
if (ia.size() != 1) {
result = l_undef; // factored polynomials not handled here
break;
}
polynomial::polynomial const *p = ia.p(0);
rational bound(0);
lp::lar_term t;
process_polynomial_check_assignment(p, bound, nl2lp, t);
nla::ineq inq(lp::lconstraint_kind::EQ, t, bound); // initial value overwritten in cases below
switch (a->get_kind()) {
case nlsat::atom::EQ:
inq = nla::ineq(l.sign() ? lp::lconstraint_kind::NE : lp::lconstraint_kind::EQ, t, bound);
break;
case nlsat::atom::LT:
inq = nla::ineq(l.sign() ? lp::lconstraint_kind::GE : lp::lconstraint_kind::LT, t, bound);
break;
case nlsat::atom::GT:
inq = nla::ineq(l.sign() ? lp::lconstraint_kind::LE : lp::lconstraint_kind::GT, t, bound);
break;
default:
UNREACHABLE();
result = l_undef;
break;
}
if (result == l_undef)
break;
if (m_nla_core.ineq_holds(inq)) {
result = l_undef;
break;
}
lemma |= inq;
}
if (result == l_false)
this->m_nla_core.m_check_feasible = true;
} // lemma_builder destructor runs here
if (result == l_undef)
m_nla_core.m_lemmas.pop_back(); // discard incomplete lemma
return result;
}
void add_monic_eq_bound(mon_eq const& m) {
@ -423,20 +635,8 @@ struct solver::imp {
add_ub(lra.get_upper_bound(v), w, lra.get_column_upper_bound_witness(v));
}
lbool r = l_undef;
statistics& st = m_nla_core.lp_settings().stats().m_st;
try {
r = m_nlsat->check();
}
catch (z3_exception&) {
if (m_limit.is_canceled()) {
r = l_undef;
}
else {
m_nlsat->collect_statistics(st);
throw;
}
}
lbool r = m_nlsat->check();
statistics &st = m_nla_core.lp_settings().stats().m_st;
m_nlsat->collect_statistics(st);
switch (r) {
@ -485,18 +685,8 @@ struct solver::imp {
add_ub(lra.get_upper_bound(v), w);
}
lbool r = l_undef;
try {
r = m_nlsat->check();
}
catch (z3_exception&) {
if (m_limit.is_canceled()) {
r = l_undef;
}
else {
throw;
}
}
lbool r = m_nlsat->check();
if (r == l_true)
return r;
@ -643,10 +833,9 @@ struct solver::imp {
unsigned w;
scoped_anum a(am());
for (unsigned v = m_values->size(); v < sz; ++v) {
if (m_nla_core.emons().is_monic_var(v)) {
if (m_nla_core.emons().is_monic_var(v)) {
am().set(a, 1);
auto &m = m_nla_core.emon(v);
for (auto x : m.vars())
am().mul(a, (*m_values)[x], a);
m_values->push_back(a);
@ -654,7 +843,7 @@ struct solver::imp {
else if (lra.column_has_term(v)) {
scoped_anum b(am());
am().set(a, 0);
for (auto const &[w, coeff] : lra.get_term(v)) {
for (auto const &[w, coeff] : lra.get_term(v)) {
am().set(b, coeff.to_mpq());
am().mul(b, (*m_values)[w], b);
am().add(a, b, a);
@ -726,15 +915,67 @@ solver::~solver() {
lbool solver::check() {
return m_imp->check();
try {
return m_imp->check();
}
catch (z3_exception &) {
statistics &st = m_imp->m_nla_core.lp_settings().stats().m_st;
m_imp->m_nlsat->collect_statistics(st);
if (m_imp->m_limit.is_canceled()) {
return l_undef;
}
else {
throw;
}
}
}
lbool solver::check(vector<dd::pdd> const& eqs) {
return m_imp->check(eqs);
try {
return m_imp->check(eqs);
}
catch (z3_exception &) {
statistics &st = m_imp->m_nla_core.lp_settings().stats().m_st;
m_imp->m_nlsat->collect_statistics(st);
if (m_imp->m_limit.is_canceled()) {
return l_undef;
}
else {
throw;
}
}
}
lbool solver::check(dd::solver::equation_vector const& eqs) {
return m_imp->check(eqs);
try {
return m_imp->check(eqs);
}
catch (z3_exception &) {
statistics &st = m_imp->m_nla_core.lp_settings().stats().m_st;
m_imp->m_nlsat->collect_statistics(st);
if (m_imp->m_limit.is_canceled()) {
return l_undef;
}
else {
throw;
}
}
}
lbool solver::check_assignment() {
try {
return m_imp->check_assignment();
}
catch (z3_exception &) {
statistics &st = m_imp->m_nla_core.lp_settings().stats().m_st;
m_imp->m_nlsat->collect_statistics(st);
if (m_imp->m_limit.is_canceled()) {
return l_undef;
}
else {
throw;
}
}
}
bool solver::need_check() {

5
src/math/lp/nra_solver.h
View file

@ -47,6 +47,11 @@ namespace nra {
*/
lbool check(dd::solver::equation_vector const& eqs);
/**
\brief Check feasibility modulo current value assignment.
*/
lbool check_assignment();
/*
\brief determine whether nra check is needed.
*/

4
src/math/lp/static_matrix.cpp
View file

@ -51,8 +51,8 @@ namespace lp {
template void static_matrix<mpq, numeric_pair<mpq> >::set(unsigned int, unsigned int, mpq const&);
template bool static_matrix<mpq, mpq>::pivot_row_to_row_given_cell(unsigned int, column_cell& , unsigned int);
template bool static_matrix<mpq, numeric_pair<mpq> >::pivot_row_to_row_given_cell(unsigned int, column_cell&, unsigned int);
template void static_matrix<mpq, mpq>::pivot_row_to_row_given_cell(unsigned int, column_cell& , unsigned int);
template void static_matrix<mpq, numeric_pair<mpq> >::pivot_row_to_row_given_cell(unsigned int, column_cell&, unsigned int);
template void static_matrix<mpq, numeric_pair<mpq> >::pivot_row_to_row_given_cell_with_sign(unsigned int, column_cell&, unsigned int, int);
template void static_matrix<mpq, mpq>::pivot_row_to_row_given_cell_with_sign(unsigned int, row_cell<empty_struct>&, unsigned int, int);
template void static_matrix<mpq, numeric_pair<mpq> >::add_rows(mpq const&, unsigned int, unsigned int);

24
src/math/lp/static_matrix.h
View file

@ -60,6 +60,14 @@ std::ostream& operator<<(std::ostream& out, const row_strip<T>& r) {
return out << "\n";
}
// Below, static_matrix has a superclass when Z3DEBUG is set, and some
// methods are overrides in that case.
#ifdef Z3DEBUG
#define DEBUG_OVERRIDE override
#else
#define DEBUG_OVERRIDE
#endif
// each assignment for this matrix should be issued only once!!!
template <typename T, typename X>
class static_matrix
@ -119,9 +127,13 @@ public:
void init_empty_matrix(unsigned m, unsigned n);
unsigned row_count() const { return static_cast<unsigned>(m_rows.size()); }
unsigned row_count() const DEBUG_OVERRIDE {
return static_cast<unsigned>(m_rows.size());
}
unsigned column_count() const { return static_cast<unsigned>(m_columns.size()); }
unsigned column_count() const DEBUG_OVERRIDE {
return static_cast<unsigned>(m_columns.size());
}
unsigned lowest_row_in_column(unsigned col);
@ -197,7 +209,7 @@ public:
void cross_out_row_from_column(unsigned col, unsigned k);
T get_elem(unsigned i, unsigned j) const;
T get_elem(unsigned i, unsigned j) const DEBUG_OVERRIDE;
unsigned number_of_non_zeroes_in_column(unsigned j) const { return static_cast<unsigned>(m_columns[j].size()); }
@ -218,8 +230,8 @@ public:
#ifdef Z3DEBUG
unsigned get_number_of_rows() const { return row_count(); }
unsigned get_number_of_columns() const { return column_count(); }
virtual void set_number_of_rows(unsigned /*m*/) { }
virtual void set_number_of_columns(unsigned /*n*/) { }
void set_number_of_rows(unsigned /*m*/) override { }
void set_number_of_columns(unsigned /*n*/) override { }
#endif
T get_balance() const;
@ -293,7 +305,7 @@ public:
// pivot row i to row ii
bool pivot_row_to_row_given_cell(unsigned i, column_cell& c, unsigned j);
void pivot_row_to_row_given_cell(unsigned i, column_cell& c, unsigned j);
void pivot_row_to_row_given_cell_with_sign(unsigned piv_row_index, column_cell& c, unsigned j, int j_sign);
void transpose_rows(unsigned i, unsigned ii) {
auto t = m_rows[i];

19
src/math/lp/static_matrix_def.h
View file

@ -48,7 +48,7 @@ namespace lp {
}
template <typename T, typename X> bool static_matrix<T, X>::pivot_row_to_row_given_cell(unsigned i,
template <typename T, typename X> void static_matrix<T, X>::pivot_row_to_row_given_cell(unsigned i,
column_cell & c, unsigned pivot_col) {
unsigned ii = c.var();
SASSERT(i < row_count() && ii < column_count() && i != ii);
@ -82,7 +82,7 @@ namespace lp {
if (is_zero(rowii[k].coeff()))
remove_element(rowii, rowii[k]);
}
return !rowii.empty();
SASSERT(!rowii.empty());
}
@ -462,12 +462,23 @@ namespace lp {
column_cell& cs = m_columns[row_el_iv.var()][column_offset];
unsigned row_offset = cs.offset();
if (column_offset != column_vals.size() - 1) {
auto & cc = column_vals[column_offset] = column_vals.back(); // copy from the tail
auto & cc = column_vals[column_offset] = column_vals.back(); // column cells are tiny, a plain copy is optimal
m_rows[cc.var()][cc.offset()].offset() = column_offset;
}
if (row_offset != row_vals.size() - 1) {
auto & rc = row_vals[row_offset] = row_vals.back(); // copy from the tail
row_cell<T> & rc = row_vals[row_offset];
row_cell<T> & tail = row_vals.back();
rc.var() = tail.var();
rc.offset() = tail.offset();
// Relocating the tail coefficient: a copy allocates a fresh bignum only when the
// source (tail) is big, so swap to steal the tail's storage exactly in that case.
// When the tail is small the copy never allocates and is cheaper than a swap. See
// Z3Prover/bench#3143.
if (tail.coeff().is_big())
rc.coeff().swap(tail.coeff());
else
rc.coeff() = tail.coeff();
m_columns[rc.var()][rc.offset()].offset() = row_offset;
}

4
src/math/polynomial/README
View file

@ -1,3 +1,5 @@
Polynomial manipulation package.
It contains support for univariate (upolynomial.*) and multivariate polynomials (polynomial.*).
Multivariate polynomial factorization does not work yet (polynomial_factorization.*), and it is disabled.
Multivariate polynomial factorization uses evaluation and bivariate Hensel lifting: evaluate away
extra variables, factor the univariate specialization, then lift to bivariate factors in Zp[x]
and verify over Z. For >2 variables, trial division checks if bivariate factors divide the original.

87
src/math/polynomial/algebraic_numbers.cpp
View file

@ -22,6 +22,7 @@ Notes:
#include "util/mpbqi.h"
#include "util/timeit.h"
#include "util/common_msgs.h"
#include "util/index_sort_with_mutations.h"
#include "math/polynomial/algebraic_numbers.h"
#include "math/polynomial/upolynomial.h"
#include "math/polynomial/sexpr2upolynomial.h"
@ -593,10 +594,57 @@ namespace algebraic_numbers {
}
}
// Sort an index permutation with a bounds-safe, mutation-aware merge
// sort. The comparator (compare/lt) is NOT pure: it MUTATES the
// algebraic numbers it compares (refining their isolating intervals) and
// may throw on the resource limit, so std::sort would be undefined
// behavior here. See util/index_sort_with_mutations.h for the rationale.
void merge_sort_roots_perm(numeral_vector & r, unsigned_vector & perm) {
unsigned n = perm.size();
if (n < 2)
return;
unsigned_vector scratch;
scratch.resize(n, 0);
// Strict, total, stable index comparator: decided sign first, then index
// tiebreak (covers the equal/limit case so the order stays deterministic).
auto idx_lt = [&](unsigned x, unsigned y) {
::sign s = compare(r[x], r[y]);
return s != sign_zero ? s == sign_neg : x < y;
};
stable_index_merge_sort(perm.data(), scratch.data(), n, idx_lt);
}
void sort_roots(numeral_vector & r) {
if (m_limit.inc()) {
// DEBUG_CODE(check_transitivity(r););
std::sort(r.begin(), r.end(), lt_proc(m_wrapper));
if (!m_limit.inc())
return;
// DEBUG_CODE(check_transitivity(r););
unsigned n = r.size();
if (n < 2)
return;
unsigned_vector perm;
perm.resize(n, 0);
for (unsigned i = 0; i < n; ++i)
perm[i] = i;
merge_sort_roots_perm(r, perm);
// Apply the permutation in place via swap cycles. anum swap is a cheap
// pointer swap (move nulls the source), so this is O(n) cheap moves.
unsigned_vector pos; // pos[v] = current position of element v
pos.resize(n, 0);
unsigned_vector at; // at[p] = element currently at position p
at.resize(n, 0);
for (unsigned i = 0; i < n; ++i) {
pos[i] = i;
at[i] = i;
}
for (unsigned target = 0; target < n; ++target) {
unsigned want = perm[target]; // element that should end up at target
unsigned cur = pos[want]; // where it currently is
if (cur == target)
continue;
unsigned other = at[target]; // element currently at target
std::swap(r[target], r[cur]);
at[target] = want; at[cur] = other;
pos[want] = target; pos[other] = cur;
}
}
@ -2018,8 +2066,15 @@ namespace algebraic_numbers {
scoped_mpbq la(bqm()), ua(bqm());
scoped_mpbq lb(bqm()), ub(bqm());
unsigned precision = 10;
if (get_interval(a, la, ua, precision) &&
get_interval(b, lb, ub, precision)) {
// Important: both intervals must be computed. Do not short-circuit with &&:
// the refined bounds la, ua, lb, ub are all used below (and beyond this
// if statement, in the interval-separation checks that compare against
// the bounds of a and b), so get_interval(b, ...) has to run even when
// get_interval(a, ...) returns false (which happens when a is rational
// and its exact root is found).
bool a_separated = get_interval(a, la, ua, precision);
bool b_separated = get_interval(b, lb, ub, precision);
if (a_separated && b_separated) {
IF_VERBOSE(9, verbose_stream() << "sturm 0\n");
if (la > ub)
return sign_pos;
@ -2028,6 +2083,20 @@ namespace algebraic_numbers {
}
IF_VERBOSE(9, verbose_stream() << "sturm 1\n");
// Check whether a can be separated from b's interval and vice versa
// this recognizes the case where the intervals overlap,
// but the anums do not lie in the intersection of the intervals.
scoped_mpq l_a(qm()), u_a(qm()), l_b(qm()), u_b(qm());
to_mpq(qm(), la, l_a);
to_mpq(qm(), ua, u_a);
to_mpq(qm(), lb, l_b);
to_mpq(qm(), ub, u_b);
if (compare(cell_a, l_b) == sign_neg) return sign_neg;
if (compare(cell_a, u_b) == sign_pos) return sign_pos;
if (compare(cell_b, l_a) == sign_neg) return sign_pos;
if (compare(cell_b, u_a) == sign_pos) return sign_neg;
//
// EXPENSIVE CASE
// Let seq be the Sturm-Tarski sequence for
@ -2620,7 +2689,8 @@ namespace algebraic_numbers {
TRACE(isolate_roots, tout << "resultant loop i: " << i << ", y: x" << y << "\np_y: " << p_y << "\n";
tout << "q: " << q << "\n";);
if (ext_pm.is_zero(q)) {
SASSERT(!nested_call);
if (nested_call)
throw algebraic_exception("resultant vanished during nested isolate_roots call");
break;
}
}
@ -2632,7 +2702,8 @@ namespace algebraic_numbers {
// until we find one that is not zero at x2v.
// In the process we will copy p_prime to the local polynomial manager, since we will need to create
// an auxiliary variable.
SASSERT(!nested_call);
if (nested_call)
throw algebraic_exception("resultant vanished during nested isolate_roots call");
unsigned n = ext_pm.degree(p_prime, x);
SASSERT(n > 0);
if (n == 1) {
@ -3447,4 +3518,4 @@ namespace algebraic_numbers {
void manager::collect_statistics(statistics & st) const {
m_imp->collect_statistics(st);
}
};
}

15
src/math/polynomial/algebraic_numbers.h
View file

@ -381,7 +381,7 @@ namespace algebraic_numbers {
class anum {
class anum {
enum anum_kind { BASIC = 0, ROOT };
void* m_cell;
public:
@ -389,6 +389,17 @@ namespace algebraic_numbers {
anum(basic_cell* cell) :m_cell(TAG(void*, cell, BASIC)) { }
anum(algebraic_cell * cell):m_cell(TAG(void*, cell, ROOT)) { }
// Move nulls the source so std::sort's inner shifts stay alias-free
// if the comparator throws between moves (avoids a later double-free).
anum(anum const &) = default;
anum & operator=(anum const &) = default;
anum(anum && other) noexcept : m_cell(other.m_cell) { other.m_cell = nullptr; }
anum & operator=(anum && other) noexcept {
m_cell = other.m_cell;
other.m_cell = nullptr;
return *this;
}
bool is_basic() const { return GET_TAG(m_cell) == BASIC; }
basic_cell * to_basic() const { SASSERT(is_basic()); return UNTAG(basic_cell*, m_cell); }
algebraic_cell * to_algebraic() const { SASSERT(!is_basic()); return UNTAG(algebraic_cell*, m_cell); }
@ -399,7 +410,7 @@ namespace algebraic_numbers {
anum& operator=(basic_cell* cell) { SASSERT(is_null()); m_cell = TAG(void*, cell, BASIC); return *this; }
anum& operator=(algebraic_cell* cell) { SASSERT(is_null()); m_cell = TAG(void*, cell, ROOT); return *this; }
};
};
}
typedef algebraic_numbers::manager anum_manager;
typedef algebraic_numbers::manager::numeral anum;

565
src/math/polynomial/polynomial.cpp
View file

@ -471,8 +471,8 @@ namespace polynomial {
void reset(monomial const * m) {
unsigned id = m->id();
SASSERT(id < m_m2pos.size());
m_m2pos[id] = UINT_MAX;
if (id < m_m2pos.size())
m_m2pos[id] = UINT_MAX;
}
void set(monomial const * m, unsigned pos) {
@ -3923,6 +3923,7 @@ namespace polynomial {
unsigned counter = 0;
while (true) {
checkpoint();
(void)counter;
SASSERT(degree(pp_u, x) >= degree(pp_v, x));
unsigned delta = degree(pp_u, x) - degree(pp_v, x);
@ -6964,16 +6965,572 @@ namespace polynomial {
}
}
// Convert Zp-lifted coefficient arrays to centered Z representatives,
// build multivariate polynomials F1(x,y) and F2(x,y), and verify q == F1*F2.
// Returns true on success. Does NOT deallocate the coefficient vectors.
bool reconstruct_lifted_factors(
polynomial const * q,
var x, var y,
unsigned deg_y,
uint64_t prime,
vector<upolynomial::scoped_numeral_vector*> const & q_coeffs,
vector<upolynomial::scoped_numeral_vector*> const & F1_coeffs,
vector<upolynomial::scoped_numeral_vector*> const & F2_coeffs,
polynomial_ref & F1_out,
polynomial_ref & F2_out) {
scoped_numeral p_num(m_manager);
m_manager.set(p_num, static_cast<int>(prime));
scoped_numeral half_p(m_manager);
m_manager.set(half_p, static_cast<int>(prime));
m_manager.div(half_p, mpz(2), half_p);
polynomial_ref F1_poly(pm()), F2_poly(pm());
F1_poly = mk_zero();
F2_poly = mk_zero();
for (unsigned j = 0; j <= deg_y; j++) {
// Center the coefficients: if coeff > p/2, subtract p
for (unsigned i = 0; i < F1_coeffs[j]->size(); i++)
if (m_manager.gt((*F1_coeffs[j])[i], half_p))
m_manager.sub((*F1_coeffs[j])[i], p_num, (*F1_coeffs[j])[i]);
for (unsigned i = 0; i < F2_coeffs[j]->size(); i++)
if (m_manager.gt((*F2_coeffs[j])[i], half_p))
m_manager.sub((*F2_coeffs[j])[i], p_num, (*F2_coeffs[j])[i]);
// Build y^j * F_coeffs[j](x)
if (!F1_coeffs[j]->empty()) {
polynomial_ref univ(pm());
univ = to_polynomial(F1_coeffs[j]->size(), F1_coeffs[j]->data(), x);
if (!is_zero(univ)) {
monomial_ref yj(pm());
yj = mk_monomial(y, j);
polynomial_ref term(pm());
term = mul(yj, univ);
F1_poly = add(F1_poly, term);
}
}
if (!F2_coeffs[j]->empty()) {
polynomial_ref univ(pm());
univ = to_polynomial(F2_coeffs[j]->size(), F2_coeffs[j]->data(), x);
if (!is_zero(univ)) {
monomial_ref yj(pm());
yj = mk_monomial(y, j);
polynomial_ref term(pm());
term = mul(yj, univ);
F2_poly = add(F2_poly, term);
}
}
}
// Verify: q == F1 * F2 over Z[x, y]
polynomial_ref product(pm());
product = mul(F1_poly, F2_poly);
if (!eq(product, q))
return false;
F1_out = F1_poly;
F2_out = F2_poly;
return true;
}
// Hensel lifting loop: for j = 1..deg_y, compute F1_coeffs[j] and F2_coeffs[j]
// using Bezout coefficients s, t such that s*f1 + t*f2 = 1 in Zp[x].
// F1_coeffs[0] and F2_coeffs[0] must already be initialized.
void hensel_lift_coefficients(
upolynomial::zp_manager & zp_upm,
unsigned deg_y,
upolynomial::scoped_numeral_vector const & f1_p,
upolynomial::scoped_numeral_vector const & f2_p,
upolynomial::scoped_numeral_vector const & s_vec,
upolynomial::scoped_numeral_vector const & t_vec,
numeral const & lc_val,
vector<upolynomial::scoped_numeral_vector*> const & q_coeffs,
vector<upolynomial::scoped_numeral_vector*> & F1_coeffs,
vector<upolynomial::scoped_numeral_vector*> & F2_coeffs) {
auto & nm = upm().m();
auto & znm = zp_upm.m();
scoped_numeral lc_inv(m_manager);
znm.set(lc_inv, lc_val);
znm.inv(lc_inv);
for (unsigned j = 1; j <= deg_y; j++) {
checkpoint();
// Compute e_j = q_coeffs[j] - sum_{a+b=j, a>0, b>0} F1_coeffs[a] * F2_coeffs[b]
upolynomial::scoped_numeral_vector e_j(nm);
if (j < q_coeffs.size() && !q_coeffs[j]->empty())
zp_upm.set(q_coeffs[j]->size(), q_coeffs[j]->data(), e_j);
for (unsigned a = 0; a <= j; a++) {
unsigned b = j - a;
if (b > deg_y) continue;
if (a == j && b == 0) continue;
if (a == 0 && b == j) continue;
if (F1_coeffs[a]->empty() || F2_coeffs[b]->empty()) continue;
upolynomial::scoped_numeral_vector prod(nm);
zp_upm.mul(F1_coeffs[a]->size(), F1_coeffs[a]->data(),
F2_coeffs[b]->size(), F2_coeffs[b]->data(), prod);
upolynomial::scoped_numeral_vector new_e(nm);
zp_upm.sub(e_j.size(), e_j.data(), prod.size(), prod.data(), new_e);
e_j.swap(new_e);
}
if (e_j.empty()) continue;
// Solve using Bezout coefficients: A = (e_j * t) mod f1
upolynomial::scoped_numeral_vector e_t(nm), A_val(nm), Q_tmp(nm);
zp_upm.mul(e_j.size(), e_j.data(), t_vec.size(), t_vec.data(), e_t);
zp_upm.div_rem(e_t.size(), e_t.data(), f1_p.size(), f1_p.data(), Q_tmp, A_val);
// B = (e_j * s) mod f2
upolynomial::scoped_numeral_vector e_s(nm), B_val(nm);
zp_upm.mul(e_j.size(), e_j.data(), s_vec.size(), s_vec.data(), e_s);
zp_upm.div_rem(e_s.size(), e_s.data(), f2_p.size(), f2_p.data(), Q_tmp, B_val);
// F1[j] = A * lc_inv, F2[j] = B
zp_upm.mul(A_val, lc_inv);
zp_upm.set(A_val.size(), A_val.data(), *F1_coeffs[j]);
zp_upm.set(B_val.size(), B_val.data(), *F2_coeffs[j]);
}
}
// Extract coefficients of q w.r.t. y as upolynomials in Zp[x], and initialize
// the lifted factor coefficient arrays with F1[0] = f1, F2[0] = lc_val * f2.
// Returns false if q is not truly bivariate in x and y.
bool extract_and_init_lift_coefficients(
upolynomial::zp_manager & zp_upm,
polynomial const * q,
var y,
unsigned deg_y,
upolynomial::scoped_numeral_vector const & f1_p,
upolynomial::scoped_numeral_vector const & f2_p,
numeral const & lc_val,
vector<upolynomial::scoped_numeral_vector*> & q_coeffs,
vector<upolynomial::scoped_numeral_vector*> & F1_coeffs,
vector<upolynomial::scoped_numeral_vector*> & F2_coeffs) {
auto & nm = upm().m();
auto & znm = zp_upm.m();
for (unsigned j = 0; j <= deg_y; j++) {
polynomial_ref cj(pm());
cj = coeff(q, y, j);
auto * vec = alloc(upolynomial::scoped_numeral_vector, nm);
if (!is_zero(cj) && is_univariate(cj))
upm().to_numeral_vector(cj, *vec);
else if (!is_zero(cj) && is_const(cj)) {
vec->push_back(numeral());
nm.set(vec->back(), cj->a(0));
}
else if (!is_zero(cj)) {
dealloc(vec);
return false;
}
for (unsigned i = 0; i < vec->size(); i++)
znm.p_normalize((*vec)[i]);
zp_upm.trim(*vec);
q_coeffs.push_back(vec);
}
// Initialize lifted factor coefficient arrays
for (unsigned j = 0; j <= deg_y; j++) {
F1_coeffs.push_back(alloc(upolynomial::scoped_numeral_vector, nm));
F2_coeffs.push_back(alloc(upolynomial::scoped_numeral_vector, nm));
}
// F1[0] = f1, F2[0] = lc_val * f2
zp_upm.set(f1_p.size(), f1_p.data(), *F1_coeffs[0]);
scoped_numeral lc_p(m_manager);
znm.set(lc_p, lc_val);
upolynomial::scoped_numeral_vector lc_f2(nm);
zp_upm.set(f2_p.size(), f2_p.data(), lc_f2);
zp_upm.mul(lc_f2, lc_p);
zp_upm.set(lc_f2.size(), lc_f2.data(), *F2_coeffs[0]);
return true;
}
// Bivariate Hensel lifting for multivariate factorization.
//
// Mathematical setup:
// We have q(x, y) in Z[x, y] with degree deg_y in y.
// Evaluating at y = 0 gives a univariate factorization
// q(x, 0) = lc_val * uf1_monic(x) * uf2_monic(x)
// where uf1_monic and uf2_monic are monic, coprime polynomials in Z[x],
// and lc_val = lc(q, x)(y=0) is an integer.
//
// Goal: lift to q(x, y) = F1(x, y) * F2(x, y) over Z[x, y].
//
// Method (linear Hensel lifting in Zp[x]):
// 1. Reduce uf1_monic, uf2_monic to f1, f2 in Zp[x] and compute
// Bezout coefficients s, t with s*f1 + t*f2 = 1 in Zp[x].
// This requires gcd(f1, f2) = 1 in Zp[x], i.e. the prime p
// must not divide the resultant of f1, f2.
// 2. Expand q, F1, F2 as polynomials in y with coefficients in Zp[x]:
// q = q_0(x) + q_1(x)*y + ... + q_{deg_y}(x)*y^{deg_y}
// F1 = F1_0(x) + F1_1(x)*y + ...
// F2 = F2_0(x) + F2_1(x)*y + ...
// The y^0 coefficients are known: F1_0 = f1, F2_0 = lc_val * f2.
// 3. For j = 1, ..., deg_y, matching the y^j coefficient of q = F1 * F2:
// q_j = sum_{a+b=j} F1_a * F2_b
// The unknowns are F1_j and F2_j. Set
// e_j = q_j - sum_{a+b=j, 0<a, 0<b} F1_a * F2_b
// so that F1_j * F2_0 + F2_j * F1_0 = e_j,
// i.e. F1_j * (lc_val * f2) + F2_j * f1 = e_j.
// From step 1 we have the Bezout identity s*f1 + t*f2 = 1 in Zp[x].
// Multiplying by e_j: (e_j*s)*f1 + (e_j*t)*f2 = e_j.
// Reducing to match degree constraints gives:
// F2_j = (e_j * s) mod f2
// F1_j = ((e_j * t) mod f1) / lc_val
// all in Zp[x].
// 4. Map Zp coefficients to centered representatives in (-p/2, p/2]
// and verify q = F1 * F2 over Z[x, y].
//
// Preconditions:
// - q is a bivariate polynomial in x and y, i.e. every coeff(q, y, j)
// is either zero, a constant, or univariate in x.
// - uf1_monic, uf2_monic are monic univariate polynomials in Z[x],
// stored as coefficient vectors, such that
// q(x, 0) = lc_val * uf1_monic(x) * uf2_monic(x).
// - gcd(uf1_monic, uf2_monic) = 1 over Q[x].
// - prime does not divide lc_val.
//
// Returns true if the lift succeeds, i.e. q = F1 * F2 holds over Z.
// Returns false if coprimality fails in Zp, if q is not bivariate,
// or if p is too small for the centered representative trick to work.
bool try_bivar_hensel_lift(
polynomial const * q,
var x, var y,
unsigned deg_y,
upolynomial::numeral_vector const & uf1_monic, // monic factor of q(x,0)/lc_val in Z[x], as coefficient vector
upolynomial::numeral_vector const & uf2_monic, // monic factor of q(x,0)/lc_val in Z[x], coprime to uf1_monic
numeral const & lc_val, // lc(q, x) evaluated at y=0, so q(x,0) = lc_val * uf1_monic * uf2_monic
uint64_t prime,
polynomial_ref & F1_out,
polynomial_ref & F2_out) {
typedef upolynomial::zp_manager zp_mgr;
typedef upolynomial::zp_numeral_manager zp_nm;
auto & nm = upm().m();
zp_mgr zp_upm(m_limit, nm.m());
scoped_numeral p_num(m_manager);
m_manager.set(p_num, static_cast<int>(prime));
zp_upm.set_zp(p_num);
zp_nm & znm = zp_upm.m();
// Convert h1, h2 to Zp
upolynomial::scoped_numeral_vector f1_p(nm), f2_p(nm);
for (unsigned i = 0; i < uf1_monic.size(); i++) {
f1_p.push_back(numeral());
znm.set(f1_p.back(), uf1_monic[i]);
}
zp_upm.trim(f1_p);
for (unsigned i = 0; i < uf2_monic.size(); i++) {
f2_p.push_back(numeral());
znm.set(f2_p.back(), uf2_monic[i]);
}
zp_upm.trim(f2_p);
// Make monic in Zp
zp_upm.mk_monic(f1_p.size(), f1_p.data());
zp_upm.mk_monic(f2_p.size(), f2_p.data());
// Extended GCD in Zp[x]: s*f1 + t*f2 = 1
upolynomial::scoped_numeral_vector s_vec(nm), t_vec(nm), d_vec(nm);
zp_upm.ext_gcd(f1_p, f2_p, s_vec, t_vec, d_vec);
// Check gcd = 1
if (d_vec.size() != 1 || !znm.is_one(d_vec[0]))
return false;
vector<upolynomial::scoped_numeral_vector*> q_coeffs, F1_coeffs, F2_coeffs;
if (!extract_and_init_lift_coefficients(zp_upm, q, y, deg_y,
f1_p, f2_p, lc_val,
q_coeffs, F1_coeffs, F2_coeffs)) {
for (auto * v : q_coeffs) dealloc(v);
return false;
}
hensel_lift_coefficients(zp_upm, deg_y, f1_p, f2_p,
s_vec, t_vec, lc_val,
q_coeffs, F1_coeffs, F2_coeffs);
bool ok = reconstruct_lifted_factors(q, x, y, deg_y, prime,
q_coeffs, F1_coeffs, F2_coeffs,
F1_out, F2_out);
for (auto * v : q_coeffs) dealloc(v);
for (auto * v : F1_coeffs) dealloc(v);
for (auto * v : F2_coeffs) dealloc(v);
return ok;
}
// Evaluation points used for multivariate factorization
static constexpr int s_factor_eval_values[] = { 0, 1, -1, 2, -2, 3, -3 };
static constexpr unsigned s_n_factor_eval_values = sizeof(s_factor_eval_values) / sizeof(s_factor_eval_values[0]);
// Primes for Hensel lifting, tried in increasing order.
// Lifting succeeds when the prime exceeds twice the largest coefficient in any factor.
// A Mignotte-style bound could automate this, but for now we try a fixed list.
static constexpr uint64_t s_factor_primes[] = { 39103, 104729, 1000003, 100000007 };
void factor_n_sqf_pp(polynomial const * p, factors & r, var x, unsigned k) {
SASSERT(degree(p, x) > 2);
SASSERT(is_primitive(p, x));
SASSERT(is_square_free(p, x));
TRACE(factor, tout << "factor square free (degree > 2):\n"; p->display(tout, m_manager); tout << "\n";);
// TODO: invoke Dejan's procedure
if (is_univariate(p)) {
factor_sqf_pp_univ(p, r, k, factor_params());
return;
}
// Try multivariate factorization. If checkpoint() throws during the
// attempt, the shared som_buffer/m_m2pos may be left dirty. Catch the
// exception, reset the buffers, return unfactored, then rethrow so
// cancellation propagates normally.
try {
if (try_multivar_factor(p, r, x, k))
return;
}
catch (...) {
m_som_buffer.reset();
m_som_buffer2.reset();
m_cheap_som_buffer.reset();
m_cheap_som_buffer2.reset();
throw;
}
// Could not factor, return p as-is
r.push_back(const_cast<polynomial*>(p), k);
}
// Returns true if factorization succeeded and factors were added to r.
bool try_multivar_factor(polynomial const * p, factors & r, var x, unsigned k) {
var_vector all_vars;
m_wrapper.vars(p, all_vars);
// Try the main variable x first (caller chose it), then others by degree
svector<var> main_vars;
main_vars.push_back(x);
svector<std::pair<unsigned, var>> var_by_deg;
for (var v : all_vars)
if (v != x)
var_by_deg.push_back(std::make_pair(degree(p, v), v));
std::sort(var_by_deg.begin(), var_by_deg.end(),
[](auto const& a, auto const& b) { return a.first > b.first; });
for (auto const& [d, v] : var_by_deg)
if (d > 1) main_vars.push_back(v);
for (var main_var : main_vars) {
unsigned deg_main = degree(p, main_var);
if (deg_main <= 1) continue;
for (var lift_var : all_vars) {
if (lift_var == main_var) continue;
checkpoint();
// Variables to evaluate away
var_vector eval_vars;
for (var v : all_vars)
if (v != main_var && v != lift_var)
eval_vars.push_back(v);
unsigned n_eval = eval_vars.size();
// Try a small number of evaluation point combos for extra variables
unsigned max_combos = (n_eval == 0) ? 1 : std::min(s_n_factor_eval_values, 5u);
for (unsigned combo = 0; combo < max_combos; combo++) {
checkpoint();
// Reduce to bivariate
polynomial_ref p_bivar(pm());
if (n_eval > 0) {
scoped_numeral_vector eval_vals(m_manager);
for (unsigned i = 0; i < n_eval; i++) {
eval_vals.push_back(numeral());
unsigned c = combo;
for (unsigned skip = 0; skip < i; skip++)
c /= s_n_factor_eval_values;
m_manager.set(eval_vals.back(), s_factor_eval_values[c % s_n_factor_eval_values]);
}
p_bivar = substitute(p, n_eval, eval_vars.data(), eval_vals.data());
}
else
p_bivar = const_cast<polynomial*>(p);
if (degree(p_bivar, main_var) != deg_main) continue;
unsigned deg_lift = degree(p_bivar, lift_var);
if (deg_lift == 0) continue;
// Find a good evaluation point a for the lift variable
for (int a : s_factor_eval_values) {
scoped_numeral val_scoped(m_manager);
m_manager.set(val_scoped, a);
numeral const & val_ref = val_scoped;
polynomial_ref p_univ(pm());
p_univ = substitute(p_bivar, 1, &lift_var, &val_ref);
if (!is_univariate(p_univ)) continue;
if (degree(p_univ, main_var) != deg_main) continue;
if (!is_square_free(p_univ, main_var)) continue;
// Factor the univariate polynomial
up_manager::scoped_numeral_vector p_univ_vec(upm().m());
polynomial_ref p_univ_ref(pm());
p_univ_ref = p_univ;
upm().to_numeral_vector(p_univ_ref, p_univ_vec);
// Make primitive before factoring
upm().get_primitive(p_univ_vec, p_univ_vec);
up_manager::factors univ_fs(upm());
upolynomial::factor_square_free(upm(), p_univ_vec, univ_fs);
unsigned nf = univ_fs.distinct_factors();
if (nf <= 1) continue;
// Translate so evaluation is at lift_var = 0
polynomial_ref q(pm());
scoped_numeral a_val(m_manager);
m_manager.set(a_val, a);
q = translate(p_bivar, lift_var, a_val);
// Get leading coefficient at evaluation point
polynomial_ref lc_poly(pm());
lc_poly = coeff(q, main_var, deg_main);
scoped_numeral lc_at_0(m_manager);
if (is_const(lc_poly))
m_manager.set(lc_at_0, lc_poly->a(0));
else {
scoped_numeral zero_val(m_manager);
m_manager.set(zero_val, 0);
numeral const & zero_ref = zero_val;
polynomial_ref lc_eval(pm());
lc_eval = substitute(lc_poly, 1, &lift_var, &zero_ref);
if (is_const(lc_eval))
m_manager.set(lc_at_0, lc_eval->a(0));
else
continue;
}
// Try splits with increasing primes.
// Only contiguous splits {0..split-1} vs {split..nf-1} are tried,
// not all subset partitions. This avoids exponential search but may
// miss some factorizations. Recursive calls on the lifted factors
// partially compensate by further splitting successful lifts.
for (uint64_t prime : s_factor_primes) {
scoped_numeral prime_num(m_manager);
m_manager.set(prime_num, static_cast<int64_t>(prime));
scoped_numeral gcd_val(m_manager);
m_manager.gcd(prime_num, lc_at_0, gcd_val);
if (!m_manager.is_one(gcd_val)) continue;
for (unsigned split = 1; split <= nf / 2; split++) {
checkpoint();
upolynomial::scoped_numeral_vector h1(upm().m()), h2(upm().m());
upm().set(univ_fs[0].size(), univ_fs[0].data(), h1);
for (unsigned i = 1; i < split; i++) {
upolynomial::scoped_numeral_vector temp(upm().m());
upm().mul(h1.size(), h1.data(), univ_fs[i].size(), univ_fs[i].data(), temp);
h1.swap(temp);
}
upm().set(univ_fs[split].size(), univ_fs[split].data(), h2);
for (unsigned i = split + 1; i < nf; i++) {
upolynomial::scoped_numeral_vector temp(upm().m());
upm().mul(h2.size(), h2.data(), univ_fs[i].size(), univ_fs[i].data(), temp);
h2.swap(temp);
}
auto & nm_ref = upm().m();
if (!h1.empty() && nm_ref.is_neg(h1.back())) {
for (unsigned i = 0; i < h1.size(); i++)
nm_ref.neg(h1[i]);
}
if (!h2.empty() && nm_ref.is_neg(h2.back())) {
for (unsigned i = 0; i < h2.size(); i++)
nm_ref.neg(h2[i]);
}
polynomial_ref F1(pm()), F2(pm());
if (!try_bivar_hensel_lift(q, main_var, lift_var, deg_lift, h1, h2, lc_at_0, prime, F1, F2))
continue;
// Translate back
scoped_numeral neg_a(m_manager);
m_manager.set(neg_a, -a);
F1 = translate(F1, lift_var, neg_a);
F2 = translate(F2, lift_var, neg_a);
if (n_eval == 0) {
// p is bivariate, factors verified by try_bivar_hensel_lift.
// Use specialized handlers for small degrees to avoid
// re-entering factor_sqf_pp which corrupts shared buffers.
polynomial_ref bivar_fs[] = { F1, F2 };
for (polynomial_ref & bf : bivar_fs) {
if (is_const(bf) || degree(bf, x) == 0) continue;
unsigned d = degree(bf, x);
if (d == 1)
factor_1_sqf_pp(bf, r, x, k);
else if (d == 2 && is_primitive(bf, x) && is_square_free(bf, x))
factor_2_sqf_pp(bf, r, x, k);
else
r.push_back(bf, k);
}
return true;
}
// Multivariate: check if bivariate factors divide original p.
// We use exact_pseudo_division, which computes Q, R with
// lc(cand)^d * p = Q * cand + R. If R=0 and cand*Q == p
// then cand is a true factor. The eq() check is needed
// because pseudo-division may introduce an lc power that
// prevents Q from being the exact quotient.
polynomial_ref cands[] = { F1, F2 };
for (polynomial_ref & cand : cands) {
if (is_const(cand)) continue;
polynomial_ref Q_div(pm()), R_div(pm());
var div_var = max_var(cand);
exact_pseudo_division(const_cast<polynomial*>(p), cand, div_var, Q_div, R_div);
if (!is_zero(R_div)) continue;
polynomial_ref check(pm());
check = mul(cand, Q_div);
if (eq(check, p)) {
// Push factors directly, using specialized handlers
// for small degrees only.
polynomial_ref parts[] = { cand, Q_div };
for (polynomial_ref & part : parts) {
if (is_const(part)) {
acc_constant(r, part->a(0));
continue;
}
if (degree(part, x) == 0) continue;
unsigned d = degree(part, x);
if (d == 1)
factor_1_sqf_pp(part, r, x, k);
else if (d == 2 && is_primitive(part, x) && is_square_free(part, x))
factor_2_sqf_pp(part, r, x, k);
else
r.push_back(part, k);
}
return true;
}
}
}
}
}
}
}
}
return false;
}
void factor_sqf_pp(polynomial const * p, factors & r, var x, unsigned k, factor_params const & params) {
SASSERT(degree(p, x) > 0);
SASSERT(is_primitive(p, x));
@ -7697,7 +8254,7 @@ namespace polynomial {
p->display_smt2(out, m_imp->m_manager, proc);
return out;
}
};
}
polynomial::polynomial * convert(polynomial::manager & sm, polynomial::polynomial * p, polynomial::manager & tm,
polynomial::var x, unsigned max_d) {

4
src/math/polynomial/polynomial.h
View file

@ -36,7 +36,7 @@ class small_object_allocator;
namespace algebraic_numbers {
class anum;
class manager;
};
}
namespace polynomial {
typedef unsigned var;
@ -1065,7 +1065,7 @@ namespace polynomial {
scoped_set_zp(manager & _m, uint64_t p):m(_m), m_modular(m.modular()), m_p(m.m()) { m_p = m.p(); m.set_zp(p); }
~scoped_set_zp() { if (m_modular) m.set_zp(m_p); else m.set_z(); }
};
};
}
typedef polynomial::polynomial_ref polynomial_ref;
typedef polynomial::polynomial_ref_vector polynomial_ref_vector;

2
src/math/polynomial/polynomial_cache.cpp
View file

@ -256,4 +256,4 @@ namespace polynomial {
dealloc(m_imp);
m_imp = alloc(imp, _m);
}
};
}

2
src/math/polynomial/polynomial_cache.h
View file

@ -40,4 +40,4 @@ namespace polynomial {
void factor(polynomial const * p, polynomial_ref_vector & distinct_factors);
void reset();
};
};
}

2
src/math/polynomial/polynomial_primes.h
View file

@ -66,5 +66,5 @@ namespace polynomial {
};
#endif
};
}

2
src/math/polynomial/polynomial_var2value.h
View file

@ -43,6 +43,6 @@ namespace polynomial {
}
};
};
}

2
src/math/polynomial/rpolynomial.cpp
View file

@ -789,4 +789,4 @@ namespace rpolynomial {
}
#endif
};
}

2
src/math/polynomial/rpolynomial.h
View file

@ -175,7 +175,7 @@ namespace rpolynomial {
return out;
}
};
};
}
typedef rpolynomial::polynomial_ref rpolynomial_ref;
typedef rpolynomial::polynomial_ref_vector rpolynomial_ref_vector;

6
src/math/polynomial/upolynomial.cpp
View file

@ -2616,6 +2616,7 @@ namespace upolynomial {
\warning This method may loop if p is not square free or if (a,b) is not an isolating interval.
*/
bool manager::isolating2refinable(unsigned sz, numeral const * p, mpbq_manager & bqm, mpbq & a, mpbq & b) {
checkpoint();
int sign_a = eval_sign_at(sz, p, a);
int sign_b = eval_sign_at(sz, p, b);
TRACE(upolynomial, tout << "sign_a: " << sign_a << ", sign_b: " << sign_b << "\n";);
@ -2631,6 +2632,7 @@ namespace upolynomial {
bqm.add(a, b, new_a);
bqm.div2(new_a);
while (true) {
checkpoint();
TRACE(upolynomial, tout << "CASE 2, a: " << bqm.to_string(a) << ", b: " << bqm.to_string(b) << ", new_a: " << bqm.to_string(new_a) << "\n";);
int sign_new_a = eval_sign_at(sz, p, new_a);
if (sign_new_a != sign_b) {
@ -2656,6 +2658,7 @@ namespace upolynomial {
bqm.add(a, b, new_b);
bqm.div2(new_b);
while (true) {
checkpoint();
TRACE(upolynomial, tout << "CASE 3, a: " << bqm.to_string(a) << ", b: " << bqm.to_string(b) << ", new_b: " << bqm.to_string(new_b) << "\n";);
int sign_new_b = eval_sign_at(sz, p, new_b);
if (sign_new_b != sign_a) {
@ -2709,6 +2712,7 @@ namespace upolynomial {
bqm.div2(new_b2);
while (true) {
checkpoint();
TRACE(upolynomial,
tout << "CASE 4\na1: " << bqm.to_string(a1) << ", b1: " << bqm.to_string(b1) << ", new_a1: " << bqm.to_string(new_a1) << "\n";
tout << "a2: " << bqm.to_string(a2) << ", b2: " << bqm.to_string(b2) << ", new_b2: " << bqm.to_string(new_b2) << "\n";);
@ -3138,4 +3142,4 @@ namespace upolynomial {
}
return out;
}
};
}

2
src/math/polynomial/upolynomial.h
View file

@ -917,4 +917,4 @@ namespace upolynomial {
std::ostream& display(std::ostream & out, upolynomial_sequence const & seq, char const * var_name = "x") const;
};
};
}

2
src/math/polynomial/upolynomial_factorization.cpp
View file

@ -1299,4 +1299,4 @@ bool factor_square_free(z_manager & upm, numeral_vector const & f, factors & fs,
return factor_square_free(upm, f, fs, 1, params);
}
}; // end upolynomial namespace
} // end upolynomial namespace

2
src/math/polynomial/upolynomial_factorization.h
View file

@ -90,5 +90,5 @@ namespace upolynomial {
That is, the factors of f are inserted as factors of degree k into fs.
*/
bool factor_square_free(z_manager & upm, numeral_vector const & f, factors & fs, unsigned k, factor_params const & ps = factor_params());
};
}

2
src/math/polynomial/upolynomial_factorization_int.h
View file

@ -416,5 +416,5 @@ namespace upolynomial {
}
}
};
};
}

37
src/math/realclosure/realclosure.cpp
View file

@ -3448,16 +3448,23 @@ namespace realclosure {
return true;
}
unsigned get_sign_condition_size(numeral const &a, unsigned i) {
algebraic * ext = to_algebraic(to_rational_function(a)->ext());
sign_condition* get_ith_sign_condition(algebraic* ext, unsigned i) {
const sign_det * sdt = ext->sdt();
if (!sdt)
return 0;
return nullptr;
sign_condition * sc = sdt->sc(ext->sc_idx());
while (i) {
if (sc) sc = sc->prev();
while (i && sc) {
sc = sc->prev();
i--;
}
return sc;
}
unsigned get_sign_condition_size(numeral const &a, unsigned i) {
algebraic * ext = to_algebraic(to_rational_function(a)->ext());
sign_condition * sc = get_ith_sign_condition(ext, i);
if (!sc)
return 0;
return ext->sdt()->qs()[sc->qidx()].size();
}
@ -3466,14 +3473,9 @@ namespace realclosure {
if (!is_algebraic(a))
return 0;
algebraic * ext = to_algebraic(to_rational_function(a)->ext());
const sign_det * sdt = ext->sdt();
if (!sdt)
sign_condition * sc = get_ith_sign_condition(ext, i);
if (!sc)
return 0;
sign_condition * sc = sdt->sc(ext->sc_idx());
while (i) {
if (sc) sc = sc->prev();
i--;
}
const polynomial & q = ext->sdt()->qs()[sc->qidx()];
return q.size();
}
@ -3483,14 +3485,9 @@ namespace realclosure {
if (!is_algebraic(a))
return numeral();
algebraic * ext = to_algebraic(to_rational_function(a)->ext());
const sign_det * sdt = ext->sdt();
if (!sdt)
sign_condition * sc = get_ith_sign_condition(ext, i);
if (!sc)
return numeral();
sign_condition * sc = sdt->sc(ext->sc_idx());
while (i) {
if (sc) sc = sc->prev();
i--;
}
const polynomial & q = ext->sdt()->qs()[sc->qidx()];
if (j >= q.size())
return numeral();
@ -6489,7 +6486,7 @@ namespace realclosure {
{
return m_imp->get_sign_condition_coefficient(a, i, j);
}
};
}
void pp(realclosure::manager::imp * imp, realclosure::polynomial const & p, realclosure::extension * ext) {
imp->display_polynomial_expr(std::cout, p, ext, false, false);

2
src/math/realclosure/realclosure.h
View file

@ -317,7 +317,7 @@ namespace realclosure {
void * data() { return m_value; }
static num mk(void * ptr) { num r; r.m_value = reinterpret_cast<value*>(ptr); return r; }
};
};
}
typedef realclosure::manager rcmanager;
typedef rcmanager::numeral rcnumeral;

2
src/math/simplex/bit_matrix.cpp
View file

@ -125,7 +125,7 @@ unsigned_vector bit_matrix::gray(unsigned n) {
auto v = gray(n-1);
auto w = v;
w.reverse();
for (auto & u : v) u |= (1 << (n-1));
for (auto & u : v) u |= (1u << (n-1));
v.append(w);
return v;
}

2
src/math/simplex/model_based_opt.cpp
View file

@ -89,7 +89,7 @@ namespace opt {
}
void model_based_opt::def::dec_ref() {
SASSERT(m_ref_count > 0);
++m_ref_count;
--m_ref_count;
if (m_ref_count == 0)
dealloc(this);
}

29
src/math/simplex/model_based_opt.h
View file

@ -85,6 +85,7 @@ namespace opt {
enum def_t { add_t, mul_t, div_t, const_t, var_t};
struct def {
def() = default;
virtual ~def() = default;
static def* from_row(row const& r, unsigned x);
def_t m_type;
unsigned m_ref_count = 0;
@ -116,9 +117,15 @@ namespace opt {
class def_ref {
def* m_def = nullptr;
public:
def_ref(def* d) {
if (d) d->inc_ref();
m_def = d;
def_ref() = default;
def_ref(def* d) : m_def(d) {
if (m_def) m_def->inc_ref();
}
def_ref(def_ref const& other) : m_def(other.m_def) {
if (m_def) m_def->inc_ref();
}
def_ref(def_ref&& other) noexcept : m_def(other.m_def) {
other.m_def = nullptr;
}
def_ref& operator=(def* d) {
if (d) d->inc_ref();
@ -136,25 +143,37 @@ namespace opt {
return *this;
}
def_ref& operator=(def_ref&& d) noexcept {
if (&d == this)
return *this;
if (m_def) m_def->dec_ref();
m_def = d.m_def;
d.m_def = nullptr;
return *this;
}
def& operator*() { return *m_def; }
def* operator->() { return m_def; }
def const& operator*() const { return *m_def; }
operator bool() const { return !!m_def; }
operator bool() const { return m_def != nullptr; }
~def_ref() { if (m_def) m_def->dec_ref(); };
~def_ref() { if (m_def) m_def->dec_ref(); }
};
struct add_def : public def {
def* x, *y;
add_def(def* x, def* y) : x(x), y(y) { m_type = add_t; x->inc_ref(); y->inc_ref(); }
~add_def() override { x->dec_ref(); y->dec_ref(); }
};
struct mul_def : public def {
def* x, *y;
mul_def(def* x, def* y) : x(x), y(y) { m_type = mul_t; x->inc_ref(); y->inc_ref(); }
~mul_def() override { x->dec_ref(); y->dec_ref(); }
};
struct div_def : public def {
def* x;
rational m_div{ 1 };
div_def(def* x, rational const& d) : x(x), m_div(d) { m_type = div_t; x->inc_ref(); }
~div_def() override { x->dec_ref(); }
};
struct var_def : public def {
var v;

2
src/math/simplex/simplex.cpp
View file

@ -74,4 +74,4 @@ namespace simplex {
}
}
}
};
}

2
src/math/simplex/simplex.h
View file

@ -204,5 +204,5 @@ namespace simplex {
void kernel(sparse_matrix<mpq_ext>& s, vector<vector<rational>>& K);
void kernel_ffe(sparse_matrix<mpq_ext> &s, vector<vector<rational>> &K);
};
}

2
src/math/simplex/simplex_def.h
View file

@ -1035,6 +1035,6 @@ namespace simplex {
}
};
}

2
src/math/simplex/sparse_matrix.h
View file

@ -355,4 +355,4 @@ namespace simplex {
typedef unsynch_mpq_inf_manager eps_manager;
};
};
}

2
src/math/simplex/sparse_matrix_def.h
View file

@ -603,5 +603,5 @@ namespace simplex {
};
}

2
src/math/subpaving/subpaving.cpp
View file

@ -271,4 +271,4 @@ namespace subpaving {
return alloc(context_mpfx_wrapper, lim, m, qm, p, a);
}
};
}

2
src/math/subpaving/subpaving.h
View file

@ -116,6 +116,6 @@ context * mk_hwf_context(reslimit& lim, f2n<hwf_manager> & m, unsynch_mpq_manage
context * mk_mpff_context(reslimit& lim, mpff_manager & m, unsynch_mpq_manager & qm, params_ref const & p = params_ref(), small_object_allocator * a = nullptr);
context * mk_mpfx_context(reslimit& lim, mpfx_manager & m, unsynch_mpq_manager & qm, params_ref const & p = params_ref(), small_object_allocator * a = nullptr);
};
}

2
src/math/subpaving/subpaving_hwf.h
View file

@ -42,5 +42,5 @@ public:
context_hwf(reslimit& lim, f2n<hwf_manager> & m, params_ref const & p, small_object_allocator * a):context_t<config_hwf>(lim, config_hwf(m), p, a) {}
};
};
}

2
src/math/subpaving/subpaving_mpf.h
View file

@ -43,5 +43,5 @@ public:
context_mpf(reslimit& lim, f2n<mpf_manager> & m, params_ref const & p, small_object_allocator * a):context_t<config_mpf>(lim, config_mpf(m), p, a) {}
};
};
}

2
src/math/subpaving/subpaving_mpff.h
View file

@ -39,5 +39,5 @@ struct config_mpff {
typedef context_t<config_mpff> context_mpff;
};
}

2
src/math/subpaving/subpaving_mpfx.h
View file

@ -39,5 +39,5 @@ struct config_mpfx {
typedef context_t<config_mpfx> context_mpfx;
};
}

2
src/math/subpaving/subpaving_mpq.h
View file

@ -37,5 +37,5 @@ struct config_mpq {
typedef context_t<config_mpq> context_mpq;
};
}

2
src/math/subpaving/subpaving_t.h
View file

@ -845,5 +845,5 @@ public:
void operator()();
};
};
}

2
src/math/subpaving/subpaving_t_def.h
View file

@ -1952,4 +1952,4 @@ bool context_t<C>::check_invariant() const {
}
};
}