diff --git a/src/math/lp/dioph_eq.cpp b/src/math/lp/dioph_eq.cpp index 1e6c25fce..6899ee1d8 100644 --- a/src/math/lp/dioph_eq.cpp +++ b/src/math/lp/dioph_eq.cpp @@ -3,12 +3,12 @@ #include "math/lp/lar_solver.h" #include "math/lp/lp_utils.h" #include +#include namespace lp { // This class represents a term with an added constant number c, in form sum {x_i*a_i} + c. class dioph_eq::imp { - class term_o:public lar_term { mpq m_c; public: @@ -34,9 +34,32 @@ namespace lp { this->c() += a * t.c(); this->m_coeffs.erase(term_column); } + + friend term_o operator*(const mpq& k, const term_o& term) { + term_o r; + for (const auto& p : term) + r.add_monomial(p.coeff()*k, p.j()); + r.c() = k*term.c(); + return r; + } +#if Z3DEBUG + friend bool operator== (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; + } +#endif + term_o& operator += (const term_o& t) { + for (const auto & p: t) { + add_monomial(p.coeff(), p.j()); + } + m_c += t.c(); + return *this; + } + }; - std::ostream& print_lar_term_L(lar_term & t, std::ostream & out) const { + std::ostream& print_lar_term_L(const lar_term & t, std::ostream & out) const { return print_linear_combination_customized(t.coeffs_as_vector(), [](int j)->std::string {return "y"+std::to_string(j);}, out ); } @@ -108,6 +131,18 @@ namespace lp { unsigned m_conflict_index = -1; // m_eprime[m_conflict_index] gives the conflict imp(int_solver& lia, lar_solver& lra): lia(lia), lra(lra) {} + + term_o row_to_term(const row_strip& row) const { + term_o t; + const auto lcm = get_denominators_lcm(row); + for (const auto & p: row) + if (lia.is_fixed(p.var())) + t.c() += lia.lower_bound(p.var()).x; + else + t.add_monomial(lcm * p.coeff(), p.var()); + return t; + } + void init() { unsigned n_of_rows = lra.A_r().row_count(); @@ -126,22 +161,16 @@ namespace lp { auto & row = lra.get_row(i); if (!all_vars_are_int(row)) continue; const auto lcm = get_denominators_lcm(row); - term_o t; - for (const auto & p: row) - if (lia.is_fixed(p.var())) - t.c() += lia.lower_bound(p.var()).x; - else - t.add_monomial(lcm * p.coeff(), p.var()); + term_o t = row_to_term(row); t.j() = i; //highjack this field to point to the original tableau row, TODO - is it used? - TRACE("dioph_eq", tout << "row = "; - lra.print_row(row, tout) << std::endl; - tout << "t:"; print_term_o(t, tout) << std::endl;); + TRACE("dioph_eq", tout << "row = "; lra.print_row(row, tout) << std::endl;); unsigned lvar = static_cast(m_f.size()); lar_term lt = lar_term(lvar); eprime_pair pair = {t, lt}; m_eprime.push_back(pair); m_f.push_back(lvar); + TRACE("dioph_eq", print_eprime_entry(lvar, tout);); } } @@ -207,6 +236,11 @@ namespace lp { } lia_move check() { + std::cout << "ddd = " << ++lp_settings::ddd << std::endl; + if (lp::lp_settings::ddd == 5) { + enable_trace("dioph_eq"); + enable_trace("dioph_eq_fresh"); + } init(); while(m_f.size()) { if (!normalize_by_gcd()) @@ -218,19 +252,19 @@ namespace lp { std::list::iterator pick_eh() { return m_f.begin(); // TODO: make a smarter joice } + void substitute_var_on_f(unsigned k, int k_sign, const term_o& k_subst, const lar_term& l_term) { - TRACE("dioph_eq", print_term_o(k_subst, tout<< "x" << k<< " -> ") << std::endl;); + TRACE("dioph_eq", print_term_o(k_subst, tout<< "x" << k<< " -> ") << std::endl; tout << "l_term:"; print_lar_term_L(l_term, tout) << std::endl;); for (unsigned e_index: m_f) { term_o& e = m_eprime[e_index].m_e; if (!e.contains(k)) continue; const mpq& k_coeff = e.get_coeff(k); - TRACE("dioph_eq", print_term_o(e, tout << "before:") << std::endl; + TRACE("dioph_eq", print_eprime_entry(e_index, tout << "before:") << std::endl; tout << "k_coeff:" << k_coeff << std::endl;); - m_eprime[e_index].m_l += mpq(k_sign)*k_coeff*l_term; + m_eprime[e_index].m_l += mpq(-k_sign)*k_coeff*l_term; e.substitute(k_subst, k); - TRACE("dioph_eq", print_term_o(e, tout << "after :") << std::endl; - print_lar_term_L(m_eprime[e_index].m_l, tout) << std::endl;); + TRACE("dioph_eq", print_eprime_entry(e_index, tout << "after :") << std::endl;); } } @@ -269,6 +303,7 @@ namespace lp { substitute_var_on_f(k, k_sign, t, e_pair.m_l) ; } +// k is the variable to substitute void fresh_var_step(eprime_pair & e_pair, mpq& ahk, unsigned k) { // step 7 from the paper auto & eh = e_pair.m_e; @@ -298,11 +333,12 @@ namespace lp { xt_term.add_monomial(q, p.j()); k_subs.add_monomial(-q, p.j()); } + eprime_pair xt_entry = {xt_term, lar_term()}; // 0 for m_l field + m_eprime.push_back(xt_entry); TRACE("dioph_eq", tout << "xt_term:"; print_term_o(xt_term, tout) << std::endl; - tout << "k_subs:"; print_term_o(k_subs, tout) << std::endl;); - eprime_pair x_t_entry = {xt_term, lar_term()}; // 0 for m_l field - m_eprime.push_back(x_t_entry); - substitute_var_on_f(k, 1, k_subs, e_pair.m_l); + tout << "k_subs:"; print_term_o(k_subs, tout) << std::endl; + print_eprime_entry(m_eprime.size()-1, tout);); + substitute_var_on_f(k, 1, k_subs, xt_entry.m_l); // the term to eliminate the fresh variable term_o xt_subs = xt_term.clone(); @@ -310,14 +346,17 @@ namespace lp { TRACE("dioph_eq", tout << "xt_subs: x~"<< fresh_index(xt) << " -> "; print_term_o(xt_subs, tout) << std::endl;); m_sigma.insert(xt, xt_subs); } + std::ostream& print_eprime_entry(unsigned i, std::ostream& out) { + out << "m_eprime[" << i << "]:{\n"; + print_term_o(m_eprime[i].m_e, out << "\tm_e:") << "," << std::endl; + print_lar_term_L(m_eprime[i].m_l, out<< "\tm_l:") << "\n}"<< std::endl; + return out; + } // this is the step 6 or 7 of the algorithm void rewrite_eqs() { auto eh_it = pick_eh(); auto eprime_entry = m_eprime[*eh_it]; - TRACE("dioph_eq", tout << "eprime_entry[" << *eh_it << "]{\n"; - print_term_o(m_eprime[*eh_it].m_e, tout << "\tm_e:") << "," << std::endl; - print_lar_term_L(m_eprime[*eh_it].m_l, tout<< "\tm_l:") << "\n}"<< std::endl;); - + TRACE("dioph_eq", print_eprime_entry(*eh_it, tout);); term_o& eh = eprime_entry.m_e; auto [ahk, k, k_sign] = find_minimal_abs_coeff(eh); TRACE("dioph_eq", tout << "ahk:" << ahk << ", k:" << k << ", k_sign:" << k_sign << std::endl;); @@ -333,10 +372,8 @@ namespace lp { void explain(lp::explanation& ex) { SASSERT(ex.empty()); + TRACE("dioph_eq", tout << "conflict:"; print_eprime_entry(m_conflict_index, tout) << std::endl;); auto & ep = m_eprime[m_conflict_index]; - for (const auto & p : ep.m_l) { - remove_fresh_variables(m_eprime[p.j()].m_e); - } u_dependency* dep = nullptr; for (const auto & pl : ep.m_l) { for (const auto & p : lra.get_row(m_eprime[pl.j()].m_e.j())) @@ -346,36 +383,26 @@ namespace lp { TRACE("dioph_eq", lra.print_expl(tout, ex);); } unsigned fresh_index(unsigned j) const {return UINT_MAX - j;} - struct fresh_lt { - bool operator()(unsigned v1, unsigned v2) const { return v1 > v2; } - }; - typedef heap fresh_queue; - + void remove_fresh_variables(term_o& t) { - fresh_queue q(fresh_index(m_last_fresh_x_var)+1); - std::set processed; + std::queue q; for (auto p : t) { - auto j = p.j(); - if (is_fresh_var(j)) { - q.insert(fresh_index(j)); - processed.insert(j); + if (is_fresh_var(p.j())) { + q.push(p.j()); } } CTRACE("dioph_eq_fresh", q.empty() == false, print_term_o(t, tout)<< std::endl); while (!q.empty()) { - unsigned j = q.erase_min(); - j = fresh_index(j); + unsigned j = q.front(); + q.pop(); const term_o& sub_t = m_sigma.find(j); TRACE("dioph_eq_fresh", tout << "x~" << fresh_index(j) << "; sub_t:"; print_term_o(sub_t, tout) << std::endl;); t.substitute(sub_t, j); TRACE("dioph_eq_fresh", tout << "sub_t:"; print_term_o(sub_t, tout) << std::endl; tout << "after sub:";print_term_o(t, tout) << std::endl;); for (auto p : sub_t) - if (is_fresh_var(p.j())) { - if (processed.find(p.j()) != processed.end()) continue; - q.insert(fresh_index(p.j())); - processed.insert(j); - } + if (is_fresh_var(p.j()) && t.contains(p.j())) + q.push(p.j()); } } diff --git a/src/math/lp/lar_term.h b/src/math/lp/lar_term.h index b316daa9f..77fcb983d 100644 --- a/src/math/lp/lar_term.h +++ b/src/math/lp/lar_term.h @@ -169,12 +169,12 @@ public: return ret; } - + friend lar_term operator*(const mpq& k, const lar_term& term) { lar_term r; - for (const auto& p : term) { + for (const auto& p : term) { r.add_monomial(p.coeff()*k, p.j()); - } + } return r; } diff --git a/src/math/lp/lp_utils.h b/src/math/lp/lp_utils.h index d6943bf50..aa41f032c 100644 --- a/src/math/lp/lp_utils.h +++ b/src/math/lp/lp_utils.h @@ -86,6 +86,10 @@ bool is_non_decreasing(const K& v) { template std::ostream& print_linear_combination_customized(const vector> & coeffs, std::function var_str, std::ostream & out) { + if (coeffs.size() == 0) { + out << "0"; + return out; + } bool first = true; for (const auto & it : coeffs) { T val = it.first;