mirror of
https://github.com/Z3Prover/z3
synced 2025-04-27 19:05:51 +00:00
patch for Sturm sequence bug #4961
This commit is contained in:
Nikolaj Bjorner
parent
2d1684bc2d
commit
7d60d8462d
4 changed files with 124 additions and 114 deletions
|
|
@ -1844,79 +1844,97 @@ namespace algebraic_numbers {
|
|||
}
|
||||
}
|
||||
|
||||
// EXPENSIVE CASE
|
||||
// Let seq be the Sturm-Tarski sequence for
|
||||
// p_a, p_a' * p_b
|
||||
// Let s_l be the number of sign variations at a_lower.
|
||||
// Let s_u be the number of sign variations at a_upper.
|
||||
// By Sturm-Tarski Theorem, we have that
|
||||
// V = s_l - s_u = #(p_b(r) > 0) - #(p_b(r) < 0) at roots r of p_a
|
||||
// Since there is only one root of p_a in (a_lower, b_lower),
|
||||
// we are evaluating the sign of p_b at a.
|
||||
// That is V is the sign of p_b at a.
|
||||
//
|
||||
// We have
|
||||
// V < 0 -> p_b(a) < 0 -> if p_b(b_lower) < 0 then b > a else b < a
|
||||
// V == 0 -> p_b(a) == 0 -> a = b
|
||||
// V > 0 -> p_b(a) > 0 -> if p_b(b_lower) > 0 then b > a else b < a
|
||||
// Simplifying we have:
|
||||
// V == 0 --> a = b
|
||||
// if (V < 0) == (p_b(b_lower) < 0) then b > a else b < a
|
||||
//
|
||||
// workaround: Sturm sequences are buggy as exemplified by several open github issues
|
||||
// instead of relying on Sturm check if a simple interval expansion allows to separate
|
||||
// a and b.
|
||||
scoped_mpbq la(bqm()), ua(bqm());
|
||||
scoped_mpbq lb(bqm()), ub(bqm());
|
||||
unsigned precision = 10;
|
||||
if (get_interval(a, la, ua, precision) &&
|
||||
get_interval(b, lb, ub, precision)) {
|
||||
IF_VERBOSE(9, verbose_stream() << "sturm 0\n");
|
||||
if (la > ub)
|
||||
return sign_pos;
|
||||
if (ua < lb)
|
||||
return sign_neg;
|
||||
}
|
||||
IF_VERBOSE(9, verbose_stream() << "sturm 1\n");
|
||||
|
||||
m_compare_sturm++;
|
||||
upolynomial::scoped_upolynomial_sequence seq(upm());
|
||||
upm().sturm_tarski_seq(cell_a->m_p_sz, cell_a->m_p, cell_b->m_p_sz, cell_b->m_p, seq);
|
||||
unsigned V1 = upm().sign_variations_at(seq, a_lower);
|
||||
unsigned V2 = upm().sign_variations_at(seq, a_upper);
|
||||
int V = V1 - V2;
|
||||
TRACE("algebraic", tout << "comparing using sturm\n";
|
||||
display_interval(tout, a) << "\n";
|
||||
display_interval(tout, b) << "\n";
|
||||
tout << "V: " << V << " V1 " << V1 << " V2 " << V2
|
||||
<< ", sign_lower(a): " << sign_lower(cell_a)
|
||||
<< ", sign_lower(b): " << sign_lower(cell_b) << "\n";
|
||||
/*upm().display(tout << "sequence: ", seq);*/
|
||||
);
|
||||
if (V == 0)
|
||||
return sign_zero;
|
||||
if ((V < 0) == (sign_lower(cell_b) < 0))
|
||||
return sign_neg;
|
||||
else
|
||||
return sign_pos;
|
||||
//
|
||||
// EXPENSIVE CASE
|
||||
// Let seq be the Sturm-Tarski sequence for
|
||||
// p_a, p_a' * p_b
|
||||
// Let s_l be the number of sign variations at a_lower.
|
||||
// Let s_u be the number of sign variations at a_upper.
|
||||
// By Sturm-Tarski Theorem, we have that
|
||||
// V = s_l - s_u = #(p_b(r) > 0) - #(p_b(r) < 0) at roots r of p_a
|
||||
// Since there is only one root of p_a in (a_lower, b_lower),
|
||||
// we are evaluating the sign of p_b at a.
|
||||
// That is V is the sign of p_b at a.
|
||||
//
|
||||
// We have
|
||||
// V < 0 -> p_b(a) < 0 -> if p_b(b_lower) < 0 then b > a else b < a
|
||||
// V == 0 -> p_b(a) == 0 -> a = b
|
||||
// V > 0 -> p_b(a) > 0 -> if p_b(b_lower) > 0 then b > a else b < a
|
||||
// Simplifying we have:
|
||||
// V == 0 --> a = b
|
||||
// if (V < 0) == (p_b(b_lower) < 0) then b > a else b < a
|
||||
//
|
||||
|
||||
// Here is an unexplored option for comparing numbers.
|
||||
//
|
||||
// The isolating intervals of a and b are still overlapping
|
||||
// Then we compute
|
||||
// r(x) = Resultant(x - y1 + y2, p1(y1), p2(y2))
|
||||
// where p1(y1) and p2(y2) are the polynomials defining a and b.
|
||||
// Remarks:
|
||||
// 1) The resultant r(x) must not be the zero polynomial,
|
||||
// since the polynomial x - y1 + y2 does not vanish in any of the roots of p1(y1) and p2(y2)
|
||||
//
|
||||
// 2) By resultant calculus, If alpha, beta1, beta2 are roots of x - y1 + y2, p1(y1), p2(y2)
|
||||
// then alpha is a root of r(x).
|
||||
// Thus, we have that a - b is a root of r(x)
|
||||
//
|
||||
// 3) If 0 is not a root of r(x), then a != b (by remark 2)
|
||||
//
|
||||
// 4) Let (l1, u1) and (l2, u2) be the intervals of a and b.
|
||||
// Then, a - b must be in (l1 - u2, u1 - l2)
|
||||
//
|
||||
// 5) Assume a != b, then if we keep refining the isolating intervals for a and b,
|
||||
// then eventually, (l1, u1) and (l2, u2) will not overlap.
|
||||
// Thus, if 0 is not a root of r(x), we can keep refining until
|
||||
// the intervals do not overlap.
|
||||
//
|
||||
// 6) If 0 is a root of r(x), we have two possibilities:
|
||||
// a) Isolate roots of r(x) in the interval (l1 - u2, u1 - l2),
|
||||
// and then keep refining (l1, u1) and (l2, u2) until they
|
||||
// (l1 - u2, u1 - l2) "convers" only one root.
|
||||
//
|
||||
// b) Compute the sturm sequence for r(x),
|
||||
// keep refining the (l1, u1) and (l2, u2) until
|
||||
// (l1 - u2, u1 - l2) contains only one root of r(x)
|
||||
|
||||
m_compare_sturm++;
|
||||
upolynomial::scoped_upolynomial_sequence seq(upm());
|
||||
upm().sturm_tarski_seq(cell_a->m_p_sz, cell_a->m_p, cell_b->m_p_sz, cell_b->m_p, seq);
|
||||
unsigned V1 = upm().sign_variations_at(seq, a_lower);
|
||||
unsigned V2 = upm().sign_variations_at(seq, a_upper);
|
||||
int V = V1 - V2;
|
||||
TRACE("algebraic", tout << "comparing using sturm\n";
|
||||
display_interval(tout, a) << "\n";
|
||||
display_interval(tout, b) << "\n";
|
||||
tout << "V: " << V << " V1 " << V1 << " V2 " << V2
|
||||
<< ", sign_lower(a): " << sign_lower(cell_a)
|
||||
<< ", sign_lower(b): " << sign_lower(cell_b) << "\n";
|
||||
/*upm().display(tout << "sequence: ", seq);*/
|
||||
);
|
||||
if (V == 0)
|
||||
return sign_zero;
|
||||
if ((V < 0) == (sign_lower(cell_b) < 0))
|
||||
return sign_neg;
|
||||
else
|
||||
return sign_pos;
|
||||
|
||||
// Here is an unexplored option for comparing numbers.
|
||||
//
|
||||
// The isolating intervals of a and b are still overlapping
|
||||
// Then we compute
|
||||
// r(x) = Resultant(x - y1 + y2, p1(y1), p2(y2))
|
||||
// where p1(y1) and p2(y2) are the polynomials defining a and b.
|
||||
// Remarks:
|
||||
// 1) The resultant r(x) must not be the zero polynomial,
|
||||
// since the polynomial x - y1 + y2 does not vanish in any of the roots of p1(y1) and p2(y2)
|
||||
//
|
||||
// 2) By resultant calculus, If alpha, beta1, beta2 are roots of x - y1 + y2, p1(y1), p2(y2)
|
||||
// then alpha is a root of r(x).
|
||||
// Thus, we have that a - b is a root of r(x)
|
||||
//
|
||||
// 3) If 0 is not a root of r(x), then a != b (by remark 2)
|
||||
//
|
||||
// 4) Let (l1, u1) and (l2, u2) be the intervals of a and b.
|
||||
// Then, a - b must be in (l1 - u2, u1 - l2)
|
||||
//
|
||||
// 5) Assume a != b, then if we keep refining the isolating intervals for a and b,
|
||||
// then eventually, (l1, u1) and (l2, u2) will not overlap.
|
||||
// Thus, if 0 is not a root of r(x), we can keep refining until
|
||||
// the intervals do not overlap.
|
||||
//
|
||||
// 6) If 0 is a root of r(x), we have two possibilities:
|
||||
// a) Isolate roots of r(x) in the interval (l1 - u2, u1 - l2),
|
||||
// and then keep refining (l1, u1) and (l2, u2) until they
|
||||
// (l1 - u2, u1 - l2) "convers" only one root.
|
||||
//
|
||||
// b) Compute the sturm sequence for r(x),
|
||||
// keep refining the (l1, u1) and (l2, u2) until
|
||||
// (l1 - u2, u1 - l2) contains only one root of r(x)
|
||||
}
|
||||
|
||||
::sign compare(numeral & a, numeral & b) {
|
||||
|
|
|
|||
|
|
@ -71,20 +71,15 @@ namespace nlsat {
|
|||
}
|
||||
|
||||
bool check_interval(anum_manager & am, interval const & i) {
|
||||
if (i.m_lower_inf) {
|
||||
SASSERT(i.m_lower_open);
|
||||
}
|
||||
if (i.m_upper_inf) {
|
||||
SASSERT(i.m_upper_open);
|
||||
}
|
||||
SASSERT(!i.m_lower_inf || i.m_lower_open);
|
||||
SASSERT(!i.m_upper_inf || i.m_upper_open);
|
||||
|
||||
if (!i.m_lower_inf && !i.m_upper_inf) {
|
||||
auto s = am.compare(i.m_lower, i.m_upper);
|
||||
TRACE("nlsat_interval", tout << "lower: "; am.display_decimal(tout, i.m_lower); tout << ", upper: "; am.display_decimal(tout, i.m_upper);
|
||||
tout << "\ns: " << s << "\n";);
|
||||
SASSERT(s <= 0);
|
||||
if (::is_zero(s)) {
|
||||
SASSERT(!i.m_lower_open && !i.m_upper_open);
|
||||
}
|
||||
SASSERT(!is_zero(s) || !i.m_lower_open && !i.m_upper_open);
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
|
@ -92,12 +87,10 @@ namespace nlsat {
|
|||
bool check_no_overlap(anum_manager & am, interval const & curr, interval const & next) {
|
||||
SASSERT(!curr.m_upper_inf);
|
||||
SASSERT(!next.m_lower_inf);
|
||||
int sign = am.compare(curr.m_upper, next.m_lower);
|
||||
CTRACE("nlsat", sign > 0, display(tout, am, curr); tout << " "; display(tout, am, next); tout << "\n";);
|
||||
SASSERT(sign <= 0);
|
||||
if (sign == 0) {
|
||||
SASSERT(curr.m_upper_open || next.m_lower_open);
|
||||
}
|
||||
sign s = am.compare(curr.m_upper, next.m_lower);
|
||||
CTRACE("nlsat", s > 0, display(tout, am, curr); tout << " "; display(tout, am, next); tout << "\n";);
|
||||
SASSERT(s <= 0);
|
||||
SASSERT(!is_zero(s) || curr.m_upper_open || next.m_lower_open);
|
||||
return true;
|
||||
}
|
||||
|
||||
|
|
@ -107,9 +100,7 @@ namespace nlsat {
|
|||
for (unsigned i = 0; i < sz; i++) {
|
||||
interval const & curr = ints[i];
|
||||
SASSERT(check_interval(am, curr));
|
||||
if (i < sz - 1) {
|
||||
SASSERT(check_no_overlap(am, curr, ints[i+1]));
|
||||
}
|
||||
SASSERT(i >= sz - 1 || check_no_overlap(am, curr, ints[i+1]));
|
||||
});
|
||||
return true;
|
||||
}
|
||||
|
|
@ -276,10 +267,9 @@ namespace nlsat {
|
|||
|
||||
interval_set * interval_set_manager::mk_union(interval_set const * s1, interval_set const * s2) {
|
||||
#if 0
|
||||
// issue #2867:
|
||||
static unsigned s_count = 0;
|
||||
s_count++;
|
||||
if (s_count == 8442) {
|
||||
if (s_count == 4470) {
|
||||
enable_trace("nlsat_interval");
|
||||
enable_trace("algebraic");
|
||||
}
|
||||
|
|
@ -416,7 +406,7 @@ namespace nlsat {
|
|||
else {
|
||||
SASSERT(l1_l2_sign > 0);
|
||||
if (u1_u2_sign == 0) {
|
||||
TRACE("nlsat_interval", tout << "l1_l2_sign > 0, u1_u2_sign == 0\n";);
|
||||
TRACE("nlsat_interval", tout << "l2 < l1 <= u1 = u2\n";);
|
||||
// Case:
|
||||
// 1) [ ]
|
||||
// [ ]
|
||||
|
|
@ -426,7 +416,7 @@ namespace nlsat {
|
|||
i2++;
|
||||
}
|
||||
else if (u1_u2_sign < 0) {
|
||||
TRACE("nlsat_interval", tout << "l1_l2_sign > 0, u1_u2_sign > 0\n";);
|
||||
TRACE("nlsat_interval", tout << "l2 < l1 <= u2 < u2\n";);
|
||||
// Case:
|
||||
// 1) [ ]
|
||||
// [ ]
|
||||
|
|
@ -436,7 +426,7 @@ namespace nlsat {
|
|||
else {
|
||||
auto u2_l1_sign = compare_upper_lower(m_am, int2, int1);
|
||||
if (u2_l1_sign < 0) {
|
||||
TRACE("nlsat_interval", tout << "l1_l2_sign > 0, u1_u2_sign > 0, u2_l1_sign < 0\n";);
|
||||
TRACE("nlsat_interval", tout << "l2 <= u2 < l1 <= u1\n";);
|
||||
// Case:
|
||||
// 1) [ ]
|
||||
// [ ]
|
||||
|
|
@ -457,7 +447,7 @@ namespace nlsat {
|
|||
i2++;
|
||||
}
|
||||
else {
|
||||
TRACE("nlsat_interval", tout << "l1_l2_sign > 0, u1_u2_sign > 0, u2_l1_sign > 0\n";);
|
||||
TRACE("nlsat_interval", tout << "l2 < l1 < u2 < u1\n";);
|
||||
SASSERT(l1_l2_sign > 0);
|
||||
SASSERT(u1_u2_sign > 0);
|
||||
SASSERT(u2_l1_sign > 0);
|
||||
|
|
@ -474,6 +464,7 @@ namespace nlsat {
|
|||
}
|
||||
SASSERT(result.size() <= 1 ||
|
||||
check_no_overlap(m_am, result[result.size() - 2], result[result.size() - 1]));
|
||||
|
||||
}
|
||||
|
||||
SASSERT(!result.empty());
|
||||
|
|
|
|||
|
|
@ -57,9 +57,8 @@ mpbq_manager::~mpbq_manager() {
|
|||
}
|
||||
|
||||
void mpbq_manager::reset(mpbq_vector & v) {
|
||||
unsigned sz = v.size();
|
||||
for (unsigned i = 0; i < sz; i++)
|
||||
reset(v[i]);
|
||||
for (auto & e : v)
|
||||
reset(e);
|
||||
v.reset();
|
||||
}
|
||||
|
||||
|
|
@ -391,23 +390,25 @@ std::string mpbq_manager::to_string(mpbq const & a) {
|
|||
return buffer.str();
|
||||
}
|
||||
|
||||
void mpbq_manager::display(std::ostream & out, mpbq const & a) {
|
||||
std::ostream& mpbq_manager::display(std::ostream & out, mpbq const & a) {
|
||||
out << m_manager.to_string(a.m_num);
|
||||
if (a.m_k > 0)
|
||||
out << "/2";
|
||||
if (a.m_k > 1)
|
||||
out << "^" << a.m_k;
|
||||
return out;
|
||||
}
|
||||
|
||||
void mpbq_manager::display_pp(std::ostream & out, mpbq const & a) {
|
||||
std::ostream& mpbq_manager::display_pp(std::ostream & out, mpbq const & a) {
|
||||
out << m_manager.to_string(a.m_num);
|
||||
if (a.m_k > 0)
|
||||
out << "/2";
|
||||
if (a.m_k > 1)
|
||||
out << "<sup>" << a.m_k << "</sup>";
|
||||
return out;
|
||||
}
|
||||
|
||||
void mpbq_manager::display_smt2(std::ostream & out, mpbq const & a, bool decimal) {
|
||||
std::ostream& mpbq_manager::display_smt2(std::ostream & out, mpbq const & a, bool decimal) {
|
||||
if (a.m_k == 0) {
|
||||
m_manager.display_smt2(out, a.m_num, decimal);
|
||||
}
|
||||
|
|
@ -421,12 +422,12 @@ void mpbq_manager::display_smt2(std::ostream & out, mpbq const & a, bool decimal
|
|||
if (decimal) out << ".0";
|
||||
out << "))";
|
||||
}
|
||||
return out;
|
||||
}
|
||||
|
||||
void mpbq_manager::display_decimal(std::ostream & out, mpbq const & a, unsigned prec) {
|
||||
std::ostream& mpbq_manager::display_decimal(std::ostream & out, mpbq const & a, unsigned prec) {
|
||||
if (is_int(a)) {
|
||||
out << m_manager.to_string(a.m_num);
|
||||
return;
|
||||
return out << m_manager.to_string(a.m_num);
|
||||
}
|
||||
mpz two(2);
|
||||
mpz ten(10);
|
||||
|
|
@ -455,17 +456,16 @@ void mpbq_manager::display_decimal(std::ostream & out, mpbq const & a, unsigned
|
|||
m_manager.del(n1);
|
||||
m_manager.del(v1);
|
||||
m_manager.del(two_k);
|
||||
return out;
|
||||
}
|
||||
|
||||
void mpbq_manager::display_decimal(std::ostream & out, mpbq const & a, mpbq const & b, unsigned prec) {
|
||||
std::ostream& mpbq_manager::display_decimal(std::ostream & out, mpbq const & a, mpbq const & b, unsigned prec) {
|
||||
mpz two(2);
|
||||
mpz ten(10);
|
||||
mpz two_k1, two_k2;
|
||||
mpz n1, v1, n2, v2;
|
||||
if (m_manager.is_neg(a.m_num) != m_manager.is_neg(b.m_num)) {
|
||||
out << "?";
|
||||
return;
|
||||
}
|
||||
if (m_manager.is_neg(a.m_num) != m_manager.is_neg(b.m_num))
|
||||
return out << "?";
|
||||
if (m_manager.is_neg(a.m_num))
|
||||
out << "-";
|
||||
m_manager.set(v1, a.m_num);
|
||||
|
|
@ -512,6 +512,7 @@ void mpbq_manager::display_decimal(std::ostream & out, mpbq const & a, mpbq cons
|
|||
m_manager.del(v2);
|
||||
m_manager.del(two_k1);
|
||||
m_manager.del(two_k2);
|
||||
return out;
|
||||
}
|
||||
|
||||
bool mpbq_manager::to_mpbq(mpq const & q, mpbq & bq) {
|
||||
|
|
|
|||
|
|
@ -258,17 +258,17 @@ public:
|
|||
void select_small_core(unsynch_mpq_manager & qm, mpq const & lower, mpq const & upper, mpbq & r);
|
||||
|
||||
|
||||
void display(std::ostream & out, mpbq const & a);
|
||||
void display_pp(std::ostream & out, mpbq const & a);
|
||||
void display_decimal(std::ostream & out, mpbq const & a, unsigned prec = 8);
|
||||
std::ostream& display(std::ostream & out, mpbq const & a);
|
||||
std::ostream& display_pp(std::ostream & out, mpbq const & a);
|
||||
std::ostream& display_decimal(std::ostream & out, mpbq const & a, unsigned prec = 8);
|
||||
/**
|
||||
\brief Display a in decimal while its digits match b digits.
|
||||
This function is useful when a and b are representing an interval [a,b] which
|
||||
contains an algebraic number
|
||||
*/
|
||||
void display_decimal(std::ostream & out, mpbq const & a, mpbq const & b, unsigned prec);
|
||||
std::ostream& display_decimal(std::ostream & out, mpbq const & a, mpbq const & b, unsigned prec);
|
||||
|
||||
void display_smt2(std::ostream & out, mpbq const & a, bool decimal);
|
||||
std::ostream& display_smt2(std::ostream & out, mpbq const & a, bool decimal);
|
||||
|
||||
/**
|
||||
\brief Approximate n as b/2^k' s.t. k' <= k.
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue