mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-06-07 14:43:23 +00:00
Code Activity

port Grobner

Browse source 
Signed-off-by: Lev Nachmanson <levnach@hotmail.com>
This commit is contained in:
Lev Nachmanson 2019-10-25 16:27:27 -07:00
parent 6009b738d6
commit a7418db611
3 changed files with 98 additions and 80 deletions
Download patch file Download diff file Expand all files Collapse all files

17
src/math/lp/nex.h
View file

@ -52,7 +52,10 @@ inline std::ostream & operator<<(std::ostream& out, expr_type t) {
class nex; class nex;
class nex_scalar; class nex_scalar;
class nex_pow; // forward definitions
// This is the class of non-linear expressions // This is the class of non-linear expressions
class nex { class nex {
public: public:
// the scalar and the variable have size 1 // the scalar and the variable have size 1
@ -69,7 +72,10 @@ public:
return true; return true;
} }
} }
virtual unsigned number_of_child_powers() const { return 0; }
virtual const nex* get_child_exp(unsigned) const { return this; }
virtual unsigned get_child_pow(unsigned) const { return 1; }
virtual bool all_factors_are_elementary() const { return true; }
bool is_sum() const { return type() == expr_type::SUM; } bool is_sum() const { return type() == expr_type::SUM; }
bool is_mul() const { return type() == expr_type::MUL; } bool is_mul() const { return type() == expr_type::MUL; }
bool is_var() const { return type() == expr_type::VAR; } bool is_var() const { return type() == expr_type::VAR; }
@ -80,7 +86,7 @@ public:
virtual bool contains(lpvar j) const { return false; } virtual bool contains(lpvar j) const { return false; }
virtual int get_degree() const = 0; virtual int get_degree() const = 0;
// simplifies the expression and also assigns the address of "this" to *e // simplifies the expression and also assigns the address of "this" to *e
virtual const rational& coeff() const { return rational::one(); }
#ifdef Z3DEBUG #ifdef Z3DEBUG
virtual void sort() {}; virtual void sort() {};
@ -93,6 +99,8 @@ std::ostream& operator<<(std::ostream& out, const nex&);
class nex_var : public nex { class nex_var : public nex {
lpvar m_j; lpvar m_j;
public: public:
unsigned number_of_child_powers() const { return 1; }
nex_pow get_nex_pow(unsigned j);
nex_var(lpvar j) : m_j(j) {} nex_var(lpvar j) : m_j(j) {}
nex_var() {} nex_var() {}
expr_type type() const { return expr_type::VAR; } expr_type type() const { return expr_type::VAR; }
@ -174,6 +182,11 @@ class nex_mul : public nex {
rational m_coeff; rational m_coeff;
vector<nex_pow> m_children; vector<nex_pow> m_children;
public: public:
virtual const nex* get_child_exp(unsigned j) const { return m_children[j].e(); }
virtual unsigned get_child_pow(unsigned j) const { return m_children[j].pow(); }
unsigned number_of_child_powers() const { return m_children.size(); }
nex_mul() : m_coeff(rational(1)) {} nex_mul() : m_coeff(rational(1)) {}
template <typename T> template <typename T>

138
src/math/lp/nla_grobner.cpp
View file

@ -244,12 +244,14 @@ nla_grobner::equation* nla_grobner::simplify_using_processed(equation* eq) {
} }
nex_mul* nla_grobner::get_highest_monomial(nex* e) const { const nex* nla_grobner::get_highest_monomial(const nex* e) const {
switch (e->type()) { switch (e->type()) {
case expr_type::MUL: case expr_type::MUL:
return to_mul(e); return to_mul(e);
case expr_type::SUM: case expr_type::SUM:
return to_mul((*to_sum(e))[0]); return *(to_sum(e)->begin());
case expr_type::VAR:
return e;
default: default:
TRACE("grobner", tout << *e << "\n";); TRACE("grobner", tout << *e << "\n";);
return nullptr; return nullptr;
@ -259,7 +261,7 @@ nex_mul* nla_grobner::get_highest_monomial(nex* e) const {
// target 2fg + 3fp + e = 0 // target 2fg + 3fp + e = 0
// target is replaced by 2(-k/3 - l/3)g + 3(-k/3 - l/3)p + e = -2/3kg -2/3lg - kp -lp + e // target is replaced by 2(-k/3 - l/3)g + 3(-k/3 - l/3)p + e = -2/3kg -2/3lg - kp -lp + e
bool nla_grobner::simplify_target_monomials(equation * source, equation * target) { bool nla_grobner::simplify_target_monomials(equation * source, equation * target) {
const nex_mul * high_mon = get_highest_monomial(source->exp()); auto * high_mon = get_highest_monomial(source->exp());
if (high_mon == nullptr) if (high_mon == nullptr)
return false; return false;
SASSERT(high_mon->all_factors_are_elementary()); SASSERT(high_mon->all_factors_are_elementary());
@ -280,7 +282,7 @@ bool nla_grobner::simplify_target_monomials(equation * source, equation * target
} }
bool nla_grobner::simplify_target_monomials_sum_check_only(nex_sum* targ_sum, bool nla_grobner::simplify_target_monomials_sum_check_only(nex_sum* targ_sum,
const nex_mul* high_mon) { const nex* high_mon) const {
for (auto t : *targ_sum) { for (auto t : *targ_sum) {
if (!t->is_mul()) continue; // what if t->is_var() if (!t->is_mul()) continue; // what if t->is_var()
if (divide_ignore_coeffs_check_only(to_mul(t), high_mon)) { if (divide_ignore_coeffs_check_only(to_mul(t), high_mon)) {
@ -295,7 +297,7 @@ bool nla_grobner::simplify_target_monomials_sum_check_only(nex_sum* targ_sum,
bool nla_grobner::simplify_target_monomials_sum(equation * source, bool nla_grobner::simplify_target_monomials_sum(equation * source,
equation * target, nex_sum* targ_sum, equation * target, nex_sum* targ_sum,
const nex_mul* high_mon) { const nex* high_mon) {
if (!simplify_target_monomials_sum_check_only(targ_sum, high_mon)) if (!simplify_target_monomials_sum_check_only(targ_sum, high_mon))
return false; return false;
unsigned targ_orig_size = targ_sum->size(); unsigned targ_orig_size = targ_sum->size();
@ -307,7 +309,7 @@ bool nla_grobner::simplify_target_monomials_sum(equation * source,
return true; return true;
} }
nex_mul* nla_grobner::divide_ignore_coeffs(nex* ej, const nex_mul* h) { nex_mul* nla_grobner::divide_ignore_coeffs(nex* ej, const nex* h) {
if (!ej->is_mul()) if (!ej->is_mul())
return nullptr; return nullptr;
nex_mul* m = to_mul(ej); nex_mul* m = to_mul(ej);
@ -317,16 +319,16 @@ nex_mul* nla_grobner::divide_ignore_coeffs(nex* ej, const nex_mul* h) {
} }
// return true if h divides t // return true if h divides t
bool nla_grobner::divide_ignore_coeffs_check_only(nex_mul* t , const nex_mul* h) { bool nla_grobner::divide_ignore_coeffs_check_only(nex_mul* t , const nex* h) const {
SASSERT(m_nex_creator.is_simplified(t) && m_nex_creator.is_simplified(h)); SASSERT(m_nex_creator.is_simplified(t) && m_nex_creator.is_simplified(h));
unsigned j = 0; // points to t unsigned j = 0; // points to t
for(auto & p : *h) { for(unsigned k = 0; k < h->number_of_child_powers(); k++) {
lpvar h_var = to_var(h->get_child_exp(k))->var();
bool p_swallowed = false; bool p_swallowed = false;
lpvar h_var = to_var(p.e())->var();
for (; j < t->size() && !p_swallowed; j++) { for (; j < t->size() && !p_swallowed; j++) {
auto &tp = (*t)[j]; auto &tp = (*t)[j];
if (to_var(tp.e())->var() == h_var) { if (to_var(tp.e())->var() == h_var) {
if (tp.pow() >= p.pow()) { if (tp.pow() >= static_cast<int>(h->get_child_pow(k))) {
p_swallowed = true; p_swallowed = true;
} }
} }
@ -341,20 +343,21 @@ bool nla_grobner::divide_ignore_coeffs_check_only(nex_mul* t , const nex_mul* h)
} }
// perform the division t / h, but ignores the coefficients // perform the division t / h, but ignores the coefficients
nex_mul * nla_grobner::divide_ignore_coeffs_perform(nex_mul* t, const nex_mul* h) { // h does not change
nex_mul * nla_grobner::divide_ignore_coeffs_perform(nex_mul* t, const nex* h) {
nex_mul * r = m_nex_creator.mk_mul(); nex_mul * r = m_nex_creator.mk_mul();
unsigned j = 0; // points to t unsigned j = 0; // points to t
for(auto & p : *h) { for(unsigned k = 0; k < h->number_of_child_powers(); k++) {
bool p_swallowed = false; lpvar h_var = to_var(h->get_child_exp(k))->var();
lpvar h_var = to_var(p.e())->var(); for (; j < t->size(); j++) {
for (; j < t->size() && !p_swallowed; j++) {
auto &tp = (*t)[j]; auto &tp = (*t)[j];
if (to_var(tp.e())->var() == h_var) { if (to_var(tp.e())->var() == h_var) {
if (tp.pow() >= p.pow()) { int h_pow = h->get_child_pow(k);
p_swallowed = true; SASSERT(tp.pow() >= h_pow);
if (tp.pow() > p.pow()) j++;
r->add_child_in_power(tp.e(), tp.pow() - p.pow()); if (tp.pow() > h_pow)
} r->add_child_in_power(tp.e(), tp.pow() - h_pow);
break;
} else { } else {
r->add_child_in_power(tp); r->add_child_in_power(tp);
} }
@ -368,7 +371,7 @@ nex_mul * nla_grobner::divide_ignore_coeffs_perform(nex_mul* t, const nex_mul* h
// and b*high_mon + e = 0, so high_mon = -e/b // and b*high_mon + e = 0, so high_mon = -e/b
// then targ_sum->children()[j] = - (c/b) * e*p // then targ_sum->children()[j] = - (c/b) * e*p
void nla_grobner::simplify_target_monomials_sum_j(equation * source, equation *target, nex_sum* targ_sum, const nex_mul* high_mon, unsigned j) { void nla_grobner::simplify_target_monomials_sum_j(equation * source, equation *target, nex_sum* targ_sum, const nex* high_mon, unsigned j) {
nex * ej = (*targ_sum)[j]; nex * ej = (*targ_sum)[j];
nex_mul * ej_over_high_mon = divide_ignore_coeffs(ej, high_mon); nex_mul * ej_over_high_mon = divide_ignore_coeffs(ej, high_mon);
@ -479,22 +482,24 @@ void nla_grobner::simplify_to_process(equation* eq) {
} }
// if e is the sum adds to r all children of e multiplied by beta, except the first one // if e is the sum then add to r all children of e multiplied by beta, except the first one
// which corresponds to the highest monomial // which corresponds to the highest monomial,
// otherwise do nothing
void nla_grobner::add_mul_skip_first(nex_sum* r, const rational& beta, nex *e, nex_mul* c) { void nla_grobner::add_mul_skip_first(nex_sum* r, const rational& beta, nex *e, nex_mul* c) {
SASSERT(e->is_sum());
if (e->is_sum()) { if (e->is_sum()) {
nex_sum *es = to_sum(e); nex_sum *es = to_sum(e);
for (unsigned j = 1; j < es->size(); j++) { for (unsigned j = 1; j < es->size(); j++) {
r->add_child(m_nex_creator.mk_mul(beta, (*es)[j], c)); r->add_child(m_nex_creator.mk_mul(beta, (*es)[j], c));
} }
TRACE("grobner", tout << "r = " << *r << "\n";); TRACE("grobner", tout << "r = " << *r << "\n";);
} else {
TRACE("grobner", tout << "e = " << *e << "\n";);
} }
} }
// let e1: alpha*ab+q=0, and e2: beta*ac+e=0, then beta*qc - alpha*eb = 0 // let e1: alpha*ab+q=0, and e2: beta*ac+e=0, then beta*qc - alpha*eb = 0
nex * nla_grobner::expr_superpose(nex* e1, nex* e2, nex_mul* ab, nex_mul* ac, nex_mul* b, nex_mul* c) { nex * nla_grobner::expr_superpose(nex* e1, nex* e2, const nex* ab, const nex* ac, nex_mul* b, nex_mul* c) {
TRACE("grobner", tout << "e1 = " << *e1 << "\ne2 = " << *e2 <<"\n";); TRACE("grobner", tout << "e1 = " << *e1 << "\ne2 = " << *e2 <<"\n";);
nex_sum * r = m_nex_creator.mk_sum(); nex_sum * r = m_nex_creator.mk_sum();
rational alpha = - ab->coeff(); rational alpha = - ab->coeff();
@ -513,8 +518,8 @@ nex * nla_grobner::expr_superpose(nex* e1, nex* e2, nex_mul* ab, nex_mul* ac, ne
// let eq1: ab+q=0, and eq2: ac+e=0, then qc - eb = 0 // let eq1: ab+q=0, and eq2: ac+e=0, then qc - eb = 0
void nla_grobner::superpose(equation * eq1, equation * eq2) { void nla_grobner::superpose(equation * eq1, equation * eq2) {
TRACE("grobner", tout << "eq1="; display_equation(tout, *eq1) << "eq2="; display_equation(tout, *eq2);); TRACE("grobner", tout << "eq1="; display_equation(tout, *eq1) << "eq2="; display_equation(tout, *eq2););
nex_mul * ab = get_highest_monomial(eq1->exp()); const nex * ab = get_highest_monomial(eq1->exp());
nex_mul * ac = get_highest_monomial(eq2->exp()); const nex * ac = get_highest_monomial(eq2->exp());
nex_mul *b, *c; nex_mul *b, *c;
TRACE("grobner", tout << "ab="; if(ab) { tout << *ab; } else { tout << "null"; }; TRACE("grobner", tout << "ab="; if(ab) { tout << *ab; } else { tout << "null"; };
tout << " , " << " ac="; if(ac) { tout << *ac; } else { tout << "null"; }; tout << "\n";); tout << " , " << " ac="; if(ac) { tout << *ac; } else { tout << "null"; }; tout << "\n";);
@ -525,66 +530,74 @@ void nla_grobner::superpose(equation * eq1, equation * eq2) {
init_equation(eq, expr_superpose( eq1->exp(), eq2->exp(), ab, ac, b, c), m_dep_manager.mk_join(eq1->dep(), eq2->dep())); init_equation(eq, expr_superpose( eq1->exp(), eq2->exp(), ab, ac, b, c), m_dep_manager.mk_join(eq1->dep(), eq2->dep()));
insert_to_simplify(eq); insert_to_simplify(eq);
} }
// Let a be the greatest common divider of ab and bc,
bool nla_grobner::find_b_c(nex_mul*ab, nex_mul* ac, nex_mul*& b, nex_mul*& c) { // then ab/a is stored in b, and ac/a is stored in c
bool nla_grobner::find_b_c(const nex* ab, const nex* ac, nex_mul*& b, nex_mul*& c) {
if (!find_b_c_check_only(ab, ac)) if (!find_b_c_check_only(ab, ac))
return false; return false;
b = m_nex_creator.mk_mul(); c = m_nex_creator.mk_mul(); b = m_nex_creator.mk_mul(); c = m_nex_creator.mk_mul();
nex_pow* bp = ab->begin(); unsigned ab_size = ab->number_of_child_powers();
nex_pow* cp = ac->begin(); unsigned ac_size = ac->number_of_child_powers();
unsigned i = 0, j = 0;
// nex_pow* bp = ab->begin();
// nex_pow* cp = ac->begin();
for (;;) { for (;;) {
if (m_nex_creator.lt(bp->e(), cp->e())) { const nex* m = ab->get_child_exp(i);
b->add_child_in_power(*bp); const nex* n = ac->get_child_exp(j);
if (++bp == ab->end()) if (m_nex_creator.lt(m, n)) {
b->add_child_in_power(const_cast<nex*>(m), ab->get_child_pow(i));
if (++i == ab_size)
break; break;
} else if (m_nex_creator.lt(cp->e(), bp->e())) { } else if (m_nex_creator.lt(n, m)) {
c->add_child_in_power(*cp); c->add_child_in_power(const_cast<nex*>(n), ac->get_child_pow(j));
if (++cp == ac->end()) if (++j == ac_size)
break; break;
} else { } else {
unsigned b_pow = bp->pow(); unsigned b_pow = ab->get_child_pow(i);
unsigned c_pow = cp->pow(); unsigned c_pow = ac->get_child_pow(j);
if (b_pow > c_pow) { if (b_pow > c_pow) {
b->add_child_in_power(bp->e(), b_pow - c_pow); b->add_child_in_power(const_cast<nex*>(m), b_pow - c_pow);
} else if (c_pow > b_pow) { } else if (c_pow > b_pow) {
c->add_child_in_power(cp->e(), c_pow - b_pow); c->add_child_in_power(const_cast<nex*>(n), c_pow - b_pow);
} // otherwise the power are equal and no child added to either b or c } // otherwise the power are equal and no child added to either b or c
bp++; cp++; i++; j++;
if (bp == ab->end() || cp == ac->end()) { if (i == ab_size || j == ac_size) {
break; break;
} }
} }
} }
while (cp != ac->end()) { while (i != ab_size) {
c->add_child_in_power(*cp); c->add_child_in_power(const_cast<nex*>(ab->get_child_exp(i)), ab->get_child_pow(i));
cp++; i++;
} }
while (bp != ab->end()) { while (j != ac_size) {
b->add_child_in_power(*bp); c->add_child_in_power(const_cast<nex*>(ac->get_child_exp(j)), ac->get_child_pow(j));
bp++; j++;
} }
TRACE("nla_grobner", tout << "b=" << *b << ", c=" <<*c << "\n";); TRACE("nla_grobner", tout << "b=" << *b << ", c=" <<*c << "\n";);
return true; return true;
} }
// Finds out if ab and bc have a non-trivial common divider
bool nla_grobner::find_b_c_check_only(const nex_mul*ab, const nex_mul* ac) const { bool nla_grobner::find_b_c_check_only(const nex* ab, const nex* ac) const {
if (ab == nullptr || ac == nullptr) if (ab == nullptr || ac == nullptr)
return false; return false;
SASSERT(m_nex_creator.is_simplified(ab) && m_nex_creator.is_simplified(ab)); SASSERT(m_nex_creator.is_simplified(ab) && m_nex_creator.is_simplified(ab));
unsigned i = 0, j = 0; // i points to ab, j points to ac unsigned i = 0, j = 0; // i points to ab, j points to ac
for (;;) { for (;;) {
if (m_nex_creator.lt((*ab)[i].e(), (*ac)[j].e())) { const nex* m = ab->get_child_exp(i);
const nex* n = ac->get_child_exp(j);
if (m_nex_creator.lt(m , n)) {
i++; i++;
if (i == ab->size()) if (i == ab->number_of_child_powers())
return false; return false;
} else if (m_nex_creator.lt((*ac)[j].e(), (*ab)[i].e())) { } else if (m_nex_creator.lt(n, m)) {
j++; j++;
if (j == ac->size()) if (j == ac->number_of_child_powers())
return false; return false;
} else { } else {
TRACE("grobner", tout << "found common " << *(*ac)[j].e() << " in " << *ab << " and " << *ac << "\n";); TRACE("grobner", tout << "found common " << *m << "\n";);
return true; return true;
} }
} }
@ -659,23 +672,16 @@ void nla_grobner::update_statistics(){
bool nla_grobner::find_conflict(ptr_vector<equation>& eqs){ bool nla_grobner::find_conflict(ptr_vector<equation>& eqs){
eqs.reset(); eqs.reset();
get_equations(eqs); get_equations(eqs);
TRACE("grobner_bug", tout << "after gb\n";);
for (equation* eq : eqs) { for (equation* eq : eqs) {
TRACE("grobner_bug", display_equation(tout, *eq);); TRACE("grobner_bug", display_equation(tout, *eq););
if (is_inconsistent(eq) || is_inconsistent2(eq)) if (is_inconsistent(eq))
return true; return true;
} }
return false; return false;
} }
bool nla_grobner::is_inconsistent(equation*) { bool nla_grobner::is_inconsistent(equation* e ) {
// NOT_IMPLEMENTED_YET(); todo implement return e->exp()->is_scalar() && (!to_scalar(e->exp())->value().is_zero());
return false;
}
bool nla_grobner::is_inconsistent2(equation*) {
// NOT_IMPLEMENTED_YET(); todo implement
return false;
} }
template <typename T, typename V> template <typename T, typename V>

21
src/math/lp/nla_grobner.h
View file

@ -107,7 +107,6 @@ private:
void update_statistics(); void update_statistics();
bool find_conflict(ptr_vector<equation>& eqs); bool find_conflict(ptr_vector<equation>& eqs);
bool is_inconsistent(equation*); bool is_inconsistent(equation*);
bool is_inconsistent2(equation*);
bool push_calculation_forward(ptr_vector<equation>& eqs, unsigned&); bool push_calculation_forward(ptr_vector<equation>& eqs, unsigned&);
void compute_basis_init(); void compute_basis_init();
bool compute_basis_loop(); bool compute_basis_loop();
@ -123,8 +122,8 @@ bool simplify_processed_with_eq(equation*);
bool canceled() { return false; } // todo, implement bool canceled() { return false; } // todo, implement
void superpose(equation * eq1, equation * eq2); void superpose(equation * eq1, equation * eq2);
void superpose(equation * eq); void superpose(equation * eq);
bool find_b_c(nex_mul*ab, nex_mul* ac, nex_mul*& b, nex_mul*& c); bool find_b_c(const nex *ab, const nex* ac, nex_mul*& b, nex_mul*& c);
bool find_b_c_check_only(const nex_mul*ab, const nex_mul* ac) const; bool find_b_c_check_only(const nex* ab, const nex* ac) const;
bool is_trivial(equation* ) const; bool is_trivial(equation* ) const;
bool is_better_choice(equation * eq1, equation * eq2); bool is_better_choice(equation * eq1, equation * eq2);
void del_equations(unsigned old_size); void del_equations(unsigned old_size);
@ -151,15 +150,15 @@ bool simplify_processed_with_eq(equation*);
m_to_superpose.insert(eq); m_to_superpose.insert(eq);
} }
void simplify_equations_to_process(); void simplify_equations_to_process();
nex_mul * get_highest_monomial(nex * e) const; const nex * get_highest_monomial(const nex * e) const;
ci_dependency* dep_from_vector( svector<lp::constraint_index> & fixed_vars_constraints); ci_dependency* dep_from_vector( svector<lp::constraint_index> & fixed_vars_constraints);
bool simplify_target_monomials_sum(equation *, equation *, nex_sum*, const nex_mul*); bool simplify_target_monomials_sum(equation *, equation *, nex_sum*, const nex*);
bool simplify_target_monomials_sum_check_only(nex_sum*, const nex_mul*); bool simplify_target_monomials_sum_check_only(nex_sum*, const nex*) const;
void simplify_target_monomials_sum_j(equation *, equation *, nex_sum*, const nex_mul*, unsigned); void simplify_target_monomials_sum_j(equation *, equation *, nex_sum*, const nex*, unsigned);
nex_mul * divide_ignore_coeffs(nex* ej, const nex_mul*); nex_mul * divide_ignore_coeffs(nex* ej, const nex*);
bool divide_ignore_coeffs_check_only(nex_mul* , const nex_mul*); bool divide_ignore_coeffs_check_only(nex_mul* , const nex*) const;
nex_mul * divide_ignore_coeffs_perform(nex_mul* , const nex_mul*); nex_mul * divide_ignore_coeffs_perform(nex_mul* , const nex*);
nex * expr_superpose(nex* e1, nex* e2, nex_mul* ab, nex_mul* ac, nex_mul* b, nex_mul* c); nex * expr_superpose(nex* e1, nex* e2, const nex* ab, const nex* ac, nex_mul* b, nex_mul* c);
void add_mul_skip_first(nex_sum* r, const rational& beta, nex *e, nex_mul* c); void add_mul_skip_first(nex_sum* r, const rational& beta, nex *e, nex_mul* c);
}; // end of grobner }; // end of grobner
} }