mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-15 13:28:47 +00:00
Code Activity

Add new rational function normalization procedure.

Browse source 
Signed-off-by: Leonardo de Moura <leonardo@microsoft.com>
This commit is contained in:
Leonardo de Moura 2013-01-15 08:39:30 -08:00
parent f0737bdf7f
commit 217c8375ce
1 changed files with 82 additions and 57 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

139
src/math/realclosure/realclosure.cpp
View file

@ -3450,6 +3450,10 @@ namespace realclosure {
// //
// --------------------------------- // ---------------------------------
bool is_monic(value_ref_buffer const & p) {
return p.size() > 0 && is_rational_one(p[p.size() - 1]);
}
/** /**
\brief Force the leading coefficient of p to be 1. \brief Force the leading coefficient of p to be 1.
*/ */
@ -4762,14 +4766,20 @@ namespace realclosure {
bool determine_sign(rational_function_value * v) { bool determine_sign(rational_function_value * v) {
if (!contains_zero(v->interval())) if (!contains_zero(v->interval()))
return true; return true;
bool r;
switch (v->ext()->knd()) { switch (v->ext()->knd()) {
case extension::TRANSCENDENTAL: determine_transcendental_sign(v); return true; // it is never zero case extension::TRANSCENDENTAL: determine_transcendental_sign(v); r = true; break; // it is never zero
case extension::INFINITESIMAL: determine_infinitesimal_sign(v); return true; // it is never zero case extension::INFINITESIMAL: determine_infinitesimal_sign(v); r = true; break; // it is never zero
case extension::ALGEBRAIC: return determine_algebraic_sign(v); case extension::ALGEBRAIC: r = determine_algebraic_sign(v); break;
default: default:
UNREACHABLE(); UNREACHABLE();
return false; r = false;
} }
TRACE("rcf_determine_sign_bug",
tout << "result: " << r << "\n";
display_compact(tout, v); tout << "\n";
tout << "sign: " << sign(v) << "\n";);
return r;
} }
bool determine_sign(value_ref & r) { bool determine_sign(value_ref & r) {
@ -4784,10 +4794,10 @@ namespace realclosure {
// --------------------------------- // ---------------------------------
/** /**
\brief Set new_p1 and new_p2 using the following normalization rules: \brief Compute polynomials new_p1 and new_p2 s.t.
- new_p1 <- p1/p2[0]; new_p2 <- one IF sz2 == 1 - p1/p2 == new_p1/new_p2, AND
- new_p1 <- one; new_p2 <- p2/p1[0]; IF sz1 == 1 - new_p2 is a Monic polynomial, AND
- new_p1 <- p1/gcd(p1, p2); new_p2 <- p2/gcd(p1, p2); Otherwise - gcd(new_p1, new_p2) == 1
*/ */
void normalize_fraction(unsigned sz1, value * const * p1, unsigned sz2, value * const * p2, value_ref_buffer & new_p1, value_ref_buffer & new_p2) { void normalize_fraction(unsigned sz1, value * const * p1, unsigned sz2, value * const * p2, value_ref_buffer & new_p1, value_ref_buffer & new_p2) {
INC_DEPTH(); INC_DEPTH();
@ -4800,47 +4810,58 @@ namespace realclosure {
div(sz1, p1, p2[0], new_p1); div(sz1, p1, p2[0], new_p1);
new_p2.reset(); new_p2.push_back(one()); new_p2.reset(); new_p2.push_back(one());
} }
else if (sz1 == 1) {
SASSERT(sz2 > 1);
// - new_p1 <- one; new_p2 <- p2/p1[0]; IF sz1 == 1
new_p1.reset(); new_p1.push_back(one());
div(sz2, p2, p1[0], new_p2);
}
else { else {
// - new_p1 <- p1/gcd(p1, p2); new_p2 <- p2/gcd(p1, p2); Otherwise value * lc = p2[sz2 - 1];
value_ref_buffer g(*this); if (is_rational_one(lc)) {
gcd(sz1, p1, sz2, p2, g); // p2 is monic
if (is_rational_one(g)) { normalize_num_monic_den(sz1, p1, sz2, p2, new_p1, new_p2);
new_p1.append(sz1, p1);
new_p2.append(sz2, p2);
}
else if (g.size() == sz1 || g.size() == sz2) {
// After dividing p1 and p2 by g, one of the quotients will have size 1.
// Thus, we have to apply the first two rules again.
value_ref_buffer tmp_p1(*this);
value_ref_buffer tmp_p2(*this);
div(sz1, p1, g.size(), g.c_ptr(), tmp_p1);
div(sz2, p2, g.size(), g.c_ptr(), tmp_p2);
if (tmp_p2.size() == 1) {
div(tmp_p1.size(), tmp_p1.c_ptr(), tmp_p2[0], new_p1);
new_p2.reset(); new_p2.push_back(one());
}
else if (tmp_p1.size() == 1) {
SASSERT(tmp_p2.size() > 1);
new_p1.reset(); new_p1.push_back(one());
div(tmp_p2.size(), tmp_p2.c_ptr(), tmp_p1[0], new_p2);
}
else {
UNREACHABLE();
}
} }
else { else {
div(sz1, p1, g.size(), g.c_ptr(), new_p1); // p2 is not monic
div(sz2, p2, g.size(), g.c_ptr(), new_p2); value_ref_buffer tmp1(*this);
SASSERT(new_p1.size() > 1); value_ref_buffer tmp2(*this);
SASSERT(new_p2.size() > 1); div(sz1, p1, lc, tmp1);
div(sz2, p2, lc, tmp2);
normalize_num_monic_den(tmp1.size(), tmp1.c_ptr(), tmp2.size(), tmp2.c_ptr(), new_p1, new_p2);
} }
} }
TRACE("normalize_fraction_bug",
display_poly(tout, sz1, p1); tout << "\n";
display_poly(tout, sz2, p2); tout << "\n";
tout << "====>\n";
display_poly(tout, new_p1.size(), new_p1.c_ptr()); tout << "\n";
display_poly(tout, new_p2.size(), new_p2.c_ptr()); tout << "\n";);
}
/**
\brief Auxiliary function for normalize_fraction.
It produces new_p1 and new_p2 s.t.
new_p1/new_p2 == p1/p2
gcd(new_p1, new_p2) == 1
Assumptions:
\pre p2 is monic
\pre sz2 > 1
*/
void normalize_num_monic_den(unsigned sz1, value * const * p1, unsigned sz2, value * const * p2,
value_ref_buffer & new_p1, value_ref_buffer & new_p2) {
SASSERT(sz2 > 1);
SASSERT(is_rational_one(p2[sz2-1]));
value_ref_buffer g(*this);
gcd(sz1, p1, sz2, p2, g);
SASSERT(is_monic(g));
if (is_rational_one(g)) {
new_p1.append(sz1, p1);
new_p2.append(sz2, p2);
}
else {
div(sz1, p1, g.size(), g.c_ptr(), new_p1);
div(sz2, p2, g.size(), g.c_ptr(), new_p2);
SASSERT(is_monic(new_p2));
}
} }
/** /**
@ -4925,10 +4946,8 @@ namespace realclosure {
value_ref_buffer new_num(*this); value_ref_buffer new_num(*this);
value_ref_buffer new_den(*this); value_ref_buffer new_den(*this);
normalize_fraction(num.size(), num.c_ptr(), ad.size(), ad.c_ptr(), new_num, new_den); normalize_fraction(num.size(), num.c_ptr(), ad.size(), ad.c_ptr(), new_num, new_den);
if (new_num.empty()) SASSERT(!new_num.empty());
r = 0; mk_add_value(a, b, new_num.size(), new_num.c_ptr(), new_den.size(), new_den.c_ptr(), r);
else
mk_add_value(a, b, new_num.size(), new_num.c_ptr(), new_den.size(), new_den.c_ptr(), r);
} }
} }
} }
@ -4986,10 +5005,8 @@ namespace realclosure {
value_ref_buffer new_num(*this); value_ref_buffer new_num(*this);
value_ref_buffer new_den(*this); value_ref_buffer new_den(*this);
normalize_fraction(num.size(), num.c_ptr(), den.size(), den.c_ptr(), new_num, new_den); normalize_fraction(num.size(), num.c_ptr(), den.size(), den.c_ptr(), new_num, new_den);
if (new_num.empty()) SASSERT(!new_num.empty());
r = 0; mk_add_value(a, b, new_num.size(), new_num.c_ptr(), new_den.size(), new_den.c_ptr(), r);
else
mk_add_value(a, b, new_num.size(), new_num.c_ptr(), new_den.size(), new_den.c_ptr(), r);
} }
} }
} }
@ -5334,7 +5351,11 @@ namespace realclosure {
polynomial const & ad = a->den(); polynomial const & ad = a->den();
scoped_mpbqi ri(bqim()); scoped_mpbqi ri(bqim());
bqim().inv(interval(a), ri); bqim().inv(interval(a), ri);
r = mk_rational_function_value_core(a->ext(), ad.size(), ad.c_ptr(), an.size(), an.c_ptr()); // The GCD of an and ad is one, we may use a simpler version of normalize
value_ref_buffer new_num(*this);
value_ref_buffer new_den(*this);
normalize_fraction(ad.size(), ad.c_ptr(), an.size(), an.c_ptr(), new_num, new_den);
r = mk_rational_function_value_core(a->ext(), new_num.size(), new_num.c_ptr(), new_den.size(), new_den.c_ptr());
swap(r->interval(), ri); swap(r->interval(), ri);
SASSERT(!contains_zero(r->interval())); SASSERT(!contains_zero(r->interval()));
} }
@ -5712,16 +5733,16 @@ namespace realclosure {
} }
} }
void display_compact(std::ostream & out, numeral const & a) const { void display_compact(std::ostream & out, value * a) const {
collect_algebraic_refs c; collect_algebraic_refs c;
c.mark(a.m_value); c.mark(a);
if (c.m_found.empty()) { if (c.m_found.empty()) {
display(out, a.m_value, true); display(out, a, true);
} }
else { else {
std::sort(c.m_found.begin(), c.m_found.end(), rank_lt_proc()); std::sort(c.m_found.begin(), c.m_found.end(), rank_lt_proc());
out << "["; out << "[";
display(out, a.m_value, true); display(out, a, true);
for (unsigned i = 0; i < c.m_found.size(); i++) { for (unsigned i = 0; i < c.m_found.size(); i++) {
algebraic * ext = c.m_found[i]; algebraic * ext = c.m_found[i];
out << "; r!" << ext->idx() << " := "; out << "; r!" << ext->idx() << " := ";
@ -5731,6 +5752,10 @@ namespace realclosure {
} }
} }
void display_compact(std::ostream & out, numeral const & a) const {
display_compact(out, a.m_value);
}
void display(std::ostream & out, numeral const & a, bool compact=false) const { void display(std::ostream & out, numeral const & a, bool compact=false) const {
if (compact) if (compact)
display_compact(out, a); display_compact(out, a);