Skip to content

Commit e041fe9

Browse files
committed
slack
Signed-off-by: Lev Nachmanson <[email protected]>
1 parent 7ca94e8 commit e041fe9

File tree

3 files changed

+213
-20
lines changed

3 files changed

+213
-20
lines changed

src/math/lp/dioph_eq.cpp

Lines changed: 206 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -291,6 +291,12 @@ namespace lp {
291291
tp, [](int j) -> std::string { return "x" + std::to_string(j); }, out);
292292
}
293293

294+
std::ostream& print_espace(std::ostream & out) const {
295+
out << "m_espace:";
296+
print_term_o(create_term_from_espace(), out) << std::endl;
297+
return out;
298+
}
299+
294300
std::ostream& print_term_o(term_o const& term, std::ostream& out) const {
295301
if (term.size() == 0 && term.c().is_zero()) {
296302
out << "0";
@@ -376,6 +382,28 @@ namespace lp {
376382
r.add_monomial(p.coeff(), p.var());
377383
}
378384
return r;
385+
}
386+
// Deep‑copy this term_with_index into 'target'
387+
void copy(term_with_index& target) const {
388+
if (&target == this)
389+
return;
390+
// Clear target's data and index
391+
target.clear();
392+
// Copy monomial data verbatim.
393+
target.m_data = m_data;
394+
395+
// Re‑create a compact m_index that tracks only the used variables.
396+
unsigned max_var = 0;
397+
for (auto const& iv : m_data)
398+
max_var = std::max(max_var, iv.var());
399+
400+
target.m_index.assign(max_var + 1, -1);
401+
for (unsigned idx = 0; idx < m_data.size(); ++idx)
402+
target.m_index[m_data[idx].var()] = static_cast<int>(idx);
403+
404+
#if Z3DEBUG
405+
SASSERT(target.invariant());
406+
#endif
379407
}
380408

381409
auto size() const { return m_data.size(); }
@@ -493,7 +521,8 @@ namespace lp {
493521
term_with_index m_lspace;
494522
// m_espace is for operations on m_e_matrix rows
495523
term_with_index m_espace;
496-
524+
term_with_index m_espace_backup;
525+
497526
bijection m_k2s;
498527
bij_map<std::pair<lar_term, unsigned>> m_fresh_k2xt_terms;
499528
// m_row2fresh_defs[i] is the set of all fresh variables xt
@@ -512,7 +541,7 @@ namespace lp {
512541
std::unordered_map<unsigned, std::unordered_set<unsigned>> m_columns_to_terms;
513542

514543
unsigned m_normalize_conflict_index = UINT_MAX; // the row index of the conflict
515-
mpq m_normalize_conflict_gcd; // the gcd of the coefficients the m_e_matrix[m_normalize_conflict_gcd].
544+
mpq m_normalize_conflict_gcd; // the gcd of the coefficients of m_e_matrix[m_normalize_conflict_gcd].
516545
void reset_conflict() { m_normalize_conflict_index = UINT_MAX; }
517546
bool has_conflict_index() const { return m_normalize_conflict_index != UINT_MAX; }
518547
void set_rewrite_conflict(unsigned idx, const mpq& gcd) {
@@ -524,7 +553,7 @@ namespace lp {
524553

525554
void undo_add_term_method(const lar_term* t) {
526555
TRACE("d_undo", tout << "t:" << t << ", t->j():" << t->j() << std::endl;);
527-
TRACE("dioph_eq", lra.print_term(*t, tout); tout << ", t->j() =" << t->j() << std::endl;);
556+
TRACE("dio", lra.print_term(*t, tout); tout << ", t->j() =" << t->j() << std::endl;);
528557
if (!contains(m_active_terms, t)) {
529558
for (auto i = m_added_terms.size(); i-- > 0; ) {
530559
if (m_added_terms[i] != t)
@@ -548,7 +577,7 @@ namespace lp {
548577
SASSERT(std::find(m_added_terms.begin(), m_added_terms.end(), t) == m_added_terms.end());
549578
SASSERT(contains(m_active_terms, t));
550579
m_active_terms.erase(t);
551-
TRACE("dioph_eq", tout << "the deleted term column in m_l_matrix" << std::endl; for (auto p : m_l_matrix.column(t->j())) { tout << "p.coeff():" << p.coeff() << ", row " << p.var() << std::endl; } tout << "m_l_matrix has " << m_l_matrix.column_count() << " columns" << std::endl; tout << "and " << m_l_matrix.row_count() << " rows" << std::endl; print_lar_term_L(*t, tout); tout << "; t->j()=" << t->j() << std::endl;);
580+
TRACE("dio", tout << "the deleted term column in m_l_matrix" << std::endl; for (auto p : m_l_matrix.column(t->j())) { tout << "p.coeff():" << p.coeff() << ", row " << p.var() << std::endl; } tout << "m_l_matrix has " << m_l_matrix.column_count() << " columns" << std::endl; tout << "and " << m_l_matrix.row_count() << " rows" << std::endl; print_lar_term_L(*t, tout); tout << "; t->j()=" << t->j() << std::endl;);
552581
shrink_matrices();
553582
}
554583

@@ -884,7 +913,7 @@ namespace lp {
884913
}
885914
void substitute_with_fresh_def(unsigned ei, unsigned j, const mpq& alpha) {
886915
const lar_term& sub_term = m_fresh_k2xt_terms.get_by_key(j).first;
887-
TRACE("dioph_eq", print_lar_term_L(sub_term, tout) << std::endl;);
916+
TRACE("dio", print_lar_term_L(sub_term, tout) << std::endl;);
888917
SASSERT(sub_term.get_coeff(j).is_one());
889918
// we need to eliminate alpha*j in ei's row
890919
add_term_to_entry(-alpha, sub_term, ei);
@@ -1192,13 +1221,13 @@ namespace lp {
11921221
// The function returns true if and only if there is no conflict.
11931222
bool normalize_e_by_gcd(unsigned ei, mpq& g) {
11941223
mpq& e = m_sum_of_fixed[ei];
1195-
TRACE("dioph_eq", print_entry(ei, tout) << std::endl;);
1224+
TRACE("dio", print_entry(ei, tout) << std::endl;);
11961225
g = gcd_of_coeffs(m_e_matrix.m_rows[ei], false);
11971226
if (g.is_zero() || g.is_one()) {
11981227
SASSERT(g.is_one() || e.is_zero());
11991228
return true;
12001229
}
1201-
TRACE("dioph_eq", tout << "g:" << g << std::endl;);
1230+
TRACE("dio", tout << "g:" << g << std::endl;);
12021231
mpq c_g = e / g;
12031232
if (c_g.is_int()) {
12041233
for (auto& p : m_e_matrix.m_rows[ei]) {
@@ -1210,7 +1239,7 @@ namespace lp {
12101239
p.coeff() /= g;
12111240
}
12121241

1213-
TRACE("dioph_eq", tout << "ep_m_e:";
1242+
TRACE("dio", tout << "ep_m_e:";
12141243
print_entry(ei, tout) << std::endl;);
12151244
SASSERT(entry_invariant(ei));
12161245
return true;
@@ -1221,7 +1250,7 @@ namespace lp {
12211250

12221251
lia_move subs_qfront_by_fresh(unsigned k, protected_queue& q, unsigned j) {
12231252
const lar_term& e = m_fresh_k2xt_terms.get_by_key(k).first;
1224-
TRACE("dioph_eq", tout << "k:" << k << ", in ";
1253+
TRACE("dio", tout << "k:" << k << ", in ";
12251254
print_term_o(create_term_from_espace(), tout) << std::endl;
12261255
tout << "subs with e:";
12271256
print_lar_term_L(e, tout) << std::endl;);
@@ -1241,7 +1270,7 @@ namespace lp {
12411270
}
12421271

12431272
// there is no change in m_l_matrix
1244-
TRACE("dioph_eq", tout << "after subs k:" << k << "\n";
1273+
TRACE("dio", tout << "after subs k:" << k << "\n";
12451274
print_term_o(create_term_from_espace(), tout) << std::endl;
12461275
tout << "m_lspace:{"; print_lar_term_L(m_lspace.m_data, tout);
12471276
tout << "}, opened:"; print_ml(m_lspace.to_term(), tout) << std::endl;);
@@ -1568,18 +1597,182 @@ namespace lp {
15681597
}
15691598
}
15701599

1600+
1601+
// returns false if all coefficients are +-1 and true otherwise
1602+
mpq find_second_smallest_coeff_in_espace() {
1603+
mpq a; // first smallest
1604+
mpq b; // second smallest
1605+
for (const auto & [c, v]: m_espace) {
1606+
if (var_is_fresh(v))
1607+
return mpq(1);
1608+
mpq ac = abs(c);
1609+
if (a.is_zero())
1610+
a = ac; // first smallest init
1611+
else if (ac < a) {
1612+
b = a; // init b
1613+
a = ac; // first smallest improved
1614+
}
1615+
else if (ac < b) {
1616+
b = ac; // second smallest improved
1617+
}
1618+
}
1619+
return b;
1620+
}
1621+
1622+
lia_move try_improve_gcd_on_espace(unsigned term_j) {
1623+
mpq second_smallest_coeff = find_second_smallest_coeff_in_espace();
1624+
TRACE("dio", tout << "second_smallest_coeff:" << second_smallest_coeff << std::endl;);
1625+
if (abs(second_smallest_coeff) <= mpq(1)) {
1626+
//can we improve here?
1627+
return lia_move::undef;
1628+
}
1629+
auto r = try_make_gcd(second_smallest_coeff, true, term_j);
1630+
if (r == lia_move::undef) {
1631+
r = try_make_gcd(second_smallest_coeff, false, term_j);
1632+
}
1633+
return r;
1634+
}
1635+
1636+
struct restore_espace {
1637+
term_with_index & m_original;
1638+
term_with_index & m_backup;
1639+
restore_espace(term_with_index & orig, term_with_index & backup): m_original(orig), m_backup(backup) {
1640+
m_original.copy(m_backup);
1641+
}
1642+
~restore_espace() {
1643+
m_backup.copy(m_original);
1644+
}
1645+
};
1646+
1647+
// g is a candidate for new gcd
1648+
lia_move try_make_gcd(const mpq& g, bool upper_bound, unsigned term_j) {
1649+
restore_espace re(m_espace, m_espace_backup);
1650+
if ((upper_bound && !lra.column_has_upper_bound(term_j)) ||
1651+
(!upper_bound && !lra.column_has_lower_bound(term_j)))
1652+
return lia_move::undef;
1653+
mpq new_bound = upper_bound? lra.get_upper_bound(term_j).x: lra.get_lower_bound(term_j).x;
1654+
TRACE("dio", tout << "upper_bound:" << upper_bound << ", new_bound:" << new_bound << std::endl;);
1655+
for (const auto &[c, v] : m_espace) {
1656+
if (abs(c) == g) continue;
1657+
if (upper_bound) {
1658+
if (!supplement_to_g_upper(c, v, g, new_bound, term_j))
1659+
return lia_move::undef;
1660+
} else {
1661+
if (!supplement_to_g_lower(c, v, g, new_bound, term_j))
1662+
return lia_move::undef;
1663+
}
1664+
}
1665+
TRACE("dio", print_espace(tout); tout << "g:" << g << std::endl;);
1666+
SASSERT(gcd_of_coeffs(m_espace.m_data, true) == g);
1667+
mpq rs_g = new_bound % g;
1668+
if (rs_g.is_neg())
1669+
rs_g += g;
1670+
SASSERT(!rs_g.is_neg());
1671+
new_bound -= rs_g;
1672+
TRACE("dio", tout << "new_bound:" << new_bound << std::endl;);
1673+
if (upper_bound) {
1674+
if (new_bound < lra.get_upper_bound(term_j).x) {
1675+
NOT_IMPLEMENTED_YET();
1676+
}
1677+
} else {
1678+
if (new_bound > lra.get_lower_bound(term_j).x) {
1679+
NOT_IMPLEMENTED_YET();
1680+
}
1681+
}
1682+
1683+
return lia_move::undef;
1684+
}
1685+
1686+
// new_bound initially is set to the original lower bound of term_j
1687+
bool supplement_to_g_lower(const mpq& c, unsigned lj, const mpq & g, mpq& new_bound, unsigned term_j) {
1688+
restore_espace re(m_espace, m_espace_backup);
1689+
auto r = c % g;
1690+
TRACE("dio", tout << "lj:" << lj << ", g:"<< g << ", new_bound:" << new_bound << ", r:" << r << std::endl;);
1691+
if (r.is_zero())
1692+
return true; // the coefficient is divisible by g
1693+
if (r.is_neg())
1694+
r += g;
1695+
SASSERT((c - r) % g == 0 && r < g && r.is_pos());
1696+
unsigned j = local_to_lar_solver(lj);
1697+
if (lra.column_is_free(j)) return false;
1698+
if (lra.column_is_bounded(j)) {
1699+
const auto& ub = lra.get_upper_bound(j).x;
1700+
const auto& lb = lra.get_lower_bound(j).x;
1701+
TRACE("dio", tout << "lb:" << lb<< ", ub:" << ub << "\n";);
1702+
/*
1703+
If lb >= 0 then we can substract r*xj from term_j and be sure that the new term does not get bigger, from the other side it cannot diminish by more than r*bu.
1704+
In this case we need to update new_bound -= r*ub.
1705+
1706+
1707+
*/
1708+
if (!lb.is_neg()) {
1709+
m_espace.add(-r, lj);
1710+
new_bound -= r * ub;
1711+
TRACE("dio", print_espace(tout) << "\n"; tout << "new_bound:" << new_bound << std::endl;);
1712+
1713+
} else {
1714+
NOT_IMPLEMENTED_YET();
1715+
}
1716+
}
1717+
NOT_IMPLEMENTED_YET();
1718+
1719+
SASSERT(r.is_pos());
1720+
// m_espace <= new_bound
1721+
r = g - r;
1722+
TRACE("dio", tout << "r:" << r << std::endl;);
1723+
// m_espace:4x2 + 2x3 + x4 - 256 >= lb
1724+
// We have something like: c = 1, lj = 4,g = 2, then r = 1.
1725+
// If we know that 0 >= x[j] >= k and
1726+
// then term = m_espace >= m_espace+ r*x_lj >= bound + r*k
1727+
m_espace.add(r, lj);
1728+
new_bound += r*lra.get_upper_bound(j).x;
1729+
TRACE("dio", print_espace(tout); tout << "new_bound:" << new_bound << std::endl; );
1730+
1731+
return true;
1732+
}
1733+
1734+
void backup_espace() {
1735+
m_espace.copy(m_espace_backup);
1736+
}
1737+
1738+
// new_bound is initially let to the original upper bound of term_j
1739+
bool supplement_to_g_upper(const mpq& c, unsigned lj, const mpq & g, mpq& new_bound, unsigned term_j) {
1740+
auto r = c % g;
1741+
TRACE("dio", tout << "r:" << r << std::endl;);
1742+
if (r.is_zero())
1743+
return true; // the coefficient is divisible by g
1744+
if (r.is_neg())
1745+
r += g;
1746+
SASSERT(r.is_pos());
1747+
unsigned j = local_to_lar_solver(lj);
1748+
// m_espace <= new_bound
1749+
r = g - r;
1750+
TRACE("dio", tout << "r:" << r << std::endl;);
1751+
if (!lra.column_is_bounded(j)) return false;
1752+
// m_espace:4x2 + 2x3 + x4 - 256
1753+
// We have something like: c = 1, lj = 4,g = 2, then r = 1.
1754+
// If we know that 0 <= x[j] <= k and
1755+
// then term = m_espace <= m_espace+ r*x_lj <= new_bound + r*k
1756+
m_espace.add(r, lj);
1757+
new_bound += r*lra.get_upper_bound(j).x;
1758+
TRACE("dio", print_espace(tout); tout << "new_bound:" << new_bound << std::endl; );
1759+
1760+
return true;
1761+
}
15711762

15721763
lia_move tighten_on_espace(unsigned j) {
15731764
mpq g = gcd_of_coeffs(m_espace.m_data, true);
1574-
if (g.is_one())
1765+
if (g.is_one()) {
15751766
return lia_move::undef;
1767+
return try_improve_gcd_on_espace(j);
1768+
}
15761769
if (g.is_zero()) {
15771770
handle_constant_term(j);
15781771
if (!m_infeas_explanation.empty())
15791772
return lia_move::conflict;
15801773
return lia_move::undef;
15811774
}
1582-
// g is not trivial, trying to tighten the bounds
1775+
// g is non-trivial, trying to tighten the bounds
15831776
auto r = tighten_bounds_for_non_trivial_gcd(g, j, true);
15841777
if (r == lia_move::undef)
15851778
r = tighten_bounds_for_non_trivial_gcd(g, j, false);
@@ -2486,7 +2679,7 @@ namespace lp {
24862679
}
24872680
SASSERT(h == f_vector[ih]);
24882681
if (min_ahk.is_one()) {
2489-
TRACE("dioph_eq", tout << "push to S:\n"; print_entry(h, tout););
2682+
TRACE("dio", tout << "push to S:\n"; print_entry(h, tout););
24902683
move_entry_from_f_to_s(kh, h);
24912684
eliminate_var_in_f(h, kh, kh_sign);
24922685
f_vector[ih] = f_vector.back();

src/math/lp/int_solver.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -56,6 +56,7 @@ class int_solver {
5656

5757
lp_settings& settings();
5858
const lp_settings& settings() const;
59+
public:
5960
bool at_bound(unsigned j) const;
6061
bool has_lower(unsigned j) const;
6162
bool has_upper(unsigned j) const;

src/math/lp/lar_solver.cpp

Lines changed: 6 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1175,7 +1175,7 @@ namespace lp {
11751175
const vector<std::pair<mpq, unsigned>>& inf_row,
11761176
int inf_sign) const {
11771177

1178-
#if 0
1178+
#if 1
11791179
impq slack(0);
11801180

11811181
for (auto& [coeff, j] : inf_row) {
@@ -1185,6 +1185,7 @@ namespace lp {
11851185

11861186
#define get_sign(_x_) (_x_.is_pos() ? 1 : (_x_.is_neg() ? -1 : 0))
11871187
int sign = get_sign(slack);
1188+
11881189
#endif
11891190

11901191
for (auto& [coeff, j] : inf_row) {
@@ -1193,15 +1194,13 @@ namespace lp {
11931194
const column& ul = m_columns[j];
11941195
u_dependency* bound_constr_i = is_upper ? ul.upper_bound_witness() : ul.lower_bound_witness();
11951196

1196-
#if 0
1197-
if (false)
1198-
;
1199-
else if(is_upper) {
1197+
#if 1
1198+
if(is_upper) {
12001199
if (ul.previous_upper() != UINT_MAX) {
12011200
auto const& [_is_upper, _j, _bound, _column] = m_column_updates[ul.previous_upper()];
12021201
auto new_slack = slack + coeff * (_bound - get_upper_bound(j));
12031202
if (sign == get_sign(new_slack)) {
1204-
//verbose_stream() << "can weaken j" << j << " " << coeff << " " << get_upper_bound(j) << " " << _bound << "\n";
1203+
// verbose_stream() << "can weaken j" << j << " " << coeff << " " << get_upper_bound(j) << " " << _bound << "\n";
12051204
slack = new_slack;
12061205
bound_constr_i = _column.upper_bound_witness();
12071206
}
@@ -1212,7 +1211,7 @@ namespace lp {
12121211
auto const& [_is_upper, _j, _bound, _column] = m_column_updates[ul.previous_lower()];
12131212
auto new_slack = slack + coeff * (_bound - get_lower_bound(j));
12141213
if (sign == get_sign(new_slack)) {
1215-
//verbose_stream() << "can weaken j" << j << " " << coeff << " " << get_lower_bound(j) << " " << _bound << "\n";
1214+
// verbose_stream() << "can weaken j" << j << " " << coeff << " " << get_lower_bound(j) << " " << _bound << "\n";
12161215
slack = new_slack;
12171216
bound_constr_i = _column.lower_bound_witness();
12181217
}

0 commit comments

Comments
 (0)