From 324396e4039d417057f466341fba45108c57dd24 Mon Sep 17 00:00:00 2001 From: Lev Date: Fri, 14 Sep 2018 17:12:49 -0700 Subject: [PATCH] separate the gomory cut functionality in a separate file Signed-off-by: Lev --- src/util/lp/CMakeLists.txt | 1 + src/util/lp/gomory.cpp | 227 +++++++++++++++++++++++++++++++++++++ src/util/lp/gomory.h | 36 ++++++ src/util/lp/int_solver.cpp | 221 +----------------------------------- src/util/lp/int_solver.h | 25 ++-- src/util/lp/lar_solver.cpp | 1 + src/util/lp/lar_solver.h | 5 +- 7 files changed, 284 insertions(+), 232 deletions(-) create mode 100644 src/util/lp/gomory.cpp create mode 100644 src/util/lp/gomory.h diff --git a/src/util/lp/CMakeLists.txt b/src/util/lp/CMakeLists.txt index bde6ed93a..539c68712 100644 --- a/src/util/lp/CMakeLists.txt +++ b/src/util/lp/CMakeLists.txt @@ -6,6 +6,7 @@ z3_add_component(lp core_solver_pretty_printer.cpp dense_matrix.cpp eta_matrix.cpp + gomory.cpp indexed_vector.cpp int_solver.cpp lar_solver.cpp diff --git a/src/util/lp/gomory.cpp b/src/util/lp/gomory.cpp new file mode 100644 index 000000000..dd3d7bbed --- /dev/null +++ b/src/util/lp/gomory.cpp @@ -0,0 +1,227 @@ +/*++ + Copyright (c) 2017 Microsoft Corporation + + Module Name: + + + + Abstract: + + + + Author: + Nikolaj Bjorner (nbjorner) + Lev Nachmanson (levnach) + + Revision History: + + + --*/ +#include "util/lp/gomory.h" +#include "util/lp/int_solver.h" +#include "util/lp/lar_solver.h" +namespace lp { + +class gomory::imp { + lar_term & m_t; // the term to return in the cut + mpq & m_k; // the right side of the cut + explanation& m_ex; // the conflict explanation + unsigned m_inf_col; // a basis column which has to be an integer but has a not integral value + const row_strip& m_row; + const int_solver& m_int_solver; + + + const impq & get_value(unsigned j) const { return m_int_solver.get_value(j); } + bool is_real(unsigned j) const { return m_int_solver.is_real(j); } + bool at_lower(unsigned j) const { return m_int_solver.at_lower(j); } + bool at_upper(unsigned j) const { return m_int_solver.at_upper(j); } + const impq & lower_bound(unsigned j) const { return m_int_solver.lower_bound(j); } + const impq & upper_bound(unsigned j) const { return m_int_solver.upper_bound(j); } + constraint_index column_lower_bound_constraint(unsigned j) const { return m_int_solver.column_lower_bound_constraint(j); } + constraint_index column_upper_bound_constraint(unsigned j) const { return m_int_solver.column_upper_bound_constraint(j); } + + void int_case_in_gomory_cut(const mpq & a, unsigned x_j, + mpq & lcm_den, const mpq& f_0, const mpq& one_minus_f_0) { + lp_assert(is_int(x_j)); + lp_assert(!a.is_int()); + mpq f_j = int_solver::fractional_part(a); + TRACE("gomory_cut_detail", + tout << a << " x_j" << x_j << " k = " << m_k << "\n"; + tout << "f_j: " << f_j << "\n"; + tout << "f_0: " << f_0 << "\n"; + tout << "1 - f_0: " << 1 - f_0 << "\n"; + tout << "at_lower(" << x_j << ") = " << at_lower(x_j) << std::endl; + ); + lp_assert (!f_j.is_zero()); + mpq new_a; + if (at_lower(x_j)) { + if (f_j <= one_minus_f_0) { + new_a = f_j / one_minus_f_0; + } + else { + new_a = (1 - f_j) / f_0; + } + m_k.addmul(new_a, lower_bound(x_j).x); + m_ex.push_justification(column_lower_bound_constraint(x_j), new_a); + } + else { + lp_assert(at_upper(x_j)); + if (f_j <= f_0) { + new_a = f_j / f_0; + } + else { + new_a = (mpq(1) - f_j) / one_minus_f_0; + } + new_a.neg(); // the upper terms are inverted + m_k.addmul(new_a, upper_bound(x_j).x); + m_ex.push_justification(column_upper_bound_constraint(x_j), new_a); + } + TRACE("gomory_cut_detail", tout << "new_a: " << new_a << " k: " << m_k << "\n";); + m_t.add_monomial(new_a, x_j); + lcm_den = lcm(lcm_den, denominator(new_a)); + } + + void real_case_in_gomory_cut(const mpq & a, unsigned x_j, const mpq& f_0, const mpq& one_minus_f_0) { + TRACE("gomory_cut_detail_real", tout << "real\n";); + mpq new_a; + if (at_lower(x_j)) { + if (a.is_pos()) { + new_a = a / one_minus_f_0; + } + else { + new_a = a / f_0; + new_a.neg(); + } + m_k.addmul(new_a, lower_bound(x_j).x); // is it a faster operation than + // k += lower_bound(x_j).x * new_a; + m_ex.push_justification(column_lower_bound_constraint(x_j), new_a); + } + else { + lp_assert(at_upper(x_j)); + if (a.is_pos()) { + new_a = a / f_0; + new_a.neg(); // the upper terms are inverted. + } + else { + new_a = a / one_minus_f_0; + } + m_k.addmul(new_a, upper_bound(x_j).x); // k += upper_bound(x_j).x * new_a; + m_ex.push_justification(column_upper_bound_constraint(x_j), new_a); + } + TRACE("gomory_cut_detail_real", tout << a << "*v" << x_j << " k: " << m_k << "\n";); + m_t.add_monomial(new_a, x_j); + } + + lia_move report_conflict_from_gomory_cut() { + lp_assert(m_k.is_pos()); + // conflict 0 >= k where k is positive + m_k.neg(); // returning 0 <= -k + return lia_move::conflict; + } + + void adjust_term_and_k_for_some_ints_case_gomory(mpq &lcm_den) { + lp_assert(!m_t.is_empty()); + auto pol = m_t.coeffs_as_vector(); + m_t.clear(); + if (pol.size() == 1) { + TRACE("gomory_cut_detail", tout << "pol.size() is 1" << std::endl;); + unsigned v = pol[0].second; + lp_assert(is_int(v)); + const mpq& a = pol[0].first; + m_k /= a; + if (a.is_pos()) { // we have av >= k + if (!m_k.is_int()) + m_k = ceil(m_k); + // switch size + m_t.add_monomial(- mpq(1), v); + m_k.neg(); + } else { + if (!m_k.is_int()) + m_k = floor(m_k); + m_t.add_monomial(mpq(1), v); + } + } else { + TRACE("gomory_cut_detail", tout << "pol.size() > 1" << std::endl;); + lcm_den = lcm(lcm_den, denominator(m_k)); + lp_assert(lcm_den.is_pos()); + if (!lcm_den.is_one()) { + // normalize coefficients of integer parameters to be integers. + for (auto & pi: pol) { + pi.first *= lcm_den; + SASSERT(!is_int(pi.second) || pi.first.is_int()); + } + m_k *= lcm_den; + } + // negate everything to return -pol <= -m_k + for (const auto & pi: pol) + m_t.add_monomial(-pi.first, pi.second); + m_k.neg(); + } + TRACE("gomory_cut_detail", tout << "k = " << m_k << std::endl;); + lp_assert(m_k.is_int()); + } +public: + lia_move create_cut() { + TRACE("gomory_cut", + tout << "applying cut at:\n"; m_int_solver.m_lar_solver->print_row(m_row, tout); tout << std::endl; + for (auto & p : m_row) { + m_int_solver.m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(p.var(), tout); + } + tout << "inf_col = " << m_inf_col << std::endl; + ); + + // gomory will be t <= k and the current solution has a property t > k + m_k = 1; + mpq lcm_den(1); + unsigned x_j; + mpq a; + bool some_int_columns = false; + mpq f_0 = int_solver::fractional_part(get_value(m_inf_col)); + mpq one_min_f_0 = 1 - f_0; + for (const auto & p : m_row) { + x_j = p.var(); + if (x_j == m_inf_col) + continue; + // make the format compatible with the format used in: Integrating Simplex with DPLL(T) + a = p.coeff(); + a.neg(); + if (is_real(x_j)) + real_case_in_gomory_cut(a, x_j, f_0, one_min_f_0); + else if (!a.is_int()) { // f_j will be zero and no monomial will be added + some_int_columns = true; + int_case_in_gomory_cut(a, x_j, lcm_den, f_0, one_min_f_0); + } + } + + if (m_t.is_empty()) + return report_conflict_from_gomory_cut(); + if (some_int_columns) + adjust_term_and_k_for_some_ints_case_gomory(lcm_den); + lp_assert(m_int_solver.current_solution_is_inf_on_cut()); + m_int_solver.m_lar_solver->subs_term_columns(m_t, m_k); + TRACE("gomory_cut", tout<<"gomory cut:"; m_int_solver.m_lar_solver->print_term(m_t, tout); tout << " <= " << m_k << std::endl;); + return lia_move::cut; + } + imp(lar_term & t, mpq & k, explanation& ex, unsigned basic_inf_int_j, const row_strip& row, const int_solver& int_slv ) : + m_t(t), + m_k(k), + m_ex(ex), + m_inf_col(basic_inf_int_j), + m_row(row), + m_int_solver(int_slv) + { + } + +}; + +lia_move gomory::create_cut() { + return m_imp->create_cut(); +} + +gomory::gomory(lar_term & t, mpq & k, explanation& ex, unsigned basic_inf_int_j, const row_strip& row, const int_solver& s) { + m_imp = alloc(imp, t, k, ex, basic_inf_int_j, row, s); +} + +gomory::~gomory() { dealloc(m_imp); } + +} diff --git a/src/util/lp/gomory.h b/src/util/lp/gomory.h new file mode 100644 index 000000000..b7946d6b0 --- /dev/null +++ b/src/util/lp/gomory.h @@ -0,0 +1,36 @@ +/*++ +Copyright (c) 2017 Microsoft Corporation + +Module Name: + + + +Abstract: + + + +Author: + Nikolaj Bjorner (nbjorner) + Lev Nachmanson (levnach) + +Revision History: + + +--*/ +#pragma once +#include "util/lp/lar_term.h" +#include "util/lp/lia_move.h" +#include "util/lp/explanation.h" +#include "util/lp/static_matrix.h" + +namespace lp { +class int_solver; +class gomory { + class imp; + imp *m_imp; +public : + gomory(lar_term & t, mpq & k, explanation& ex, unsigned basic_inf_int_j, const row_strip& row, const int_solver& s); + lia_move create_cut(); + ~gomory(); +}; +} diff --git a/src/util/lp/int_solver.cpp b/src/util/lp/int_solver.cpp index f12e93103..cc01a0038 100644 --- a/src/util/lp/int_solver.cpp +++ b/src/util/lp/int_solver.cpp @@ -8,6 +8,7 @@ #include "util/lp/lp_utils.h" #include #include "util/lp/monomial.h" +#include "util/lp/gomory.h" namespace lp { @@ -101,12 +102,7 @@ bool int_solver::is_gomory_cut_target(const row_strip& row) { unsigned j; for (const auto & p : row) { j = p.var(); - if (is_base(j)) continue; - if (is_free(j)) - return false; - if (!at_bound(j)) - return false; - if (!is_zero(get_value(j).y)) { + if (!is_base(j) && (!at_bound(j) || !is_zero(get_value(j).y))) { TRACE("gomory_cut", tout << "row is not gomory cut target:\n"; display_column(tout, j); tout << "infinitesimal: " << !is_zero(get_value(j).y) << "\n";); @@ -117,36 +113,6 @@ bool int_solver::is_gomory_cut_target(const row_strip& row) { } -void int_solver::real_case_in_gomory_cut(const mpq & a, unsigned x_j, const mpq& f_0, const mpq& one_minus_f_0) { - TRACE("gomory_cut_detail_real", tout << "real\n";); - mpq new_a; - if (at_low(x_j)) { - if (a.is_pos()) { - new_a = a / one_minus_f_0; - } - else { - new_a = a / f_0; - new_a.neg(); - } - m_k->addmul(new_a, lower_bound(x_j).x); // is it a faster operation than - // k += lower_bound(x_j).x * new_a; - m_ex->push_justification(column_lower_bound_constraint(x_j), new_a); - } - else { - lp_assert(at_upper(x_j)); - if (a.is_pos()) { - new_a = a / f_0; - new_a.neg(); // the upper terms are inverted. - } - else { - new_a = a / one_minus_f_0; - } - m_k->addmul(new_a, upper_bound(x_j).x); // k += upper_bound(x_j).x * new_a; - m_ex->push_justification(column_upper_bound_constraint(x_j), new_a); - } - TRACE("gomory_cut_detail_real", tout << a << "*v" << x_j << " k: " << *m_k << "\n";); - m_t->add_monomial(new_a, x_j); -} constraint_index int_solver::column_upper_bound_constraint(unsigned j) const { return m_lar_solver->get_column_upper_bound_witness(j); @@ -157,99 +123,6 @@ constraint_index int_solver::column_lower_bound_constraint(unsigned j) const { } -void int_solver::int_case_in_gomory_cut(const mpq & a, unsigned x_j, - mpq & lcm_den, const mpq& f_0, const mpq& one_minus_f_0) { - lp_assert(is_int(x_j)); - lp_assert(!a.is_int()); - mpq f_j = fractional_part(a); - TRACE("gomory_cut_detail", - tout << a << " x_j" << x_j << " k = " << *m_k << "\n"; - tout << "f_j: " << f_j << "\n"; - tout << "f_0: " << f_0 << "\n"; - tout << "1 - f_0: " << 1 - f_0 << "\n"; - tout << "at_low(" << x_j << ") = " << at_low(x_j) << std::endl; - ); - lp_assert (!f_j.is_zero()); - mpq new_a; - if (at_low(x_j)) { - if (f_j <= one_minus_f_0) { - new_a = f_j / one_minus_f_0; - } - else { - new_a = (1 - f_j) / f_0; - } - m_k->addmul(new_a, lower_bound(x_j).x); - m_ex->push_justification(column_lower_bound_constraint(x_j), new_a); - } - else { - lp_assert(at_upper(x_j)); - if (f_j <= f_0) { - new_a = f_j / f_0; - } - else { - new_a = (mpq(1) - f_j) / one_minus_f_0; - } - new_a.neg(); // the upper terms are inverted - m_k->addmul(new_a, upper_bound(x_j).x); - m_ex->push_justification(column_upper_bound_constraint(x_j), new_a); - } - TRACE("gomory_cut_detail", tout << "new_a: " << new_a << " k: " << *m_k << "\n";); - m_t->add_monomial(new_a, x_j); - lcm_den = lcm(lcm_den, denominator(new_a)); -} - -lia_move int_solver::report_conflict_from_gomory_cut() { - TRACE("empty_pol",); - lp_assert(m_k->is_pos()); - // conflict 0 >= k where k is positive - m_k->neg(); // returning 0 <= -k - return lia_move::conflict; -} - -void int_solver::gomory_cut_adjust_t_and_k(vector> & pol, - lar_term & t, - mpq &k, - bool some_ints, - mpq & lcm_den) { - if (!some_ints) - return; - - t.clear(); - if (pol.size() == 1) { - unsigned v = pol[0].second; - lp_assert(is_int(v)); - bool k_is_int = k.is_int(); - const mpq& a = pol[0].first; - k /= a; - if (a.is_pos()) { // we have av >= k - if (!k_is_int) - k = ceil(k); - // switch size - t.add_monomial(- mpq(1), v); - k.neg(); - } else { - if (!k_is_int) - k = floor(k); - t.add_monomial(mpq(1), v); - } - } else if (some_ints) { - lcm_den = lcm(lcm_den, denominator(k)); - lp_assert(lcm_den.is_pos()); - if (!lcm_den.is_one()) { - // normalize coefficients of integer parameters to be integers. - for (auto & pi: pol) { - pi.first *= lcm_den; - SASSERT(!is_int(pi.second) || pi.first.is_int()); - } - k *= lcm_den; - } - // negate everything to return -pol <= -k - for (const auto & pi: pol) - t.add_monomial(-pi.first, pi.second); - k.neg(); - } -} - bool int_solver::current_solution_is_inf_on_cut() const { const auto & x = m_lar_solver->m_mpq_lar_core_solver.m_r_x; impq v = m_t->apply(x); @@ -261,95 +134,11 @@ bool int_solver::current_solution_is_inf_on_cut() const { return v * sign > (*m_k) * sign; } -void int_solver::adjust_term_and_k_for_some_ints_case_gomory(mpq &lcm_den) { - lp_assert(!m_t->is_empty()); - auto pol = m_t->coeffs_as_vector(); - m_t->clear(); - if (pol.size() == 1) { - TRACE("gomory_cut_detail", tout << "pol.size() is 1" << std::endl;); - unsigned v = pol[0].second; - lp_assert(is_int(v)); - const mpq& a = pol[0].first; - (*m_k) /= a; - if (a.is_pos()) { // we have av >= k - if (!(*m_k).is_int()) - (*m_k) = ceil((*m_k)); - // switch size - m_t->add_monomial(- mpq(1), v); - (*m_k).neg(); - } else { - if (!(*m_k).is_int()) - (*m_k) = floor((*m_k)); - m_t->add_monomial(mpq(1), v); - } - } else { - TRACE("gomory_cut_detail", tout << "pol.size() > 1" << std::endl;); - lcm_den = lcm(lcm_den, denominator((*m_k))); - lp_assert(lcm_den.is_pos()); - if (!lcm_den.is_one()) { - // normalize coefficients of integer parameters to be integers. - for (auto & pi: pol) { - pi.first *= lcm_den; - SASSERT(!is_int(pi.second) || pi.first.is_int()); - } - (*m_k) *= lcm_den; - } - // negate everything to return -pol <= -(*m_k) - for (const auto & pi: pol) - m_t->add_monomial(-pi.first, pi.second); - (*m_k).neg(); - } - TRACE("gomory_cut_detail", tout << "k = " << (*m_k) << std::endl;); - lp_assert((*m_k).is_int()); -} - - - - lia_move int_solver::mk_gomory_cut( unsigned inf_col, const row_strip & row) { - lp_assert(column_is_int_inf(inf_col)); - TRACE("gomory_cut", - tout << "applying cut at:\n"; m_lar_solver->print_row(row, tout); tout << std::endl; - for (auto & p : row) { - m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(p.var(), tout); - } - tout << "inf_col = " << inf_col << std::endl; - ); - - // gomory will be t <= k and the current solution has a property t > k - *m_k = 1; - mpq lcm_den(1); - unsigned x_j; - mpq a; - bool some_int_columns = false; - mpq f_0 = int_solver::fractional_part(get_value(inf_col)); - mpq one_min_f_0 = 1 - f_0; - for (const auto & p : row) { - x_j = p.var(); - if (x_j == inf_col) - continue; - // make the format compatible with the format used in: Integrating Simplex with DPLL(T) - a = p.coeff(); - a.neg(); - if (is_real(x_j)) - real_case_in_gomory_cut(a, x_j, f_0, one_min_f_0); - else if (!a.is_int()) { // f_j will be zero and no monomial will be added - some_int_columns = true; - int_case_in_gomory_cut(a, x_j, lcm_den, f_0, one_min_f_0); - } - } - - if (m_t->is_empty()) - return report_conflict_from_gomory_cut(); - if (some_int_columns) - adjust_term_and_k_for_some_ints_case_gomory(lcm_den); - - lp_assert(current_solution_is_inf_on_cut()); - m_lar_solver->subs_term_columns(*m_t); - TRACE("gomory_cut", tout<<"precut:"; m_lar_solver->print_term(*m_t, tout); tout << " <= " << *m_k << std::endl;); - return lia_move::cut; + gomory gc(*m_t, *m_k, *m_ex, inf_col, row, *this); + return gc.create_cut(); } lia_move int_solver::proceed_with_gomory_cut(unsigned j) { @@ -1121,7 +910,7 @@ bool int_solver::at_bound(unsigned j) const { } } -bool int_solver::at_low(unsigned j) const { +bool int_solver::at_lower(unsigned j) const { auto & mpq_solver = m_lar_solver->m_mpq_lar_core_solver.m_r_solver; switch (mpq_solver.m_column_types[j] ) { case column_type::fixed: diff --git a/src/util/lp/int_solver.h b/src/util/lp/int_solver.h index 414ca3006..dfe51711c 100644 --- a/src/util/lp/int_solver.h +++ b/src/util/lp/int_solver.h @@ -53,6 +53,13 @@ public: bool move_non_basic_column_to_bounds(unsigned j); lia_move check_wrapper(lar_term& t, mpq& k, explanation& ex); bool is_base(unsigned j) const; + bool is_real(unsigned j) const; + const impq & lower_bound(unsigned j) const; + const impq & upper_bound(unsigned j) const; + bool is_int(unsigned j) const; + const impq & get_value(unsigned j) const; + bool at_lower(unsigned j) const; + bool at_upper(unsigned j) const; private: @@ -79,10 +86,7 @@ private: void add_to_explanation_from_fixed_or_boxed_column(unsigned j); lia_move patch_nbasic_columns(); bool get_freedom_interval_for_column(unsigned j, bool & inf_l, impq & l, bool & inf_u, impq & u, mpq & m); - const impq & lower_bound(unsigned j) const; - const impq & upper_bound(unsigned j) const; - bool is_int(unsigned j) const; - bool is_real(unsigned j) const; +private: bool is_boxed(unsigned j) const; bool is_fixed(unsigned j) const; bool is_free(unsigned j) const; @@ -91,7 +95,6 @@ private: void set_value_for_nbasic_column_ignore_old_values(unsigned j, const impq & new_val); bool non_basic_columns_are_at_bounds() const; bool is_feasible() const; - const impq & get_value(unsigned j) const; bool column_is_int_inf(unsigned j) const; void trace_inf_rows() const; lia_move branch_or_sat(); @@ -104,13 +107,9 @@ private: bool move_non_basic_columns_to_bounds(); void branch_infeasible_int_var(unsigned); lia_move mk_gomory_cut(unsigned inf_col, const row_strip& row); - lia_move report_conflict_from_gomory_cut(); - void adjust_term_and_k_for_some_ints_case_gomory(mpq& lcm_den); lia_move proceed_with_gomory_cut(unsigned j); bool is_gomory_cut_target(const row_strip&); bool at_bound(unsigned j) const; - bool at_low(unsigned j) const; - bool at_upper(unsigned j) const; bool has_low(unsigned j) const; bool has_upper(unsigned j) const; unsigned row_of_basic_column(unsigned j) const; @@ -125,17 +124,13 @@ public: lp_assert(is_rational(n)); return n.x - floor(n.x); } -private: - void real_case_in_gomory_cut(const mpq & a, unsigned x_j, const mpq& f_0, const mpq& one_minus_f_0); - void int_case_in_gomory_cut(const mpq & a, unsigned x_j, mpq & lcm_den, const mpq& f_0, const mpq& one_minus_f_0); constraint_index column_upper_bound_constraint(unsigned j) const; constraint_index column_lower_bound_constraint(unsigned j) const; - void display_row_info(std::ostream & out, unsigned row_index) const; - void gomory_cut_adjust_t_and_k(vector> & pol, lar_term & t, mpq &k, bool num_ints, mpq &lcm_den); bool current_solution_is_inf_on_cut() const; -public: + bool shift_var(unsigned j, unsigned range); private: + void display_row_info(std::ostream & out, unsigned row_index) const; unsigned random(); bool has_inf_int() const; lia_move create_branch_on_column(int j); diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index 89044c7d6..3a53b6068 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -1645,6 +1645,7 @@ void lar_solver::push_and_register_term(lar_term* t) { // terms var_index lar_solver::add_term(const vector> & coeffs, const mpq &m_v) { + TRACE("add_term_lar_solver", print_linear_combination_of_column_indices(coeffs, tout);); if (strategy_is_undecided()) return add_term_undecided(coeffs, m_v); diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h index 6ef0ea596..3c0ed4fbf 100644 --- a/src/util/lp/lar_solver.h +++ b/src/util/lp/lar_solver.h @@ -535,7 +535,7 @@ public: return m_columns_to_ul_pairs()[j].lower_bound_witness(); } - void subs_term_columns(lar_term& t) { + void subs_term_columns(lar_term& t, mpq & rs) { vector> columns_to_subs; for (const auto & m : t.m_coeffs) { unsigned tj = adjust_column_index_to_term_index(m.first); @@ -545,9 +545,12 @@ public: for (const auto & p : columns_to_subs) { auto it = t.m_coeffs.find(p.first); lp_assert(it != t.m_coeffs.end()); + const lar_term& lt = get_term(p.second); mpq v = it->second; t.m_coeffs.erase(it); t.m_coeffs[p.second] = v; + if (lt.m_v.is_zero()) continue; + rs -= v * lt.m_v; } }