lean2/src/util/numerics/mpbq.cpp
Leonardo de Moura 26097475fd Use fullpath in #include directives.
Signed-off-by: Leonardo de Moura <leonardo@microsoft.com>
2013-09-12 20:04:10 -07:00

353 lines
8.2 KiB
C++

/*
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 <limits.h>
#include "util/numerics/mpbq.h"
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;
if (is_zero()) {
m_k = 0;
return;
}
unsigned k = m_num.power_of_two_multiple();
if (k > m_k)
k = m_k;
div2k(m_num, m_num, k);
m_k -= k;
}
int cmp(mpbq const & a, mpbq const & b) {
static thread_local mpz tmp;
if (a.m_k == b.m_k) {
return cmp(a.m_num, b.m_num);
} else if (a.m_k < b.m_k) {
mul2k(tmp, a.m_num, b.m_k - a.m_k);
return cmp(tmp, b.m_num);
} else {
lean_assert(a.m_k > b.m_k);
mul2k(tmp, b.m_num, a.m_k - b.m_k);
return cmp(a.m_num, tmp);
}
}
int cmp(mpbq const & a, mpz const & b) {
static thread_local mpz tmp;
if (a.m_k == 0) {
return cmp(a.m_num, b);
} else {
mul2k(tmp, b, a.m_k);
return cmp(a.m_num, tmp);
}
}
int cmp(mpbq const & a, mpq const & b) {
if (a.is_integer() && b.is_integer()) {
return -cmp(b, a.m_num);
} else {
static thread_local mpz tmp1;
static thread_local mpz tmp2;
// tmp1 <- numerator(a)*denominator(b)
denominator(tmp1, b); tmp1 *= a.m_num;
// tmp2 <- numerator(b)*denominator(a)
numerator(tmp2, b); mul2k(tmp2, tmp2, a.m_k);
return cmp(tmp1, tmp2);
}
}
mpbq & mpbq::operator+=(mpbq const & a) {
if (m_k == a.m_k) {
m_num += a.m_num;
} else if (m_k < a.m_k) {
mul2k(m_num, m_num, a.m_k - m_k);
m_k = a.m_k;
m_num += a.m_num;
} else {
lean_assert(m_k > a.m_k);
static thread_local mpz tmp;
mul2k(tmp, a.m_num, m_k - a.m_k);
m_num += tmp;
}
normalize();
return *this;
}
template<typename T>
mpbq & mpbq::add_int(T const & a) {
if (m_k == 0) {
m_num += a;
} else {
lean_assert(m_k > 0);
static thread_local mpz tmp;
tmp = a;
mul2k(tmp, tmp, m_k);
m_num += tmp;
}
normalize();
return *this;
}
mpbq & mpbq::operator+=(unsigned a) { return add_int<unsigned>(a); }
mpbq & mpbq::operator+=(int a) { return add_int<int>(a); }
mpbq & mpbq::operator-=(mpbq const & a) {
if (m_k == a.m_k) {
m_num -= a.m_num;
} else if (m_k < a.m_k) {
mul2k(m_num, m_num, a.m_k - m_k);
m_k = a.m_k;
m_num -= a.m_num;
} else {
lean_assert(m_k > a.m_k);
static thread_local mpz tmp;
mul2k(tmp, a.m_num, m_k - a.m_k);
m_num -= tmp;
}
normalize();
return *this;
}
template<typename T>
mpbq & mpbq::sub_int(T const & a) {
if (m_k == 0) {
m_num -= a;
} else {
lean_assert(m_k > 0);
static thread_local mpz tmp;
tmp = a;
mul2k(tmp, tmp, m_k);
m_num -= tmp;
}
normalize();
return *this;
}
mpbq & mpbq::operator-=(unsigned a) { return sub_int<unsigned>(a); }
mpbq & mpbq::operator-=(int a) { return sub_int<int>(a); }
mpbq & mpbq::operator*=(mpbq const & a) {
m_num *= a.m_num;
if (m_k == 0 || a.m_k == 0) {
m_k += a.m_k;
normalize();
} else {
m_k += a.m_k;
}
return *this;
}
template<typename T>
mpbq & mpbq::mul_int(T const & a) {
m_num *= a;
normalize();
return *this;
}
mpbq & mpbq::operator*=(unsigned a) { return mul_int<unsigned>(a); }
mpbq & mpbq::operator*=(int a) { return mul_int<int>(a); }
void power(mpbq & a, mpbq const & b, unsigned k) {
lean_assert(static_cast<unsigned long long>(k) * static_cast<unsigned long long>(b.m_k) <= static_cast<unsigned long long>(UINT_MAX));
// We don't need to normalize because:
// If b.m_k == 0, then b is an integer, and the result be an integer
// If b.m_k > 0, then b.m_num must be odd, and the (b.m_num)^k will also be odd
a.m_k = b.m_k * k;
power(a.m_num, b.m_num, k);
}
int mpbq::magnitude_lb() const {
int s = m_num.sgn();
if (s < 0) {
return m_num.mlog2() - m_k + 1;
} else if (s == 0) {
return 0;
} else {
lean_assert(s > 0);
return m_num.log2() - m_k;
}
}
int mpbq::magnitude_ub() const {
int s = m_num.sgn();
if (s < 0) {
return m_num.mlog2() - m_k;
} else if (s == 0) {
return 0;
} else {
lean_assert(s > 0);
return m_num.log2() - m_k + 1;
}
}
void mul2(mpbq & a) {
if (a.m_k == 0) {
mul2k(a.m_num, a.m_num, 1);
} else {
a.m_k--;
}
}
void mul2k(mpbq & a, unsigned k) {
if (k == 0)
return;
if (a.m_k < k) {
mul2k(a.m_num, a.m_num, k - a.m_k);
a.m_k = 0;
} else {
lean_assert(a.m_k >= k);
a.m_k -= k;
}
}
bool root_lower(mpbq & a, mpbq const & b, unsigned n) {
bool r = root(a.m_num, b.m_num, n);
if (!r)
--a.m_num;
if (b.m_k % n == 0) {
a.m_k = b.m_k / n;
a.normalize();
return r;
} else if (a.m_num.is_neg()) {
a.m_k = b.m_k / n;
a.normalize();
return false;
} else {
a.m_k = b.m_k / n;
a.m_k++;
a.normalize();
return false;
}
}
bool root_upper(mpbq & a, mpbq const & b, unsigned n) {
bool r = root(a.m_num, b.m_num, n);
if (b.m_k % n == 0) {
a.m_k = b.m_k / n;
a.normalize();
return r;
} else if (a.m_num.is_neg()) {
a.m_k = b.m_k / n;
a.m_k++;
a.normalize();
return false;
} else {
a.m_k = b.m_k / n;
a.normalize();
return false;
}
}
void refine_upper(mpq const & q, mpbq & l, mpbq & u) {
lean_assert(l < q && q < u);
lean_assert(!q.get_denominator().is_power_of_two());
mpbq mid;
while (true) {
mid = l + u;
div2(mid);
if (mid > q) {
swap(u, mid);
lean_assert(l < q && q < u);
return;
}
swap(l, mid);
}
}
void refine_lower(mpq const & q, mpbq & l, mpbq & u) {
lean_assert(l < q && q < u);
lean_assert(!q.get_denominator().is_power_of_two());
mpbq mid;
while (true) {
mid = l + u;
div2(mid);
if (mid < q) {
swap(l, mid);
lean_assert(l < q && q < u);
return;
}
swap(u, mid);
}
}
bool lt_1div2k(mpbq const & a, unsigned k) {
if (a.m_num.is_nonpos())
return true;
if (a.m_k <= k) {
// since a.m_num >= 1
return false;
} else {
lean_assert(a.m_k > k);
static thread_local mpz tmp;
tmp = 1;
mul2k(tmp, tmp, a.m_k - k);
return a.m_num < tmp;
}
}
std::ostream & operator<<(std::ostream & out, mpbq const & v) {
if (v.m_k == 0) {
out << v.m_num;
} else if (v.m_k == 1) {
out << v.m_num << "/2";
} else {
out << v.m_num << "/2^" << v.m_k;
}
return out;
}
void display_decimal(std::ostream & out, mpbq const & a, unsigned prec) {
if (a.is_integer()) {
out << a.m_num;
return;
} else {
mpz two_k;
mpz n1, v1;
if (a.is_neg())
out << "-";
v1 = abs(a.m_num);
power(two_k, mpz(2), a.m_k);
n1 = rem(v1, two_k);
v1 = v1/two_k;
lean_assert(!n1.is_zero());
out << v1;
out << ".";
for (unsigned i = 0; i < prec; i++) {
n1 *= 10;
v1 = n1/two_k;
n1 = rem(n1, two_k);
out << v1;
if (n1.is_zero())
return;
}
out << "?";
}
}
}
void print(lean::mpbq const & n) { std::cout << n << std::endl; }