Fix bugs in mpbq.

Signed-off-by: Leonardo de Moura <leonardo@microsoft.com>

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);
     }
 }
 

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<unsigned>(a); }
+
     friend mpbq operator+(mpbq a, mpbq const & b) { return a += b; }
     template<typename T>
     friend mpbq operator+(mpbq a, T const & b) { return a += mpbq(b); }
@@ -150,6 +161,9 @@ public:
     template<typename T>
     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; }
 

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));

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; }

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)

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 <iostream>
+#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;
+
+}

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)

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;
+}
+}

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);
+}