From 217c8375ce0fc9de72446e79039cdf710492d50b Mon Sep 17 00:00:00 2001 From: Leonardo de Moura Date: Tue, 15 Jan 2013 08:39:30 -0800 Subject: [PATCH] Add new rational function normalization procedure. Signed-off-by: Leonardo de Moura --- src/math/realclosure/realclosure.cpp | 139 ++++++++++++++++----------- 1 file changed, 82 insertions(+), 57 deletions(-) diff --git a/src/math/realclosure/realclosure.cpp b/src/math/realclosure/realclosure.cpp index 7982e28b4..8f9697d98 100644 --- a/src/math/realclosure/realclosure.cpp +++ b/src/math/realclosure/realclosure.cpp @@ -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. */ @@ -4762,14 +4766,20 @@ namespace realclosure { bool determine_sign(rational_function_value * v) { if (!contains_zero(v->interval())) return true; + bool r; switch (v->ext()->knd()) { - case extension::TRANSCENDENTAL: determine_transcendental_sign(v); return true; // it is never zero - case extension::INFINITESIMAL: determine_infinitesimal_sign(v); return true; // it is never zero - case extension::ALGEBRAIC: return determine_algebraic_sign(v); + case extension::TRANSCENDENTAL: determine_transcendental_sign(v); r = true; break; // it is never zero + case extension::INFINITESIMAL: determine_infinitesimal_sign(v); r = true; break; // it is never zero + case extension::ALGEBRAIC: r = determine_algebraic_sign(v); break; default: 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) { @@ -4784,10 +4794,10 @@ namespace realclosure { // --------------------------------- /** - \brief Set new_p1 and new_p2 using the following normalization rules: - - new_p1 <- p1/p2[0]; new_p2 <- one IF sz2 == 1 - - new_p1 <- one; new_p2 <- p2/p1[0]; IF sz1 == 1 - - new_p1 <- p1/gcd(p1, p2); new_p2 <- p2/gcd(p1, p2); Otherwise + \brief Compute polynomials new_p1 and new_p2 s.t. + - p1/p2 == new_p1/new_p2, AND + - new_p2 is a Monic polynomial, AND + - 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) { INC_DEPTH(); @@ -4800,47 +4810,58 @@ namespace realclosure { div(sz1, p1, p2[0], new_p1); 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 { - // - new_p1 <- p1/gcd(p1, p2); new_p2 <- p2/gcd(p1, p2); Otherwise - value_ref_buffer g(*this); - gcd(sz1, p1, sz2, p2, g); - if (is_rational_one(g)) { - 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(); - } + value * lc = p2[sz2 - 1]; + if (is_rational_one(lc)) { + // p2 is monic + normalize_num_monic_den(sz1, p1, sz2, p2, new_p1, new_p2); } else { - div(sz1, p1, g.size(), g.c_ptr(), new_p1); - div(sz2, p2, g.size(), g.c_ptr(), new_p2); - SASSERT(new_p1.size() > 1); - SASSERT(new_p2.size() > 1); + // p2 is not monic + value_ref_buffer tmp1(*this); + value_ref_buffer tmp2(*this); + 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_den(*this); normalize_fraction(num.size(), num.c_ptr(), ad.size(), ad.c_ptr(), new_num, new_den); - if (new_num.empty()) - r = 0; - else - mk_add_value(a, b, new_num.size(), new_num.c_ptr(), new_den.size(), new_den.c_ptr(), r); + SASSERT(!new_num.empty()); + 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_den(*this); normalize_fraction(num.size(), num.c_ptr(), den.size(), den.c_ptr(), new_num, new_den); - if (new_num.empty()) - r = 0; - else - mk_add_value(a, b, new_num.size(), new_num.c_ptr(), new_den.size(), new_den.c_ptr(), r); + SASSERT(!new_num.empty()); + 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(); scoped_mpbqi ri(bqim()); 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); 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; - c.mark(a.m_value); + c.mark(a); if (c.m_found.empty()) { - display(out, a.m_value, true); + display(out, a, true); } else { std::sort(c.m_found.begin(), c.m_found.end(), rank_lt_proc()); out << "["; - display(out, a.m_value, true); + display(out, a, true); for (unsigned i = 0; i < c.m_found.size(); i++) { algebraic * ext = c.m_found[i]; 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 { if (compact) display_compact(out, a);