diff --git a/src/math/polynomial/algebraic_numbers.cpp b/src/math/polynomial/algebraic_numbers.cpp index 779a07d2d..b4763555f 100644 --- a/src/math/polynomial/algebraic_numbers.cpp +++ b/src/math/polynomial/algebraic_numbers.cpp @@ -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) { diff --git a/src/nlsat/nlsat_interval_set.cpp b/src/nlsat/nlsat_interval_set.cpp index 69b09e600..d07dd84a3 100644 --- a/src/nlsat/nlsat_interval_set.cpp +++ b/src/nlsat/nlsat_interval_set.cpp @@ -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()); diff --git a/src/util/mpbq.cpp b/src/util/mpbq.cpp index 3edbdab67..74b8deb49 100644 --- a/src/util/mpbq.cpp +++ b/src/util/mpbq.cpp @@ -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 << "" << a.m_k << ""; + 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) { diff --git a/src/util/mpbq.h b/src/util/mpbq.h index cfd6971bc..13cd6ae0c 100644 --- a/src/util/mpbq.h +++ b/src/util/mpbq.h @@ -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.