diff --git a/src/math/dd/dd_pdd.cpp b/src/math/dd/dd_pdd.cpp index 300a8350a..858552ffb 100644 --- a/src/math/dd/dd_pdd.cpp +++ b/src/math/dd/dd_pdd.cpp @@ -338,6 +338,56 @@ namespace dd { } } + /* + * Compare leading monomials. + * The pdd format makes lexicographic comparison easy: compare based on + * the top variable and descend depending on whether hi(x) == hi(y) + */ + bool pdd_manager::lt(pdd const& a, pdd const& b) { + PDD x = a.root; + PDD y = b.root; + if (x == y) return false; + while (true) { + SASSERT(x != y); + if (is_val(x)) + return !is_val(y) || val(x) < val(y); + if (is_val(y)) + return false; + if (level(x) == level(y)) { + if (hi(x) == hi(y)) { + x = lo(x); + y = lo(y); + } + else { + x = hi(x); + y = hi(y); + } + } + else { + return level(x) > level(y); + } + } + } + + /** + Compare leading terms of pdds + */ + bool pdd_manager::different_leading_term(pdd const& a, pdd const& b) { + PDD x = a.root; + PDD y = b.root; + while (true) { + if (x == y) return false; + if (is_val(x) || is_val(y)) return true; + if (level(x) == level(y)) { + x = hi(x); + y = hi(y); + } + else { + return false; + } + } + } + void pdd_manager::push(PDD b) { m_pdd_stack.push_back(b); } diff --git a/src/math/dd/dd_pdd.h b/src/math/dd/dd_pdd.h index e06d6b3a2..b148a7b75 100644 --- a/src/math/dd/dd_pdd.h +++ b/src/math/dd/dd_pdd.h @@ -141,7 +141,6 @@ namespace dd { small_object_allocator m_alloc; mutable svector m_mark; mutable unsigned m_mark_level; - mutable svector m_count; mutable svector m_todo; bool m_disable_gc; bool m_is_new_node; @@ -230,6 +229,8 @@ namespace dd { pdd mul(rational const& c, pdd const& b); pdd reduce(pdd const& a, pdd const& b); bool try_spoly(pdd const& a, pdd const& b, pdd& r); + bool lt(pdd const& a, pdd const& b); + bool different_leading_term(pdd const& a, pdd const& b); std::ostream& display(std::ostream& out); std::ostream& display(std::ostream& out, pdd const& b); @@ -243,7 +244,8 @@ namespace dd { pdd_manager* m; pdd(unsigned root, pdd_manager* m): root(root), m(m) { m->inc_ref(root); } public: - pdd(pdd & other): root(other.root), m(other.m) { m->inc_ref(root); } + pdd(pdd_manager& pm): root(0), m(&pm) { SASSERT(is_zero()); } + pdd(pdd const& other): root(other.root), m(other.m) { m->inc_ref(root); } pdd(pdd && other): root(0), m(other.m) { std::swap(root, other.root); } pdd& operator=(pdd const& other); ~pdd() { m->dec_ref(root); } @@ -252,6 +254,7 @@ namespace dd { unsigned var() const { return m->var(root); } rational const& val() const { SASSERT(is_val()); return m->val(root); } bool is_val() const { return m->is_val(root); } + bool is_zero() const { return m->is_zero(root); } pdd operator+(pdd const& other) const { return m->add(*this, other); } pdd operator-(pdd const& other) const { return m->sub(*this, other); } @@ -259,10 +262,13 @@ namespace dd { pdd operator*(rational const& other) const { return m->mul(other, *this); } pdd operator+(rational const& other) const { return m->add(other, *this); } pdd reduce(pdd const& other) const { return m->reduce(*this, other); } + bool different_leading_term(pdd const& other) const { return m->different_leading_term(*this, other); } std::ostream& display(std::ostream& out) const { return m->display(out, *this); } bool operator==(pdd const& other) const { return root == other.root; } bool operator!=(pdd const& other) const { return root != other.root; } + bool operator<(pdd const& b) const { return m->lt(*this, b); } + unsigned size() const { return m->pdd_size(*this); } }; @@ -276,6 +282,7 @@ namespace dd { inline pdd operator-(pdd const& b, int x) { return b + (-rational(x)); } + std::ostream& operator<<(std::ostream& out, pdd const& b); } diff --git a/src/math/grobner/CMakeLists.txt b/src/math/grobner/CMakeLists.txt index 1b56c775f..0ec842f75 100644 --- a/src/math/grobner/CMakeLists.txt +++ b/src/math/grobner/CMakeLists.txt @@ -1,6 +1,7 @@ z3_add_component(grobner SOURCES grobner.cpp + pdd_grobner.cpp COMPONENT_DEPENDENCIES ast ) diff --git a/src/math/grobner/pdd_grobner.cpp b/src/math/grobner/pdd_grobner.cpp new file mode 100644 index 000000000..3b961cd96 --- /dev/null +++ b/src/math/grobner/pdd_grobner.cpp @@ -0,0 +1,253 @@ +/*++ + Copyright (c) 2019 Microsoft Corporation + + Abstract: + + Grobner core based on pdd representation of polynomials + + Author: + Nikolaj Bjorner (nbjorner) + Lev Nachmanson (levnach) + + --*/ + +#include "math/grobner/pdd_grobner.h" + +namespace dd { + + /*** + The main algorithm maintains two sets (S, A), + where S is m_processed, and A is m_to_simplify. + Initially S is empty and A contains the initial equations. + + Each step proceeds as follows: + - pick a in A, and remove a from A + - simplify a using S + - simplify S using a + - for s in S: + b = superpose(a, s) + add b to A + - add a to S + - simplify A using a + */ + + grobner::grobner(reslimit& lim, pdd_manager& m) : + m(m), + m_limit(lim), + m_conflict(false) + {} + + + grobner::~grobner() { + del_equations(0); + } + + void grobner::compute_basis() { + while (!done() && compute_basis_step()) { + TRACE("grobner", display(tout);); + } + } + + bool grobner::compute_basis_step() { + m_stats.m_compute_steps++; + equation* e = pick_next(); + if (!e) return false; + equation& eq = *e; + SASSERT(!eq.is_processed()); + if (!simplify_using(eq, m_processed)) return false; + if (eliminate(eq)) return true; + if (check_conflict(eq)) return false; + if (!simplify_using(m_processed, eq)) return false; + TRACE("grobner", tout << "eq = "; display_equation(tout, eq);); + superpose(eq); + toggle_processed(eq); + return simplify_using(m_to_simplify, eq); + } + + grobner::equation* grobner::pick_next() { + equation* eq = nullptr; + ptr_buffer to_delete; + for (auto* curr : m_to_simplify) { + if (!eq|| is_simpler(*curr, *eq)) { + eq = curr; + } + } + for (auto* e : to_delete) + del_equation(e); + if (eq) + m_to_simplify.erase(eq); + return eq; + } + + void grobner::superpose(equation const & eq) { + for (equation * target : m_processed) { + superpose(eq, *target); + } + } + + bool grobner::simplify_using(equation& eq, equation_set const& eqs) { + bool simplified, changed_leading_term; + TRACE("grobner", tout << "simplifying: "; display_equation(tout, eq); tout << "using equalities of size " << eqs.size() << "\n";); + do { + simplified = false; + for (equation* p : eqs) { + if (simplify_source_target(*p, eq, changed_leading_term)) { + simplified = true; + } + if (canceled() || eq.poly().is_val()) { + break; + } + } + } + while (simplified && !eq.poly().is_val()); + + TRACE("grobner", tout << "simplification result: "; display_equation(tout, eq);); + + check_conflict(eq); + return !done(); + } + + /* + Use the given equation to simplify equations in set + */ + bool grobner::simplify_using(equation_set& set, equation const& eq) { + TRACE("grobner", tout << "poly " << eq.poly() << "\n";); + ptr_buffer to_delete; + ptr_buffer to_move; + bool changed_leading_term = false; + for (equation* target : set) { + if (canceled()) + return false; + if (simplify_source_target(eq, *target, changed_leading_term)) { + if (is_trivial(*target)) + to_delete.push_back(target); + else if (check_conflict(*target)) + return false; + else if (changed_leading_term && target->is_processed()) { + to_move.push_back(target); + } + } + } + for (equation* eq : to_delete) + del_equation(eq); + for (equation* eq : to_move) + toggle_processed(*eq); + return true; + } + + // return true iff simplified + bool grobner::simplify_source_target(equation const& source, equation& target, bool& changed_leading_term) { + TRACE("grobner", tout << "simplifying: " << target.poly() << "\nusing: " << source.poly() << "\n";); + m_stats.m_simplified++; + pdd r = source.poly().reduce(target.poly()); + if (r == source.poly()) { + return false; + } + changed_leading_term = target.is_processed() && m.different_leading_term(r, source.poly()); + target = r; + target = m_dep_manager.mk_join(target.dep(), source.dep()); + update_stats_max_degree_and_size(target); + + TRACE("grobner", tout << "simplified " << target.poly() << "\n";); + return true; + } + + // let eq1: ab+q=0, and eq2: ac+e=0, then qc - eb = 0 + void grobner::superpose(equation const& eq1, equation const& eq2) { + TRACE("grobner_d", tout << "eq1="; display_equation(tout, eq1) << "eq2="; display_equation(tout, eq2);); + pdd r(m); + if (!m.try_spoly(eq1.poly(), eq2.poly(), r)) return; + m_stats.m_superposed++; + if (r.is_zero()) return; + unsigned idx = m_equations.size(); + equation* eq = alloc(equation, r, m_dep_manager.mk_join(eq1.dep(), eq2.dep()), idx); + update_stats_max_degree_and_size(*eq); + if (r.is_val()) set_conflict(); + m_to_simplify.insert(eq); + } + + void grobner::toggle_processed(equation& eq) { + if (eq.is_processed()) { + m_to_simplify.insert(&eq); + m_processed.erase(&eq); + } + else { + m_processed.insert(&eq); + m_to_simplify.erase(&eq); + } + eq.set_processed(!eq.is_processed()); + } + + grobner::equation_set const& grobner::equations() { + m_all_eqs.reset(); + for (equation* eq : m_equations) if (eq) m_all_eqs.insert(eq); + return m_all_eqs; + } + + void grobner::reset() { + del_equations(0); + SASSERT(m_equations.empty()); + m_processed.reset(); + m_to_simplify.reset(); + m_stats.reset(); + m_conflict = false; + } + + void grobner::add(pdd const& p, u_dependency * dep) { + if (p.is_zero()) return; + if (p.is_val()) set_conflict(); + equation * eq = alloc(equation, p, dep, m_equations.size()); + m_equations.push_back(eq); + m_to_simplify.insert(eq); + } + + void grobner::del_equations(unsigned old_size) { + for (unsigned i = m_equations.size(); i-- > old_size; ) { + del_equation(m_equations[i]); + } + m_equations.shrink(old_size); + } + + void grobner::del_equation(equation* eq) { + if (!eq) return; + unsigned idx = eq->idx(); + m_processed.erase(eq); + m_to_simplify.erase(eq); + dealloc(m_equations[idx]); + m_equations[idx] = nullptr; + } + + bool grobner::canceled() { + return m_limit.get_cancel_flag(); + } + + bool grobner::done() { + return m_to_simplify.size() + m_processed.size() >= m_config.m_eqs_threshold || canceled() || m_conflict; + } + + void grobner::update_stats_max_degree_and_size(const equation& e) { + // TBD m_stats.m_max_expr_size = std::max(m_stats.m_max_expr_size, e.poly().size()); + // TBD m_stats.m_max_expr_degree = std::max(m_stats.m_max_expr_degree, e.poly.degree()); + } + + std::ostream& grobner::print_stats(std::ostream & out) const { + return out << "stats:\nsteps = " << m_stats.m_compute_steps << "\nsimplified: " << + m_stats.m_simplified << "\nsuperposed: " << + m_stats.m_superposed << "\nexpr degree: " << m_stats.m_max_expr_degree << + "\nexpr size: " << m_stats.m_max_expr_size << "\n"; + } + + std::ostream& grobner::display_equation(std::ostream & out, const equation & eq) const { + out << "expr = " << eq.poly() << "\n"; + m_print_dep(eq.dep(), out); + return out; + } + + std::ostream& grobner::display(std::ostream& out) const { + out << "processed\n"; for (auto e : m_processed) display_equation(out, *e); + out << "to_simplify\n"; for (auto e : m_to_simplify) display_equation(out, *e); + return out; + } + +} + diff --git a/src/math/grobner/pdd_grobner.h b/src/math/grobner/pdd_grobner.h new file mode 100644 index 000000000..cf07ca548 --- /dev/null +++ b/src/math/grobner/pdd_grobner.h @@ -0,0 +1,126 @@ +/*++ + Copyright (c) 2017 Microsoft Corporation + + Module Name: + + + + Abstract: + + + + Author: + Nikolaj Bjorner (nbjorner) + Lev Nachmanson (levnach) + + Revision History: + + + --*/ +#pragma once + +#include "math/dd/dd_pdd.h" +#include "util/dependency.h" +#include "util/obj_hashtable.h" +#include "util/region.h" +#include "util/rlimit.h" + +namespace dd { + +class grobner { +public: + struct stats { + unsigned m_simplified; + unsigned m_max_expr_size; + unsigned m_max_expr_degree; + unsigned m_superposed; + unsigned m_compute_steps; + void reset() { memset(this, 0, sizeof(*this)); } + }; + + struct config { + unsigned m_eqs_threshold; + }; + +private: + class equation { + bool m_processed; //!< state + unsigned m_idx; //!< position at m_equations + pdd m_poly; // simplified expressionted monomials + u_dependency * m_dep; //!< justification for the equality + public: + equation(pdd const& p, u_dependency* d, unsigned idx): + m_processed(false), + m_idx(idx), + m_poly(p), + m_dep(d) + {} + + const pdd& poly() const { return m_poly; } + u_dependency * dep() const { return m_dep; } + unsigned hash() const { return m_idx; } + unsigned idx() const { return m_idx; } + void operator=(pdd& p) { m_poly = p; } + void operator=(u_dependency* d) { m_dep = d; } + bool is_processed() const { return m_processed; } + void set_processed(bool p) { m_processed = p; } + }; + + typedef obj_hashtable equation_set; + typedef ptr_vector equation_vector; + typedef std::function print_dep_t; + + pdd_manager& m; + reslimit& m_limit; + stats m_stats; + config m_config; + print_dep_t m_print_dep; + equation_vector m_equations; + equation_set m_processed; + equation_set m_to_simplify; + mutable u_dependency_manager m_dep_manager; + equation_set m_all_eqs; + bool m_conflict; +public: + grobner(reslimit& lim, pdd_manager& m); + ~grobner(); + + void operator=(print_dep_t& pd) { m_print_dep = pd; } + void operator=(config const& c) { m_config = c; } + + void reset(); + void compute_basis(); + void add(pdd const&, u_dependency * dep); + equation_set const& equations(); + u_dependency_manager& dep() const { return m_dep_manager; } + std::ostream& display_equation(std::ostream& out, const equation& eq) const; + std::ostream& display(std::ostream& out) const; + +private: + bool compute_basis_step(); + equation* pick_next(); + bool canceled(); + bool done(); + void superpose(equation const& eq1, equation const& eq2); + void superpose(equation const& eq); + bool simplify_source_target(equation const& source, equation& target, bool& changed_leading_term); + bool simplify_using(equation_set& set, equation const& eq); + bool simplify_using(equation& eq, equation_set const& eqs); + void toggle_processed(equation& eq); + + bool eliminate(equation& eq) { return is_trivial(eq) && (del_equation(&eq), true); } + bool is_trivial(equation const& eq) const { return eq.poly().is_zero(); } + bool is_simpler(equation const& eq1, equation const& eq2) { return eq1.poly() < eq2.poly(); } + bool check_conflict(equation const& eq) { return eq.poly().is_val() && !is_trivial(eq) && (set_conflict(), true); } + void set_conflict() { m_conflict = true; } + + void del_equations(unsigned old_size); + void del_equation(equation* eq); + + void update_stats_max_degree_and_size(const equation& e); + + std::ostream& print_stats(std::ostream&) const; + +}; + +} diff --git a/src/util/dependency.h b/src/util/dependency.h index 0b6772852..c1aaac724 100644 --- a/src/util/dependency.h +++ b/src/util/dependency.h @@ -321,6 +321,8 @@ public: // Implement old dependency manager used by interval and groebner typedef scoped_dependency_manager v_dependency_manager; typedef scoped_dependency_manager::dependency v_dependency; +typedef scoped_dependency_manager u_dependency_manager; +typedef scoped_dependency_manager::dependency u_dependency; #endif /* DEPENDENCY_H_ */