// Copyright 2018 vDaniel Parker // Distributed under the Boost license, Version 1.0. // (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) // See https://github.com/danielaparker/jsoncons for latest version #ifndef JSONCONS_BIGINT_HPP #define JSONCONS_BIGINT_HPP #include #include // std::vector #include #include #include // assert #include // std::numeric_limits #include // std::max, std::min, std::reverse #include // std::string #include // std::memcpy #include // std::fmod #include // std::allocator #include // std::initializer_list #include // std::enable_if #include namespace jsoncons { /* This implementation is based on Chapter 2 and Appendix A of Ammeraal, L. (1996) Algorithms and Data Structures in C++, Chichester: John Wiley. */ namespace detail { template class basic_bigint_base { public: using allocator_type = Allocator; using basic_type_allocator_type = typename std::allocator_traits:: template rebind_alloc; private: basic_type_allocator_type alloc_; public: using allocator_traits_type = std::allocator_traits; using stored_allocator_type = allocator_type; using pointer = typename allocator_traits_type::pointer; using value_type = typename allocator_traits_type::value_type; using size_type = typename allocator_traits_type::size_type; using pointer_traits = std::pointer_traits; basic_bigint_base() : alloc_() { } explicit basic_bigint_base(const allocator_type& alloc) : alloc_(basic_type_allocator_type(alloc)) { } basic_type_allocator_type get_allocator() const { return alloc_; } }; } // namespace detail template > class basic_bigint : protected detail::basic_bigint_base { using base_t = detail::basic_bigint_base; public: using size_type = typename base_t::size_type; using value_type = typename base_t::value_type; using base_t::get_allocator; using bigint_type = basic_bigint; static constexpr uint64_t max_basic_type = (std::numeric_limits::max)(); static constexpr uint64_t basic_type_bits = sizeof(uint64_t) * 8; // Number of bits static constexpr uint64_t basic_type_halfBits = basic_type_bits/2; static constexpr uint16_t word_length = 4; // Use multiples of word_length words static constexpr uint64_t r_mask = (uint64_t(1) << basic_type_halfBits) - 1; static constexpr uint64_t l_mask = max_basic_type - r_mask; static constexpr uint64_t l_bit = max_basic_type - (max_basic_type >> 1); private: struct common_storage { uint8_t is_dynamic_:1; uint8_t is_negative_:1; size_type length_; }; struct short_storage { uint8_t is_dynamic_:1; uint8_t is_negative_:1; size_type length_; uint64_t values_[2]; short_storage() : is_dynamic_(false), is_negative_(false), length_(0), values_{0,0} { } template short_storage(T n, typename std::enable_if::value && sizeof(T) <= sizeof(int64_t) && std::is_signed::value>::type* = 0) : is_dynamic_(false), is_negative_(n < 0), length_(n == 0 ? 0 : 1) { values_[0] = n < 0 ? (uint64_t(0)-static_cast(n)) : static_cast(n); values_[1] = 0; } template short_storage(T n, typename std::enable_if::value && sizeof(T) <= sizeof(int64_t) && !std::is_signed::value>::type* = 0) : is_dynamic_(false), is_negative_(false), length_(n == 0 ? 0 : 1) { values_[0] = n; values_[1] = 0; } template short_storage(T n, typename std::enable_if::value && sizeof(int64_t) < sizeof(T) && std::is_signed::value>::type* = 0) : is_dynamic_(false), is_negative_(n < 0), length_(n == 0 ? 0 : 2) { using unsigned_type = typename std::make_unsigned::type; auto u = n < 0 ? (unsigned_type(0)-static_cast(n)) : static_cast(n); values_[0] = uint64_t(u & max_basic_type);; u >>= basic_type_bits; values_[1] = uint64_t(u & max_basic_type);; } template short_storage(T n, typename std::enable_if::value && sizeof(int64_t) < sizeof(T) && !std::is_signed::value>::type* = 0) : is_dynamic_(false), is_negative_(false), length_(n == 0 ? 0 : 2) { values_[0] = uint64_t(n & max_basic_type);; n >>= basic_type_bits; values_[1] = uint64_t(n & max_basic_type);; } short_storage(const short_storage& stor) : is_dynamic_(false), is_negative_(stor.is_negative_), length_(stor.length_) { values_[0] = stor.values_[0]; values_[1] = stor.values_[1]; } short_storage& operator=(const short_storage& stor) = delete; }; struct dynamic_storage { using real_allocator_type = typename std::allocator_traits:: template rebind_alloc; using pointer = typename std::allocator_traits::pointer; uint8_t is_dynamic_:1; uint8_t is_negative_:1; size_type length_; size_type capacity_; pointer data_; dynamic_storage() : is_dynamic_(true), is_negative_(false), length_(0), capacity_(0), data_(nullptr) { } dynamic_storage(const dynamic_storage& stor, const real_allocator_type& alloc) : is_dynamic_(true), is_negative_(stor.is_negative_), length_(stor.length_), capacity_(0), data_(nullptr) { create(stor.length_, alloc); std::memcpy(data_, stor.data_, stor.length_*sizeof(uint64_t)); } dynamic_storage(dynamic_storage&& stor) noexcept : is_dynamic_(true), is_negative_(stor.is_negative_), length_(stor.length_), capacity_(stor.capacity_), data_(stor.data_) { stor.length_ = 0; stor.capacity_ = 0; stor.data_ = nullptr; } void create(size_type length, real_allocator_type alloc) { capacity_ = round_up(length); data_ = std::allocator_traits::allocate(alloc, capacity_); JSONCONS_TRY { std::allocator_traits::construct(alloc, type_traits::to_plain_pointer(data_)); } JSONCONS_CATCH(...) { std::allocator_traits::deallocate(alloc, data_, capacity_); JSONCONS_RETHROW; } } void destroy(const real_allocator_type& a) noexcept { if (data_ != nullptr) { real_allocator_type alloc(a); std::allocator_traits::destroy(alloc, type_traits::to_plain_pointer(data_)); std::allocator_traits::deallocate(alloc, data_,capacity_); } } void reserve(size_type n, const real_allocator_type& a) { real_allocator_type alloc(a); size_type capacity_new = round_up(n); uint64_t* data_old = data_; data_ = std::allocator_traits::allocate(alloc, capacity_new); if (length_ > 0) { std::memcpy( data_, data_old, length_*sizeof(uint64_t)); } if (capacity_ > 0 && data_ != nullptr) { std::allocator_traits::deallocate(alloc, data_old, capacity_); } capacity_ = capacity_new; } // Find suitable new block size constexpr size_type round_up(size_type i) const { return (i/word_length + 1) * word_length; } }; union { common_storage common_stor_; short_storage short_stor_; dynamic_storage dynamic_stor_; }; public: basic_bigint() { ::new (&short_stor_) short_storage(); } explicit basic_bigint(const Allocator& alloc) : base_t(alloc) { ::new (&short_stor_) short_storage(); } basic_bigint(const basic_bigint& n) : base_t(n.get_allocator()) { if (!n.is_dynamic()) { ::new (&short_stor_) short_storage(n.short_stor_); } else { ::new (&dynamic_stor_) dynamic_storage(n.dynamic_stor_, get_allocator()); } } basic_bigint(basic_bigint&& other) noexcept : base_t(other.get_allocator()) { if (!other.is_dynamic()) { ::new (&short_stor_) short_storage(other.short_stor_); } else { ::new (&dynamic_stor_) dynamic_storage(std::move(other.dynamic_stor_)); } } template basic_bigint(Integer n, typename std::enable_if::value>::type* = 0) { ::new (&short_stor_) short_storage(n); } ~basic_bigint() noexcept { destroy(); } constexpr bool is_dynamic() const { return common_stor_.is_dynamic_; } constexpr size_type length() const { return common_stor_.length_; } constexpr size_type capacity() const { return is_dynamic() ? dynamic_stor_.capacity_ : 2; } bool is_negative() const { return common_stor_.is_negative_; } void is_negative(bool value) { common_stor_.is_negative_ = value; } constexpr const uint64_t* data() const { return is_dynamic() ? dynamic_stor_.data_ : short_stor_.values_; } uint64_t* data() { return is_dynamic() ? dynamic_stor_.data_ : short_stor_.values_; } template static basic_bigint from_string(const std::basic_string& s) { return from_string(s.data(), s.length()); } template static basic_bigint from_string(const CharT* s) { return from_string(s, std::char_traits::length(s)); } template static basic_bigint from_string(const CharT* data, size_type length) { bool neg; if (*data == '-') { neg = true; data++; --length; } else { neg = false; } basic_bigint v = 0; for (size_type i = 0; i < length; i++) { CharT c = data[i]; switch (c) { case '0':case '1':case '2':case '3':case '4':case '5':case '6':case '7':case '8': case '9': v = (v * 10) + (uint64_t)(c - '0'); break; default: JSONCONS_THROW(std::runtime_error(std::string("Invalid digit ") + "\'" + (char)c + "\'")); } } if (neg) { v.common_stor_.is_negative_ = true; } return v; } template static basic_bigint from_string_radix(const CharT* data, size_type length, uint8_t radix) { if (!(radix >= 2 && radix <= 16)) { JSONCONS_THROW(std::runtime_error("Unsupported radix")); } bool neg; if (*data == '-') { neg = true; data++; --length; } else { neg = false; } basic_bigint v = 0; for (size_type i = 0; i < length; i++) { CharT c = data[i]; uint64_t d; switch (c) { case '0':case '1':case '2':case '3':case '4':case '5':case '6':case '7':case '8': case '9': d = (uint64_t)(c - '0'); break; case 'a':case 'b':case 'c':case 'd':case 'e':case 'f': d = (uint64_t)(c - ('a' - 10)); break; case 'A':case 'B':case 'C':case 'D':case 'E':case 'F': d = (uint64_t)(c - ('A' - 10)); break; default: JSONCONS_THROW(std::runtime_error(std::string("Invalid digit in radix ") + std::to_string(radix) + ": \'" + (char)c + "\'")); } if (d >= radix) { JSONCONS_THROW(std::runtime_error(std::string("Invalid digit in radix ") + std::to_string(radix) + ": \'" + (char)c + "\'")); } v = (v * radix) + d; } if ( neg ) { v.common_stor_.is_negative_ = true; } return v; } static basic_bigint from_bytes_be(int signum, const uint8_t* str, std::size_t n) { static const double radix_log2 = std::log2(next_power_of_two(256)); // Estimate how big the result will be, so we can pre-allocate it. double bits = radix_log2 * n; double big_digits = std::ceil(bits / 64.0); //std::cout << "ESTIMATED: " << big_digits << "\n"; bigint_type v = 0; v.reserve(static_cast(big_digits)); if (n > 0) { for (std::size_t i = 0; i < n; i++) { v = (v * 256) + (uint64_t)(str[i]); } } //std::cout << "ACTUAL: " << v.length() << "\n"; if (signum < 0) { v.common_stor_.is_negative_ = true; } return v; } uint64_t* begin() { return is_dynamic() ? dynamic_stor_.data_ : short_stor_.values_; } const uint64_t* begin() const { return is_dynamic() ? dynamic_stor_.data_ : short_stor_.values_; } uint64_t* end() { return begin() + length(); } const uint64_t* end() const { return begin() + length(); } void resize(size_type n) { size_type len_old = common_stor_.length_; reserve(n); common_stor_.length_ = n; if ( common_stor_.length_ > len_old ) { memset( data()+len_old, 0, (common_stor_.length_ - len_old)*sizeof(uint64_t) ); } } void reserve(size_type n) { if (capacity() < n) { if (!is_dynamic()) { size_type size = short_stor_.length_; size_type is_neg = short_stor_.is_negative_; uint64_t values[2] = {short_stor_.values_[0], short_stor_.values_[1]}; ::new (&dynamic_stor_) dynamic_storage(); dynamic_stor_.reserve(n, get_allocator()); dynamic_stor_.length_ = size; dynamic_stor_.is_negative_ = is_neg; dynamic_stor_.data_[0] = values[0]; dynamic_stor_.data_[1] = values[1]; } else { dynamic_stor_.reserve(n, get_allocator()); } } } // operators bool operator!() const { return length() == 0 ? true : false; } basic_bigint operator-() const { basic_bigint v(*this); v.common_stor_.is_negative_ = !v.is_negative(); return v; } basic_bigint& operator=( const basic_bigint& y ) { if ( this != &y ) { resize( y.length() ); common_stor_.is_negative_ = y.is_negative(); if ( y.length() > 0 ) { std::memcpy( data(), y.data(), y.length()*sizeof(uint64_t) ); } } return *this; } basic_bigint& operator+=( const basic_bigint& y ) { if ( is_negative() != y.is_negative() ) return *this -= -y; uint64_t d; uint64_t carry = 0; resize( (std::max)(y.length(), length()) + 1 ); for (size_type i = 0; i < length(); i++ ) { if ( i >= y.length() && carry == 0 ) break; d = data()[i] + carry; carry = d < carry; if ( i < y.length() ) { data()[i] = d + y.data()[i]; if ( data()[i] < d ) carry = 1; } else data()[i] = d; } reduce(); return *this; } basic_bigint& operator-=( const basic_bigint& y ) { if ( is_negative() != y.is_negative() ) return *this += -y; if ( (!is_negative() && y > *this) || (is_negative() && y < *this) ) return *this = -(y - *this); uint64_t borrow = 0; uint64_t d; for (size_type i = 0; i < length(); i++ ) { if ( i >= y.length() && borrow == 0 ) break; d = data()[i] - borrow; borrow = d > data()[i]; if ( i < y.length()) { data()[i] = d - y.data()[i]; if ( data()[i] > d ) borrow = 1; } else data()[i] = d; } reduce(); return *this; } basic_bigint& operator*=( int64_t y ) { *this *= uint64_t(y < 0 ? -y : y); if ( y < 0 ) common_stor_.is_negative_ = !is_negative(); return *this; } basic_bigint& operator*=( uint64_t y ) { size_type len0 = length(); uint64_t hi; uint64_t lo; uint64_t dig = data()[0]; uint64_t carry = 0; resize( length() + 1 ); size_type i = 0; for (i = 0; i < len0; i++ ) { DDproduct( dig, y, hi, lo ); data()[i] = lo + carry; dig = data()[i+1]; carry = hi + (data()[i] < lo); } data()[i] = carry; reduce(); return *this; } basic_bigint& operator*=( basic_bigint y ) { if ( length() == 0 || y.length() == 0 ) return *this = 0; bool difSigns = is_negative() != y.is_negative(); if ( length() + y.length() == 2 ) // length() = y.length() = 1 { uint64_t a = data()[0], b = y.data()[0]; data()[0] = a * b; if ( data()[0] / a != b ) { resize( 2 ); DDproduct( a, b, data()[1], data()[0] ); } common_stor_.is_negative_ = difSigns; return *this; } if ( length() == 1 ) // && y.length() > 1 { uint64_t digit = data()[0]; *this = y; *this *= digit; } else { if ( y.length() == 1 ) *this *= y.data()[0]; else { size_type lenProd = length() + y.length(), jA, jB; uint64_t sumHi = 0, sumLo, hi, lo, sumLo_old, sumHi_old, carry=0; basic_bigint x = *this; resize( lenProd ); // Give *this length lenProd for (size_type i = 0; i < lenProd; i++ ) { sumLo = sumHi; sumHi = carry; carry = 0; for ( jA=0; jA < x.length(); jA++ ) { jB = i - jA; if ( jB >= 0 && jB < y.length() ) { DDproduct( x.data()[jA], y.data()[jB], hi, lo ); sumLo_old = sumLo; sumHi_old = sumHi; sumLo += lo; if ( sumLo < sumLo_old ) sumHi++; sumHi += hi; carry += (sumHi < sumHi_old); } } data()[i] = sumLo; } } } reduce(); common_stor_.is_negative_ = difSigns; return *this; } basic_bigint& operator/=( const basic_bigint& divisor ) { basic_bigint r; divide( divisor, *this, r, false ); return *this; } basic_bigint& operator%=( const basic_bigint& divisor ) { basic_bigint q; divide( divisor, q, *this, true ); return *this; } basic_bigint& operator<<=( uint64_t k ) { size_type q = (size_type)(k / basic_type_bits); if ( q ) // Increase common_stor_.length_ by q: { resize(length() + q); for (size_type i = length(); i-- > 0; ) data()[i] = ( i < q ? 0 : data()[i - q]); k %= basic_type_bits; } if ( k ) // 0 < k < basic_type_bits: { uint64_t k1 = basic_type_bits - k; uint64_t mask = (uint64_t(1) << k) - uint64_t(1); resize( length() + 1 ); for (size_type i = length(); i-- > 0; ) { data()[i] <<= k; if ( i > 0 ) data()[i] |= (data()[i-1] >> k1) & mask; } } reduce(); return *this; } basic_bigint& operator>>=(uint64_t k) { size_type q = (size_type)(k / basic_type_bits); if ( q >= length() ) { resize( 0 ); return *this; } if (q > 0) { memmove( data(), data()+q, (size_type)((length() - q)*sizeof(uint64_t)) ); resize( length() - q ); k %= basic_type_bits; if ( k == 0 ) { reduce(); return *this; } } size_type n = (size_type)(length() - 1); int64_t k1 = basic_type_bits - k; uint64_t mask = (uint64_t(1) << k) - 1; for (size_type i = 0; i <= n; i++) { data()[i] >>= k; if ( i < n ) data()[i] |= ((data()[i+1] & mask) << k1); } reduce(); return *this; } basic_bigint& operator++() { *this += 1; return *this; } basic_bigint operator++(int) { basic_bigint old = *this; ++(*this); return old; } basic_bigint& operator--() { *this -= 1; return *this; } basic_bigint operator--(int) { basic_bigint old = *this; --(*this); return old; } basic_bigint& operator|=( const basic_bigint& a ) { if ( length() < a.length() ) { resize( a.length() ); } const uint64_t* qBegin = a.begin(); const uint64_t* q = a.end() - 1; uint64_t* p = begin() + a.length() - 1; while ( q >= qBegin ) { *p-- |= *q--; } reduce(); return *this; } basic_bigint& operator^=( const basic_bigint& a ) { if ( length() < a.length() ) { resize( a.length() ); } const uint64_t* qBegin = a.begin(); const uint64_t* q = a.end() - 1; uint64_t* p = begin() + a.length() - 1; while ( q >= qBegin ) { *p-- ^= *q--; } reduce(); return *this; } basic_bigint& operator&=( const basic_bigint& a ) { size_type old_length = length(); resize( (std::min)( length(), a.length() ) ); const uint64_t* pBegin = begin(); uint64_t* p = end() - 1; const uint64_t* q = a.begin() + length() - 1; while ( p >= pBegin ) { *p-- &= *q--; } if ( old_length > length() ) { memset( data() + length(), 0, old_length - length() ); } reduce(); return *this; } explicit operator bool() const { return length() != 0 ? true : false; } explicit operator int64_t() const { int64_t x = 0; if ( length() > 0 ) { x = data() [0]; } return is_negative() ? -x : x; } explicit operator uint64_t() const { uint64_t u = 0; if ( length() > 0 ) { u = data() [0]; } return u; } explicit operator double() const { double x = 0.0; double factor = 1.0; double values = (double)max_basic_type + 1.0; const uint64_t* p = begin(); const uint64_t* pEnd = end(); while ( p < pEnd ) { x += *p*factor; factor *= values; ++p; } return is_negative() ? -x : x; } explicit operator long double() const { long double x = 0.0; long double factor = 1.0; long double values = (long double)max_basic_type + 1.0; const uint64_t* p = begin(); const uint64_t* pEnd = end(); while ( p < pEnd ) { x += *p*factor; factor *= values; ++p; } return is_negative() ? -x : x; } template void write_bytes_be(int& signum, std::vector& data) const { basic_bigint n(*this); signum = (n < 0) ? -1 : (n > 0 ? 1 : 0); basic_bigint divisor(256); while (n >= 256) { basic_bigint q; basic_bigint r; n.divide(divisor, q, r, true); n = q; data.push_back((uint8_t)(uint64_t)r); } if (n >= 0) { data.push_back((uint8_t)(uint64_t)n); } std::reverse(data.begin(),data.end()); } std::string to_string() const { std::string s; write_string(s); return s; } template void write_string(std::basic_string& data) const { basic_bigint v(*this); size_t len = (v.length() * basic_type_bits / 3) + 2; data.reserve(len); static uint64_t p10 = 1; static uint64_t ip10 = 0; if ( v.length() == 0 ) { data.push_back('0'); } else { uint64_t r; if ( p10 == 1 ) { while ( p10 <= (std::numeric_limits::max)()/10 ) { p10 *= 10; ip10++; } } // p10 is max unsigned power of 10 basic_bigint R; basic_bigint LP10 = p10; // LP10 = p10 = ::pow(10, ip10) do { v.divide( LP10, v, R, true ); r = (R.length() ? R.data()[0] : 0); for ( size_type j=0; j < ip10; j++ ) { data.push_back(char(r % 10 + '0')); r /= 10; if ( r + v.length() == 0 ) break; } } while ( v.length() ); if (is_negative()) { data.push_back('-'); } std::reverse(data.begin(),data.end()); } } std::string to_string_hex() const { std::string s; write_string_hex(s); return s; } template void write_string_hex(std::basic_string& data) const { basic_bigint v(*this); std::size_t len = (v.length() * basic_bigint::basic_type_bits / 3) + 2; data.reserve(len); // 1/3 > ln(2)/ln(10) static uint64_t p10 = 1; static uint64_t ip10 = 0; if ( v.length() == 0 ) { data.push_back('0'); } else { uint64_t r; if ( p10 == 1 ) { while ( p10 <= (std::numeric_limits::max)()/16 ) { p10 *= 16; ip10++; } } // p10 is max unsigned power of 16 basic_bigint R; basic_bigint LP10 = p10; // LP10 = p10 = ::pow(16, ip10) do { v.divide( LP10, v, R, true ); r = (R.length() ? R.data()[0] : 0); for ( size_type j=0; j < ip10; j++ ) { uint8_t c = r % 16; data.push_back((c < 10) ? ('0' + c) : ('A' - 10 + c)); r /= 16; if ( r + v.length() == 0 ) break; } } while (v.length()); if (is_negative()) { data.push_back('-'); } std::reverse(data.begin(),data.end()); } } // Global Operators friend bool operator==( const basic_bigint& x, const basic_bigint& y ) noexcept { return x.compare(y) == 0 ? true : false; } friend bool operator==( const basic_bigint& x, int y ) noexcept { return x.compare(y) == 0 ? true : false; } friend bool operator!=( const basic_bigint& x, const basic_bigint& y ) noexcept { return x.compare(y) != 0 ? true : false; } friend bool operator!=( const basic_bigint& x, int y ) noexcept { return x.compare(basic_bigint(y)) != 0 ? true : false; } friend bool operator<( const basic_bigint& x, const basic_bigint& y ) noexcept { return x.compare(y) < 0 ? true : false; } friend bool operator<( const basic_bigint& x, int64_t y ) noexcept { return x.compare(y) < 0 ? true : false; } friend bool operator>( const basic_bigint& x, const basic_bigint& y ) noexcept { return x.compare(y) > 0 ? true : false; } friend bool operator>( const basic_bigint& x, int y ) noexcept { return x.compare(basic_bigint(y)) > 0 ? true : false; } friend bool operator<=( const basic_bigint& x, const basic_bigint& y ) noexcept { return x.compare(y) <= 0 ? true : false; } friend bool operator<=( const basic_bigint& x, int y ) noexcept { return x.compare(y) <= 0 ? true : false; } friend bool operator>=( const basic_bigint& x, const basic_bigint& y ) noexcept { return x.compare(y) >= 0 ? true : false; } friend bool operator>=( const basic_bigint& x, int y ) noexcept { return x.compare(y) >= 0 ? true : false; } friend basic_bigint operator+( basic_bigint x, const basic_bigint& y ) { return x += y; } friend basic_bigint operator+( basic_bigint x, int64_t y ) { return x += y; } friend basic_bigint operator-( basic_bigint x, const basic_bigint& y ) { return x -= y; } friend basic_bigint operator-( basic_bigint x, int64_t y ) { return x -= y; } friend basic_bigint operator*( int64_t x, const basic_bigint& y ) { return basic_bigint(y) *= x; } friend basic_bigint operator*( basic_bigint x, const basic_bigint& y ) { return x *= y; } friend basic_bigint operator*( basic_bigint x, int64_t y ) { return x *= y; } friend basic_bigint operator/( basic_bigint x, const basic_bigint& y ) { return x /= y; } friend basic_bigint operator/( basic_bigint x, int y ) { return x /= y; } friend basic_bigint operator%( basic_bigint x, const basic_bigint& y ) { return x %= y; } friend basic_bigint operator<<( basic_bigint u, unsigned k ) { return u <<= k; } friend basic_bigint operator<<( basic_bigint u, int k ) { return u <<= k; } friend basic_bigint operator>>( basic_bigint u, unsigned k ) { return u >>= k; } friend basic_bigint operator>>( basic_bigint u, int k ) { return u >>= k; } friend basic_bigint operator|( basic_bigint x, const basic_bigint& y ) { return x |= y; } friend basic_bigint operator|( basic_bigint x, int y ) { return x |= y; } friend basic_bigint operator|( basic_bigint x, unsigned y ) { return x |= y; } friend basic_bigint operator^( basic_bigint x, const basic_bigint& y ) { return x ^= y; } friend basic_bigint operator^( basic_bigint x, int y ) { return x ^= y; } friend basic_bigint operator^( basic_bigint x, unsigned y ) { return x ^= y; } friend basic_bigint operator&( basic_bigint x, const basic_bigint& y ) { return x &= y; } friend basic_bigint operator&( basic_bigint x, int y ) { return x &= y; } friend basic_bigint operator&( basic_bigint x, unsigned y ) { return x &= y; } friend basic_bigint abs( const basic_bigint& a ) { if ( a.is_negative() ) { return -a; } return a; } friend basic_bigint power( basic_bigint x, unsigned n ) { basic_bigint y = 1; while ( n ) { if ( n & 1 ) { y *= x; } x *= x; n >>= 1; } return y; } friend basic_bigint sqrt( const basic_bigint& a ) { basic_bigint x = a; basic_bigint b = a; basic_bigint q; b <<= 1; while ( b >>= 2, b > 0 ) { x >>= 1; } while ( x > (q = a/x) + 1 || x < q - 1 ) { x += q; x >>= 1; } return x < q ? x : q; } template friend std::basic_ostream& operator<<(std::basic_ostream& os, const basic_bigint& v) { std::basic_string s; v.write_string(s); os << s; return os; } int compare( const basic_bigint& y ) const noexcept { if ( is_negative() != y.is_negative() ) return y.is_negative() - is_negative(); int code = 0; if ( length() == 0 && y.length() == 0 ) code = 0; else if ( length() < y.length() ) code = -1; else if ( length() > y.length() ) code = +1; else { for (size_type i = length(); i-- > 0; ) { if (data()[i] > y.data()[i]) { code = 1; break; } else if (data()[i] < y.data()[i]) { code = -1; break; } } } return is_negative() ? -code : code; } void divide( basic_bigint denom, basic_bigint& quot, basic_bigint& rem, bool remDesired ) const { if ( denom.length() == 0 ) { JSONCONS_THROW(std::runtime_error( "Zero divide." )); } bool quot_neg = is_negative() ^ denom.is_negative(); bool rem_neg = is_negative(); int x = 0; basic_bigint num = *this; num.common_stor_.is_negative_ = denom.common_stor_.is_negative_ = false; if ( num < denom ) { quot = uint64_t(0); rem = num; rem.common_stor_.is_negative_ = rem_neg; return; } if ( denom.length() == 1 && num.length() == 1 ) { quot = uint64_t( num.data()[0]/denom.data()[0] ); rem = uint64_t( num.data()[0]%denom.data()[0] ); quot.common_stor_.is_negative_ = quot_neg; rem.common_stor_.is_negative_ = rem_neg; return; } else if (denom.length() == 1 && (denom.data()[0] & l_mask) == 0 ) { // Denominator fits into a half word uint64_t divisor = denom.data()[0], dHi = 0, q1, r, q2, dividend; quot.resize(length()); for (size_type i=length(); i-- > 0; ) { dividend = (dHi << basic_type_halfBits) | (data()[i] >> basic_type_halfBits); q1 = dividend/divisor; r = dividend % divisor; dividend = (r << basic_type_halfBits) | (data()[i] & r_mask); q2 = dividend/divisor; dHi = dividend % divisor; quot.data()[i] = (q1 << basic_type_halfBits) | q2; } quot.reduce(); rem = dHi; quot.common_stor_.is_negative_ = quot_neg; rem.common_stor_.is_negative_ = rem_neg; return; } basic_bigint num0 = num, denom0 = denom; int second_done = normalize(denom, num, x); size_type l = denom.length() - 1; size_type n = num.length() - 1; quot.resize(n - l); for (size_type i=quot.length(); i-- > 0; ) quot.data()[i] = 0; rem = num; if ( rem.data()[n] >= denom.data()[l] ) { rem.resize(rem.length() + 1); n++; quot.resize(quot.length() + 1); } uint64_t d = denom.data()[l]; for ( size_type k = n; k > l; k-- ) { uint64_t q = DDquotient(rem.data()[k], rem.data()[k-1], d); subtractmul( rem.data() + k - l - 1, denom.data(), l + 1, q ); quot.data()[k - l - 1] = q; } quot.reduce(); quot.common_stor_.is_negative_ = quot_neg; if ( remDesired ) { unnormalize(rem, x, second_done); rem.common_stor_.is_negative_ = rem_neg; } } private: void destroy() noexcept { if (is_dynamic()) { dynamic_stor_.destroy(get_allocator()); } } void DDproduct( uint64_t A, uint64_t B, uint64_t& hi, uint64_t& lo ) const // Multiplying two digits: (hi, lo) = A * B { uint64_t hiA = A >> basic_type_halfBits, loA = A & r_mask, hiB = B >> basic_type_halfBits, loB = B & r_mask, mid1, mid2, old; lo = loA * loB; hi = hiA * hiB; mid1 = loA * hiB; mid2 = hiA * loB; old = lo; lo += mid1 << basic_type_halfBits; hi += (lo < old) + (mid1 >> basic_type_halfBits); old = lo; lo += mid2 << basic_type_halfBits; hi += (lo < old) + (mid2 >> basic_type_halfBits); } uint64_t DDquotient( uint64_t A, uint64_t B, uint64_t d ) const // Divide double word (A, B) by d. Quotient = (qHi, qLo) { uint64_t left, middle, right, qHi, qLo, x, dLo1, dHi = d >> basic_type_halfBits, dLo = d & r_mask; qHi = A/(dHi + 1); // This initial guess of qHi may be too small. middle = qHi * dLo; left = qHi * dHi; x = B - (middle << basic_type_halfBits); A -= (middle >> basic_type_halfBits) + left + (x > B); B = x; dLo1 = dLo << basic_type_halfBits; // Increase qHi if necessary: while ( A > dHi || (A == dHi && B >= dLo1) ) { x = B - dLo1; A -= dHi + (x > B); B = x; qHi++; } qLo = ((A << basic_type_halfBits) | (B >> basic_type_halfBits))/(dHi + 1); // This initial guess of qLo may be too small. right = qLo * dLo; middle = qLo * dHi; x = B - right; A -= (x > B); B = x; x = B - (middle << basic_type_halfBits); A -= (middle >> basic_type_halfBits) + (x > B); B = x; // Increase qLo if necessary: while ( A || B >= d ) { x = B - d; A -= (x > B); B = x; qLo++; } return (qHi << basic_type_halfBits) + qLo; } void subtractmul( uint64_t* a, uint64_t* b, size_type n, uint64_t& q ) const // a -= q * b: b in n positions; correct q if necessary { uint64_t hi, lo, d, carry = 0; size_type i; for ( i = 0; i < n; i++ ) { DDproduct( b[i], q, hi, lo ); d = a[i]; a[i] -= lo; if ( a[i] > d ) carry++; d = a[i + 1]; a[i + 1] -= hi + carry; carry = a[i + 1] > d; } if ( carry ) // q was too large { q--; carry = 0; for ( i = 0; i < n; i++ ) { d = a[i] + carry; carry = d < carry; a[i] = d + b[i]; if ( a[i] < d ) carry = 1; } a[n] = 0; } } int normalize( basic_bigint& denom, basic_bigint& num, int& x ) const { size_type r = denom.length() - 1; uint64_t y = denom.data()[r]; x = 0; while ( (y & l_bit) == 0 ) { y <<= 1; x++; } denom <<= x; num <<= x; if ( r > 0 && denom.data()[r] < denom.data()[r-1] ) { denom *= max_basic_type; num *= max_basic_type; return 1; } return 0; } void unnormalize( basic_bigint& rem, int x, int secondDone ) const { if ( secondDone ) { rem /= max_basic_type; } if ( x > 0 ) { rem >>= x; } else { rem.reduce(); } } size_type round_up(size_type i) const // Find suitable new block size { return (i/word_length + 1) * word_length; } void reduce() { uint64_t* p = end() - 1; uint64_t* pBegin = begin(); while ( p >= pBegin ) { if ( *p ) { break; } --common_stor_.length_; --p; } if ( length() == 0 ) { common_stor_.is_negative_ = false; } } static uint64_t next_power_of_two(uint64_t n) { n = n - 1; n |= n >> 1; n |= n >> 2; n |= n >> 4; n |= n >> 8; n |= n >> 16; n |= n >> 32; return n + 1; } }; using bigint = basic_bigint>; #if !defined(JSONCONS_NO_DEPRECATED) JSONCONS_DEPRECATED_MSG("Instead, use bigint") typedef bigint bignum; #endif } #endif