From 1d761ea9a599171dfb1c73a83072b023303ba41b Mon Sep 17 00:00:00 2001 From: Leonardo de Moura Date: Sat, 12 Jan 2013 15:45:43 -0800 Subject: [PATCH] Add clean_denominators procedure Signed-off-by: Leonardo de Moura --- src/math/realclosure/realclosure.cpp | 376 +++++++++++++++++++++++++-- src/math/realclosure/realclosure.h | 2 + src/test/rcf.cpp | 23 +- 3 files changed, 373 insertions(+), 28 deletions(-) diff --git a/src/math/realclosure/realclosure.cpp b/src/math/realclosure/realclosure.cpp index d8babe671..998e5abb2 100644 --- a/src/math/realclosure/realclosure.cpp +++ b/src/math/realclosure/realclosure.cpp @@ -2416,13 +2416,32 @@ namespace realclosure { return new (allocator()) rational_value(); } - rational_value * mk_rational(mpq & v) { + /** + \brief Make a rational and swap its value with v + */ + rational_value * mk_rational_and_swap(mpq & v) { + SASSERT(!qm().is_zero(v)); rational_value * r = mk_rational(); ::swap(r->m_value, v); return r; } + rational_value * mk_rational(mpq const & v) { + SASSERT(!qm().is_zero(v)); + rational_value * r = mk_rational(); + qm().set(r->m_value, v); + return r; + } + + rational_value * mk_rational(mpz const & v) { + SASSERT(!qm().is_zero(v)); + rational_value * r = mk_rational(); + qm().set(r->m_value, v); + return r; + } + rational_value * mk_rational(mpbq const & v) { + SASSERT(!bqm().is_zero(v)); scoped_mpq v_q(qm()); // v as a rational ::to_mpq(qm(), v, v_q); return mk_rational(v_q); @@ -2921,7 +2940,7 @@ namespace realclosure { for (unsigned i = 1; i < sz; i++) { mpq i_mpq(i); value_ref a_i(*this); - a_i = mk_rational(i_mpq); + a_i = mk_rational_and_swap(i_mpq); mul(a_i, p[i], a_i); r.push_back(a_i); } @@ -3086,30 +3105,42 @@ namespace realclosure { return true; } + /** + \brief r <- p(b) + */ + void mk_polynomial_value(unsigned n, value * const * p, value * b, value_ref & r) { + SASSERT(n > 0); + if (n == 1 || b == 0) { + r = p[0]; + } + else { + SASSERT(n >= 2); + // We compute the result using the Horner Sequence + // ((a_{n-1}*b + a_{n-2})*b + a_{n-3})*b + a_{n-4} ... + // where a_i's are the coefficients of p. + mul(p[n - 1], b, r); // r <- a_{n-1} * b + unsigned i = n - 1; + while (i > 0) { + --i; + if (p[i] != 0) + add(r, p[i], r); // r <- r + a_i + if (i > 0) + mul(r, b, r); // r <- r * b + } + } + } + /** \brief Evaluate the sign of p(b) by computing a value object. */ int expensive_eval_sign_at(unsigned n, value * const * p, mpbq const & b) { SASSERT(n > 0); SASSERT(p[n - 1] != 0); - value_ref bv(*this); - bv = mk_rational(b); - // We compute the result using the Horner Sequence - // ((a_{n-1}*bv + a_{n-2})*bv + a_{n-3})*bv + a_{n-4} ... - // where a_i's are the coefficients of p. - value_ref r(*this); - // r <- a_{n-1} * bv - mul(p[n - 1], bv, r); - unsigned i = n - 1; - while (i > 0) { - checkpoint(); - --i; - if (p[i] != 0) - add(r, p[i], r); // r <- r + a_i - if (i > 0) - mul(r, bv, r); // r <- r * bv - } - return sign(r); + value_ref _b(*this); + _b = mk_rational(b); + value_ref pb(*this); + mk_polynomial_value(n, p, _b, pb); + return sign(pb); } /** @@ -4290,7 +4321,7 @@ namespace realclosure { if (qm().is_zero(v)) r = 0; else - r = mk_rational(v); + r = mk_rational_and_swap(v); } else { INC_DEPTH(); @@ -4319,7 +4350,7 @@ namespace realclosure { if (qm().is_zero(v)) r = 0; else - r = mk_rational(v); + r = mk_rational_and_swap(v); } else { value_ref neg_b(*this); @@ -4353,7 +4384,7 @@ namespace realclosure { scoped_mpq v(qm()); qm().set(v, to_mpq(a)); qm().neg(v); - r = mk_rational(v); + r = mk_rational_and_swap(v); } else { neg_rf(to_rational_function(a), r); @@ -4495,7 +4526,7 @@ namespace realclosure { else if (is_nz_rational(a) && is_nz_rational(b)) { scoped_mpq v(qm()); qm().mul(to_mpq(a), to_mpq(b), v); - r = mk_rational(v); + r = mk_rational_and_swap(v); } else { INC_DEPTH(); @@ -4530,7 +4561,7 @@ namespace realclosure { else if (is_nz_rational(a) && is_nz_rational(b)) { scoped_mpq v(qm()); qm().div(to_mpq(a), to_mpq(b), v); - r = mk_rational(v); + r = mk_rational_and_swap(v); } else { value_ref inv_b(*this); @@ -4561,7 +4592,7 @@ namespace realclosure { if (is_nz_rational(a)) { scoped_mpq v(qm()); qm().inv(to_mpq(a), v); - r = mk_rational(v); + r = mk_rational_and_swap(v); } else { inv_rf(to_rational_function(a), r); @@ -4651,6 +4682,280 @@ namespace realclosure { return compare(a.m_value, b.m_value); } + // --------------------------------- + // + // Structural equality + // + // --------------------------------- + + /** + \brief Values a and b are said to be "structurally" equal if: + - a and b are 0. + - a and b are rationals and compare(a, b) == 0 + - a and b are rational function values p_a(x)/q_a(x) and p_b(y)/q_b(y) where x and y are field extensions, and + * x == y (pointer equality, i.e., they are the same field extension object). + * Every coefficient of p_a is structurally equal to every coefficient of p_b + * Every coefficient of q_a is structurally equal to every coefficient of q_b + Clearly structural equality implies equality, but the reverse is not true. + */ + bool struct_eq(value * a, value * b) const { + if (a == b) + return true; + else if (a == 0 || b == 0) + return false; + else if (is_nz_rational(a) && is_nz_rational(b)) + return qm().eq(to_mpq(a), to_mpq(b)); + else if (is_nz_rational(a) || is_nz_rational(b)) + return false; + else { + SASSERT(is_rational_function(a)); + SASSERT(is_rational_function(b)); + rational_function_value * rf_a = to_rational_function(a); + rational_function_value * rf_b = to_rational_function(b); + if (rf_a->ext() != rf_b->ext()) + return false; + return struct_eq(rf_a->num(), rf_b->num()) && struct_eq(rf_a->den(), rf_b->den()); + } + } + + /** + Auxiliary method for + bool struct_eq(value * a, value * b) + */ + bool struct_eq(unsigned sz_a, value * const * p_a, unsigned sz_b, value * const * p_b) const { + if (sz_a != sz_b) + return false; + for (unsigned i = 0; i < sz_a; i++) { + if (!struct_eq(p_a[i], p_b[i])) + return false; + } + return true; + } + + /** + Auxiliary method for + bool struct_eq(value * a, value * b) + */ + bool struct_eq(polynomial const & p_a, polynomial const & p_b) const { + return struct_eq(p_a.size(), p_a.c_ptr(), p_b.size(), p_b.c_ptr()); + } + + // --------------------------------- + // + // Clean denominators + // + // --------------------------------- + + /** + \brief We say 'a' has "clean" denominators if + - a is 0 + - a is a rational_value that is an integer + - a is a rational_function_value of the form p_a(x)/1 where the coefficients of p_a also have clean denominators. + */ + bool has_clean_denominators(value * a) const { + if (a == 0) + return true; + else if (is_nz_rational(a)) + return qm().is_int(to_mpq(a)); + else { + rational_function_value * rf_a = to_rational_function(a); + return is_rational_one(rf_a->den()) && has_clean_denominators(rf_a->num()); + } + } + + /** + \brief See comment at has_clean_denominators(value * a) + */ + bool has_clean_denominators(polynomial const & p) const { + unsigned sz = p.size(); + for (unsigned i = 0; i < sz; i++) { + if (!has_clean_denominators(p[i])) + return false; + } + return true; + } + + /** + \brief "Clean" the denominators of 'a'. That is, return p and q s.t. + a == p/q + and + has_clean_denominators(p) and has_clean_denominators(q) + */ + void clean_denominators_core(value * a, value_ref & p, value_ref & q) { + INC_DEPTH(); + TRACE("rcf_clean", tout << "clean_denominators_core [" << m_exec_depth << "]\na: "; display(tout, a, false); tout << "\n";); + p.reset(); q.reset(); + if (a == 0) { + p = a; + q = one(); + } + else if (is_nz_rational(a)) { + p = mk_rational(to_mpq(a).numerator()); + q = mk_rational(to_mpq(a).denominator()); + } + else { + rational_function_value * rf_a = to_rational_function(a); + value_ref_buffer p_num(*this), p_den(*this); + value_ref d_num(*this), d_den(*this); + clean_denominators_core(rf_a->num(), p_num, d_num); + clean_denominators_core(rf_a->den(), p_den, d_den); + value_ref x(*this); + x = mk_rational_function_value(rf_a->ext()); + mk_polynomial_value(p_num.size(), p_num.c_ptr(), x, p); + mk_polynomial_value(p_den.size(), p_den.c_ptr(), x, q); + if (!struct_eq(d_den, d_num)) { + mul(p, d_den, p); + mul(q, d_num, q); + } + } + } + + /** + \brief Clean the denominators of the polynomial p, it returns clean_p and d s.t. + p = clean_p/d + and has_clean_denominators(clean_p) && has_clean_denominators(d) + */ + void clean_denominators_core(polynomial const & p, value_ref_buffer & norm_p, value_ref & d) { + value_ref_buffer nums(*this), dens(*this); + value_ref a_n(*this), a_d(*this); + bool all_one = true; + for (unsigned i = 0; i < p.size(); i++) { + if (p[i]) { + clean_denominators_core(p[i], a_n, a_d); + nums.push_back(a_n); + if (!is_rational_one(a_d)) + all_one = false; + dens.push_back(a_d); + } + else { + nums.push_back(0); + dens.push_back(0); + } + } + if (all_one) { + norm_p = nums; + d = one(); + } + else { + // Compute lcm of the integer elements in dens. + // This is a little trick to control the coefficient growth. + // We don't compute lcm of the other elements of dens because it is too expensive. + scoped_mpq lcm_z(qm()); + bool found_z = false; + SASSERT(nums.size() == p.size()); + SASSERT(dens.size() == p.size()); + for (unsigned i = 0; i < p.size(); i++) { + if (!dens[i]) + continue; + if (is_nz_rational(dens[i])) { + mpq const & _d = to_mpq(dens[i]); + SASSERT(qm().is_int(_d)); + if (!found_z) { + found_z = true; + qm().set(lcm_z, _d); + } + else { + qm().lcm(lcm_z, _d, lcm_z); + } + } + } + + value_ref lcm(*this); + if (found_z) { + lcm = mk_rational(lcm_z); + } + else { + lcm = one(); + } + + // Compute the multipliers for nums. + // Compute norm_p and d + // + // We do NOT use GCD to compute the LCM of the denominators of non-rational values. + // However, we detect structurally equivalent denominators. + // + // Thus a/(b+1) + c/(b+1) is converted into a*c/(b+1) instead of (a*(b+1) + c*(b+1))/(b+1)^2 + norm_p.reset(); + d = lcm; + value_ref_buffer multipliers(*this); + value_ref m(*this); + for (unsigned i = 0; i < p.size(); i++) { + if (!nums[i]) { + norm_p.push_back(0); + } + else { + SASSERT(dens[i]); + bool is_z; + if (!is_nz_rational(dens[i])) { + m = lcm; + is_z = false; + } + else { + scoped_mpq num_z(qm()); + qm().div(lcm_z, to_mpq(dens[i]), num_z); + SASSERT(qm().is_int(num_z)); + m = mk_rational_and_swap(num_z); + is_z = true; + } + bool found_lt_eq = false; + for (unsigned j = 0; j < p.size(); j++) { + TRACE("rcf_clean_bug", tout << "j: " << j << " "; display(tout, m, false); tout << "\n";); + if (!dens[j]) + continue; + if (i != j && !is_nz_rational(dens[j])) { + if (struct_eq(dens[i], dens[j])) { + if (j < i) + found_lt_eq = true; + } + else { + mul(m, dens[j], m); + } + } + } + if (!is_z && !found_lt_eq) { + mul(dens[i], d, d); + } + mul(m, nums[i], m); + norm_p.push_back(m); + } + } + } + SASSERT(norm_p.size() == p.size()); + } + + void clean_denominators(value * a, value_ref & p, value_ref & q) { + if (has_clean_denominators(a)) { + p = a; + q = one(); + } + else { + clean_denominators_core(a, p, q); + } + } + + void clean_denominators(polynomial const & p, value_ref_buffer & norm_p, value_ref & d) { + if (has_clean_denominators(p)) { + norm_p.append(p.size(), p.c_ptr()); + d = one(); + } + else { + clean_denominators_core(p, norm_p, d); + } + } + + void clean_denominators(numeral const & a, numeral & p, numeral & q) { + value_ref _p(*this), _q(*this); + clean_denominators(a.m_value, _p, _q); + set(p, _p); + set(q, _q); + } + + // --------------------------------- + // + // "Pretty printing" + // + // --------------------------------- + struct collect_algebraic_refs { char_vector m_visited; // Set of visited algebraic extensions. ptr_vector m_found; // vector/list of visited algebraic extensions. @@ -4683,11 +4988,20 @@ namespace realclosure { } }; + static unsigned num_nz_coeffs(polynomial const & p) { + unsigned r = 0; + for (unsigned i = 0; i < p.size(); i++) { + if (p[i]) + r++; + } + return r; + } + bool use_parenthesis(value * v) const { if (is_zero(v) || is_nz_rational(v)) return false; rational_function_value * rf = to_rational_function(v); - return rf->num().size() > 1 || !is_rational_one(rf->den()); + return num_nz_coeffs(rf->num()) > 1 || !is_rational_one(rf->den()); } template @@ -5140,6 +5454,10 @@ namespace realclosure { m_imp->display_interval(out, a); } + void manager::clean_denominators(numeral const & a, numeral & p, numeral & q) { + save_interval_ctx ctx(this); + m_imp->clean_denominators(a, p, q); + } }; void pp(realclosure::manager::imp * imp, realclosure::polynomial const & p, realclosure::extension * ext) { @@ -5162,6 +5480,10 @@ void pp(realclosure::manager::imp * imp, realclosure::manager::imp::value_ref_bu pp(imp, p[i]); } +void pp(realclosure::manager::imp * imp, realclosure::manager::imp::value_ref const & v) { + pp(imp, v.get()); +} + void pp(realclosure::manager::imp * imp, realclosure::mpbqi const & i) { imp->bqim().display(std::cout, i); std::cout << std::endl; diff --git a/src/math/realclosure/realclosure.h b/src/math/realclosure/realclosure.h index a106ebc26..1e0ac4ba5 100644 --- a/src/math/realclosure/realclosure.h +++ b/src/math/realclosure/realclosure.h @@ -263,6 +263,8 @@ namespace realclosure { void display_interval(std::ostream & out, numeral const & a) const; + + void clean_denominators(numeral const & a, numeral & p, numeral & q); }; class value; diff --git a/src/test/rcf.cpp b/src/test/rcf.cpp index ac240fa9d..5bef0dfef 100644 --- a/src/test/rcf.cpp +++ b/src/test/rcf.cpp @@ -135,15 +135,36 @@ static void tst_lin_indep(unsigned m, unsigned n, int _A[], unsigned ex_sz, unsi A.set(i, j, _A[i*n + j]); unsigned_vector r; r.resize(A.n()); - mm.linear_independent_rows(A, r.c_ptr()); + scoped_mpz_matrix B(mm); + mm.linear_independent_rows(A, r.c_ptr(), B); SASSERT(r.size() == ex_sz); for (unsigned i = 0; i < ex_sz; i++) { SASSERT(r[i] == ex_r[i]); } } +static void tst_denominators() { + unsynch_mpq_manager qm; + rcmanager m(qm); + scoped_rcnumeral a(m); + scoped_rcnumeral t(m); + scoped_rcnumeral eps(m); + m.mk_pi(a); + m.inv(a); + m.mk_infinitesimal("eps", eps); + t = (a - eps*2) / (a*eps + 1); + // t = t + a * 2; + scoped_rcnumeral n(m), d(m); + std::cout << t << "\n"; + m.clean_denominators(t, n, d); + std::cout << "---->\n" << n << "\n" << d << "\n"; +} void tst_rcf() { + enable_trace("rcf_clean"); + enable_trace("rcf_clean_bug"); + tst_denominators(); + return; 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); }