diff --git a/src/math/lp/dioph_eq.cpp b/src/math/lp/dioph_eq.cpp index 2783b83c7..363c3c6ca 100644 --- a/src/math/lp/dioph_eq.cpp +++ b/src/math/lp/dioph_eq.cpp @@ -35,11 +35,11 @@ namespace lp { } 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; + 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) { @@ -58,7 +58,7 @@ namespace lp { }; std::ostream& print_S(std::ostream & out) const { - out << "S:\n"; + out << "S:\n"; for (unsigned i : m_s) { out << "x" << m_eprime[i].m_e.j() << " ->\n"; print_eprime_entry(i, out); @@ -103,9 +103,9 @@ namespace lp { if (!term.c().is_zero()) { if (term.c().is_pos()) - out << " + " << term.c(); + out << " + " << term.c(); else - out << " - " << -term.c(); + out << " - " << -term.c(); } out << ", j():"<< term.j(); return out; @@ -137,6 +137,7 @@ namespace lp { std::list m_f; // F = {λ(t):t in m_f} // set S std::list m_s; // S = {λ(t): t in m_s} + u_map m_k2s; // k is substituted by using equation in m_eprime[m_k2s[k]] unsigned m_conflict_index = -1; // m_eprime[m_conflict_index] gives the conflict @@ -179,9 +180,9 @@ namespace lp { continue; } eprime_pair pair = {t, lar_term(i)}; - m_f.push_back(m_f.size()); + m_f.push_back(static_cast(m_f.size())); m_eprime.push_back(pair); - TRACE("dioph_eq", print_eprime_entry(m_f.size() - 1, tout);); + TRACE("dioph_eq", print_eprime_entry(static_cast(m_f.size()) - 1, tout);); } } @@ -189,11 +190,11 @@ namespace lp { mpq gcd_of_coeffs(const term_o& t) { mpq g(0); for (const auto & p : t) { - if (g.is_zero()) { + if (g.is_zero()) g = abs(p.coeff()); - } else { + else g = gcd(g, p.coeff()); - } + if (g.is_one()) break; } return g; } @@ -214,7 +215,7 @@ namespace lp { if (new_c.is_int() == false) { TRACE("dioph_eq", print_term_o(ep.m_e, tout << "conflict m_e:") << std::endl; - tout << "g:" << g << std::endl; + tout << "g:" << g << std::endl; print_lar_term_L(ep.m_l, tout << "m_l:") << std::endl; for (const auto & p: ep.m_l) { tout << p.coeff() << "("; print_term_o(m_eprime[p.j()].m_e, tout) << ")"<< std::endl; @@ -222,7 +223,7 @@ namespace lp { tout << "S:\n"; for (const auto& t : m_sigma) { tout << "x" << t.m_key << " -> "; - print_term_o(t.m_value, tout) << std::endl; + print_term_o(t.m_value, tout) << std::endl; } ); return false; @@ -233,9 +234,9 @@ namespace lp { ep.m_e.c() = new_c; ep.m_l *= (1/g); TRACE("dioph_eq", - tout << "ep_m_e:"; print_term_o(ep.m_e, tout) << std::endl; - tout << "ep.ml:"; - print_lar_term_L(ep.m_l, tout ) << std::endl;); + tout << "ep_m_e:"; print_term_o(ep.m_e, tout) << std::endl; + tout << "ep.ml:"; + print_lar_term_L(ep.m_l, tout ) << std::endl;); } return true; @@ -250,6 +251,91 @@ namespace lp { } return true; } + void init_term_from_constraint(term_o& t, const lar_base_constraint& c) { + for (const auto& p: c.coeffs()) { + t.add_monomial(p.first, p.second); + } + t.c() = - c.rhs(); + } + void subs_k(const eprime_pair& e, unsigned k, term_o& t, std::queue & q) { + const mpq& k_coeff = e.m_e.get_coeff(k); + bool is_one = k_coeff.is_one(); + SASSERT(is_one || k_coeff.is_minus_one()); + t.erase_var(k); + if (is_one) { + for (const auto& p: e.m_e) { + if (p.j() == k) continue; + t.add_monomial(-p.coeff(), p.j()); + if (m_k2s.contains(p.j())) + q.push(p.j()); + } + t.c() -= e.m_e.c(); + } else { // k_coeff is -1 + for (const auto& p: e.m_e) { + if (p.j() == k) continue; + t.add_monomial(p.coeff(), p.j()); + if (m_k2s.contains(p.j())) + q.push(p.j()); + } + t.c() += e.m_e.c(); + } + } + void process_q_with_S(std::queue &q, term_o& t, const lar_base_constraint& c) { + while (!q.empty()) { + unsigned k = q.front(); + TRACE("dioph_eq", tout << "k:" << k << std::endl;); + q.pop(); + const eprime_pair& e = m_eprime[m_k2s[k]]; + TRACE("dioph_eq", tout << "k:" << k << " in: "; print_eprime_entry(m_k2s[k], tout) << std::endl;); + subs_k(e, k, t, q); + } + TRACE("dioph_eq", tout << "substituted t:"; print_term_o(t, tout) << std::endl;); + } + + void tighten_with_S(const lar_base_constraint& c) { + TRACE("dioph_eq", tout << "constraint:"; lra.constraints().display(tout, c) << std::endl;); + SASSERT(c.is_active()); + term_o t; + std::queue q; + for (const auto& p: c.coeffs()) { + if (m_k2s.contains(p.second)) + q.push(p.second); + } + if (q.empty()) return; + init_term_from_constraint(t, c); + TRACE("dioph_eq", tout << "t:"; print_term_o(t, tout) << std::endl;); + process_q_with_S(q, t, c); + TRACE("dioph_eq", tout << "tighten:"; print_term_o(t, tout) << std::endl;); + mpq g = gcd_of_coeffs(t); + if (g.is_one()) return; + mpq new_c = t.c() / g; + if (new_c.is_int()) return; + std::cout << "."; return; + if (c.kind() == lconstraint_kind::EQ) { + TRACE("dioph_eq", tout << "conflict:" <::iterator pick_eh() { return m_f.begin(); // TODO: make a smarter joice @@ -281,9 +367,12 @@ namespace lp { const mpq& k_coeff = e.get_coeff(k); TRACE("dioph_eq", print_eprime_entry(e_index, tout << "before:") << std::endl; - tout << "k_coeff:" << k_coeff << std::endl;); + tout << "k_coeff:" << k_coeff << std::endl;); if (!l_term.is_empty()) { - add_operator(m_eprime[e_index].m_l, -k_sign*k_coeff, l_term); + if (k_sign == 1) + add_operator(m_eprime[e_index].m_l, -k_coeff, l_term); + else + add_operator(m_eprime[e_index].m_l, k_coeff, l_term); } e.substitute(k_subst, k); TRACE("dioph_eq", print_eprime_entry(e_index, tout << "after :") << std::endl;); @@ -369,7 +458,7 @@ namespace lp { print_eprime_entry(m_eprime.size()-1, tout);); substitute_var_on_f(k, 1, k_subs, xt_entry.m_l, e_index); - // the term to eliminate the fresh variable + // the term to eliminate the fresh variable term_o xt_subs = xt_term.clone(); xt_subs.add_monomial(mpq(1), xt); TRACE("dioph_eq", tout << "xt_subs: x~"<< fresh_index(xt) << " -> "; print_term_o(xt_subs, tout) << std::endl;); @@ -377,11 +466,19 @@ namespace lp { } std::ostream& print_eprime_entry(unsigned i, std::ostream& out) const { 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; + 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 + + // k is the index of the index of the variable with the coefficient +-1 that is being substituted + void move_from_f_to_s(unsigned k, std::list::iterator it) { + m_s.push_back(*it); + m_k2s.insert(k, *it); + m_f.erase(it); + } + + // 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]; @@ -392,8 +489,7 @@ namespace lp { if (ahk.is_one()) { eprime_entry.m_e.j() = k; TRACE("dioph_eq", tout << "push to S:\n"; print_eprime_entry(*eh_it, tout);); - m_s.push_back(*eh_it); - m_f.erase(eh_it); + move_from_f_to_s(k, eh_it); eliminate_var_in_f(eprime_entry, k , k_sign); } else { fresh_var_step(*eh_it, k); @@ -416,25 +512,25 @@ namespace lp { unsigned fresh_index(unsigned j) const {return UINT_MAX - j;} void remove_fresh_variables(term_o& t) { - std::queue q; - for (auto p : t) { - 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.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()) && t.contains(p.j())) + std::queue q; + for (auto p : t) { + 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.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()) && t.contains(p.j())) q.push(p.j()); - } + } } bool is_fresh_var(unsigned j) const { diff --git a/src/math/lp/lar_term.h b/src/math/lp/lar_term.h index acd268122..f98fc4868 100644 --- a/src/math/lp/lar_term.h +++ b/src/math/lp/lar_term.h @@ -59,7 +59,6 @@ public: } } - void add_var(lpvar j) { rational c(1); add_monomial(c, j); @@ -71,10 +70,15 @@ public: unsigned size() const { return static_cast(m_coeffs.size()); } - u_map const & coeffs() const { + u_map const &coeffs() const { return m_coeffs; } + u_map &coeffs() { + return m_coeffs; + } + + void subst_by_term(const lar_term& t, unsigned term_column) { auto it = this->m_coeffs.find_core(term_column); if (it == nullptr) return; @@ -138,6 +142,16 @@ public: SASSERT(it != nullptr); return it->get_data().m_value; } + + mpq & get_coeff(unsigned j){ + auto* it = m_coeffs.find_core(j); + SASSERT(it != nullptr); + return it->get_data().m_value; + } + + void erase_var(unsigned j) { + m_coeffs.erase(j); + } // the monomial ax[j] is substituted by ax[k] void subst_index(unsigned j, unsigned k) { auto* it = m_coeffs.find_core(j);