From 14827e94f0348712d4a2cc11475737eb6188aec4 Mon Sep 17 00:00:00 2001 From: Leonardo de Moura Date: Fri, 4 Jan 2013 15:01:27 -0800 Subject: [PATCH] Fix typos and bugs. Add tests. Signed-off-by: Leonardo de Moura --- src/math/realclosure/realclosure.cpp | 159 ++++++++++++++++++++++----- src/math/realclosure/realclosure.h | 6 +- src/test/algebraic.cpp | 2 +- src/test/rcf.cpp | 7 ++ 4 files changed, 144 insertions(+), 30 deletions(-) diff --git a/src/math/realclosure/realclosure.cpp b/src/math/realclosure/realclosure.cpp index e426bc0cf..eae2f6087 100644 --- a/src/math/realclosure/realclosure.cpp +++ b/src/math/realclosure/realclosure.cpp @@ -76,11 +76,23 @@ namespace realclosure { numeral m_upper; unsigned m_lower_inf:1; unsigned m_upper_inf:1; - interval():m_lower_inf(true), m_upper_inf(true) {} - interval(numeral & l, numeral & u):m_lower_inf(false), m_upper_inf(false) { + unsigned m_lower_open:1; + unsigned m_upper_open:1; + interval():m_lower_inf(true), m_upper_inf(true), m_lower_open(true), m_upper_open(true) {} + interval(numeral & l, numeral & u):m_lower_inf(false), m_upper_inf(false), m_lower_open(true), m_upper_open(true) { swap(m_lower, l); swap(m_upper, u); } + numeral & lower() { return m_lower; } + numeral & upper() { return m_upper; } + void set_lower_is_inf(bool f) { m_lower_inf = f; } + void set_upper_is_inf(bool f) { m_upper_inf = f; } + void set_lower_is_open(bool f) { m_lower_open = f; } + void set_upper_is_open(bool f) { m_upper_open = f; } + numeral const & lower() const { return m_lower; } + numeral const & upper() const { return m_upper; } + bool lower_is_inf() const { return m_lower_inf; } + bool upper_is_inf() const { return m_upper_inf; } }; void set_rounding(bool to_plus_inf) { m_manager.m_to_plus_inf = to_plus_inf; } @@ -92,16 +104,16 @@ namespace realclosure { numeral const & upper(interval const & a) const { return a.m_upper; } numeral & lower(interval & a) { return a.m_lower; } numeral & upper(interval & a) { return a.m_upper; } - bool lower_is_open(interval const & a) const { return true; } - bool upper_is_open(interval const & a) const { return true; } + bool lower_is_open(interval const & a) const { return a.m_lower_open; } + bool upper_is_open(interval const & a) const { return a.m_upper_open; } bool lower_is_inf(interval const & a) const { return a.m_lower_inf; } bool upper_is_inf(interval const & a) const { return a.m_upper_inf; } // Setters void set_lower(interval & a, numeral const & n) { m_manager.set(a.m_lower, n); } void set_upper(interval & a, numeral const & n) { m_manager.set(a.m_upper, n); } - void set_lower_is_open(interval & a, bool v) { SASSERT(v); } - void set_upper_is_open(interval & a, bool v) { SASSERT(v); } + void set_lower_is_open(interval & a, bool v) { a.m_lower_open = v; } + void set_upper_is_open(interval & a, bool v) { a.m_upper_open = v; } void set_lower_is_inf(interval & a, bool v) { a.m_lower_inf = v; } void set_upper_is_inf(interval & a, bool v) { a.m_upper_inf = v; } @@ -124,14 +136,16 @@ namespace realclosure { // --------------------------------- // - // Values: rational and composite + // Values are represented as + // - arbitrary precision rationals (mpq) + // - rational functions on field extensions // // --------------------------------- struct value { - unsigned m_ref_count; - bool m_rational; - mpbqi m_interval; // approximation as a binary rational + unsigned m_ref_count; //!< Reference counter + bool m_rational; //!< True if the value is represented as an abitrary precision rational value. + mpbqi m_interval; //!< approximation as an interval with binary rational end-points value(bool rat):m_ref_count(0), m_rational(rat) {} bool is_rational() const { return m_rational; } mpbqi const & interval() const { return m_interval; } @@ -153,8 +167,8 @@ namespace realclosure { // it contains a nonzero coefficient at position k > 0. It is not a constant polynomial on m_ext. polynomial m_p; extension * m_ext; - bool m_real; - mpbqi m_interval; // approximation as a binary rational + bool m_real; //!< True if the polynomial expression does not depend on infinitesimal values + mpbqi m_interval; //!< approximation as an interval with binary rational end-points polynomial_expr():m_ext(0), m_real(false) {} extension * ext() const { return m_ext; } bool is_real() const { return m_real; } @@ -204,7 +218,7 @@ namespace realclosure { typedef int sign; - typedef std::pair p2s; + typedef std::pair p2s; typedef sarray signs; @@ -249,7 +263,7 @@ namespace realclosure { struct algebraic : public extension { polynomial m_p; signs m_signs; - bool m_real; + bool m_real; //!< True if the polynomial p does not depend on infinitesimal extensions. algebraic(unsigned idx):extension(ALGEBRAIC, idx), m_real(false) {} @@ -271,8 +285,7 @@ namespace realclosure { struct infinitesimal : public extension { symbol m_name; - mpbqi m_interval; - + infinitesimal(unsigned idx, symbol const & n):extension(INFINITESIMAL, idx), m_name(n) {} void display(std::ostream & out) const { @@ -296,7 +309,7 @@ namespace realclosure { mpbqi_manager m_bqim; ptr_vector m_extensions[3]; value * m_one; - unsigned m_eps_prec; + unsigned m_ini_precision; //!< initial precision for transcendentals, infinitesimals, etc. volatile bool m_cancel; struct scoped_polynomial_seq { @@ -352,6 +365,8 @@ namespace realclosure { m_one = mk_rational(one); inc_ref(m_one); m_cancel = false; + + updt_params(p); } ~imp() { @@ -398,7 +413,7 @@ namespace realclosure { } void updt_params(params_ref const & p) { - m_eps_prec = p.get_uint("eps_prec", 24); + m_ini_precision = p.get_uint("initial_precision", 24); // TODO } @@ -502,24 +517,48 @@ namespace realclosure { a.m_value = 0; } + /** + \brief Return true if the given interval is smaller than 1/2^k + */ + bool check_precision(mpbqi const & interval, unsigned k) { + if (interval.lower_is_inf() || interval.upper_is_inf()) + return false; + scoped_mpbq w(bqm()); + bqm().sub(interval.upper(), interval.lower(), w); + return bqm().lt_1div2k(w, k); + } + + /** + \brief Return true if v is zero. + */ static bool is_zero(value * v) { return v == 0; } + /** + \brief Return true if v is represented using a nonzero arbitrary precision rational value. + */ static bool is_nz_rational(value * v) { SASSERT(v != 0); return v->is_rational(); } + /** + \brief Return true if v is one. + */ bool is_one(value * v) const { + // TODO: make the check complete return !is_zero(v) && is_nz_rational(v) && qm().is_one(to_mpq(v)); } + /** + \brief Return true if v is a represented as a rational function of the set of field extensions. + */ static bool is_rational_function(value * v) { SASSERT(v != 0); return !(v->is_rational()); } - + static rational_value * to_nz_rational(value * v) { SASSERT(is_nz_rational(v)); return static_cast(v); @@ -659,7 +698,7 @@ namespace realclosure { m_extensions[extension::INFINITESIMAL].push_back(eps); value * p[2] = { 0, one() }; mpbq zero(0); - mpbq tiny(1, m_eps_prec); + mpbq tiny(1, m_ini_precision); mpbqi interval(zero, tiny); polynomial_expr * numerator = mk_polynomial_expr(2, p, eps, interval); rational_function_value * v = alloc(rational_function_value, numerator, 0); @@ -747,6 +786,45 @@ namespace realclosure { return to_rational_function(a)->is_real(); } + void mpq_to_mpbqi(mpq const & q, mpbqi & interval) { + interval.set_lower_is_inf(false); + interval.set_upper_is_inf(false); + if (bqm().to_mpbq(q, interval.lower())) { + bqm().set(interval.upper(), interval.lower()); + interval.set_lower_is_open(false); + interval.set_upper_is_open(false); + } + else { + bqm().set(interval.upper(), interval.lower()); + bqm().mul2(interval.upper()); + interval.set_lower_is_open(true); + interval.set_upper_is_open(true); + if (qm().is_neg(q)) { + ::swap(interval.lower(), interval.upper()); + } + while (bqim().contains_zero(interval) || !check_precision(interval, m_ini_precision)) { + checkpoint(); + bqm().refine_lower(q, interval.lower(), interval.upper()); + bqm().refine_upper(q, interval.lower(), interval.upper()); + } + } + } + + void initialize_rational_value_interval(value * a) { + // For rational values, we only compute the binary intervals if needed. + SASSERT(is_nz_rational(a)); + mpq_to_mpbqi(to_mpq(a), a->m_interval); + } + + mpbqi const & interval(value * a) const { + SASSERT(a != 0); + if (bqim().contains_zero(a->m_interval)) { + SASSERT(is_nz_rational(a)); + const_cast(this)->initialize_rational_value_interval(a); + } + return a->m_interval; + } + rational_value * mk_rational() { return new (allocator()) rational_value(); } @@ -757,6 +835,18 @@ namespace realclosure { return r; } + template + void update_mpq_value(value * a, T & v) { + SASSERT(is_nz_rational(a)); + qm().set(to_mpq(a), v); + bqim().reset(a->m_interval); + } + + template + void update_mpq_value(numeral & a, T & v) { + update_mpq_value(a.m_value, v); + } + void set(numeral & a, int n) { if (n == 0) { reset(a); @@ -769,7 +859,7 @@ namespace realclosure { inc_ref(a.m_value); } SASSERT(is_unique_nz_rational(a)); - qm().set(to_mpq(a), n); + update_mpq_value(a, n); } void set(numeral & a, mpz const & n) { @@ -784,7 +874,7 @@ namespace realclosure { inc_ref(a.m_value); } SASSERT(is_unique_nz_rational(a)); - qm().set(to_mpq(a), n); + update_mpq_value(a, n); } void set(numeral & a, mpq const & n) { @@ -799,7 +889,7 @@ namespace realclosure { inc_ref(a.m_value); } SASSERT(is_unique_nz_rational(a)); - qm().set(to_mpq(a), n); + update_mpq_value(a, n); } void set(numeral & a, numeral const & n) { @@ -1448,10 +1538,24 @@ namespace realclosure { void display(std::ostream & out, numeral const & a) const { display(out, a.m_value, false); } + + void display_non_rational_in_decimal(std::ostream & out, numeral const & a, unsigned precision) { + SASSERT(!is_zero(a)); + SASSERT(!is_nz_rational(a)); + mpbqi const & i = interval(a.m_value); + if (check_precision(i, precision*4)) { + // hack + if (bqm().is_int(i.lower())) + bqm().display_decimal(out, i.upper(), precision); + else + bqm().display_decimal(out, i.lower(), precision); + } + else { + out << "?"; + } + } void display_decimal(std::ostream & out, numeral const & a, unsigned precision) const { - if (!is_real(a)) - throw exception("number cannot be printed in decimal notation because it is not a real"); if (is_zero(a)) { out << "0"; } @@ -1459,8 +1563,7 @@ namespace realclosure { qm().display_decimal(out, to_mpq(a), precision); } else { - - // TODO + const_cast(this)->display_non_rational_in_decimal(out, a, precision); } } @@ -1468,7 +1571,7 @@ namespace realclosure { if (is_zero(a)) out << "[0, 0]"; else - bqim().display(out, a.m_value->interval()); + bqim().display(out, interval(a.m_value)); } }; diff --git a/src/math/realclosure/realclosure.h b/src/math/realclosure/realclosure.h index 6d1907a63..0c2fe8560 100644 --- a/src/math/realclosure/realclosure.h +++ b/src/math/realclosure/realclosure.h @@ -362,6 +362,10 @@ inline std::ostream & operator<<(std::ostream & out, rc_decimal_pp const & n) { return out; } +inline rc_decimal_pp decimal_pp(scoped_rcnumeral const & n, unsigned prec = 10) { + return rc_decimal_pp(n, prec); +} + struct rc_interval_pp { rcmanager & m; rcnumeral const & n; @@ -374,7 +378,7 @@ inline std::ostream & operator<<(std::ostream & out, rc_interval_pp const & n) { return out; } -inline rc_interval_pp interval_pp(rc_interval_pp const & n) { +inline rc_interval_pp interval_pp(scoped_rcnumeral const & n) { return rc_interval_pp(n); } diff --git a/src/test/algebraic.cpp b/src/test/algebraic.cpp index d758ee5d8..7c0fb4a95 100644 --- a/src/test/algebraic.cpp +++ b/src/test/algebraic.cpp @@ -208,7 +208,7 @@ void tst_refine_mpbq(int n, int d) { } void tst_refine_mpbq() { - tst_refine_mpbq(-5, 7); + tst_refine_mpbq(5, 7); } void tst_mpbq_root() { diff --git a/src/test/rcf.cpp b/src/test/rcf.cpp index e71a48dc5..87743e138 100644 --- a/src/test/rcf.cpp +++ b/src/test/rcf.cpp @@ -27,7 +27,14 @@ static void tst1() { scoped_rcnumeral eps(m); m.mk_infinitesimal("eps", eps); std::cout << sym_pp(eps) << std::endl; + std::cout << interval_pp(a) << std::endl; std::cout << interval_pp(eps) << std::endl; + mpq aux; + qm.set(aux, 1, 3); + m.set(a, aux); + std::cout << interval_pp(a) << std::endl; + std::cout << decimal_pp(eps, 4) << std::endl; + std::cout << decimal_pp(a) << std::endl; } void tst_rcf() {