From b662bc8dc7996ef47ac30bf52f395c0d45f812f0 Mon Sep 17 00:00:00 2001 From: Leonardo de Moura Date: Wed, 9 Jan 2013 11:16:04 -0800 Subject: [PATCH] Add lower and upper bounds for negative and positive roots Signed-off-by: Leonardo de Moura --- src/math/realclosure/realclosure.cpp | 140 +++++++++++++++++++-------- src/test/rcf.cpp | 2 +- 2 files changed, 98 insertions(+), 44 deletions(-) diff --git a/src/math/realclosure/realclosure.cpp b/src/math/realclosure/realclosure.cpp index 434f933b6..200a93f58 100644 --- a/src/math/realclosure/realclosure.cpp +++ b/src/math/realclosure/realclosure.cpp @@ -1015,7 +1015,7 @@ namespace realclosure { t->m_k++; t->m_proc(t->m_k, qim(), i); int m = magnitude(i); - TRACE("rcf", + TRACE("rcf_transcendental", tout << "refine_transcendental_interval rational: " << m << "\nrational interval: "; qim().display(tout, i); tout << std::endl;); unsigned k; @@ -1033,7 +1033,7 @@ namespace realclosure { void refine_transcendental_interval(transcendental * t, unsigned prec) { while (!check_precision(t->interval(), prec)) { - TRACE("rcf", tout << "refine_transcendental_interval: " << magnitude(t->interval()) << std::endl;); + TRACE("rcf_transcendental", tout << "refine_transcendental_interval: " << magnitude(t->interval()) << std::endl;); checkpoint(); refine_transcendental_interval(t); } @@ -1165,7 +1165,7 @@ namespace realclosure { Similarly, the procedure will fail if one of the a_{n-k} has an interval of the form (l, oo) or (-oo, u). Both cases can only happen if the values of the coefficients depend on infinitesimal values. */ - bool pos_root_upper_bound(unsigned n, value * const * p, unsigned & N) { + bool pos_root_upper_bound(unsigned n, value * const * p, int & N) { SASSERT(n > 1); SASSERT(!is_zero(p[n-1])); int lc_sign = sign(p[n-1]); @@ -1173,15 +1173,15 @@ namespace realclosure { int lc_mag; if (!abs_lower_magnitude(interval(p[n-1]), lc_mag)) return false; - N = 0; + N = -static_cast(m_ini_precision); for (unsigned k = 2; k <= n; k++) { value * a = p[n - k]; if (!is_zero(a) && sign(a) != lc_sign) { int a_mag; if (!abs_upper_magnitude(interval(a), a_mag)) return false; - int C = (a_mag - lc_mag)/k + 1 /* compensate imprecision on division */ + 1 /* Knuth's bound is 2 x Max {... } */; - if (C > 0 && static_cast(C) > N) + int C = (a_mag - lc_mag)/static_cast(k) + 1 /* compensate imprecision on division */ + 1 /* Knuth's bound is 2 x Max {... } */; + if (C > N) N = C; } } @@ -1204,13 +1204,13 @@ namespace realclosure { } /** - \brief Find negative root upper bound using Knuth's approach. + \brief Find negative root lower bound using Knuth's approach. This is similar to pos_root_upper_bound. In principle, we can use the same algorithm. We just have to adjust the coefficients by using the transformation p(-x). */ - bool neg_root_upper_bound(unsigned n, value * const * as, unsigned & N) { + bool neg_root_lower_bound(unsigned n, value * const * as, int & N) { SASSERT(n > 1); SASSERT(!is_zero(as[n-1])); scoped_mpbqi aux(bqim()); @@ -1219,7 +1219,7 @@ namespace realclosure { int lc_mag; if (!abs_lower_magnitude(aux, lc_mag)) return false; - N = 0; + N = -static_cast(m_ini_precision); for (unsigned k = 2; k <= n; k++) { value * a = as[n - k]; if (!is_zero(a)) { @@ -1229,8 +1229,8 @@ namespace realclosure { int a_mag; if (!abs_upper_magnitude(aux, a_mag)) return false; - int C = (a_mag - lc_mag)/k + 1 /* compensate imprecision on division */ + 1 /* Knuth's bound is 2 x Max {... } */; - if (C > 0 && static_cast(C) > N) + int C = (a_mag - lc_mag)/static_cast(k) + 1 /* compensate imprecision on division */ + 1 /* Knuth's bound is 2 x Max {... } */; + if (C > N) N = C; } } @@ -1238,6 +1238,53 @@ namespace realclosure { return true; } + /** + \brief q <- x^{n-1}*p(1/x) + + Given p(x) a_{n-1} * x^{n-1} + ... + a_0, this method stores + a_0 * x^{n-1} + ... + a_{n-1} into q. + */ + void reverse(unsigned n, value * const * p, value_ref_buffer & q) { + unsigned i = n; + while (i > 0) { + --i; + q.push_back(p[i]); + } + } + + /** + \brief To compute the lower bound for positive roots we computer the upper bound for the polynomial q(x) = x^{n-1}*p(1/x). + Assume U is an upper bound for roots of q(x), i.e., (r > 0 and q(r) = 0) implies r < U. + Note that if r is a root for q(x), then 1/r is a root for p(x) and 1/U is a lower bound for positive roots of p(x). + The polynomial q(x) is just p(x) "reversed". + */ + bool pos_root_lower_bound(unsigned n, value * const * p, int & N) { + value_ref_buffer q(*this); + reverse(n, p, q); + if (pos_root_upper_bound(n, q.c_ptr(), N)) { + N = -N; + return true; + } + else { + return false; + } + } + + /** + \brief See comment on pos_root_lower_bound. + */ + bool neg_root_upper_bound(unsigned n, value * const * p, int & N) { + value_ref_buffer q(*this); + reverse(n, p, q); + if (pos_root_lower_bound(n, q.c_ptr(), N)) { + N = -N; + return true; + } + else { + return false; + } + } + /** \brief Root isolation for polynomials that are - nonlinear (degree > 2) @@ -1245,18 +1292,25 @@ namespace realclosure { - square free - nonconstant */ - void nl_nz_sqf_isolate_roots(unsigned n, value * const * as, numeral_vector & roots) { + void nl_nz_sqf_isolate_roots(unsigned n, value * const * p, numeral_vector & roots) { SASSERT(n > 2); - SASSERT(!is_zero(as[0])); - SASSERT(!is_zero(as[n-1])); - unsigned pos_N, neg_N; - bool has_neg_upper = neg_root_upper_bound(n, as, neg_N); - bool has_pos_upper = pos_root_upper_bound(n, as, pos_N); - TRACE("rcf", - display_poly(tout, n, as); tout << "\n"; - tout << "has_neg_upper: " << has_neg_upper << " neg_N: " << neg_N << "\n"; - tout << "has_pos_upper: " << has_pos_upper << " pos_N: " << pos_N << "\n";); - + SASSERT(!is_zero(p[0])); + SASSERT(!is_zero(p[n-1])); + int pos_lower_N, pos_upper_N, neg_lower_N, neg_upper_N; + bool has_neg_lower = neg_root_lower_bound(n, p, neg_lower_N); + bool has_neg_upper = neg_root_upper_bound(n, p, neg_upper_N); + bool has_pos_lower = pos_root_lower_bound(n, p, pos_lower_N); + bool has_pos_upper = pos_root_upper_bound(n, p, pos_upper_N); + TRACE("rcf_isolate", + display_poly(tout, n, p); tout << "\n"; + if (has_neg_lower) tout << "-2^" << neg_lower_N; else tout << "-oo"; + tout << " < neg-roots < "; + if (has_neg_upper) tout << "-2^" << neg_upper_N; else tout << "0"; + tout << "\n"; + if (has_pos_lower) tout << "2^" << pos_lower_N; else tout << "0"; + tout << " < pos-roots < "; + if (has_pos_upper) tout << "2^" << pos_upper_N; else tout << "oo"; + tout << "\n";); // TODO } @@ -1266,62 +1320,62 @@ namespace realclosure { - square free - nonconstant */ - void nz_sqf_isolate_roots(unsigned n, value * const * as, numeral_vector & roots) { + void nz_sqf_isolate_roots(unsigned n, value * const * p, numeral_vector & roots) { SASSERT(n > 1); - SASSERT(!is_zero(as[0])); - SASSERT(!is_zero(as[n-1])); + SASSERT(!is_zero(p[0])); + SASSERT(!is_zero(p[n-1])); if (n == 2) { // we don't need a field extension for linear polynomials. numeral r; value_ref v(*this); - neg(as[0], v); - div(v, as[1], v); + neg(p[0], v); + div(v, p[1], v); set(r, v); roots.push_back(r); } else { - nl_nz_sqf_isolate_roots(n, as, roots); + nl_nz_sqf_isolate_roots(n, p, roots); } } /** \brief Root isolation for polynomials where 0 is not a root. */ - void nz_isolate_roots(unsigned n, value * const * as, numeral_vector & roots) { + void nz_isolate_roots(unsigned n, value * const * p, numeral_vector & roots) { SASSERT(n > 0); - SASSERT(!is_zero(as[0])); - SASSERT(!is_zero(as[n-1])); + SASSERT(!is_zero(p[0])); + SASSERT(!is_zero(p[n-1])); if (n == 1) { // constant polynomial return; } value_ref_buffer sqf(*this); - square_free(n, as, sqf); + square_free(n, p, sqf); nz_sqf_isolate_roots(sqf.size(), sqf.c_ptr(), roots); } /** \brief Root isolation entry point. */ - void isolate_roots(unsigned n, numeral const * as, numeral_vector & roots) { + void isolate_roots(unsigned n, numeral const * p, numeral_vector & roots) { SASSERT(n > 0); - SASSERT(!is_zero(as[n-1])); + SASSERT(!is_zero(p[n-1])); if (n == 1) { // constant polynomial return; } unsigned i = 0; for (; i < n; i++) { - if (!is_zero(as[i])) + if (!is_zero(p[i])) break; } SASSERT(i < n); - SASSERT(!is_zero(as[i])); - ptr_buffer as_values; + SASSERT(!is_zero(p[i])); + ptr_buffer nz_p; for (; i < n; i++) - as_values.push_back(as[i].m_value); - nz_isolate_roots(as_values.size(), as_values.c_ptr(), roots); - if (as_values.size() < n) { + nz_p.push_back(p[i].m_value); + nz_isolate_roots(nz_p.size(), nz_p.c_ptr(), roots); + if (nz_p.size() < n) { // zero is a root roots.push_back(numeral()); } @@ -1725,7 +1779,7 @@ namespace realclosure { \brief r <- rem(p1, p2) */ void rem(unsigned sz1, value * const * p1, unsigned sz2, value * const * p2, value_ref_buffer & r) { - TRACE("rcf", + TRACE("rcf_rem", tout << "rem\n"; display_poly(tout, sz1, p1); tout << "\n"; display_poly(tout, sz2, p2); tout << "\n";); @@ -1744,7 +1798,7 @@ namespace realclosure { checkpoint(); sz1 = r.size(); if (sz1 < sz2) { - TRACE("rcf", tout << "rem result\n"; display_poly(tout, r.size(), r.c_ptr()); tout << "\n";); + TRACE("rcf_rem", tout << "rem result\n"; display_poly(tout, r.size(), r.c_ptr()); tout << "\n";); return; } unsigned m_n = sz1 - sz2; @@ -2027,7 +2081,7 @@ namespace realclosure { refine_transcendental_interval(to_transcendental(v->ext()), _prec); update_rf_interval(v, prec); - TRACE("rcf", tout << "after update_rf_interval: " << magnitude(v->interval()) << " "; + TRACE("rcf_transcendental", tout << "after update_rf_interval: " << magnitude(v->interval()) << " "; bqim().display(tout, v->interval()); tout << std::endl;); if (check_precision(v->interval(), prec)) diff --git a/src/test/rcf.cpp b/src/test/rcf.cpp index 0583c54a0..40c3e838f 100644 --- a/src/test/rcf.cpp +++ b/src/test/rcf.cpp @@ -143,7 +143,7 @@ static void tst_lin_indep(unsigned m, unsigned n, int _A[], unsigned ex_sz, unsi void tst_rcf() { - // tst1(); + tst1(); tst2(); { int A[] = {0, 1, 1, 1, 0, 1, 1, 1, -1}; int c[] = {10, 4, -4}; int b[] = {-2, 4, 6}; tst_solve(3, A, b, c, true); } { int A[] = {1, 1, 1, 0, 1, 1, 0, 1, 1}; int c[] = {3, 2, 2}; int b[] = {1, 1, 1}; tst_solve(3, A, b, c, false); }