From 03cc3739d4017ed43866a9623ca014f5a06f456d Mon Sep 17 00:00:00 2001 From: Leonardo de Moura Date: Sun, 21 Jul 2013 20:12:04 -0700 Subject: [PATCH] Fix bugs in mpbq. Signed-off-by: Leonardo de Moura --- src/numerics/mpbq.cpp | 28 +++++++++++++++++-- src/numerics/mpbq.h | 14 ++++++++++ src/numerics/mpq.cpp | 9 ++++++ src/numerics/mpq.h | 2 ++ src/tests/numerics/CMakeLists.txt | 3 ++ src/tests/numerics/mpbq.cpp | 46 +++++++++++++++++++++++++++++++ src/util/CMakeLists.txt | 2 +- src/util/bit_tricks.cpp | 32 +++++++++++++++++++++ src/util/bit_tricks.h | 12 ++++++++ 9 files changed, 145 insertions(+), 3 deletions(-) create mode 100644 src/tests/numerics/mpbq.cpp create mode 100644 src/util/bit_tricks.cpp create mode 100644 src/util/bit_tricks.h diff --git a/src/numerics/mpbq.cpp b/src/numerics/mpbq.cpp index 2cec33cb2..4c3c0f9c7 100644 --- a/src/numerics/mpbq.cpp +++ b/src/numerics/mpbq.cpp @@ -9,6 +9,30 @@ Author: Leonardo de Moura namespace lean { +bool set(mpbq & a, mpq const & b) { + if (b.is_integer()) { + numerator(a.m_num, b); + a.m_k = 0; + return true; + } + else { + static thread_local mpz d; + denominator(d, b); + unsigned shift; + if (d.is_power_of_two(shift)) { + numerator(a.m_num, b); + a.m_k = shift; + lean_assert(a == b); + return true; + } + else { + numerator(a.m_num, b); + a.m_k = d.log2() + 1; + return false; + } + } +} + void mpbq::normalize() { if (m_k == 0) return; @@ -51,7 +75,7 @@ int cmp(mpbq const & a, mpz const & b) { int cmp(mpbq const & a, mpq const & b) { if (a.is_integer() && b.is_integer()) { - return a.m_num < b; + return -cmp(b, a.m_num); } else { static thread_local mpz tmp1; @@ -60,7 +84,7 @@ int cmp(mpbq const & a, mpq const & b) { denominator(tmp1, b); tmp1 *= a.m_num; // tmp2 <- numerator(b)*denominator(a) numerator(tmp2, b); mul2k(tmp2, tmp2, a.m_k); - return tmp1 < tmp2; + return cmp(tmp1, tmp2); } } diff --git a/src/numerics/mpbq.h b/src/numerics/mpbq.h index e2d94b75c..2078e6254 100644 --- a/src/numerics/mpbq.h +++ b/src/numerics/mpbq.h @@ -7,6 +7,7 @@ Author: Leonardo de Moura #pragma once #include "mpz.h" #include "mpq.h" +#include "bit_tricks.h" namespace lean { @@ -23,12 +24,19 @@ public: mpbq(mpbq const & v):m_num(v.m_num), m_k(v.m_k) {} mpbq(mpbq && v); explicit mpbq(mpz const & v):m_num(v), m_k(0) {} + explicit mpbq(mpq const & v) { set(*this, v); } explicit mpbq(int n):m_num(n), m_k(0) {} mpbq(int n, unsigned k):m_num(n), m_k(k) { normalize(); } ~mpbq() {} + // a <- b + // Return true iff b is a binary rational. + friend bool set(mpbq & a, mpq const & b); + mpbq & operator=(mpbq const & v) { m_num = v.m_num; m_k = v.m_k; return *this; } mpbq & operator=(mpbq && v) { swap(*this, v); return *this; } + mpbq & operator=(mpz const & v) { m_num = v; m_k = 0; return *this; } + mpbq & operator=(mpq const & v) { set(*this, v); return *this; } mpbq & operator=(unsigned int v) { m_num = v; m_k = 0; return *this; } mpbq & operator=(int v) { m_num = v; m_k = 0; return *this; } @@ -132,6 +140,9 @@ public: mpbq & operator*=(unsigned a); mpbq & operator*=(int a); + mpbq & operator/=(unsigned a) { lean_assert(is_power_of_two(a)); div2k(*this, log2(a)); return *this; } + mpbq & operator/=(int a) { if (a < 0) { a = -a; neg(); } return *this /= static_cast(a); } + friend mpbq operator+(mpbq a, mpbq const & b) { return a += b; } template friend mpbq operator+(mpbq a, T const & b) { return a += mpbq(b); } @@ -150,6 +161,9 @@ public: template friend mpbq operator*(T const & a, mpbq b) { return b *= mpbq(a); } + friend mpbq operator/(mpbq a, unsigned b) { return a /= b; } + friend mpbq operator/(mpbq a, int b) { return a /= b; } + mpbq & operator++() { return operator+=(1); } mpbq operator++(int) { mpbq r(*this); ++(*this); return r; } diff --git a/src/numerics/mpq.cpp b/src/numerics/mpq.cpp index 90dad4f9a..96249c29b 100644 --- a/src/numerics/mpq.cpp +++ b/src/numerics/mpq.cpp @@ -5,9 +5,18 @@ Released under Apache 2.0 license as described in the file LICENSE. Author: Leonardo de Moura */ #include "mpq.h" +#include "mpbq.h" namespace lean { +mpq & mpq::operator=(mpbq const & b) { + *this = 2; + power(*this, *this, b.get_k()); + inv(); + *this *= b.get_numerator(); + return *this; +} + int cmp(mpq const & a, mpz const & b) { if (a.is_integer()) { return mpz_cmp(mpq_numref(a.m_val), mpq::zval(b)); diff --git a/src/numerics/mpq.h b/src/numerics/mpq.h index 893b3b171..b68ac8cda 100644 --- a/src/numerics/mpq.h +++ b/src/numerics/mpq.h @@ -8,6 +8,7 @@ Author: Leonardo de Moura #include "mpz.h" namespace lean { +class mpbq; /** \brief Wrapper for GMP rationals @@ -24,6 +25,7 @@ public: mpq & operator=(mpz const & v) { mpq_set_z(m_val, v.m_val); return *this; } mpq & operator=(mpq const & v) { mpq_set(m_val, v.m_val); return *this; } mpq & operator=(mpq && v) { swap(*this, v); return *this; } + mpq & operator=(mpbq const & b); mpq & operator=(char const * v) { mpq_set_str(m_val, v, 10); return *this; } mpq & operator=(unsigned long int v) { mpq_set_ui(m_val, v, 1u); return *this; } mpq & operator=(long int v) { mpq_set_si(m_val, v, 1); return *this; } diff --git a/src/tests/numerics/CMakeLists.txt b/src/tests/numerics/CMakeLists.txt index 19fd29e0e..8d34b39f7 100644 --- a/src/tests/numerics/CMakeLists.txt +++ b/src/tests/numerics/CMakeLists.txt @@ -1,3 +1,6 @@ add_executable(mpq mpq.cpp) target_link_libraries(mpq ${EXTRA_LIBS}) add_test(mpq ${CMAKE_CURRENT_BINARY_DIR}/mpq) +add_executable(mpbq mpbq.cpp) +target_link_libraries(mpbq ${EXTRA_LIBS}) +add_test(mpbq ${CMAKE_CURRENT_BINARY_DIR}/mpbq) diff --git a/src/tests/numerics/mpbq.cpp b/src/tests/numerics/mpbq.cpp new file mode 100644 index 000000000..46b899d64 --- /dev/null +++ b/src/tests/numerics/mpbq.cpp @@ -0,0 +1,46 @@ +/* +Copyright (c) 2013 Microsoft Corporation. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. + +Author: Leonardo de Moura +*/ +#include +#include "test.h" +#include "mpbq.h" +#include "mpq.h" +using namespace lean; + +static void tst1() { + mpbq a(1,1); + std::cout << a << "\n"; + lean_assert(a == mpq(1,2)); + std::cout << mpbq(mpq(1,3)) << "\n"; + lean_assert(!set(a, mpq(1,3))); + lean_assert(a == mpbq(1,2)); + mpq b; + b = a; + lean_assert(b == mpq(1,4)); + lean_assert(inv(b) == 4); + lean_assert(a/2 == mpbq(1,3)); + lean_assert(a/1 == a); + lean_assert(a/8 == mpbq(1,5)); + lean_assert((3*a)/8 == mpbq(3,5)); + mpbq l(0), u(1); + mpq v(1,3); + while (true) { + lean_assert(l < v); + lean_assert(v < u); + std::cout << mpbq::decimal(l,20) << " " << v << " " << mpbq::decimal(u, 20) << "\n"; + if (lt_1div2k((u - l)/2, 50)) + break; + refine_lower(v, l, u); + refine_upper(v, l, u); + } +} + +int main() { + continue_on_violation(true); + tst1(); + return has_violations() ? 1 : 0; + +} diff --git a/src/util/CMakeLists.txt b/src/util/CMakeLists.txt index 410441219..cc94ad9a8 100644 --- a/src/util/CMakeLists.txt +++ b/src/util/CMakeLists.txt @@ -1 +1 @@ -add_library(util trace.cpp debug.cpp name.cpp exception.cpp verbosity.cpp interrupt.cpp hash.cpp escaped.cpp) +add_library(util trace.cpp debug.cpp name.cpp exception.cpp verbosity.cpp interrupt.cpp hash.cpp escaped.cpp bit_tricks.cpp) diff --git a/src/util/bit_tricks.cpp b/src/util/bit_tricks.cpp new file mode 100644 index 000000000..00b01cbe6 --- /dev/null +++ b/src/util/bit_tricks.cpp @@ -0,0 +1,32 @@ +/* +Copyright (c) 2013 Microsoft Corporation. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. + +Author: Leonardo de Moura +*/ +namespace lean { +unsigned log2(unsigned v) { + unsigned r = 0; + if (v & 0xFFFF0000) { + v >>= 16; + r |= 16; + } + if (v & 0xFF00) { + v >>= 8; + r |= 8; + } + if (v & 0xF0) { + v >>= 4; + r |= 4; + } + if (v & 0xC) { + v >>= 2; + r |= 2; + } + if (v & 0x2) { + v >>= 1; + r |= 1; + } + return r; +} +} diff --git a/src/util/bit_tricks.h b/src/util/bit_tricks.h new file mode 100644 index 000000000..0bffaff41 --- /dev/null +++ b/src/util/bit_tricks.h @@ -0,0 +1,12 @@ +/* +Copyright (c) 2013 Microsoft Corporation. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. + +Author: Leonardo de Moura +*/ +#pragma once + +namespace lean { +inline bool is_power_of_two(unsigned v) { return !(v & (v - 1)) && v; } +unsigned log2(unsigned v); +}