z3-z3-4.13.0.src.util.mpz.cpp Maven / Gradle / Ivy
The newest version!
/*++
Copyright (c) 2006 Microsoft Corporation
Module Name:
mpz.cpp
Abstract:
Author:
Leonardo de Moura (leonardo) 2010-06-17.
Revision History:
--*/
#include
#include
#include
#include "util/mpz.h"
#include "util/buffer.h"
#include "util/trace.h"
#include "util/hash.h"
#include "util/bit_util.h"
#if defined(_MP_INTERNAL)
#include "util/mpn.h"
#elif defined(_MP_GMP)
#include
#else
#error No multi-precision library selected.
#endif
// Available GCD algorithms
// #define EUCLID_GCD
// #define BINARY_GCD
// #define LS_BINARY_GCD
// #define LEHMER_GCD
#if defined(_MP_GMP)
// Use LEHMER only if not using GMP
#define EUCLID_GCD
#else
#define LEHMER_GCD
#endif
#if defined(__GNUC__)
#define _trailing_zeros32(X) __builtin_ctz(X)
#elif defined(_WINDOWS) && (defined(_M_X86) || (defined(_M_X64) && !defined(_M_ARM64EC)))
// This is needed for _tzcnt_u32 and friends.
#include
#define _trailing_zeros32(X) _tzcnt_u32(X)
#else
static uint32_t _trailing_zeros32(uint32_t x) {
uint32_t r = 0;
for (; 0 == (x & 1) && r < 32; ++r, x >>= 1);
return r;
}
#endif
#if (defined(__LP64__) || defined(_WIN64)) && defined(_M_X64) && !defined(_M_ARM64EC)
#if defined(__GNUC__)
#define _trailing_zeros64(X) __builtin_ctzll(X)
#else
#define _trailing_zeros64(X) _tzcnt_u64(X)
#endif
#else
static uint64_t _trailing_zeros64(uint64_t x) {
uint64_t r = 0;
for (; 0 == (x & 1) && r < 64; ++r, x >>= 1);
return r;
}
#endif
unsigned trailing_zeros(uint32_t x) {
return static_cast(_trailing_zeros32(x));
}
unsigned trailing_zeros(uint64_t x) {
return static_cast(_trailing_zeros64(x));
}
#define _bit_min(x, y) (y + ((x - y) & ((int)(x - y) >> 31)))
#define _bit_max(x, y) (x - ((x - y) & ((int)(x - y) >> 31)))
unsigned u_gcd(unsigned u, unsigned v) {
if (u == 0) return v;
if (v == 0) return u;
unsigned shift = _trailing_zeros32(u | v);
u >>= _trailing_zeros32(u);
if (u == 1 || v == 1) return 1 << shift;
if (u == v) return u << shift;
do {
v >>= _trailing_zeros32(v);
unsigned diff = u - v;
unsigned mdiff = diff & (unsigned)((int)diff >> 31);
u = v + mdiff; // min
v = diff - 2 * mdiff; // if v <= u: u - v, if v > u: v - u = u - v - 2 * (u - v)
}
while (v != 0);
return u << shift;
}
uint64_t u64_gcd(uint64_t u, uint64_t v) {
if (u == 0) return v;
if (v == 0) return u;
if (u == 1 || v == 1) return 1;
auto shift = _trailing_zeros64(u | v);
u >>= _trailing_zeros64(u);
do {
v >>= _trailing_zeros64(v);
if (u > v) std::swap(u, v);
v -= u;
}
while (v != 0);
return u << shift;
}
template
mpz_manager::mpz_manager():
m_allocator("mpz_manager") {
#ifndef _MP_GMP
if (sizeof(digit_t) == sizeof(uint64_t)) {
// 64-bit machine
m_init_cell_capacity = 4;
}
else {
m_init_cell_capacity = 6;
}
set(m_int_min, -static_cast(INT_MIN));
#else
// GMP
mpz_init(m_tmp);
mpz_init(m_tmp2);
mpz_init(m_two32);
mpz_set_ui(m_two32, UINT_MAX);
mpz_set_ui(m_tmp, 1);
mpz_add(m_two32, m_two32, m_tmp);
mpz_init(m_uint64_max);
unsigned max_l = static_cast(UINT64_MAX);
unsigned max_h = static_cast(UINT64_MAX >> 32);
mpz_set_ui(m_uint64_max, max_h);
mpz_mul(m_uint64_max, m_two32, m_uint64_max);
mpz_set_ui(m_tmp, max_l);
mpz_add(m_uint64_max, m_uint64_max, m_tmp);
mpz_init(m_int64_max);
mpz_init(m_int64_min);
max_l = static_cast(INT64_MAX % static_cast(UINT_MAX));
max_h = static_cast(INT64_MAX / static_cast(UINT_MAX));
mpz_set_ui(m_int64_max, max_h);
mpz_set_ui(m_tmp, UINT_MAX);
mpz_mul(m_int64_max, m_tmp, m_int64_max);
mpz_set_ui(m_tmp, max_l);
mpz_add(m_int64_max, m_tmp, m_int64_max);
mpz_neg(m_int64_min, m_int64_max);
mpz_sub_ui(m_int64_min, m_int64_min, 1);
#endif
mpz one(1);
set(m_two64, (uint64_t)UINT64_MAX);
add(m_two64, one, m_two64);
}
template
mpz_manager::~mpz_manager() {
del(m_two64);
#ifndef _MP_GMP
del(m_int_min);
#else
mpz_clear(m_tmp);
mpz_clear(m_tmp2);
mpz_clear(m_two32);
mpz_clear(m_uint64_max);
mpz_clear(m_int64_max);
mpz_clear(m_int64_min);
#endif
}
#ifndef _MP_GMP
template
mpz_cell * mpz_manager::allocate(unsigned capacity) {
SASSERT(capacity >= m_init_cell_capacity);
mpz_cell * cell;
#ifdef SINGLE_THREAD
cell = reinterpret_cast(m_allocator.allocate(cell_size(capacity)));
#else
if (SYNCH) {
cell = reinterpret_cast(memory::allocate(cell_size(capacity)));
}
else {
cell = reinterpret_cast(m_allocator.allocate(cell_size(capacity)));
}
#endif
cell->m_capacity = capacity;
return cell;
}
template
void mpz_manager::deallocate(bool is_heap, mpz_cell * ptr) {
if (is_heap) {
#ifdef SINGLE_THREAD
m_allocator.deallocate(cell_size(ptr->m_capacity), ptr);
#else
if (SYNCH) {
memory::deallocate(ptr);
}
else {
m_allocator.deallocate(cell_size(ptr->m_capacity), ptr);
}
#endif
}
}
template
mpz_manager::sign_cell::sign_cell(mpz_manager& m, mpz const& a):
m_local(reinterpret_cast(m_bytes)), m_a(a) {
m_local.m_ptr->m_capacity = capacity;
m.get_sign_cell(a, m_sign, m_cell, m_local.m_ptr);
}
#endif
template
void mpz_manager::del(mpz_manager* m, mpz & a) {
if (a.m_ptr) {
SASSERT(m);
m->deallocate(a.m_owner == mpz_self, a.m_ptr);
a.m_ptr = nullptr;
a.m_kind = mpz_small;
a.m_owner = mpz_self;
}
}
template
void mpz_manager::add(mpz const & a, mpz const & b, mpz & c) {
STRACE("mpz", tout << "[mpz] " << to_string(a) << " + " << to_string(b) << " == ";);
if (is_small(a) && is_small(b)) {
set_i64(c, i64(a) + i64(b));
}
else {
big_add(a, b, c);
}
STRACE("mpz", tout << to_string(c) << "\n";);
}
template
void mpz_manager::sub(mpz const & a, mpz const & b, mpz & c) {
STRACE("mpz", tout << "[mpz] " << to_string(a) << " - " << to_string(b) << " == ";);
if (is_small(a) && is_small(b)) {
set_i64(c, i64(a) - i64(b));
}
else {
big_sub(a, b, c);
}
STRACE("mpz", tout << to_string(c) << "\n";);
}
template
void mpz_manager::set_big_i64(mpz & c, int64_t v) {
#ifndef _MP_GMP
if (c.m_ptr == nullptr) {
c.m_ptr = allocate(m_init_cell_capacity);
c.m_owner = mpz_self;
}
c.m_kind = mpz_large;
SASSERT(capacity(c) >= m_init_cell_capacity);
uint64_t _v;
if (v == std::numeric_limits::min()) {
// min-int is even
_v = -(v/2);
c.m_val = -1;
}
else if (v < 0) {
_v = -v;
c.m_val = -1;
}
else {
_v = v;
c.m_val = 1;
}
if (sizeof(digit_t) == sizeof(uint64_t)) {
// 64-bit machine
digits(c)[0] = static_cast(_v);
c.m_ptr->m_size = 1;
}
else {
// 32-bit machine
digits(c)[0] = static_cast(_v);
digits(c)[1] = static_cast(_v >> 32);
c.m_ptr->m_size = digits(c)[1] == 0 ? 1 : 2;
}
#else
if (c.m_ptr == nullptr) {
c.m_ptr = allocate();
c.m_owner = mpz_self;
}
c.m_kind = mpz_large;
uint64_t _v;
bool sign = v < 0;
if (v == std::numeric_limits::min()) {
_v = -(v/2);
}
else if (v < 0) {
_v = -v;
}
else {
_v = v;
}
mpz_set_ui(*c.m_ptr, static_cast(_v));
MPZ_BEGIN_CRITICAL();
mpz_set_ui(m_tmp, static_cast(_v >> 32));
mpz_mul(m_tmp, m_tmp, m_two32);
mpz_add(*c.m_ptr, *c.m_ptr, m_tmp);
MPZ_END_CRITICAL();
if (sign)
mpz_neg(*c.m_ptr, *c.m_ptr);
#endif
if (v == std::numeric_limits::min()) {
big_add(c, c, c);
}
}
template
void mpz_manager::set_big_ui64(mpz & c, uint64_t v) {
#ifndef _MP_GMP
if (c.m_ptr == nullptr) {
c.m_ptr = allocate(m_init_cell_capacity);
c.m_owner = mpz_self;
}
c.m_kind = mpz_large;
SASSERT(capacity(c) >= m_init_cell_capacity);
c.m_val = 1;
if (sizeof(digit_t) == sizeof(uint64_t)) {
// 64-bit machine
digits(c)[0] = static_cast(v);
c.m_ptr->m_size = 1;
}
else {
// 32-bit machine
digits(c)[0] = static_cast(v);
digits(c)[1] = static_cast(v >> 32);
c.m_ptr->m_size = digits(c)[1] == 0 ? 1 : 2;
}
#else
if (c.m_ptr == nullptr) {
c.m_ptr = allocate();
c.m_owner = mpz_self;
}
c.m_kind = mpz_large;
mpz_set_ui(*c.m_ptr, static_cast(v));
MPZ_BEGIN_CRITICAL();
mpz_set_ui(m_tmp, static_cast(v >> 32));
mpz_mul(m_tmp, m_tmp, m_two32);
mpz_add(*c.m_ptr, *c.m_ptr, m_tmp);
MPZ_END_CRITICAL();
#endif
}
#ifdef _MP_GMP
template
mpz_manager::ensure_mpz_t::ensure_mpz_t(mpz const& a) {
if (is_small(a)) {
m_result = &m_local;
mpz_init(m_local);
mpz_set_si(m_local, a.m_val);
}
else {
m_result = a.m_ptr;
}
}
template
mpz_manager::ensure_mpz_t::~ensure_mpz_t() {
if (m_result == &m_local) {
mpz_clear(m_local);
}
}
#endif
#ifndef _MP_GMP
template
void mpz_manager::set(mpz_cell& src, mpz & a, int sign, unsigned sz) {
unsigned i = sz;
for (; i > 0 && src.m_digits[i-1] == 0; --i) ;
if (i == 0) {
// src is zero
set(a, 0);
return;
}
unsigned d = src.m_digits[0];
if (i == 1 && d <= INT_MAX) {
// src fits is a fixnum
a.m_val = sign < 0 ? -static_cast(d) : static_cast(d);
a.m_kind = mpz_small;
return;
}
set_digits(a, i, src.m_digits);
a.m_val = sign;
SASSERT(a.m_kind == mpz_large);
}
#endif
template
void mpz_manager::set(mpz & a, char const * val) {
set(a, 0);
mpz ten(10);
mpz tmp;
char const * str = val;
bool sign = false;
while (str[0] == ' ') ++str;
if (str[0] == '-')
sign = true;
while (str[0]) {
if ('0' <= str[0] && str[0] <= '9') {
SASSERT(str[0] - '0' <= 9);
mul(a, ten, tmp);
add(tmp, mk_z(str[0] - '0'), a);
}
++str;
}
del(tmp);
if (sign)
neg(a);
}
template
void mpz_manager::set_digits(mpz & target, unsigned sz, digit_t const * digits) {
// remove zero digits
while (sz > 0 && digits[sz - 1] == 0)
sz--;
if (sz == 0)
set(target, 0);
else if (sz == 1)
set(target, digits[0]);
else {
#ifndef _MP_GMP
target.m_val = 1; // number is positive.
if (target.m_ptr == nullptr) {
unsigned c = sz < m_init_cell_capacity ? m_init_cell_capacity : sz;
target.m_ptr = allocate(c);
target.m_ptr->m_size = sz;
target.m_ptr->m_capacity = c;
target.m_kind = mpz_large;
target.m_owner = mpz_self;
memcpy(target.m_ptr->m_digits, digits, sizeof(digit_t) * sz);
}
else if (capacity(target) < sz) {
SASSERT(sz > m_init_cell_capacity);
mpz_cell* ptr = allocate(sz);
memcpy(ptr->m_digits, digits, sizeof(digit_t) * sz);
ptr->m_size = sz;
ptr->m_capacity = sz;
deallocate(target);
target.m_val = 1;
target.m_ptr = ptr;
target.m_kind = mpz_large;
target.m_owner = mpz_self;
}
else {
target.m_ptr->m_size = sz;
if (target.m_ptr->m_digits != digits)
memcpy(target.m_ptr->m_digits, digits, sizeof(digit_t) * sz);
target.m_kind = mpz_large;
}
#else
mk_big(target);
// reset
mpz_set_ui(*target.m_ptr, digits[sz - 1]);
SASSERT(sz > 0);
unsigned i = sz - 1;
MPZ_BEGIN_CRITICAL();
while (i > 0) {
--i;
mpz_mul_2exp(*target.m_ptr, *target.m_ptr, 32);
mpz_set_ui(m_tmp, digits[i]);
mpz_add(*target.m_ptr, *target.m_ptr, m_tmp);
}
MPZ_END_CRITICAL();
#endif
}
}
template
void mpz_manager::mul(mpz const & a, mpz const & b, mpz & c) {
STRACE("mpz", tout << "[mpz] " << to_string(a) << " * " << to_string(b) << " == ";);
if (is_small(a) && is_small(b)) {
set_i64(c, i64(a) * i64(b));
}
else {
big_mul(a, b, c);
}
STRACE("mpz", tout << to_string(c) << "\n";);
}
// d <- a + b*c
template
void mpz_manager::addmul(mpz const & a, mpz const & b, mpz const & c, mpz & d) {
if (is_one(b)) {
add(a, c, d);
}
else if (is_minus_one(b)) {
sub(a, c, d);
}
else {
mpz tmp;
mul(b,c,tmp);
add(a,tmp,d);
del(tmp);
}
}
// d <- a - b*c
template
void mpz_manager::submul(mpz const & a, mpz const & b, mpz const & c, mpz & d) {
if (is_one(b)) {
sub(a, c, d);
}
else if (is_minus_one(b)) {
add(a, c, d);
}
else {
mpz tmp;
mul(b,c,tmp);
sub(a,tmp,d);
del(tmp);
}
}
template
void mpz_manager::machine_div_rem(mpz const & a, mpz const & b, mpz & q, mpz & r) {
STRACE("mpz", tout << "[mpz-ext] divrem(" << to_string(a) << ", " << to_string(b) << ") == ";);
if (is_small(a) && is_small(b)) {
int64_t _a = i64(a);
int64_t _b = i64(b);
set_i64(q, _a / _b);
set_i64(r, _a % _b);
}
else {
big_div_rem(a, b, q, r);
}
STRACE("mpz", tout << "(" << to_string(q) << ", " << to_string(r) << ")\n";);
}
template
void mpz_manager::machine_div(mpz const & a, mpz const & b, mpz & c) {
STRACE("mpz", tout << "[mpz-ext] machine-div(" << to_string(a) << ", " << to_string(b) << ") == ";);
if (is_small(b) && i64(b) == 0)
throw default_exception("division by 0");
if (is_small(a) && is_small(b))
set_i64(c, i64(a) / i64(b));
else
big_div(a, b, c);
STRACE("mpz", tout << to_string(c) << "\n";);
}
template
void mpz_manager::reset(mpz & a) {
deallocate(a);
set(a, 0);
}
template
void mpz_manager::rem(mpz const & a, mpz const & b, mpz & c) {
STRACE("mpz", tout << "[mpz-ext] rem(" << to_string(a) << ", " << to_string(b) << ") == ";);
if (is_small(a) && is_small(b)) {
set_i64(c, i64(a) % i64(b));
}
else {
big_rem(a, b, c);
}
STRACE("mpz", tout << to_string(c) << "\n";);
}
template
void mpz_manager::div_gcd(mpz const& a, mpz const& b, mpz & c) {
STRACE("mpz", tout << "[mpz-ext] div(" << to_string(a) << ", " << to_string(b) << ") == ";);
if (is_one(b)) {
set(c, a);
}
else {
machine_div(a, b, c);
}
STRACE("mpz", tout << to_string(c) << "\n";);
}
template
void mpz_manager::div(mpz const & a, mpz const & b, mpz & c) {
STRACE("mpz", tout << "[mpz-ext] div(" << to_string(a) << ", " << to_string(b) << ") == ";);
if (is_one(b)) {
set(c, a);
}
else if (is_neg(a)) {
mpz tmp;
machine_div_rem(a, b, c, tmp);
if (!is_zero(tmp)) {
if (is_neg(b))
add(c, mk_z(1), c);
else
sub(c, mk_z(1), c);
}
del(tmp);
}
else {
machine_div(a, b, c);
}
STRACE("mpz", tout << to_string(c) << "\n";);
}
template
void mpz_manager::mod(mpz const & a, mpz const & b, mpz & c) {
STRACE("mpz", tout << "[mpz-ext] mod(" << to_string(a) << ", " << to_string(b) << ") == ";);
rem(a, b, c);
if (is_neg(c)) {
if (is_pos(b))
add(c, b, c);
else
sub(c, b, c);
}
STRACE("mpz", tout << to_string(c) << "\n";);
}
template
void mpz_manager::neg(mpz & a) {
STRACE("mpz", tout << "[mpz] 0 - " << to_string(a) << " == ";);
if (is_small(a) && a.m_val == INT_MIN) {
// neg(INT_MIN) is not a small int
set_big_i64(a, - static_cast(INT_MIN));
return;
}
#ifndef _MP_GMP
a.m_val = -a.m_val;
#else
if (is_small(a)) {
a.m_val = -a.m_val;
}
else {
mpz_neg(*a.m_ptr, *a.m_ptr);
}
#endif
STRACE("mpz", tout << to_string(a) << "\n";);
}
template
void mpz_manager::abs(mpz & a) {
if (is_small(a)) {
if (a.m_val < 0) {
if (a.m_val == INT_MIN) {
// abs(INT_MIN) is not a small int
set_big_i64(a, - static_cast(INT_MIN));
}
else
a.m_val = -a.m_val;
}
}
else {
#ifndef _MP_GMP
a.m_val = 1;
#else
mpz_abs(*a.m_ptr, *a.m_ptr);
#endif
}
}
// TBD: replace use of 'tmp' by 'c'.
#ifndef _MP_GMP
template
template
void mpz_manager::big_add_sub(mpz const & a, mpz const & b, mpz & c) {
sign_cell ca(*this, a), cb(*this, b);
int sign_b = cb.sign();
mpz_stack tmp;
if (SUB)
sign_b = -sign_b;
unsigned real_sz;
if (ca.sign() == sign_b) {
unsigned sz = std::max(ca.cell()->m_size, cb.cell()->m_size)+1;
allocate_if_needed(tmp, sz);
m_mpn_manager.add(ca.cell()->m_digits, ca.cell()->m_size,
cb.cell()->m_digits, cb.cell()->m_size,
tmp.m_ptr->m_digits, sz, &real_sz);
SASSERT(real_sz <= sz);
set(*tmp.m_ptr, c, ca.sign(), real_sz);
}
else {
digit_t borrow;
int r = m_mpn_manager.compare(ca.cell()->m_digits, ca.cell()->m_size,
cb.cell()->m_digits, cb.cell()->m_size);
if (r == 0) {
set(c, 0);
}
else if (r < 0) {
// a < b
unsigned sz = cb.cell()->m_size;
allocate_if_needed(tmp, sz);
m_mpn_manager.sub(cb.cell()->m_digits,
cb.cell()->m_size,
ca.cell()->m_digits,
ca.cell()->m_size,
tmp.m_ptr->m_digits,
&borrow);
SASSERT(borrow == 0);
set(*tmp.m_ptr, c, sign_b, sz);
}
else {
// a > b
unsigned sz = ca.cell()->m_size;
allocate_if_needed(tmp, sz);
m_mpn_manager.sub(ca.cell()->m_digits,
ca.cell()->m_size,
cb.cell()->m_digits,
cb.cell()->m_size,
tmp.m_ptr->m_digits,
&borrow);
SASSERT(borrow == 0);
set(*tmp.m_ptr, c, ca.sign(), sz);
}
}
del(tmp);
}
#endif
template
void mpz_manager::big_add(mpz const & a, mpz const & b, mpz & c) {
#ifndef _MP_GMP
big_add_sub(a, b, c);
#else
// GMP version
ensure_mpz_t a1(a), b1(b);
mk_big(c);
mpz_add(*c.m_ptr, a1(), b1());
#endif
}
template
void mpz_manager::big_sub(mpz const & a, mpz const & b, mpz & c) {
#ifndef _MP_GMP
big_add_sub(a, b, c);
#else
// GMP version
ensure_mpz_t a1(a), b1(b);
mk_big(c);
mpz_sub(*c.m_ptr, a1(), b1());
#endif
}
template
void mpz_manager::big_mul(mpz const & a, mpz const & b, mpz & c) {
#ifndef _MP_GMP
// TBD replace tmp by c.
mpz_stack tmp;
sign_cell ca(*this, a), cb(*this, b);
unsigned sz = ca.cell()->m_size + cb.cell()->m_size;
allocate_if_needed(tmp, sz);
m_mpn_manager.mul(ca.cell()->m_digits,
ca.cell()->m_size,
cb.cell()->m_digits,
cb.cell()->m_size,
tmp.m_ptr->m_digits);
set(*tmp.m_ptr, c, ca.sign() == cb.sign() ? 1 : -1, sz);
del(tmp);
#else
// GMP version
ensure_mpz_t a1(a), b1(b);
mk_big(c);
mpz_mul(*c.m_ptr, a1(), b1());
#endif
}
template
void mpz_manager::big_div_rem(mpz const & a, mpz const & b, mpz & q, mpz & r) {
#ifndef _MP_GMP
quot_rem_core(a, b, q, r);
#else
// GMP version
ensure_mpz_t a1(a), b1(b);
mk_big(q);
mk_big(r);
mpz_tdiv_qr(*q.m_ptr, *r.m_ptr, a1(), b1());
#endif
}
#ifndef _MP_GMP
template
template
void mpz_manager::quot_rem_core(mpz const & a, mpz const & b, mpz & q, mpz & r)
{
/*
+26 / +7 = +3, remainder is +5
-26 / +7 = -3, remainder is -5
+26 / -7 = -3, remainder is +5
-26 / -7 = +3, remainder is -5
*/
mpz_stack q1, r1;
sign_cell ca(*this, a), cb(*this, b);
if (cb.cell()->m_size > ca.cell()->m_size) {
if (MODE == REM_ONLY || MODE == QUOT_AND_REM)
set(r, a);
if (MODE == QUOT_ONLY || MODE == QUOT_AND_REM)
set(q, 0);
return;
}
unsigned q_sz = ca.cell()->m_size - cb.cell()->m_size + 1;
unsigned r_sz = cb.cell()->m_size;
allocate_if_needed(q1, q_sz);
allocate_if_needed(r1, r_sz);
m_mpn_manager.div(ca.cell()->m_digits, ca.cell()->m_size,
cb.cell()->m_digits, cb.cell()->m_size,
q1.m_ptr->m_digits,
r1.m_ptr->m_digits);
if (MODE == QUOT_ONLY || MODE == QUOT_AND_REM)
set(*q1.m_ptr, q, ca.sign() == cb.sign() ? 1 : -1, q_sz);
if (MODE == REM_ONLY || MODE == QUOT_AND_REM)
set(*r1.m_ptr, r, ca.sign(), r_sz);
del(q1);
del(r1);
}
#endif
template
void mpz_manager::big_div(mpz const & a, mpz const & b, mpz & c) {
#ifndef _MP_GMP
mpz dummy;
quot_rem_core(a, b, c, dummy);
SASSERT(is_zero(dummy));
del(dummy);
#else
// GMP version
ensure_mpz_t a1(a), b1(b);
mk_big(c);
mpz_tdiv_q(*c.m_ptr, a1(), b1());
#endif
}
template
void mpz_manager::big_rem(mpz const & a, mpz const & b, mpz & c) {
#ifndef _MP_GMP
mpz dummy;
quot_rem_core(a, b, dummy, c);
SASSERT(is_zero(dummy));
del(dummy);
#else
// GMP version
ensure_mpz_t a1(a), b1(b);
mk_big(c);
mpz_tdiv_r(*c.m_ptr, a1(), b1());
#endif
}
template
void mpz_manager::gcd(mpz const & a, mpz const & b, mpz & c) {
static_assert(sizeof(a.m_val) == sizeof(int), "size mismatch");
static_assert(sizeof(mpz) <= 16, "mpz size overflow");
if (is_small(a) && is_small(b) && a.m_val != INT_MIN && b.m_val != INT_MIN) {
int _a = a.m_val;
int _b = b.m_val;
if (_a < 0) _a = -_a;
if (_b < 0) _b = -_b;
unsigned r = u_gcd(_a, _b);
set(c, r);
}
else {
#ifdef _MP_GMP
ensure_mpz_t a1(a), b1(b);
mk_big(c);
mpz_gcd(*c.m_ptr, a1(), b1());
return;
#endif
if (is_zero(a)) {
set(c, b);
abs(c);
return;
}
if (is_zero(b)) {
set(c, a);
abs(c);
return;
}
#ifdef BINARY_GCD
// Binary GCD for big numbers
// - It doesn't use division
// - The initial experiments, don't show any performance improvement
// - It only works with _MP_INTERNAL
mpz u, v, diff;
set(u, a);
set(v, b);
abs(u);
abs(v);
unsigned k_u = power_of_two_multiple(u);
unsigned k_v = power_of_two_multiple(v);
unsigned k = k_u < k_v ? k_u : k_v;
machine_div2k(u, k_u);
while (true) {
machine_div2k(v, k_v);
if (lt(u, v)) {
sub(v, u, v);
}
else {
sub(u, v, diff);
swap(u, v);
swap(v, diff);
}
if (is_zero(v) || is_one(v))
break;
// reset least significant bit
if (is_small(v))
v.m_val &= ~1;
else
v.m_ptr->m_digits[0] &= ~static_cast(1);
k_v = power_of_two_multiple(v);
}
mul2k(u, k, c);
del(u); del(v); del(diff);
#endif // BINARY_GCD
#ifdef EUCLID_GCD
mpz tmp1;
mpz tmp2;
mpz aux;
set(tmp1, a);
set(tmp2, b);
abs(tmp1);
abs(tmp2);
if (lt(tmp1, tmp2))
swap(tmp1, tmp2);
if (is_zero(tmp2)) {
swap(c, tmp1);
}
else {
while (true) {
if (is_uint64(tmp1) && is_uint64(tmp2)) {
set(c, u64_gcd(get_uint64(tmp1), get_uint64(tmp2)));
break;
}
rem(tmp1, tmp2, aux);
if (is_zero(aux)) {
swap(c, tmp2);
break;
}
swap(tmp1, tmp2);
swap(tmp2, aux);
}
}
del(tmp1); del(tmp2); del(aux);
#endif // EUCLID_GCD
#ifdef LS_BINARY_GCD
mpz u, v, t, u1, u2;
set(u, a);
set(v, b);
abs(u);
abs(v);
if (lt(u, v))
swap(u, v);
while (!is_zero(v)) {
// Basic idea:
// compute t = 2^e*v such that t <= u < 2t
// u := min{u - t, 2t - u}
//
// The assignment u := min{u - t, 2t - u}
// can be replaced with u := u - t
//
// Since u and v are positive, we have:
// 2^{log2(u)} <= u < 2^{(log2(u) + 1)}
// 2^{log2(v)} <= v < 2^{(log2(v) + 1)}
// -->
// 2^{log2(v)}*2^{log2(u)-log2(v)} <= v*2^{log2(u)-log2(v)} < 2^{log2(v) + 1}*2^{log2(u)-log2(v)}
// -->
// 2^{log2(u)} <= v*2^{log2(u)-log2(v)} < 2^{log2(u) + 1}
//
// Now, let t be v*2^{log2(u)-log2(v)}
// If t <= u, then we found t
// Otherwise t = t div 2
unsigned k_u = log2(u);
unsigned k_v = log2(v);
SASSERT(k_v <= k_u);
unsigned e = k_u - k_v;
mul2k(v, e, t);
sub(u, t, u1);
if (is_neg(u1)) {
// t is too big
machine_div2k(t, 1);
// Now, u1 contains u - 2t
neg(u1);
// Now, u1 contains 2t - u
sub(u, t, u2); // u2 := u - t
}
else {
// u1 contains u - t
mul2k(t, 1);
sub(t, u, u2);
// u2 contains 2t - u
}
SASSERT(is_nonneg(u1));
SASSERT(is_nonneg(u2));
if (lt(u1, u2))
swap(u, u1);
else
swap(u, u2);
if (lt(u, v))
swap(u,v);
}
swap(u, c);
del(u); del(v); del(t); del(u1); del(u2);
#endif // LS_BINARY_GCD
#ifdef LEHMER_GCD
// For now, it only works if sizeof(digit_t) == sizeof(unsigned)
static_assert(sizeof(digit_t) == sizeof(unsigned), "");
int64_t a_hat, b_hat, A, B, C, D, T, q, a_sz, b_sz;
mpz a1, b1, t, r, tmp;
set(a1, a);
set(b1, b);
abs(a1);
abs(b1);
if (lt(a1, b1))
swap(a1, b1);
while (true) {
SASSERT(ge(a1, b1));
if (is_small(b1)) {
if (is_small(a1)) {
unsigned r = u_gcd(a1.m_val, b1.m_val);
set(c, r);
break;
}
else {
while (!is_zero(b1)) {
SASSERT(ge(a1, b1));
rem(a1, b1, tmp);
swap(a1, b1);
swap(b1, tmp);
}
swap(c, a1);
break;
}
}
SASSERT(!is_small(a1));
SASSERT(!is_small(b1));
a_sz = a1.m_ptr->m_size;
b_sz = b1.m_ptr->m_size;
SASSERT(b_sz <= a_sz);
a_hat = a1.m_ptr->m_digits[a_sz - 1];
b_hat = (b_sz == a_sz) ? b1.m_ptr->m_digits[b_sz - 1] : 0;
A = 1;
B = 0;
C = 0;
D = 1;
while (true) {
// Loop invariants
SASSERT(a_hat + A <= static_cast(UINT_MAX) + 1);
SASSERT(a_hat + B < static_cast(UINT_MAX) + 1);
SASSERT(b_hat + C < static_cast(UINT_MAX) + 1);
SASSERT(b_hat + D <= static_cast(UINT_MAX) + 1);
// overflows can't happen since I'm using int64
if (b_hat + C == 0 || b_hat + D == 0)
break;
q = (a_hat + A)/(b_hat + C);
if (q != (a_hat + B)/(b_hat + D))
break;
T = A - q*C;
A = C;
C = T;
T = B - q*D;
B = D;
D = T;
T = a_hat - q*b_hat;
a_hat = b_hat;
b_hat = T;
}
SASSERT(ge(a1, b1));
if (B == 0) {
rem(a1, b1, t);
swap(a1, b1);
swap(b1, t);
SASSERT(ge(a1, b1));
}
else {
// t <- A*a1
set(tmp, A);
mul(a1, tmp, t);
// t <- t + B*b1
set(tmp, B);
addmul(t, tmp, b1, t);
// r <- C*a1
set(tmp, C);
mul(a1, tmp, r);
// r <- r + D*b1
set(tmp, D);
addmul(r, tmp, b1, r);
// a <- t
swap(a1, t);
// b <- r
swap(b1, r);
SASSERT(ge(a1, b1));
}
}
del(a1); del(b1); del(r); del(t); del(tmp);
#endif // LEHMER_GCD
}
}
template
unsigned mpz_manager::size_info(mpz const & a) {
if (is_small(a))
return 1;
#ifndef _MP_GMP
return a.m_ptr->m_size + 1;
#else
return mpz_size(*a.m_ptr);
#endif
}
template
struct mpz_manager::sz_lt {
mpz_manager & m;
mpz const * m_as;
sz_lt(mpz_manager & _m, mpz const * as):m(_m), m_as(as) {}
bool operator()(unsigned p1, unsigned p2) {
return m.size_info(m_as[p1]) < m.size_info(m_as[p2]);
}
};
template
void mpz_manager::gcd(unsigned sz, mpz const * as, mpz & g) {
#if 0
// Optimization: sort numbers by size. Motivation: compute the gcd of the small ones first.
// The optimization did not really help.
switch (sz) {
case 0:
set(g, 0);
return;
case 1:
set(g, as[0]);
abs(g);
return;
case 2:
gcd(as[0], as[1], g);
return;
default:
break;
}
unsigned i;
for (i = 0; i < sz; i++) {
if (!is_small(as[i]))
break;
}
if (i != sz) {
// array has big numbers
sbuffer p;
for (i = 0; i < sz; i++)
p.push_back(i);
sz_lt lt(*this, as);
std::sort(p.begin(), p.end(), lt);
TRACE("mpz_gcd", for (unsigned i = 0; i < sz; i++) tout << p[i] << ":" << size_info(as[p[i]]) << " "; tout << "\n";);
gcd(as[p[0]], as[p[1]], g);
for (i = 2; i < sz; i++) {
if (is_one(g))
return;
gcd(g, as[p[i]], g);
}
return;
}
else {
gcd(as[0], as[1], g);
for (unsigned i = 2; i < sz; i++) {
if (is_one(g))
return;
gcd(g, as[i], g);
}
}
#else
// Vanilla implementation
switch (sz) {
case 0:
set(g, 0);
return;
case 1:
set(g, as[0]);
abs(g);
return;
default:
break;
}
gcd(as[0], as[1], g);
for (unsigned i = 2; i < sz; i++) {
if (is_one(g))
return;
gcd(g, as[i], g);
}
#endif
}
template
void mpz_manager::gcd(mpz const & r1, mpz const & r2, mpz & a, mpz & b, mpz & r) {
mpz tmp1, tmp2;
mpz aux, quot;
set(tmp1, r1);
set(tmp2, r2);
set(a, 1);
set(b, 0);
mpz nexta, nextb;
set(nexta, 0);
set(nextb, 1);
abs(tmp1);
abs(tmp2);
if (lt(tmp1, tmp2)) {
swap(tmp1, tmp2);
swap(nexta, nextb);
swap(a, b);
}
// tmp1 >= tmp2 >= 0
// quot_rem in one function would be faster.
while (is_pos(tmp2)) {
SASSERT(ge(tmp1, tmp2));
// aux = tmp2
set(aux, tmp2);
// quot = div(tmp1, tmp2);
machine_div(tmp1, tmp2, quot);
// tmp2 = tmp1 % tmp2
rem(tmp1, tmp2, tmp2);
// tmp1 = aux
set(tmp1, aux);
// aux = nexta
set(aux, nexta);
// nexta = a - (quot*nexta)
mul(quot, nexta, nexta);
sub(a, nexta, nexta);
// a = axu
set(a, aux);
// aux = nextb
set(aux, nextb);
// nextb = b - (quot*nextb)
mul(nextb, quot, nextb);
sub(b, nextb, nextb);
// b = aux
set(b, aux);
}
if (is_neg(r1))
neg(a);
if (is_neg(r2))
neg(b);
// SASSERT((a*r1) + (b*r2) == tmp1);
#ifdef Z3DEBUG
mul(a, r1, nexta);
mul(b, r2, nextb);
add(nexta, nextb, nexta);
SASSERT(eq(nexta, tmp1));
#endif
set(r, tmp1);
del(tmp1);
del(tmp2);
del(aux);
del(quot);
del(nexta);
del(nextb);
}
template
void mpz_manager::lcm(mpz const & a, mpz const & b, mpz & c) {
if (is_one(b)) {
set(c, a);
TRACE("lcm_bug", tout << "1. lcm(" << to_string(a) << ", " << to_string(b) << ") = " << to_string(c) << "\n";);
}
else if (is_one(a) || eq(a, b)) {
set(c, b);
TRACE("lcm_bug", tout << "2. lcm(" << to_string(a) << ", " << to_string(b) << ") = " << to_string(c) << "\n";);
}
else {
mpz r;
gcd(a, b, r);
TRACE("lcm_bug", tout << "gcd(" << to_string(a) << ", " << to_string(b) << ") = " << to_string(r) << "\n";);
if (eq(r, a)) {
set(c, b);
TRACE("lcm_bug", tout << "3. lcm(" << to_string(a) << ", " << to_string(b) << ") = " << to_string(c) << "\n";);
}
else if (eq(r, b)) {
set(c, a);
TRACE("lcm_bug", tout << "4. lcm(" << to_string(a) << ", " << to_string(b) << ") = " << to_string(c) << "\n";);
}
else {
// c contains gcd(a, b)
// so c divides a, and machine_div(a, c) is equal to div(a, c)
machine_div(a, r, r);
mul(r, b, c);
TRACE("lcm_bug", tout << "5. lcm(" << to_string(a) << ", " << to_string(b) << ") = " << to_string(c) << "\n";);
}
del(r);
}
}
template
void mpz_manager::bitwise_or(mpz const & a, mpz const & b, mpz & c) {
SASSERT(is_nonneg(a));
SASSERT(is_nonneg(b));
TRACE("mpz", tout << "is_small(a): " << is_small(a) << ", is_small(b): " << is_small(b) << "\n";);
if (is_small(a) && is_small(b)) {
c.m_val = a.m_val | b.m_val;
c.m_kind = mpz_small;
}
else {
#ifndef _MP_GMP
mpz a1, b1, a2, b2, m, tmp;
set(a1, a);
set(b1, b);
set(m, 1);
set(c, 0);
while (!is_zero(a1) && !is_zero(b1)) {
TRACE("mpz", tout << "a1: " << to_string(a1) << ", b1: " << to_string(b1) << "\n";);
mod(a1, m_two64, a2);
mod(b1, m_two64, b2);
TRACE("mpz", tout << "a2: " << to_string(a2) << ", b2: " << to_string(b2) << "\n";);
uint64_t v = get_uint64(a2) | get_uint64(b2);
TRACE("mpz", tout << "uint(a2): " << get_uint64(a2) << ", uint(b2): " << get_uint64(b2) << "\n";);
set(tmp, v);
mul(tmp, m, tmp);
add(c, tmp, c); // c += m * v
mul(m, m_two64, m);
div(a1, m_two64, a1);
div(b1, m_two64, b1);
}
if (!is_zero(a1)) {
mul(a1, m, a1);
add(c, a1, c);
}
if (!is_zero(b1)) {
mul(b1, m, b1);
add(c, b1, c);
}
del(a1); del(b1); del(a2); del(b2); del(m); del(tmp);
#else
ensure_mpz_t a1(a), b1(b);
mk_big(c);
mpz_ior(*c.m_ptr, a1(), b1());
#endif
}
}
template
void mpz_manager::bitwise_and(mpz const & a, mpz const & b, mpz & c) {
SASSERT(is_nonneg(a));
SASSERT(is_nonneg(b));
if (is_small(a) && is_small(b)) {
c.m_val = a.m_val & b.m_val;
c.m_kind = mpz_small;
}
else {
#ifndef _MP_GMP
mpz a1, b1, a2, b2, m, tmp;
set(a1, a);
set(b1, b);
set(m, 1);
set(c, 0);
while (!is_zero(a1) && !is_zero(b1)) {
mod(a1, m_two64, a2);
mod(b1, m_two64, b2);
uint64_t v = get_uint64(a2) & get_uint64(b2);
set(tmp, v);
mul(tmp, m, tmp);
add(c, tmp, c); // c += m * v
mul(m, m_two64, m);
div(a1, m_two64, a1);
div(b1, m_two64, b1);
}
del(a1); del(b1); del(a2); del(b2); del(m); del(tmp);
#else
ensure_mpz_t a1(a), b1(b);
mk_big(c);
mpz_and(*c.m_ptr, a1(), b1());
#endif
}
}
template
void mpz_manager::bitwise_xor(mpz const & a, mpz const & b, mpz & c) {
SASSERT(is_nonneg(a));
SASSERT(is_nonneg(b));
if (is_small(a) && is_small(b)) {
set_i64(c, i64(a) ^ i64(b));
}
else {
#ifndef _MP_GMP
mpz a1, b1, a2, b2, m, tmp;
set(a1, a);
set(b1, b);
set(m, 1);
set(c, 0);
while (!is_zero(a1) && !is_zero(b1)) {
mod(a1, m_two64, a2);
mod(b1, m_two64, b2);
uint64_t v = get_uint64(a2) ^ get_uint64(b2);
set(tmp, v);
mul(tmp, m, tmp);
add(c, tmp, c); // c += m * v
mul(m, m_two64, m);
div(a1, m_two64, a1);
div(b1, m_two64, b1);
}
if (!is_zero(a1)) {
mul(a1, m, a1);
add(c, a1, c);
}
if (!is_zero(b1)) {
mul(b1, m, b1);
add(c, b1, c);
}
del(a1); del(b1); del(a2); del(b2); del(m); del(tmp);
#else
ensure_mpz_t a1(a), b1(b);
mk_big(c);
mpz_xor(*c.m_ptr, a1(), b1());
#endif
}
}
template
void mpz_manager::bitwise_not(unsigned sz, mpz const & a, mpz & c) {
SASSERT(is_nonneg(a));
if (is_small(a) && sz <= 64) {
uint64_t v = ~get_uint64(a);
unsigned zero_out = 64 - sz;
v = (v << zero_out) >> zero_out;
set(c, v);
}
else {
mpz a1, a2, m, tmp;
set(a1, a);
set(m, 1);
set(c, 0);
while (sz > 0) {
mod(a1, m_two64, a2);
uint64_t n = get_uint64(a2);
uint64_t v = ~n;
SASSERT(~v == n);
if (sz < 64) {
uint64_t mask = (1ull << static_cast(sz)) - 1ull;
v = mask & v;
}
TRACE("bitwise_not", tout << "sz: " << sz << ", v: " << v << ", n: " << n << "\n";);
set(tmp, v);
SASSERT(get_uint64(tmp) == v);
mul(tmp, m, tmp);
add(c, tmp, c); // c += m * v
mul(m, m_two64, m);
div(a1, m_two64, a1);
sz -= (sz<64) ? sz : 64;
}
del(a1); del(a2); del(m); del(tmp);
TRACE("bitwise_not", tout << "sz: " << sz << " a: " << to_string(a) << " c: " << to_string(c) << "\n";);
}
}
template
void mpz_manager::big_set(mpz & target, mpz const & source) {
#ifndef _MP_GMP
if (&target == &source)
return;
target.m_val = source.m_val;
if (target.m_ptr == nullptr) {
target.m_ptr = allocate(capacity(source));
target.m_ptr->m_size = size(source);
target.m_ptr->m_capacity = capacity(source);
target.m_kind = mpz_large;
target.m_owner = mpz_self;
memcpy(target.m_ptr->m_digits, source.m_ptr->m_digits, sizeof(digit_t) * size(source));
}
else if (capacity(target) < size(source)) {
deallocate(target);
target.m_ptr = allocate(capacity(source));
target.m_ptr->m_size = size(source);
target.m_ptr->m_capacity = capacity(source);
target.m_kind = mpz_large;
target.m_owner = mpz_self;
memcpy(target.m_ptr->m_digits, source.m_ptr->m_digits, sizeof(digit_t) * size(source));
}
else {
target.m_ptr->m_size = size(source);
memcpy(target.m_ptr->m_digits, source.m_ptr->m_digits, sizeof(digit_t) * size(source));
target.m_kind = mpz_large;
}
#else
// GMP version
mk_big(target);
mpz_set(*target.m_ptr, *source.m_ptr);
#endif
}
template
int mpz_manager::big_compare(mpz const & a, mpz const & b) {
#ifndef _MP_GMP
if (sign(a) > 0) {
// a is positive
if (sign(b) > 0) {
// a & b are positive
sign_cell ca(*this, a), cb(*this, b);
return m_mpn_manager.compare(ca.cell()->m_digits, ca.cell()->m_size,
cb.cell()->m_digits, cb.cell()->m_size);
}
else {
// b is negative
return 1; // a > b
}
}
else {
// a is negative
if (sign(b) > 0) {
// b is positive
return -1; // a < b
}
else {
// a & b are negative
sign_cell ca(*this, a), cb(*this, b);
return m_mpn_manager.compare(cb.cell()->m_digits, cb.cell()->m_size,
ca.cell()->m_digits, ca.cell()->m_size);
}
}
#else
// GMP version
ensure_mpz_t a1(a), b1(b);
return mpz_cmp(a1(), b1());
#endif
}
template
bool mpz_manager::is_uint64(mpz const & a) const {
#ifndef _MP_GMP
if (a.m_val < 0)
return false;
if (is_small(a))
return true;
if (sizeof(digit_t) == sizeof(uint64_t)) {
return size(a) <= 1;
}
else {
return size(a) <= 2;
}
#else
// GMP version
if (is_small(a))
return a.m_val >= 0;
return is_nonneg(a) && mpz_cmp(*a.m_ptr, m_uint64_max) <= 0;
#endif
}
template
bool mpz_manager::is_int64(mpz const & a) const {
if (is_small(a))
return true;
#ifndef _MP_GMP
if (!is_abs_uint64(a))
return false;
uint64_t num = big_abs_to_uint64(a);
uint64_t msb = static_cast(1) << 63;
uint64_t msb_val = msb & num;
if (a.m_val >= 0) {
// non-negative number.
return (0 == msb_val);
}
else {
// negative number.
// either the high bit is 0, or
// the number is 2^64 which can be represented.
//
return 0 == msb_val || (msb_val == num);
}
#else
// GMP version
return mpz_cmp(m_int64_min, *a.m_ptr) <= 0 && mpz_cmp(*a.m_ptr, m_int64_max) <= 0;
#endif
}
template
uint64_t mpz_manager::get_uint64(mpz const & a) const {
if (is_small(a))
return static_cast(a.m_val);
#ifndef _MP_GMP
SASSERT(a.m_ptr->m_size > 0);
return big_abs_to_uint64(a);
#else
// GMP version
if (sizeof(uint64_t) == sizeof(unsigned long)) {
return mpz_get_ui(*a.m_ptr);
}
else {
MPZ_BEGIN_CRITICAL();
mpz_manager * _this = const_cast(this);
mpz_set(_this->m_tmp, *a.m_ptr);
mpz_mod(_this->m_tmp, m_tmp, m_two32);
uint64_t r = static_cast(mpz_get_ui(m_tmp));
mpz_set(_this->m_tmp, *a.m_ptr);
mpz_div(_this->m_tmp, m_tmp, m_two32);
r += static_cast(mpz_get_ui(m_tmp)) << static_cast(32);
MPZ_END_CRITICAL();
return r;
}
#endif
}
template
int64_t mpz_manager::get_int64(mpz const & a) const {
if (is_small(a))
return static_cast(a.m_val);
#ifndef _MP_GMP
SASSERT(is_int64(a));
uint64_t num = big_abs_to_uint64(a);
if (a.m_val < 0) {
if (num != 0 && (num << 1) == 0)
return INT64_MIN;
return -static_cast(num);
}
return static_cast(num);
#else
// GMP
if (sizeof(int64_t) == sizeof(long) || mpz_fits_slong_p(*a.m_ptr)) {
return mpz_get_si(*a.m_ptr);
}
else {
MPZ_BEGIN_CRITICAL();
mpz_manager * _this = const_cast(this);
mpz_mod(_this->m_tmp, *a.m_ptr, m_two32);
int64_t r = static_cast(mpz_get_ui(m_tmp));
mpz_div(_this->m_tmp, *a.m_ptr, m_two32);
r += static_cast(mpz_get_si(m_tmp)) << static_cast(32);
MPZ_END_CRITICAL();
return r;
}
#endif
}
template
double mpz_manager::get_double(mpz const & a) const {
if (is_small(a))
return static_cast(a.m_val);
#ifndef _MP_GMP
double r = 0.0;
double d = 1.0;
unsigned sz = size(a);
for (unsigned i = 0; i < sz; i++) {
r += d * static_cast(digits(a)[i]);
if (sizeof(digit_t) == sizeof(uint64_t))
d *= (1.0 + static_cast(UINT64_MAX)); // 64-bit version, multiply by 2^64
else
d *= (1.0 + static_cast(UINT_MAX)); // 32-bit version, multiply by 2^32
}
if (!(r >= 0.0)) {
r = static_cast(UINT64_MAX); // some large number
}
return a.m_val < 0 ? -r : r;
#else
return mpz_get_d(*a.m_ptr);
#endif
}
template
void mpz_manager::display(std::ostream & out, mpz const & a) const {
if (is_small(a)) {
out << a.m_val;
}
else {
#ifndef _MP_GMP
if (a.m_val < 0)
out << "-";
if (sizeof(digit_t) == 4) {
sbuffer buffer(11*size(a), 0);
out << m_mpn_manager.to_string(digits(a), size(a), buffer.begin(), buffer.size());
}
else {
sbuffer buffer(21*size(a), 0);
out << m_mpn_manager.to_string(digits(a), size(a), buffer.begin(), buffer.size());
}
#else
// GMP version
size_t sz = mpz_sizeinbase(*a.m_ptr, 10) + 2;
sbuffer buffer(sz, 0);
mpz_get_str(buffer.data(), 10, *a.m_ptr);
out << buffer.data();
#endif
}
}
template
void mpz_manager::display_smt2(std::ostream & out, mpz const & a, bool decimal) const {
if (is_neg(a)) {
mpz_manager * _this = const_cast*>(this);
_scoped_numeral > tmp(*_this);
_this->set(tmp, a);
_this->neg(tmp);
out << "(- ";
display(out, tmp);
if (decimal)
out << ".0";
out << ")";
}
else {
display(out, a);
if (decimal)
out << ".0";
}
}
template
void mpz_manager::display_hex(std::ostream & out, mpz const & a, unsigned num_bits) const {
SASSERT(num_bits % 4 == 0);
std::ios fmt(nullptr);
fmt.copyfmt(out);
out << std::hex;
if (is_small(a)) {
out << std::setw(num_bits/4) << std::setfill('0') << get_uint64(a);
} else {
#ifndef _MP_GMP
digit_t *ds = digits(a);
unsigned sz = size(a);
unsigned bitSize = sz * sizeof(digit_t) * 8;
unsigned firstDigitSize;
if (num_bits >= bitSize) {
firstDigitSize = sizeof(digit_t) * 2;
for (unsigned i = 0; i < (num_bits - bitSize)/4; ++i) {
out << "0";
}
} else {
firstDigitSize = num_bits % (sizeof(digit_t) * 8) / 4;
}
out << std::setfill('0') << std::setw(firstDigitSize) << ds[sz-1] << std::setw(sizeof(digit_t)*2);
for (unsigned i = 1; i < sz; ++i) {
out << ds[sz-i-1];
}
#else
// GMP version
size_t sz = mpz_sizeinbase(*(a.m_ptr), 16);
unsigned requiredLength = num_bits / 4;
unsigned padding = requiredLength > sz ? requiredLength - sz : 0;
sbuffer buffer(sz, 0);
mpz_get_str(buffer.data(), 16, *(a.m_ptr));
for (unsigned i = 0; i < padding; ++i) {
out << "0";
}
out << buffer.data() + (sz > requiredLength ? sz - requiredLength : 0);
#endif
}
out.copyfmt(fmt);
}
static void display_binary_data(std::ostream &out, uint64_t val, uint64_t numBits) {
for (uint64_t shift = numBits; shift-- > 64ull; ) out << "0";
if (numBits > 64) numBits = 64;
for (uint64_t shift = numBits; shift-- > 0; ) {
if (val & (1ull << shift)) {
out << "1";
} else {
out << "0";
}
}
}
template
void mpz_manager::display_bin(std::ostream & out, mpz const & a, unsigned num_bits) const {
if (is_small(a)) {
display_binary_data(out, get_uint64(a), num_bits);
}
else {
#ifndef _MP_GMP
digit_t *ds = digits(a);
unsigned sz = size(a);
const unsigned digitBitSize = sizeof(digit_t) * 8;
unsigned bitSize = sz * digitBitSize;
unsigned firstDigitLength;
if (num_bits > bitSize) {
firstDigitLength = 0;
for (unsigned i = 0; i < (num_bits - bitSize); ++i) {
out << "0";
}
} else {
firstDigitLength = num_bits % digitBitSize;
}
for (unsigned i = 0; i < sz; ++i) {
if (i == 0 && firstDigitLength != 0) {
display_binary_data(out, ds[sz-1], firstDigitLength);
} else {
display_binary_data(out, ds[sz-i-1], digitBitSize);
}
}
#else
// GMP version
size_t sz = mpz_sizeinbase(*(a.m_ptr), 2);
unsigned padding = num_bits > sz ? num_bits - sz : 0;
sbuffer buffer(sz, 0);
mpz_get_str(buffer.data(), 2, *(a.m_ptr));
for (unsigned i = 0; i < padding; ++i) {
out << "0";
}
out << buffer.data() + (sz > num_bits ? sz - num_bits : 0);
#endif
}
}
template
std::string mpz_manager::to_string(mpz const & a) const {
std::ostringstream buffer;
display(buffer, a);
return buffer.str();
}
template
unsigned mpz_manager::hash(mpz const & a) {
if (is_small(a))
return ::abs(a.m_val);
#ifndef _MP_GMP
unsigned sz = size(a);
if (sz == 1)
return static_cast(digits(a)[0]);
return string_hash(reinterpret_cast(digits(a)), sz * sizeof(digit_t), 17);
#else
return mpz_get_si(*a.m_ptr);
#endif
}
template
void mpz_manager::power(mpz const & a, unsigned p, mpz & b) {
#ifdef _MP_GMP
if (!is_small(a)) {
mk_big(b);
mpz_pow_ui(*b.m_ptr, *a.m_ptr, p);
return;
}
#endif
#ifndef _MP_GMP
if (is_small(a)) {
if (a.m_val == 2) {
if (p < 8 * sizeof(int) - 1) {
b.m_val = 1 << p;
b.m_kind = mpz_small;
}
else {
unsigned sz = p/(8 * sizeof(digit_t)) + 1;
unsigned shift = p%(8 * sizeof(digit_t));
SASSERT(sz > 0);
allocate_if_needed(b, sz);
SASSERT(b.m_ptr->m_capacity >= sz);
b.m_ptr->m_size = sz;
for (unsigned i = 0; i < sz - 1; i++)
b.m_ptr->m_digits[i] = 0;
b.m_ptr->m_digits[sz-1] = 1 << shift;
b.m_val = 1;
b.m_kind = mpz_large;
}
return;
}
if (a.m_val == 0) {
SASSERT(p != 0);
set(b, 0);
return;
}
if (a.m_val == 1) {
set(b, 1);
return;
}
}
#endif
// general purpose
unsigned mask = 1;
mpz power;
set(power, a);
set(b, 1);
while (mask <= p) {
if (mask & p)
mul(b, power, b);
mul(power, power, power);
mask = mask << 1;
}
del(power);
}
template
bool mpz_manager::is_power_of_two(mpz const & a) {
unsigned shift;
return is_power_of_two(a, shift);
}
template
bool mpz_manager::is_power_of_two(mpz const & a, unsigned & shift) {
if (is_nonpos(a))
return false;
if (is_small(a)) {
if (::is_power_of_two(a.m_val)) {
shift = ::log2((unsigned)a.m_val);
return true;
}
else {
return false;
}
}
#ifndef _MP_GMP
mpz_cell * c = a.m_ptr;
unsigned sz = c->m_size;
digit_t * ds = c->m_digits;
for (unsigned i = 0; i < sz - 1; i++) {
if (ds[i] != 0)
return false;
}
digit_t v = ds[sz-1];
if (!(v & (v - 1)) && v) {
shift = log2(a);
return true;
}
else {
return false;
}
#else
if (mpz_popcount(*a.m_ptr) == 1) {
shift = log2(a);
return true;
}
else {
return false;
}
#endif
}
// Expand capacity of a
#ifndef _MP_GMP
template
void mpz_manager::ensure_capacity(mpz & a, unsigned capacity) {
if (capacity <= 1)
return;
if (capacity < m_init_cell_capacity)
capacity = m_init_cell_capacity;
if (is_small(a)) {
int val = a.m_val;
allocate_if_needed(a, capacity);
a.m_kind = mpz_large;
SASSERT(a.m_ptr->m_capacity >= capacity);
if (val == INT_MIN) {
unsigned intmin_sz = m_int_min.m_ptr->m_size;
for (unsigned i = 0; i < intmin_sz; i++)
a.m_ptr->m_digits[i] = m_int_min.m_ptr->m_digits[i];
a.m_val = -1;
a.m_ptr->m_size = m_int_min.m_ptr->m_size;
}
else if (val < 0) {
a.m_ptr->m_digits[0] = -val;
a.m_val = -1;
a.m_ptr->m_size = 1;
}
else {
a.m_ptr->m_digits[0] = val;
a.m_val = 1;
a.m_ptr->m_size = 1;
}
}
else if (a.m_ptr->m_capacity < capacity) {
mpz_cell * new_cell = allocate(capacity);
SASSERT(new_cell->m_capacity == capacity);
unsigned old_sz = a.m_ptr->m_size;
new_cell->m_size = old_sz;
for (unsigned i = 0; i < old_sz; i++)
new_cell->m_digits[i] = a.m_ptr->m_digits[i];
deallocate(a);
a.m_ptr = new_cell;
a.m_owner = mpz_self;
a.m_kind = mpz_large;
}
}
template
void mpz_manager::normalize(mpz & a) {
mpz_cell * c = a.m_ptr;
digit_t * ds = c->m_digits;
unsigned i = c->m_size;
for (; i > 0; --i) {
if (ds[i-1] != 0)
break;
}
if (i == 0) {
// a is zero...
set(a, 0);
return;
}
if (i == 1 && ds[0] <= INT_MAX) {
// a is small
int val = a.m_val < 0 ? -static_cast(ds[0]) : static_cast(ds[0]);
a.m_val = val;
a.m_kind = mpz_small;
return;
}
// adjust size
c->m_size = i;
}
#endif
template
void mpz_manager::machine_div2k(mpz & a, unsigned k) {
if (k == 0 || is_zero(a))
return;
if (is_small(a)) {
if (k < 32) {
int64_t twok = 1ull << ((int64_t)k);
int64_t val = a.m_val;
a.m_val = (int)(val/twok);
}
else {
a.m_val = 0;
}
return;
}
#ifndef _MP_GMP
unsigned digit_shift = k / (8 * sizeof(digit_t));
mpz_cell * c = a.m_ptr;
unsigned sz = c->m_size;
if (digit_shift >= sz) {
set(a, 0);
return;
}
unsigned bit_shift = k % (8 * sizeof(digit_t));
unsigned comp_shift = (8 * sizeof(digit_t)) - bit_shift;
unsigned new_sz = sz - digit_shift;
SASSERT(new_sz >= 1);
digit_t * ds = c->m_digits;
TRACE("mpz_2k", tout << "bit_shift: " << bit_shift << ", comp_shift: " << comp_shift << ", new_sz: " << new_sz << ", sz: " << sz << "\n";);
if (new_sz < sz) {
unsigned i = 0;
unsigned j = digit_shift;
if (bit_shift != 0) {
for (; i < new_sz - 1; i++, j++) {
ds[i] = ds[j];
ds[i] >>= bit_shift;
ds[i] |= (ds[j+1] << comp_shift);
}
ds[i] = ds[j];
ds[i] >>= bit_shift;
}
else {
for (; i < new_sz; i++, j++) {
ds[i] = ds[j];
}
}
}
else {
SASSERT(new_sz == sz);
SASSERT(bit_shift != 0);
unsigned i = 0;
for (; i < new_sz - 1; i++) {
ds[i] >>= bit_shift;
ds[i] |= (ds[i+1] << comp_shift);
}
ds[i] >>= bit_shift;
}
c->m_size = new_sz;
normalize(a);
#else
ensure_mpz_t a1(a);
MPZ_BEGIN_CRITICAL();
mpz_tdiv_q_2exp(m_tmp, a1(), k);
mk_big(a);
mpz_swap(*a.m_ptr, m_tmp);
MPZ_END_CRITICAL();
#endif
}
template
void mpz_manager::mul2k(mpz & a, unsigned k) {
if (k == 0 || is_zero(a))
return;
if (is_small(a) && k < 32) {
set_i64(a, i64(a) * (static_cast(1) << k));
return;
}
#ifndef _MP_GMP
TRACE("mpz_mul2k", tout << "mul2k\na: " << to_string(a) << "\nk: " << k << "\n";);
unsigned word_shift = k / (8 * sizeof(digit_t));
unsigned bit_shift = k % (8 * sizeof(digit_t));
unsigned old_sz = is_small(a) ? 1 : a.m_ptr->m_size;
unsigned new_sz = old_sz + word_shift + 1;
ensure_capacity(a, new_sz);
TRACE("mpz_mul2k", tout << "word_shift: " << word_shift << "\nbit_shift: " << bit_shift << "\nold_sz: " << old_sz << "\nnew_sz: " << new_sz
<< "\na after ensure capacity:\n" << to_string(a) << "\n";
tout << a.m_kind << "\n";);
SASSERT(!is_small(a));
mpz_cell * cell_a = a.m_ptr;
old_sz = cell_a->m_size;
digit_t * ds = cell_a->m_digits;
for (unsigned i = old_sz; i < new_sz; i++)
ds[i] = 0;
cell_a->m_size = new_sz;
if (word_shift > 0) {
unsigned j = old_sz;
unsigned i = old_sz + word_shift;
while (j > 0) {
--j; --i;
ds[i] = ds[j];
}
while (i > 0) {
--i;
ds[i] = 0;
}
}
if (bit_shift > 0) {
DEBUG_CODE({
for (unsigned i = 0; i < word_shift; i++) {
SASSERT(ds[i] == 0);
}
});
unsigned comp_shift = (8 * sizeof(digit_t)) - bit_shift;
digit_t prev = 0;
for (unsigned i = word_shift; i < new_sz; i++) {
digit_t new_prev = (ds[i] >> comp_shift);
ds[i] <<= bit_shift;
ds[i] |= prev;
prev = new_prev;
}
}
normalize(a);
TRACE("mpz_mul2k", tout << "mul2k result:\n" << to_string(a) << "\n";);
#else
ensure_mpz_t a1(a);
mk_big(a);
mpz_mul_2exp(*a.m_ptr, a1(), k);
#endif
}
#ifndef _MP_GMP
static_assert(sizeof(digit_t) == 4 || sizeof(digit_t) == 8, "");
#endif
template
unsigned mpz_manager::power_of_two_multiple(mpz const & a) {
if (is_zero(a))
return 0;
if (is_small(a)) {
unsigned r = 0;
int v = a.m_val;
#define COUNT_DIGIT_RIGHT_ZEROS() \
if (v % (1 << 16) == 0) { \
r += 16; \
v /= (1 << 16); \
} \
if (v % (1 << 8) == 0) { \
r += 8; \
v /= (1 << 8); \
} \
if (v % (1 << 4) == 0) { \
r += 4; \
v /= (1 << 4); \
} \
if (v % (1 << 2) == 0) { \
r += 2; \
v /= (1 << 2); \
} \
if (v % 2 == 0) { \
r++; \
}
COUNT_DIGIT_RIGHT_ZEROS();
return r;
}
#ifndef _MP_GMP
mpz_cell * c = a.m_ptr;
unsigned sz = c->m_size;
unsigned r = 0;
digit_t * source = c->m_digits;
for (unsigned i = 0; i < sz; i++) {
if (source[i] != 0) {
digit_t v = source[i];
if (sizeof(digit_t) == 8) {
// TODO: we can remove this if after we move to MPN
// In MPN the digit_t is always an unsigned integer
if (static_cast(v) % (static_cast(1) << 32) == 0) {
r += 32;
v = static_cast(static_cast(v) / (static_cast(1) << 32));
}
}
COUNT_DIGIT_RIGHT_ZEROS();
return r;
}
r += (8 * sizeof(digit_t));
}
return r;
#else
return mpz_scan1(*a.m_ptr, 0);
#endif
}
template
unsigned mpz_manager::log2(mpz const & a) {
if (is_nonpos(a))
return 0;
if (is_small(a))
return ::log2((unsigned)a.m_val);
#ifndef _MP_GMP
static_assert(sizeof(digit_t) == 8 || sizeof(digit_t) == 4, "");
mpz_cell * c = a.m_ptr;
unsigned sz = c->m_size;
digit_t * ds = c->m_digits;
if (sizeof(digit_t) == 8)
return (sz - 1)*64 + uint64_log2(ds[sz-1]);
else
return (sz - 1)*32 + ::log2(static_cast(ds[sz-1]));
#else
unsigned r = mpz_sizeinbase(*a.m_ptr, 2);
SASSERT(r > 0);
return r - 1;
#endif
}
template
unsigned mpz_manager::mlog2(mpz const & a) {
if (is_nonneg(a))
return 0;
if (is_small(a) && a.m_val == INT_MIN)
return ::log2((unsigned)a.m_val);
if (is_small(a))
return ::log2((unsigned)-a.m_val);
#ifndef _MP_GMP
static_assert(sizeof(digit_t) == 8 || sizeof(digit_t) == 4, "");
mpz_cell * c = a.m_ptr;
unsigned sz = c->m_size;
digit_t * ds = c->m_digits;
if (sizeof(digit_t) == 8)
return (sz - 1)*64 + uint64_log2(ds[sz-1]);
else
return (sz - 1)*32 + ::log2(static_cast(ds[sz-1]));
#else
MPZ_BEGIN_CRITICAL();
mpz_neg(m_tmp, *a.m_ptr);
unsigned r = mpz_sizeinbase(m_tmp, 2);
MPZ_END_CRITICAL();
SASSERT(r > 0);
return r - 1;
#endif
}
template
unsigned mpz_manager::bitsize(mpz const & a) {
if (is_nonneg(a))
return log2(a) + 1;
else
return mlog2(a) + 1;
}
template
unsigned mpz_manager::next_power_of_two(mpz const & a) {
if (is_nonpos(a))
return 0;
if (is_one(a))
return 0;
unsigned shift;
if (is_power_of_two(a, shift))
return shift;
else
return log2(a) + 1;
}
template
bool mpz_manager::is_perfect_square(mpz const & a, mpz & root) {
if (is_neg(a))
return false;
set(root, 0);
if (is_zero(a)) {
return true;
}
if (is_one(a)) {
set(root, 1);
return true;
}
mpz lo, hi, mid, sq_lo, sq_hi, sq_mid;
set(lo, 1);
set(hi, a);
set(sq_lo, 1);
mul(hi, hi, sq_hi);
bool result;
// lo*lo <= *this < hi*hi
while (true) {
SASSERT(lt(lo, hi));
SASSERT(le(sq_lo, a) && lt(a, sq_hi));
if (eq(sq_lo, a)) {
set(root, lo);
result = true;
break;
}
mpz & tmp = mid;
add(lo, mpz(1), tmp);
if (eq(tmp, hi)) {
set(root, hi);
result = false;
break;
}
add(hi, lo, tmp);
div(tmp, mpz(2), mid);
SASSERT(lt(lo, mid) && lt(mid, hi));
mul(mid, mid, sq_mid);
if (gt(sq_mid, a)) {
set(hi, mid);
set(sq_hi, sq_mid);
}
else {
set(lo, mid);
set(sq_lo, sq_mid);
}
}
del(lo);
del(hi);
del(mid);
del(sq_lo);
del(sq_hi);
del(sq_mid);
return result;
}
static unsigned div_l(unsigned k, unsigned n) {
return k/n;
}
static unsigned div_u(unsigned k, unsigned n) {
return k%n == 0 ? k/n : k/n + 1;
}
template
bool mpz_manager::root(mpz & a, unsigned n) {
SASSERT(n % 2 != 0 || is_nonneg(a));
if (is_zero(a)) {
return true; // precise
}
// Initial approximation
//
// We have that:
// a > 0 -> 2^{log2(a)} <= a <= 2^{(log2(a) + 1)}
// a < 0 -> -2^{log2(a) + 1} <= a <= -2^{log2(a)}
//
// Thus
// a > 0 -> 2^{div_l(log2(a), n)} <= a^{1/n} <= 2^{div_u(log2(a) + 1, n)}
// a < 0 -> -2^{div_u(log2(a) + 1, n)} <= a^{1/n} <= -2^{div_l(log2(a), n)}
//
mpz lower;
mpz upper;
mpz mid;
mpz mid_n;
if (is_pos(a)) {
unsigned k = log2(a);
power(mpz(2), div_l(k, n), lower);
power(mpz(2), div_u(k + 1, n), upper);
}
else {
unsigned k = mlog2(a);
power(mpz(2), div_u(k + 1, n), lower);
power(mpz(2), div_l(k, n), upper);
neg(lower);
neg(upper);
}
bool result;
SASSERT(le(lower, upper));
if (eq(lower, upper)) {
swap(a, lower);
result = true;
}
else {
// Refine using bisection. TODO: use Newton's method if this is a bottleneck
while (true) {
add(upper, lower, mid);
machine_div2k(mid, 1);
TRACE("mpz", tout << "upper: "; display(tout, upper); tout << "\nlower: "; display(tout, lower); tout << "\nmid: "; display(tout, mid); tout << "\n";);
power(mid, n, mid_n);
if (eq(mid_n, a)) {
swap(a, mid);
result = true;
break;
}
if (eq(mid, lower) || eq(mid, upper)) {
swap(a, upper);
result = false;
break;
}
if (lt(mid_n, a)) {
// new lower bound
swap(mid, lower);
}
else {
SASSERT(lt(a, mid_n));
// new upper bound
swap(mid, upper);
}
}
}
del(lower);
del(upper);
del(mid);
del(mid_n);
return result;
}
template
bool mpz_manager::decompose(mpz const & a, svector & digits) {
digits.reset();
if (is_small(a)) {
if (a.m_val < 0) {
digits.push_back(-a.m_val);
return true;
}
else {
digits.push_back(a.m_val);
return false;
}
}
else {
#ifndef _MP_GMP
mpz_cell * cell_a = a.m_ptr;
unsigned sz = cell_a->m_size;
for (unsigned i = 0; i < sz; i++) {
digits.push_back(cell_a->m_digits[i]);
}
return a.m_val < 0;
#else
bool r = is_neg(a);
MPZ_BEGIN_CRITICAL();
mpz_set(m_tmp, *a.m_ptr);
mpz_abs(m_tmp, m_tmp);
while (mpz_sgn(m_tmp) != 0) {
mpz_tdiv_r_2exp(m_tmp2, m_tmp, 32);
unsigned v = mpz_get_ui(m_tmp2);
digits.push_back(v);
mpz_tdiv_q_2exp(m_tmp, m_tmp, 32);
}
MPZ_END_CRITICAL();
return r;
#endif
}
}
template
bool mpz_manager::get_bit(mpz const & a, unsigned index) {
if (is_small(a)) {
SASSERT(a.m_val >= 0);
if (index >= 8*sizeof(digit_t))
return false;
return 0 != (a.m_val & (1ull << (digit_t)index));
}
unsigned i = index / (sizeof(digit_t)*8);
unsigned o = index % (sizeof(digit_t)*8);
#ifndef _MP_GMP
mpz_cell * cell_a = a.m_ptr;
unsigned sz = cell_a->m_size;
if (sz*sizeof(digit_t)*8 <= index)
return false;
return 0 != (cell_a->m_digits[i] & (1ull << (digit_t)o));
#else
SASSERT(!is_neg(a));
svector digits;
decompose(a, digits);
if (digits.size()*sizeof(digit_t)*8 <= index)
return false;
return 0 != (digits[i] & (1ull << (digit_t)o));
#endif
}
template
bool mpz_manager::divides(mpz const & a, mpz const & b) {
_scoped_numeral > tmp(*this);
bool r;
if (is_zero(a)) {
// I assume 0 | 0.
// Remark a|b is a shorthand for (exists x. a x = b)
// If b is zero, any x will do. If b != 0, then a does not divide b
r = is_zero(b);
}
else {
rem(b, a, tmp);
r = is_zero(tmp);
}
STRACE("divides", tout << "[mpz] Divisible["; display(tout, b); tout << ", "; display(tout, a); tout << "] == " << (r?"True":"False") << "\n";);
TRACE("divides_bug", tout << "tmp: "; display(tout, tmp); tout << "\n";);
return r;
}
#ifndef SINGLE_THREAD
template class mpz_manager;
#endif
template class mpz_manager;