boost/multiprecision/gmp.hpp
/////////////////////////////////////////////////////////////////////////////// // Copyright 2011 John Maddock. // Copyright 2021 Matt Borland. Distributed under the Boost // Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) #ifndef BOOST_MP_GMP_HPP #define BOOST_MP_GMP_HPP #include <boost/multiprecision/detail/standalone_config.hpp> #include <boost/multiprecision/number.hpp> #include <boost/multiprecision/debug_adaptor.hpp> #include <boost/multiprecision/detail/integer_ops.hpp> #include <boost/multiprecision/detail/float128_functions.hpp> #include <boost/multiprecision/detail/digits.hpp> #include <boost/multiprecision/detail/atomic.hpp> #include <boost/multiprecision/detail/hash.hpp> #include <boost/multiprecision/detail/no_exceptions_support.hpp> #include <boost/multiprecision/detail/assert.hpp> #include <boost/multiprecision/detail/fpclassify.hpp> #include <boost/multiprecision/detail/string_helpers.hpp> #include <algorithm> #include <cctype> #include <cfloat> #include <climits> #include <clocale> #include <cmath> #include <cstdint> #include <cstdlib> #include <cstring> #include <iomanip> #include <iostream> #include <limits> #include <memory> #include <type_traits> #include <utility> // // Some includes we need from Boost.Math, since we rely on that library to provide these functions: // #ifdef BOOST_MP_MATH_AVAILABLE #include <boost/math/special_functions/asinh.hpp> #include <boost/math/special_functions/acosh.hpp> #include <boost/math/special_functions/atanh.hpp> #include <boost/math/special_functions/cbrt.hpp> #include <boost/math/special_functions/expm1.hpp> #include <boost/math/special_functions/gamma.hpp> #endif #ifdef BOOST_MSVC #pragma warning(push) #pragma warning(disable : 4127) #endif #include <gmp.h> #ifdef BOOST_MSVC #pragma warning(pop) #endif #if defined(__MPIR_VERSION) && defined(__MPIR_VERSION_MINOR) && defined(__MPIR_VERSION_PATCHLEVEL) #define BOOST_MP_MPIR_VERSION (__MPIR_VERSION * 10000 + __MPIR_VERSION_MINOR * 100 + __MPIR_VERSION_PATCHLEVEL) #else #define BOOST_MP_MPIR_VERSION 0 #endif namespace boost { namespace multiprecision { namespace backends { #ifdef BOOST_MSVC // warning C4127: conditional expression is constant #pragma warning(push) //#pragma warning(disable : 4127) #endif template <unsigned digits10> struct gmp_float; struct gmp_int; struct gmp_rational; } // namespace backends template <> struct number_category<backends::gmp_int> : public std::integral_constant<int, number_kind_integer> {}; template <> struct number_category<backends::gmp_rational> : public std::integral_constant<int, number_kind_rational> {}; template <unsigned digits10> struct number_category<backends::gmp_float<digits10> > : public std::integral_constant<int, number_kind_floating_point> {}; namespace backends { // // Within this file, the only functions we mark as noexcept are those that manipulate // (but don't create) an mpf_t. All other types may allocate at pretty much any time // via a user-supplied allocator, and therefore throw. // namespace detail { template <unsigned digits10> struct gmp_float_imp { #ifdef BOOST_HAS_LONG_LONG using signed_types = std::tuple<long, long long> ; using unsigned_types = std::tuple<unsigned long, unsigned long long>; #else using signed_types = std::tuple<long> ; using unsigned_types = std::tuple<unsigned long>; #endif using float_types = std::tuple<double, long double>; using exponent_type = long ; gmp_float_imp() noexcept { m_data[0]._mp_d = nullptr; // uninitialized m_data m_data[0]._mp_prec = 1; } gmp_float_imp(const gmp_float_imp& o) { // // We have to do an init followed by a set here, otherwise *this may be at // a lower precision than o: seems like mpf_init_set copies just enough bits // to get the right value, but if it's then used in further calculations // things go badly wrong!! // mpf_init2(m_data, preserve_source_precision() ? mpf_get_prec(o.data()) : boost::multiprecision::detail::digits10_2_2(get_default_precision())); if (o.m_data[0]._mp_d) mpf_set(m_data, o.m_data); } // rvalue copy gmp_float_imp(gmp_float_imp&& o) noexcept { if ((this->get_default_options() == variable_precision_options::preserve_target_precision) && (mpf_get_prec(o.data()) != boost::multiprecision::detail::digits10_2_2(get_default_precision()))) { mpf_init2(m_data, boost::multiprecision::detail::digits10_2_2(get_default_precision())); *this = static_cast<const gmp_float_imp&>(o); } else { m_data[0] = o.m_data[0]; o.m_data[0]._mp_d = nullptr; } } gmp_float_imp& operator=(const gmp_float_imp& o) { if (m_data[0]._mp_d == nullptr) { mpf_init2(m_data, preserve_source_precision() ? mpf_get_prec(o.data()) : boost::multiprecision::detail::digits10_2_2(get_default_precision())); mpf_set(m_data, o.m_data); } else if (preserve_source_precision() && (mpf_get_prec(data()) != mpf_get_prec(o.data()))) { mpf_t t; mpf_init2(t, mpf_get_prec(o.data())); mpf_set(t, o.data()); mpf_swap(data(), t); mpf_clear(t); } else { mpf_set(m_data, o.m_data); } return *this; } // rvalue assign gmp_float_imp& operator=(gmp_float_imp&& o) noexcept { if ((this->get_default_options() == variable_precision_options::preserve_target_precision) && (mpf_get_prec(o.data()) != mpf_get_prec(data()))) *this = static_cast<const gmp_float_imp&>(o); else { mpf_swap(m_data, o.m_data); } return *this; } #ifdef BOOST_HAS_LONG_LONG #if defined(ULLONG_MAX) && (ULLONG_MAX == ULONG_MAX) gmp_float_imp& operator=(unsigned long long i) { *this = static_cast<unsigned long>(i); return *this; } #else gmp_float_imp& operator=(unsigned long long i) { if (m_data[0]._mp_d == nullptr) { mpf_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); } unsigned long long mask = ((((1uLL << (std::numeric_limits<unsigned long>::digits - 1)) - 1) << 1) | 1uLL); unsigned shift = 0; mpf_t t; mpf_init2(t, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); mpf_set_ui(m_data, 0); while (i) { mpf_set_ui(t, static_cast<unsigned long>(i & mask)); if (shift) mpf_mul_2exp(t, t, shift); mpf_add(m_data, m_data, t); shift += std::numeric_limits<unsigned long>::digits; i >>= std::numeric_limits<unsigned long>::digits; } mpf_clear(t); return *this; } #endif gmp_float_imp& operator=(long long i) { if (m_data[0]._mp_d == nullptr) { mpf_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); } bool neg = i < 0; *this = static_cast<unsigned long long>(boost::multiprecision::detail::unsigned_abs(i)); if (neg) mpf_neg(m_data, m_data); return *this; } #endif gmp_float_imp& operator=(unsigned long i) { if (m_data[0]._mp_d == nullptr) { mpf_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); } mpf_set_ui(m_data, i); return *this; } gmp_float_imp& operator=(long i) { if (m_data[0]._mp_d == nullptr) { mpf_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); } mpf_set_si(m_data, i); return *this; } #ifdef BOOST_HAS_INT128 gmp_float_imp& operator=(uint128_type i) { if (m_data[0]._mp_d == nullptr) { mpf_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); } unsigned long mask = ((((1uLL << (std::numeric_limits<unsigned long>::digits - 1)) - 1) << 1) | 1uLL); unsigned shift = 0; mpf_t t; mpf_init2(t, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); mpf_set_ui(m_data, 0); while (i) { mpf_set_ui(t, static_cast<unsigned long>(i & mask)); if (shift) mpf_mul_2exp(t, t, shift); mpf_add(m_data, m_data, t); shift += std::numeric_limits<unsigned long>::digits; i >>= std::numeric_limits<unsigned long>::digits; } mpf_clear(t); return *this; } gmp_float_imp& operator=(int128_type i) { if (m_data[0]._mp_d == nullptr) { mpf_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); } bool neg = i < 0; *this = static_cast<uint128_type>(boost::multiprecision::detail::unsigned_abs(i)); if (neg) mpf_neg(m_data, m_data); return *this; } #endif gmp_float_imp& operator=(double d) { if (m_data[0]._mp_d == nullptr) { mpf_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); } mpf_set_d(m_data, d); return *this; } template <class F> gmp_float_imp& assign_float(F a) { BOOST_MP_FLOAT128_USING using std::floor; using std::frexp; using std::ldexp; if (m_data[0]._mp_d == nullptr) { mpf_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); } if (a == 0) { mpf_set_si(m_data, 0); return *this; } if (a == 1) { mpf_set_si(m_data, 1); return *this; } BOOST_MP_ASSERT(!BOOST_MP_ISINF(a)); BOOST_MP_ASSERT(!BOOST_MP_ISNAN(a)); int e; F f, term; mpf_set_ui(m_data, 0u); f = frexp(a, &e); constexpr int shift = std::numeric_limits<int>::digits - 1; while (f) { // extract int sized bits from f: f = ldexp(f, shift); term = floor(f); e -= shift; mpf_mul_2exp(m_data, m_data, shift); if (term > 0) mpf_add_ui(m_data, m_data, static_cast<unsigned>(term)); else mpf_sub_ui(m_data, m_data, static_cast<unsigned>(-term)); f -= term; } if (e > 0) mpf_mul_2exp(m_data, m_data, e); else if (e < 0) mpf_div_2exp(m_data, m_data, -e); return *this; } gmp_float_imp& operator=(long double a) { return assign_float(a); } #ifdef BOOST_HAS_FLOAT128 gmp_float_imp& operator= (float128_type a) { return assign_float(a); } #endif gmp_float_imp& operator=(const char* s) { if (m_data[0]._mp_d == nullptr) { mpf_init2(m_data, multiprecision::detail::digits10_2_2(digits10 ? digits10 : (unsigned)get_default_precision())); } if (s && (*s == '+')) ++s; // Leading "+" sign not supported by mpf_set_str: // // Validate the string as mpf_set_str does a poor job of this: // static const char* digits = "0123456789"; const char* p = s; if (*s == '-') ++s; s += boost::multiprecision::detail::find_first_not_of(s, s + std::strlen(s), digits); std::lconv const* l = std::localeconv(); std::size_t len = strlen(l->decimal_point); if (std::find(l->decimal_point, l->decimal_point + len, *s) != l->decimal_point + len) { ++s; s += boost::multiprecision::detail::find_first_not_of(s, s + std::strlen(s), digits); } if ((*s == 'e') || (*s == 'E')) { ++s; if ((*s == '+') || (*s == '-')) ++s; s += boost::multiprecision::detail::find_first_not_of(s, s + std::strlen(s), digits); } if(*s) BOOST_MP_THROW_EXCEPTION(std::runtime_error(std::string("The string \"") + s + std::string("\"could not be interpreted as a valid floating point number."))); s = p; if (0 != mpf_set_str(m_data, s, 10)) BOOST_MP_THROW_EXCEPTION(std::runtime_error(std::string("The string \"") + s + std::string("\"could not be interpreted as a valid floating point number."))); return *this; } void swap(gmp_float_imp& o) noexcept { mpf_swap(m_data, o.m_data); } std::string str(std::streamsize digits, std::ios_base::fmtflags f) const { BOOST_MP_ASSERT(m_data[0]._mp_d); bool scientific = (f & std::ios_base::scientific) == std::ios_base::scientific; bool fixed = (f & std::ios_base::fixed) == std::ios_base::fixed; std::streamsize org_digits(digits); if (scientific && digits) ++digits; std::string result; mp_exp_t e; void* (*alloc_func_ptr)(size_t); void* (*realloc_func_ptr)(void*, size_t, size_t); void (*free_func_ptr)(void*, size_t); mp_get_memory_functions(&alloc_func_ptr, &realloc_func_ptr, &free_func_ptr); if (mpf_sgn(m_data) == 0) { e = 0; result = "0"; if (fixed && digits) ++digits; } else { char* ps = mpf_get_str(nullptr, &e, 10, static_cast<std::size_t>(digits), m_data); --e; // To match with what our formatter expects. if (fixed) { // Oops we actually need a different number of digits to what we asked for: (*free_func_ptr)((void*)ps, std::strlen(ps) + 1); digits += e + 1; if (digits == 0) { // We need to get *all* the digits and then possibly round up, // we end up with either "0" or "1" as the result. ps = mpf_get_str(nullptr, &e, 10, 0, m_data); --e; unsigned offset = *ps == '-' ? 1 : 0; if (ps[offset] > '5') { ++e; ps[offset] = '1'; ps[offset + 1] = 0; } else if (ps[offset] == '5') { unsigned i = offset + 1; bool round_up = false; while (ps[i] != 0) { if (ps[i] != '0') { round_up = true; break; } ++i; } if (round_up) { ++e; ps[offset] = '1'; ps[offset + 1] = 0; } else { ps[offset] = '0'; ps[offset + 1] = 0; } } else { ps[offset] = '0'; ps[offset + 1] = 0; } } else if (digits > 0) { mp_exp_t old_e = e; ps = mpf_get_str(nullptr, &e, 10, static_cast<std::size_t>(digits), m_data); --e; // To match with what our formatter expects. if (old_e > e) { // in some cases, when we ask for more digits of precision, it will // change the number of digits to the left of the decimal, if that // happens, account for it here. // example: cout << fixed << setprecision(3) << mpf_float_50("99.9809") digits -= old_e - e; (*free_func_ptr)((void*)ps, std::strlen(ps) + 1); ps = mpf_get_str(nullptr, &e, 10, static_cast<std::size_t>(digits), m_data); --e; // To match with what our formatter expects. } } else { ps = mpf_get_str(nullptr, &e, 10, 1, m_data); --e; unsigned offset = *ps == '-' ? 1 : 0; ps[offset] = '0'; ps[offset + 1] = 0; } } result = ps; (*free_func_ptr)((void*)ps, std::strlen(ps) + 1); } boost::multiprecision::detail::format_float_string(result, e, org_digits, f, mpf_sgn(m_data) == 0); return result; } ~gmp_float_imp() noexcept { if (m_data[0]._mp_d) { mpf_clear(m_data); } } void negate() noexcept { BOOST_MP_ASSERT(m_data[0]._mp_d); mpf_neg(m_data, m_data); } int compare(const gmp_float<digits10>& o) const noexcept { BOOST_MP_ASSERT(m_data[0]._mp_d && o.m_data[0]._mp_d); return mpf_cmp(m_data, o.m_data); } int compare(long i) const noexcept { BOOST_MP_ASSERT(m_data[0]._mp_d); return mpf_cmp_si(m_data, i); } int compare(unsigned long i) const noexcept { BOOST_MP_ASSERT(m_data[0]._mp_d); return mpf_cmp_ui(m_data, i); } template <class V> typename std::enable_if<boost::multiprecision::detail::is_arithmetic<V>::value, int>::type compare(V v) const { gmp_float<digits10> d; d = v; return compare(d); } mpf_t& data() noexcept { BOOST_MP_ASSERT(m_data[0]._mp_d); return m_data; } const mpf_t& data() const noexcept { BOOST_MP_ASSERT(m_data[0]._mp_d); return m_data; } protected: mpf_t m_data; static unsigned& get_default_precision() noexcept { static BOOST_MP_THREAD_LOCAL unsigned val(get_global_default_precision()); return val; } static boost::multiprecision::detail::precision_type& get_global_default_precision() noexcept { static boost::multiprecision::detail::precision_type val(50); return val; } #ifndef BOOST_MT_NO_ATOMIC_INT static std::atomic<variable_precision_options>& get_global_default_options() noexcept #else static variable_precision_options& get_global_default_options() noexcept #endif { #ifndef BOOST_MT_NO_ATOMIC_INT static std::atomic<variable_precision_options> val{variable_precision_options::preserve_related_precision}; #else static variable_precision_options val{variable_precision_options::preserve_related_precision}; #endif return val; } static variable_precision_options& get_default_options()noexcept { static BOOST_MP_THREAD_LOCAL variable_precision_options val(get_global_default_options()); return val; } static bool preserve_source_precision() noexcept { return get_default_options() >= variable_precision_options::preserve_source_precision; } }; class gmp_char_ptr { private: char* ptr_val; void* (*alloc_func_ptr)(size_t); void* (*realloc_func_ptr)(void*, size_t, size_t); void (*free_func_ptr)(void*, size_t); public: gmp_char_ptr() = delete; explicit gmp_char_ptr(char* val_) : ptr_val {val_} { mp_get_memory_functions(&alloc_func_ptr, &realloc_func_ptr, &free_func_ptr); } ~gmp_char_ptr() noexcept { (*free_func_ptr)((void*)ptr_val, sizeof(*ptr_val)); ptr_val = nullptr; } inline char* get() noexcept { return ptr_val; } }; } // namespace detail template <unsigned digits10> struct gmp_float : public detail::gmp_float_imp<digits10> { gmp_float() { mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(digits10)); } gmp_float(const gmp_float& o) : detail::gmp_float_imp<digits10>(o) {} template <unsigned D> gmp_float(const gmp_float<D>& o, typename std::enable_if<D <= digits10>::type* = nullptr); template <unsigned D> explicit gmp_float(const gmp_float<D>& o, typename std::enable_if<!(D <= digits10)>::type* = nullptr); gmp_float(const gmp_int& o); gmp_float(const gmp_rational& o); gmp_float(const mpf_t val) { mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(digits10)); mpf_set(this->m_data, val); } gmp_float(const mpz_t val) { mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(digits10)); mpf_set_z(this->m_data, val); } gmp_float(const mpq_t val) { mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(digits10)); mpf_set_q(this->m_data, val); } // rvalue copy gmp_float(gmp_float&& o) noexcept : detail::gmp_float_imp<digits10>(static_cast<detail::gmp_float_imp<digits10>&&>(o)) {} gmp_float& operator=(const gmp_float& o) { *static_cast<detail::gmp_float_imp<digits10>*>(this) = static_cast<detail::gmp_float_imp<digits10> const&>(o); return *this; } gmp_float& operator=(gmp_float&& o) noexcept { *static_cast<detail::gmp_float_imp<digits10>*>(this) = static_cast<detail::gmp_float_imp<digits10>&&>(o); return *this; } template <unsigned D> gmp_float& operator=(const gmp_float<D>& o); gmp_float& operator=(const gmp_int& o); gmp_float& operator=(const gmp_rational& o); gmp_float& operator=(const mpf_t val) { if (this->m_data[0]._mp_d == nullptr) mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(digits10)); mpf_set(this->m_data, val); return *this; } gmp_float& operator=(const mpz_t val) { if (this->m_data[0]._mp_d == nullptr) mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(digits10)); mpf_set_z(this->m_data, val); return *this; } gmp_float& operator=(const mpq_t val) { if (this->m_data[0]._mp_d == nullptr) mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(digits10)); mpf_set_q(this->m_data, val); return *this; } template <class V> typename std::enable_if<std::is_assignable<detail::gmp_float_imp<digits10>, V>::value, gmp_float&>::type operator=(const V& v) { *static_cast<detail::gmp_float_imp<digits10>*>(this) = v; return *this; } }; template <> struct gmp_float<0> : public detail::gmp_float_imp<0> { // // We have a problem with mpf_t in that the precision we request isn't what we get. // As a result the front end can end up chasing it's tail trying to create a variable // with the the correct precision to hold the result of an expression. // See: https://github.com/boostorg/multiprecision/issues/164 // The problem is made worse by the fact that our conversions from base10 to 2 and // vice-versa do not exactly round trip (and probably never will). // The workaround is to keep track of the precision requested, and always return // that as the current actual precision. // private: unsigned requested_precision; public: gmp_float() : requested_precision(get_default_precision()) { mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(requested_precision)); } gmp_float(const mpf_t val) : requested_precision(get_default_precision()) { mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(requested_precision)); mpf_set(this->m_data, val); } gmp_float(const mpz_t val) : requested_precision(get_default_precision()) { mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(requested_precision)); mpf_set_z(this->m_data, val); } gmp_float(const mpq_t val) : requested_precision(get_default_precision()) { mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(requested_precision)); mpf_set_q(this->m_data, val); } gmp_float(const gmp_float& o) : detail::gmp_float_imp<0>(o), requested_precision(preserve_source_precision() ? o.requested_precision : get_default_precision()) {} template <unsigned D> gmp_float(const gmp_float<D>& o) { mpf_init2(this->m_data, preserve_related_precision() ? mpf_get_prec(o.data()) : multiprecision::detail::digits10_2_2(get_default_precision())); mpf_set(this->m_data, o.data()); requested_precision = preserve_related_precision() ? D : get_default_precision(); } // rvalue copy gmp_float(gmp_float&& o) noexcept : detail::gmp_float_imp<0>(static_cast<detail::gmp_float_imp<0>&&>(o)), requested_precision((this->get_default_options() != variable_precision_options::preserve_target_precision) ? o.requested_precision : get_default_precision()) {} gmp_float(const gmp_int& o); gmp_float(const gmp_rational& o); gmp_float(const gmp_float& o, unsigned digits10) : requested_precision(digits10) { mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(digits10)); mpf_set(this->m_data, o.data()); } template <class V> gmp_float(const V& o, unsigned digits10) : requested_precision(digits10) { mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(digits10)); *this = o; } #ifndef BOOST_NO_CXX17_HDR_STRING_VIEW // // Support for new types in C++17 // template <class Traits> gmp_float(const std::basic_string_view<char, Traits>& o, unsigned digits10) : requested_precision(digits10) { using default_ops::assign_from_string_view; mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(digits10)); assign_from_string_view(*this, o); } #endif gmp_float& operator=(const gmp_float& o) { *static_cast<detail::gmp_float_imp<0>*>(this) = static_cast<detail::gmp_float_imp<0> const&>(o); if(preserve_source_precision()) requested_precision = o.requested_precision; return *this; } // rvalue copy gmp_float& operator=(gmp_float&& o) noexcept { *static_cast<detail::gmp_float_imp<0>*>(this) = static_cast<detail::gmp_float_imp<0>&&>(o); if ((this->get_default_options() != variable_precision_options::preserve_target_precision)) requested_precision = o.requested_precision; return *this; } template <unsigned D> gmp_float& operator=(const gmp_float<D>& o) { if (this->m_data[0]._mp_d == nullptr) { mpf_init2(this->m_data, preserve_related_precision() ? mpf_get_prec(o.data()) : multiprecision::detail::digits10_2_2(get_default_precision())); } else if(preserve_related_precision()) { mpf_set_prec(this->m_data, mpf_get_prec(o.data())); } mpf_set(this->m_data, o.data()); if (preserve_related_precision()) requested_precision = D; return *this; } gmp_float& operator=(const gmp_int& o); gmp_float& operator=(const gmp_rational& o); gmp_float& operator=(const mpf_t val) { if (this->m_data[0]._mp_d == nullptr) { requested_precision = get_default_precision(); mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(requested_precision)); } mpf_set(this->m_data, val); return *this; } gmp_float& operator=(const mpz_t val) { if (this->m_data[0]._mp_d == nullptr) { requested_precision = get_default_precision(); mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(requested_precision)); } mpf_set_z(this->m_data, val); return *this; } gmp_float& operator=(const mpq_t val) { if (this->m_data[0]._mp_d == nullptr) { requested_precision = get_default_precision(); mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(requested_precision)); } mpf_set_q(this->m_data, val); return *this; } template <class V> typename std::enable_if<std::is_assignable<detail::gmp_float_imp<0>, V>::value, gmp_float&>::type operator=(const V& v) { constexpr unsigned d10 = std::is_floating_point<V>::value ? std::numeric_limits<V>::digits10 : std::numeric_limits<V>::digits10 ? 1 + std::numeric_limits<V>::digits10 : 1 + boost::multiprecision::detail::digits2_2_10(std::numeric_limits<V>::digits); if((thread_default_variable_precision_options() >= variable_precision_options::preserve_all_precision) && (precision() < d10)) this->precision(d10); *static_cast<detail::gmp_float_imp<0>*>(this) = v; return *this; } static unsigned default_precision() noexcept { return get_global_default_precision(); } static void default_precision(unsigned v) noexcept { get_global_default_precision() = v; } static unsigned thread_default_precision() noexcept { return get_default_precision(); } static void thread_default_precision(unsigned v) noexcept { get_default_precision() = v; } unsigned precision() const noexcept { return requested_precision; } void precision(unsigned digits10) noexcept { requested_precision = digits10; mpf_set_prec(this->m_data, multiprecision::detail::digits10_2_2(requested_precision)); } // // Variable precision options: // static variable_precision_options default_variable_precision_options()noexcept { return get_global_default_options(); } static variable_precision_options thread_default_variable_precision_options()noexcept { return get_default_options(); } static void default_variable_precision_options(variable_precision_options opts) { get_global_default_options() = opts; } static void thread_default_variable_precision_options(variable_precision_options opts) { get_default_options() = opts; } static bool preserve_source_precision() { return get_default_options() >= variable_precision_options::preserve_source_precision; } static bool preserve_related_precision() { return get_default_options() >= variable_precision_options::preserve_related_precision; } static bool preserve_all_precision() { return get_default_options() >= variable_precision_options::preserve_all_precision; } // // swap: // void swap(gmp_float& o) { std::swap(requested_precision, o.requested_precision); gmp_float_imp<0>::swap(o); } }; template <unsigned digits10, class T> inline typename std::enable_if<boost::multiprecision::detail::is_arithmetic<T>::value, bool>::type eval_eq(const gmp_float<digits10>& a, const T& b) noexcept { return a.compare(b) == 0; } template <unsigned digits10, class T> inline typename std::enable_if<boost::multiprecision::detail::is_arithmetic<T>::value, bool>::type eval_lt(const gmp_float<digits10>& a, const T& b) noexcept { return a.compare(b) < 0; } template <unsigned digits10, class T> inline typename std::enable_if<boost::multiprecision::detail::is_arithmetic<T>::value, bool>::type eval_gt(const gmp_float<digits10>& a, const T& b) noexcept { return a.compare(b) > 0; } template <unsigned D1, unsigned D2> inline void eval_add(gmp_float<D1>& result, const gmp_float<D2>& o) { mpf_add(result.data(), result.data(), o.data()); } template <unsigned D1, unsigned D2> inline void eval_subtract(gmp_float<D1>& result, const gmp_float<D2>& o) { mpf_sub(result.data(), result.data(), o.data()); } template <unsigned D1, unsigned D2> inline void eval_multiply(gmp_float<D1>& result, const gmp_float<D2>& o) { mpf_mul(result.data(), result.data(), o.data()); } template <unsigned digits10> inline bool eval_is_zero(const gmp_float<digits10>& val) noexcept { return mpf_sgn(val.data()) == 0; } template <unsigned D1, unsigned D2> inline void eval_divide(gmp_float<D1>& result, const gmp_float<D2>& o) { if (eval_is_zero(o)) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); mpf_div(result.data(), result.data(), o.data()); } template <unsigned digits10> inline void eval_add(gmp_float<digits10>& result, unsigned long i) { mpf_add_ui(result.data(), result.data(), i); } template <unsigned digits10> inline void eval_subtract(gmp_float<digits10>& result, unsigned long i) { mpf_sub_ui(result.data(), result.data(), i); } template <unsigned digits10> inline void eval_multiply(gmp_float<digits10>& result, unsigned long i) { mpf_mul_ui(result.data(), result.data(), i); } template <unsigned digits10> inline void eval_divide(gmp_float<digits10>& result, unsigned long i) { if (i == 0) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); mpf_div_ui(result.data(), result.data(), i); } template <unsigned digits10> inline void eval_add(gmp_float<digits10>& result, long i) { using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; if (i > 0) mpf_add_ui(result.data(), result.data(), static_cast<local_uint_type>(i)); else if (i < 0) mpf_sub_ui(result.data(), result.data(), static_cast<local_uint_type>(-i)); } template <unsigned digits10> inline void eval_subtract(gmp_float<digits10>& result, long i) { using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; if (i > 0) mpf_sub_ui(result.data(), result.data(), static_cast<local_uint_type>(i)); else if (i < 0) mpf_add_ui(result.data(), result.data(), static_cast<local_uint_type>(-i)); } template <unsigned digits10> inline void eval_multiply(gmp_float<digits10>& result, long i) { using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; mpf_mul_ui(result.data(), result.data(), static_cast<local_uint_type>(boost::multiprecision::detail::unsigned_abs(i))); if (i < 0) mpf_neg(result.data(), result.data()); } template <unsigned digits10> inline void eval_divide(gmp_float<digits10>& result, long i) { if (i == 0) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; mpf_div_ui(result.data(), result.data(), static_cast<local_uint_type>(boost::multiprecision::detail::unsigned_abs(i))); if (i < 0) mpf_neg(result.data(), result.data()); } // // Specialised 3 arg versions of the basic operators: // template <unsigned D1, unsigned D2, unsigned D3> inline void eval_add(gmp_float<D1>& a, const gmp_float<D2>& x, const gmp_float<D3>& y) { mpf_add(a.data(), x.data(), y.data()); } template <unsigned D1, unsigned D2> inline void eval_add(gmp_float<D1>& a, const gmp_float<D2>& x, unsigned long y) { mpf_add_ui(a.data(), x.data(), y); } template <unsigned D1, unsigned D2> inline void eval_add(gmp_float<D1>& a, const gmp_float<D2>& x, long y) { if (y < 0) mpf_sub_ui(a.data(), x.data(), boost::multiprecision::detail::unsigned_abs(y)); else mpf_add_ui(a.data(), x.data(), y); } template <unsigned D1, unsigned D2> inline void eval_add(gmp_float<D1>& a, unsigned long x, const gmp_float<D2>& y) { mpf_add_ui(a.data(), y.data(), x); } template <unsigned D1, unsigned D2> inline void eval_add(gmp_float<D1>& a, long x, const gmp_float<D2>& y) { using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; if (x < 0) { mpf_ui_sub(a.data(), static_cast<local_uint_type>(-x), y.data()); mpf_neg(a.data(), a.data()); } else mpf_add_ui(a.data(), y.data(), static_cast<local_uint_type>(x)); } template <unsigned D1, unsigned D2, unsigned D3> inline void eval_subtract(gmp_float<D1>& a, const gmp_float<D2>& x, const gmp_float<D3>& y) { mpf_sub(a.data(), x.data(), y.data()); } template <unsigned D1, unsigned D2> inline void eval_subtract(gmp_float<D1>& a, const gmp_float<D2>& x, unsigned long y) { mpf_sub_ui(a.data(), x.data(), y); } template <unsigned D1, unsigned D2> inline void eval_subtract(gmp_float<D1>& a, const gmp_float<D2>& x, long y) { using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; if (y < 0) mpf_add_ui(a.data(), x.data(), static_cast<local_uint_type>(-y)); else mpf_sub_ui(a.data(), x.data(), static_cast<local_uint_type>(y)); } template <unsigned D1, unsigned D2> inline void eval_subtract(gmp_float<D1>& a, unsigned long x, const gmp_float<D2>& y) { mpf_ui_sub(a.data(), x, y.data()); } template <unsigned D1, unsigned D2> inline void eval_subtract(gmp_float<D1>& a, long x, const gmp_float<D2>& y) { using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; if (x < 0) { mpf_add_ui(a.data(), y.data(), static_cast<local_uint_type>(-x)); mpf_neg(a.data(), a.data()); } else mpf_ui_sub(a.data(), static_cast<local_uint_type>(x), y.data()); } template <unsigned D1, unsigned D2, unsigned D3> inline void eval_multiply(gmp_float<D1>& a, const gmp_float<D2>& x, const gmp_float<D3>& y) { mpf_mul(a.data(), x.data(), y.data()); } template <unsigned D1, unsigned D2> inline void eval_multiply(gmp_float<D1>& a, const gmp_float<D2>& x, unsigned long y) { mpf_mul_ui(a.data(), x.data(), y); } template <unsigned D1, unsigned D2> inline void eval_multiply(gmp_float<D1>& a, const gmp_float<D2>& x, long y) { using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; if (y < 0) { mpf_mul_ui(a.data(), x.data(), static_cast<local_uint_type>(-y)); a.negate(); } else mpf_mul_ui(a.data(), x.data(), static_cast<local_uint_type>(y)); } template <unsigned D1, unsigned D2> inline void eval_multiply(gmp_float<D1>& a, unsigned long x, const gmp_float<D2>& y) { mpf_mul_ui(a.data(), y.data(), x); } template <unsigned D1, unsigned D2> inline void eval_multiply(gmp_float<D1>& a, long x, const gmp_float<D2>& y) { using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; if (x < 0) { mpf_mul_ui(a.data(), y.data(), static_cast<local_uint_type>(-x)); mpf_neg(a.data(), a.data()); } else mpf_mul_ui(a.data(), y.data(), static_cast<local_uint_type>(x)); } template <unsigned D1, unsigned D2, unsigned D3> inline void eval_divide(gmp_float<D1>& a, const gmp_float<D2>& x, const gmp_float<D3>& y) { if (eval_is_zero(y)) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); mpf_div(a.data(), x.data(), y.data()); } template <unsigned D1, unsigned D2> inline void eval_divide(gmp_float<D1>& a, const gmp_float<D2>& x, unsigned long y) { if (y == 0) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); mpf_div_ui(a.data(), x.data(), y); } template <unsigned D1, unsigned D2> inline void eval_divide(gmp_float<D1>& a, const gmp_float<D2>& x, long y) { if (y == 0) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; if (y < 0) { mpf_div_ui(a.data(), x.data(), static_cast<local_uint_type>(-y)); a.negate(); } else mpf_div_ui(a.data(), x.data(), static_cast<local_uint_type>(y)); } template <unsigned D1, unsigned D2> inline void eval_divide(gmp_float<D1>& a, unsigned long x, const gmp_float<D2>& y) { if (eval_is_zero(y)) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); mpf_ui_div(a.data(), x, y.data()); } template <unsigned D1, unsigned D2> inline void eval_divide(gmp_float<D1>& a, long x, const gmp_float<D2>& y) { if (eval_is_zero(y)) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); if (x < 0) { mpf_ui_div(a.data(), boost::multiprecision::detail::unsigned_abs(x), y.data()); mpf_neg(a.data(), a.data()); } else { using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; mpf_ui_div(a.data(), static_cast<local_uint_type>(x), y.data()); } } template <unsigned digits10> inline int eval_get_sign(const gmp_float<digits10>& val) noexcept { return mpf_sgn(val.data()); } template <unsigned digits10> inline void eval_convert_to(unsigned long* result, const gmp_float<digits10>& val) noexcept { if (0 == mpf_fits_ulong_p(val.data())) *result = (std::numeric_limits<unsigned long>::max)(); else *result = static_cast<unsigned long>(mpf_get_ui(val.data())); } template <unsigned digits10> inline void eval_convert_to(long* result, const gmp_float<digits10>& val) noexcept { if (0 == mpf_fits_slong_p(val.data())) { *result = (std::numeric_limits<long>::max)(); *result *= mpf_sgn(val.data()); } else *result = static_cast<long>(mpf_get_si(val.data())); } #ifdef BOOST_MP_STANDALONE template <unsigned digits10> inline void eval_convert_to(long double* result, const gmp_float<digits10>& val) noexcept { mp_exp_t exp = 0; detail::gmp_char_ptr val_char_ptr {mpf_get_str(nullptr, &exp, 10, LDBL_DIG, val.data())}; auto temp_string = std::string(val_char_ptr.get()); if(exp > 0 && static_cast<std::size_t>(exp) < temp_string.size()) { if(temp_string.front() == '-') { ++exp; } temp_string.insert(static_cast<std::size_t>(exp), static_cast<std::size_t>(1u), '.'); } *result = std::strtold(temp_string.c_str(), nullptr); if((temp_string.size() == 2ul && *result < 0.0l) || (static_cast<std::size_t>(exp) > temp_string.size())) { *result *= std::pow(10l, exp-1); } } #endif // BOOST_MP_STANDALONE template <unsigned digits10> inline void eval_convert_to(double* result, const gmp_float<digits10>& val) noexcept { *result = mpf_get_d(val.data()); } #ifdef BOOST_HAS_LONG_LONG template <unsigned digits10> inline void eval_convert_to(long long* result, const gmp_float<digits10>& val) { gmp_float<digits10> t(val); if (eval_get_sign(t) < 0) t.negate(); long digits = std::numeric_limits<long long>::digits - std::numeric_limits<long>::digits; if (digits > 0) mpf_div_2exp(t.data(), t.data(), digits); if (!mpf_fits_slong_p(t.data())) { if (eval_get_sign(val) < 0) *result = (std::numeric_limits<long long>::min)(); else *result = (std::numeric_limits<long long>::max)(); return; }; *result = mpf_get_si(t.data()); while (digits > 0) { *result <<= digits; digits -= std::numeric_limits<unsigned long>::digits; mpf_mul_2exp(t.data(), t.data(), digits >= 0 ? std::numeric_limits<unsigned long>::digits : std::numeric_limits<unsigned long>::digits + digits); unsigned long l = static_cast<unsigned long>(mpf_get_ui(t.data())); if (digits < 0) l >>= -digits; *result |= l; } if (eval_get_sign(val) < 0) *result = -*result; } template <unsigned digits10> inline void eval_convert_to(unsigned long long* result, const gmp_float<digits10>& val) { gmp_float<digits10> t(val); long digits = std::numeric_limits<long long>::digits - std::numeric_limits<long>::digits; if (digits > 0) mpf_div_2exp(t.data(), t.data(), digits); if (!mpf_fits_ulong_p(t.data())) { *result = (std::numeric_limits<long long>::max)(); return; } *result = mpf_get_ui(t.data()); while (digits > 0) { *result <<= digits; digits -= std::numeric_limits<unsigned long>::digits; mpf_mul_2exp(t.data(), t.data(), digits >= 0 ? std::numeric_limits<unsigned long>::digits : std::numeric_limits<unsigned long>::digits + digits); unsigned long l = static_cast<unsigned long>(mpf_get_ui(t.data())); if (digits < 0) l >>= -digits; *result |= l; } } #endif #ifdef BOOST_HAS_FLOAT128 template <unsigned digits10> inline void eval_convert_to(float128_type* result, const gmp_float<digits10>& val) { *result = float128_procs::strtoflt128(val.str(0, std::ios_base::scientific).c_str(), nullptr); } #endif // // Native non-member operations: // template <unsigned Digits10> inline void eval_sqrt(gmp_float<Digits10>& result, const gmp_float<Digits10>& val) { mpf_sqrt(result.data(), val.data()); } template <unsigned Digits10> inline void eval_abs(gmp_float<Digits10>& result, const gmp_float<Digits10>& val) { mpf_abs(result.data(), val.data()); } template <unsigned Digits10> inline void eval_fabs(gmp_float<Digits10>& result, const gmp_float<Digits10>& val) { mpf_abs(result.data(), val.data()); } template <unsigned Digits10> inline void eval_ceil(gmp_float<Digits10>& result, const gmp_float<Digits10>& val) { mpf_ceil(result.data(), val.data()); } template <unsigned Digits10> inline void eval_floor(gmp_float<Digits10>& result, const gmp_float<Digits10>& val) { mpf_floor(result.data(), val.data()); } template <unsigned Digits10> inline void eval_trunc(gmp_float<Digits10>& result, const gmp_float<Digits10>& val) { mpf_trunc(result.data(), val.data()); } template <unsigned Digits10> inline void eval_ldexp(gmp_float<Digits10>& result, const gmp_float<Digits10>& val, long e) { if (e > 0) mpf_mul_2exp(result.data(), val.data(), static_cast<mp_bitcnt_t>(e)); else if (e < 0) mpf_div_2exp(result.data(), val.data(), static_cast<mp_bitcnt_t>(-e)); else result = val; } template <unsigned Digits10> inline void eval_frexp(gmp_float<Digits10>& result, const gmp_float<Digits10>& val, int* e) { #if (BOOST_MP_MPIR_VERSION >= 20600) && (BOOST_MP_MPIR_VERSION < 30000) mpir_si v; mpf_get_d_2exp(&v, val.data()); #else long v; mpf_get_d_2exp(&v, val.data()); #endif *e = static_cast<int>(v); eval_ldexp(result, val, -v); } template <unsigned Digits10> inline void eval_frexp(gmp_float<Digits10>& result, const gmp_float<Digits10>& val, long* e) { #if (BOOST_MP_MPIR_VERSION >= 20600) && (BOOST_MP_MPIR_VERSION < 30000) mpir_si v; mpf_get_d_2exp(&v, val.data()); *e = v; eval_ldexp(result, val, -v); #else mpf_get_d_2exp(e, val.data()); eval_ldexp(result, val, -*e); #endif } template <unsigned Digits10> inline std::size_t hash_value(const gmp_float<Digits10>& val) { std::size_t result = 0; for (int i = 0; i < std::abs(val.data()[0]._mp_size); ++i) boost::multiprecision::detail::hash_combine(result, val.data()[0]._mp_d[i]); boost::multiprecision::detail::hash_combine(result, val.data()[0]._mp_exp, val.data()[0]._mp_size); return result; } struct gmp_int { #ifdef BOOST_HAS_LONG_LONG using signed_types = std::tuple<long, long long> ; using unsigned_types = std::tuple<unsigned long, unsigned long long>; #else using signed_types = std::tuple<long> ; using unsigned_types = std::tuple<unsigned long>; #endif using float_types = std::tuple<double, long double>; gmp_int() { mpz_init(this->m_data); } gmp_int(const gmp_int& o) { if (o.m_data[0]._mp_d) mpz_init_set(m_data, o.m_data); else mpz_init(this->m_data); } // rvalue gmp_int(gmp_int&& o) noexcept { m_data[0] = o.m_data[0]; o.m_data[0]._mp_d = nullptr; } explicit gmp_int(const mpf_t val) { mpz_init(this->m_data); mpz_set_f(this->m_data, val); } gmp_int(const mpz_t val) { mpz_init_set(this->m_data, val); } gmp_int(long i) { mpz_init_set_si(this->m_data, i); } gmp_int(unsigned long i) { mpz_init_set_ui(this->m_data, i); } explicit gmp_int(const mpq_t val) { mpz_init(this->m_data); mpz_set_q(this->m_data, val); } template <unsigned Digits10> explicit gmp_int(const gmp_float<Digits10>& o) { mpz_init(this->m_data); mpz_set_f(this->m_data, o.data()); } explicit gmp_int(const gmp_rational& o); gmp_int& operator=(const gmp_int& o) { if (m_data[0]._mp_d == nullptr) mpz_init(this->m_data); mpz_set(m_data, o.m_data); return *this; } // rvalue copy gmp_int& operator=(gmp_int&& o) noexcept { mpz_swap(m_data, o.m_data); return *this; } #ifdef BOOST_HAS_LONG_LONG #if defined(ULLONG_MAX) && (ULLONG_MAX == ULONG_MAX) gmp_int& operator=(unsigned long long i) { *this = static_cast<unsigned long>(i); return *this; } #else gmp_int& operator=(unsigned long long i) { if (m_data[0]._mp_d == nullptr) mpz_init(this->m_data); unsigned long long mask = ((((1uLL << (std::numeric_limits<unsigned long>::digits - 1)) - 1) << 1) | 1uLL); unsigned shift = 0; mpz_t t; mpz_set_ui(m_data, 0); mpz_init_set_ui(t, 0); while (i) { mpz_set_ui(t, static_cast<unsigned long>(i & mask)); if (shift) mpz_mul_2exp(t, t, shift); mpz_add(m_data, m_data, t); shift += std::numeric_limits<unsigned long>::digits; i >>= std::numeric_limits<unsigned long>::digits; } mpz_clear(t); return *this; } #endif gmp_int& operator=(long long i) { if (m_data[0]._mp_d == nullptr) mpz_init(this->m_data); bool neg = i < 0; *this = boost::multiprecision::detail::unsigned_abs(i); if (neg) mpz_neg(m_data, m_data); return *this; } #endif #ifdef BOOST_HAS_INT128 gmp_int& operator=(uint128_type i) { if (m_data[0]._mp_d == nullptr) mpz_init(this->m_data); uint128_type mask = ((((1uLL << (std::numeric_limits<unsigned long>::digits - 1)) - 1) << 1) | 1uLL); unsigned shift = 0; mpz_t t; mpz_set_ui(m_data, 0); mpz_init_set_ui(t, 0); while (i) { mpz_set_ui(t, static_cast<unsigned long>(i & mask)); if (shift) mpz_mul_2exp(t, t, shift); mpz_add(m_data, m_data, t); shift += std::numeric_limits<unsigned long>::digits; i >>= std::numeric_limits<unsigned long>::digits; } mpz_clear(t); return *this; } gmp_int& operator=(int128_type i) { if (m_data[0]._mp_d == nullptr) mpz_init(this->m_data); bool neg = i < 0; *this = boost::multiprecision::detail::unsigned_abs(i); if (neg) mpz_neg(m_data, m_data); return *this; } #endif gmp_int& operator=(unsigned long i) { if (m_data[0]._mp_d == nullptr) mpz_init(this->m_data); mpz_set_ui(m_data, i); return *this; } gmp_int& operator=(long i) { if (m_data[0]._mp_d == nullptr) mpz_init(this->m_data); mpz_set_si(m_data, i); return *this; } gmp_int& operator=(double d) { if (m_data[0]._mp_d == nullptr) mpz_init(this->m_data); mpz_set_d(m_data, d); return *this; } template <class F> gmp_int& assign_float(F a) { BOOST_MP_FLOAT128_USING using std::floor; using std::frexp; using std::ldexp; if (m_data[0]._mp_d == nullptr) mpz_init(this->m_data); if (a == 0) { mpz_set_si(m_data, 0); return *this; } if (a == 1) { mpz_set_si(m_data, 1); return *this; } BOOST_MP_ASSERT(!BOOST_MP_ISINF(a)); BOOST_MP_ASSERT(!BOOST_MP_ISNAN(a)); int e; F f, term; mpz_set_ui(m_data, 0u); f = frexp(a, &e); constexpr int shift = std::numeric_limits<int>::digits - 1; while (f != static_cast<F>(0.0f)) { // extract int sized bits from f: f = ldexp(f, shift); term = floor(f); e -= shift; mpz_mul_2exp(m_data, m_data, shift); if (term > 0) mpz_add_ui(m_data, m_data, static_cast<unsigned>(term)); else mpz_sub_ui(m_data, m_data, static_cast<unsigned>(-term)); f -= term; } if (e > 0) mpz_mul_2exp(m_data, m_data, static_cast<mp_bitcnt_t>(e)); else if (e < 0) mpz_div_2exp(m_data, m_data, static_cast<mp_bitcnt_t>(-e)); return *this; } gmp_int& operator=(long double a) { return assign_float(a); } gmp_int& operator=(const char* s) { if (m_data[0]._mp_d == nullptr) mpz_init(this->m_data); std::size_t n = s ? std::strlen(s) : 0; int radix = 10; if (n && (*s == '0')) { if ((n > 1) && ((s[1] == 'x') || (s[1] == 'X'))) { radix = 16; s += 2; n -= 2; } else { radix = 8; n -= 1; } } if (n) { if (0 != mpz_set_str(m_data, s, radix)) BOOST_MP_THROW_EXCEPTION(std::runtime_error(std::string("The string \"") + s + std::string("\"could not be interpreted as a valid integer."))); } else mpz_set_ui(m_data, 0); return *this; } #ifdef BOOST_HAS_FLOAT128 gmp_int& operator=(float128_type a) { return assign_float(a); } #endif gmp_int& operator=(const mpf_t val) { if (m_data[0]._mp_d == nullptr) mpz_init(this->m_data); mpz_set_f(this->m_data, val); return *this; } gmp_int& operator=(const mpz_t val) { if (m_data[0]._mp_d == nullptr) mpz_init(this->m_data); mpz_set(this->m_data, val); return *this; } gmp_int& operator=(const mpq_t val) { if (m_data[0]._mp_d == nullptr) mpz_init(this->m_data); mpz_set_q(this->m_data, val); return *this; } template <unsigned Digits10> gmp_int& operator=(const gmp_float<Digits10>& o) { if (m_data[0]._mp_d == nullptr) mpz_init(this->m_data); mpz_set_f(this->m_data, o.data()); return *this; } gmp_int& operator=(const gmp_rational& o); void swap(gmp_int& o) { mpz_swap(m_data, o.m_data); } std::string str(std::streamsize /*digits*/, std::ios_base::fmtflags f) const { BOOST_MP_ASSERT(m_data[0]._mp_d); int base = 10; if ((f & std::ios_base::oct) == std::ios_base::oct) base = 8; else if ((f & std::ios_base::hex) == std::ios_base::hex) base = 16; // // sanity check, bases 8 and 16 are only available for positive numbers: // if ((base != 10) && (mpz_sgn(m_data) < 0)) BOOST_MP_THROW_EXCEPTION(std::runtime_error("Formatted output in bases 8 or 16 is only available for positive numbers")); void* (*alloc_func_ptr)(size_t); void* (*realloc_func_ptr)(void*, size_t, size_t); void (*free_func_ptr)(void*, size_t); const char* ps = mpz_get_str(nullptr, base, m_data); std::string s = ps; mp_get_memory_functions(&alloc_func_ptr, &realloc_func_ptr, &free_func_ptr); (*free_func_ptr)((void*)ps, std::strlen(ps) + 1); if (f & std::ios_base::uppercase) for (size_t i = 0; i < s.length(); ++i) s[i] = static_cast<char>(std::toupper(s[i])); if ((base != 10) && (f & std::ios_base::showbase)) { int pos = s[0] == '-' ? 1 : 0; const char* pp = base == 8 ? "0" : (f & std::ios_base::uppercase) ? "0X" : "0x"; s.insert(static_cast<std::string::size_type>(pos), pp); } if ((f & std::ios_base::showpos) && (s[0] != '-')) s.insert(static_cast<std::string::size_type>(0), 1, '+'); return s; } ~gmp_int() noexcept { if (m_data[0]._mp_d) mpz_clear(m_data); } void negate() noexcept { BOOST_MP_ASSERT(m_data[0]._mp_d); mpz_neg(m_data, m_data); } int compare(const gmp_int& o) const noexcept { BOOST_MP_ASSERT(m_data[0]._mp_d && o.m_data[0]._mp_d); return mpz_cmp(m_data, o.m_data); } int compare(long i) const noexcept { BOOST_MP_ASSERT(m_data[0]._mp_d); return mpz_cmp_si(m_data, i); } int compare(unsigned long i) const noexcept { BOOST_MP_ASSERT(m_data[0]._mp_d); return mpz_cmp_ui(m_data, i); } template <class V> int compare(V v) const { gmp_int d; d = v; return compare(d); } mpz_t& data() noexcept { BOOST_MP_ASSERT(m_data[0]._mp_d); return m_data; } const mpz_t& data() const noexcept { BOOST_MP_ASSERT(m_data[0]._mp_d); return m_data; } protected: mpz_t m_data; }; template <class T> inline typename std::enable_if<boost::multiprecision::detail::is_arithmetic<T>::value, bool>::type eval_eq(const gmp_int& a, const T& b) { return a.compare(b) == 0; } template <class T> inline typename std::enable_if<boost::multiprecision::detail::is_arithmetic<T>::value, bool>::type eval_lt(const gmp_int& a, const T& b) { return a.compare(b) < 0; } template <class T> inline typename std::enable_if<boost::multiprecision::detail::is_arithmetic<T>::value, bool>::type eval_gt(const gmp_int& a, const T& b) { return a.compare(b) > 0; } inline bool eval_is_zero(const gmp_int& val) { return mpz_sgn(val.data()) == 0; } inline void eval_add(gmp_int& t, const gmp_int& o) { mpz_add(t.data(), t.data(), o.data()); } inline void eval_multiply_add(gmp_int& t, const gmp_int& a, const gmp_int& b) { mpz_addmul(t.data(), a.data(), b.data()); } inline void eval_multiply_subtract(gmp_int& t, const gmp_int& a, const gmp_int& b) { mpz_submul(t.data(), a.data(), b.data()); } inline void eval_subtract(gmp_int& t, const gmp_int& o) { mpz_sub(t.data(), t.data(), o.data()); } inline void eval_multiply(gmp_int& t, const gmp_int& o) { mpz_mul(t.data(), t.data(), o.data()); } inline void eval_divide(gmp_int& t, const gmp_int& o) { if (eval_is_zero(o)) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); mpz_tdiv_q(t.data(), t.data(), o.data()); } inline void eval_modulus(gmp_int& t, const gmp_int& o) { mpz_tdiv_r(t.data(), t.data(), o.data()); } inline void eval_add(gmp_int& t, unsigned long i) { mpz_add_ui(t.data(), t.data(), i); } inline void eval_multiply_add(gmp_int& t, const gmp_int& a, unsigned long i) { mpz_addmul_ui(t.data(), a.data(), i); } inline void eval_multiply_subtract(gmp_int& t, const gmp_int& a, unsigned long i) { mpz_submul_ui(t.data(), a.data(), i); } inline void eval_subtract(gmp_int& t, unsigned long i) { mpz_sub_ui(t.data(), t.data(), i); } inline void eval_multiply(gmp_int& t, unsigned long i) { mpz_mul_ui(t.data(), t.data(), i); } inline void eval_modulus(gmp_int& t, unsigned long i) { mpz_tdiv_r_ui(t.data(), t.data(), i); } inline void eval_divide(gmp_int& t, unsigned long i) { if (i == 0) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); mpz_tdiv_q_ui(t.data(), t.data(), i); } inline void eval_add(gmp_int& t, long i) { using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; if (i > 0) mpz_add_ui(t.data(), t.data(), static_cast<local_uint_type>(i)); else if (i < 0) mpz_sub_ui(t.data(), t.data(), static_cast<local_uint_type>(-i)); } inline void eval_multiply_add(gmp_int& t, const gmp_int& a, long i) { using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; if (i > 0) mpz_addmul_ui(t.data(), a.data(), static_cast<local_uint_type>(i)); else mpz_submul_ui(t.data(), a.data(), static_cast<local_uint_type>(-i)); } inline void eval_multiply_subtract(gmp_int& t, const gmp_int& a, long i) { using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; if (i > 0) mpz_submul_ui(t.data(), a.data(), static_cast<local_uint_type>(i)); else mpz_addmul_ui(t.data(), a.data(), static_cast<local_uint_type>(-i)); } inline void eval_subtract(gmp_int& t, long i) { using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; if (i > 0) mpz_sub_ui(t.data(), t.data(), static_cast<local_uint_type>(i)); else if (i < 0) mpz_add_ui(t.data(), t.data(), static_cast<local_uint_type>(-i)); } inline void eval_multiply(gmp_int& t, long i) { using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; mpz_mul_ui(t.data(), t.data(), static_cast<local_uint_type>(boost::multiprecision::detail::unsigned_abs(i))); if (i < 0) mpz_neg(t.data(), t.data()); } inline void eval_modulus(gmp_int& t, long i) { using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; mpz_tdiv_r_ui(t.data(), t.data(), static_cast<local_uint_type>(boost::multiprecision::detail::unsigned_abs(i))); } inline void eval_divide(gmp_int& t, long i) { if (i == 0) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; mpz_tdiv_q_ui(t.data(), t.data(), static_cast<local_uint_type>(boost::multiprecision::detail::unsigned_abs(i))); if (i < 0) mpz_neg(t.data(), t.data()); } template <class UI> inline void eval_left_shift(gmp_int& t, UI i) { mpz_mul_2exp(t.data(), t.data(), static_cast<unsigned long>(i)); } template <class UI> inline void eval_right_shift(gmp_int& t, UI i) { mpz_fdiv_q_2exp(t.data(), t.data(), static_cast<unsigned long>(i)); } template <class UI> inline void eval_left_shift(gmp_int& t, const gmp_int& v, UI i) { mpz_mul_2exp(t.data(), v.data(), static_cast<unsigned long>(i)); } template <class UI> inline void eval_right_shift(gmp_int& t, const gmp_int& v, UI i) { mpz_fdiv_q_2exp(t.data(), v.data(), static_cast<unsigned long>(i)); } inline void eval_bitwise_and(gmp_int& result, const gmp_int& v) { mpz_and(result.data(), result.data(), v.data()); } inline void eval_bitwise_or(gmp_int& result, const gmp_int& v) { mpz_ior(result.data(), result.data(), v.data()); } inline void eval_bitwise_xor(gmp_int& result, const gmp_int& v) { mpz_xor(result.data(), result.data(), v.data()); } inline void eval_add(gmp_int& t, const gmp_int& p, const gmp_int& o) { mpz_add(t.data(), p.data(), o.data()); } inline void eval_subtract(gmp_int& t, const gmp_int& p, const gmp_int& o) { mpz_sub(t.data(), p.data(), o.data()); } inline void eval_multiply(gmp_int& t, const gmp_int& p, const gmp_int& o) { mpz_mul(t.data(), p.data(), o.data()); } inline void eval_divide(gmp_int& t, const gmp_int& p, const gmp_int& o) { if (eval_is_zero(o)) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); mpz_tdiv_q(t.data(), p.data(), o.data()); } inline void eval_modulus(gmp_int& t, const gmp_int& p, const gmp_int& o) { mpz_tdiv_r(t.data(), p.data(), o.data()); } inline void eval_add(gmp_int& t, const gmp_int& p, unsigned long i) { mpz_add_ui(t.data(), p.data(), i); } inline void eval_subtract(gmp_int& t, const gmp_int& p, unsigned long i) { mpz_sub_ui(t.data(), p.data(), i); } inline void eval_multiply(gmp_int& t, const gmp_int& p, unsigned long i) { mpz_mul_ui(t.data(), p.data(), i); } inline void eval_modulus(gmp_int& t, const gmp_int& p, unsigned long i) { mpz_tdiv_r_ui(t.data(), p.data(), i); } inline void eval_divide(gmp_int& t, const gmp_int& p, unsigned long i) { if (i == 0) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); mpz_tdiv_q_ui(t.data(), p.data(), i); } inline void eval_add(gmp_int& t, const gmp_int& p, long i) { using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; if (i > 0) mpz_add_ui(t.data(), p.data(), static_cast<local_uint_type>(i)); else mpz_sub_ui(t.data(), p.data(), static_cast<local_uint_type>(-i)); } inline void eval_subtract(gmp_int& t, const gmp_int& p, long i) { using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; if (i > 0) mpz_sub_ui(t.data(), p.data(), static_cast<local_uint_type>(i)); else mpz_add_ui(t.data(), p.data(), static_cast<local_uint_type>(-i)); } inline void eval_multiply(gmp_int& t, const gmp_int& p, long i) { mpz_mul_ui(t.data(), p.data(), boost::multiprecision::detail::unsigned_abs(i)); if (i < 0) mpz_neg(t.data(), t.data()); } inline void eval_modulus(gmp_int& t, const gmp_int& p, long i) { mpz_tdiv_r_ui(t.data(), p.data(), boost::multiprecision::detail::unsigned_abs(i)); } inline void eval_divide(gmp_int& t, const gmp_int& p, long i) { if (i == 0) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); mpz_tdiv_q_ui(t.data(), p.data(), boost::multiprecision::detail::unsigned_abs(i)); if (i < 0) mpz_neg(t.data(), t.data()); } inline void eval_bitwise_and(gmp_int& result, const gmp_int& u, const gmp_int& v) { mpz_and(result.data(), u.data(), v.data()); } inline void eval_bitwise_or(gmp_int& result, const gmp_int& u, const gmp_int& v) { mpz_ior(result.data(), u.data(), v.data()); } inline void eval_bitwise_xor(gmp_int& result, const gmp_int& u, const gmp_int& v) { mpz_xor(result.data(), u.data(), v.data()); } inline void eval_complement(gmp_int& result, const gmp_int& u) { mpz_com(result.data(), u.data()); } inline int eval_get_sign(const gmp_int& val) { return mpz_sgn(val.data()); } inline void eval_convert_to(unsigned long* result, const gmp_int& val) { if (mpz_sgn(val.data()) < 0) { BOOST_MP_THROW_EXCEPTION(std::range_error("Conversion from negative integer to an unsigned type results in undefined behaviour")); } else *result = static_cast<unsigned long>(mpz_get_ui(val.data())); } inline void eval_convert_to(long* result, const gmp_int& val) { if (0 == mpz_fits_slong_p(val.data())) { *result = mpz_sgn(val.data()) < 0 ? (std::numeric_limits<long>::min)() : (std::numeric_limits<long>::max)(); } else *result = static_cast<long>(mpz_get_si(val.data())); } inline void eval_convert_to(long double* result, const gmp_int& val) { detail::gmp_char_ptr val_char_ptr {mpz_get_str(nullptr, 10, val.data())}; *result = std::strtold(val_char_ptr.get(), nullptr); } inline void eval_convert_to(double* result, const gmp_int& val) { *result = mpz_get_d(val.data()); } #ifdef BOOST_HAS_LONG_LONG inline void eval_convert_to(unsigned long long* result, const gmp_int& val) { if (mpz_sgn(val.data()) < 0) { BOOST_MP_THROW_EXCEPTION(std::range_error("Conversion from negative integer to an unsigned type results in undefined behaviour")); } *result = 0; gmp_int t(val); unsigned parts = sizeof(unsigned long long) / sizeof(unsigned long); for (unsigned i = 0; i < parts; ++i) { unsigned long long part = mpz_get_ui(t.data()); if (i) *result |= part << (i * sizeof(unsigned long) * CHAR_BIT); else *result = part; mpz_tdiv_q_2exp(t.data(), t.data(), sizeof(unsigned long) * CHAR_BIT); } } inline void eval_convert_to(long long* result, const gmp_int& val) { int s = mpz_sgn(val.data()); *result = 0; gmp_int t(val); unsigned parts = sizeof(unsigned long long) / sizeof(unsigned long); unsigned long long unsigned_result = 0; for (unsigned i = 0; i < parts; ++i) { unsigned long long part = mpz_get_ui(t.data()); if (i) unsigned_result |= part << (i * sizeof(unsigned long) * CHAR_BIT); else unsigned_result = part; mpz_tdiv_q_2exp(t.data(), t.data(), sizeof(unsigned long) * CHAR_BIT); } // // Overflow check: // bool overflow = false; if (mpz_sgn(t.data())) { overflow = true; } if ((s > 0) && (unsigned_result > static_cast<unsigned long long>((std::numeric_limits<long long>::max)()))) overflow = true; if((s < 0) && (unsigned_result > 1u - static_cast<unsigned long long>((std::numeric_limits<long long>::min)() + 1))) overflow = true; if(overflow) *result = s < 0 ? (std::numeric_limits<long long>::min)() : (std::numeric_limits<long long>::max)(); else *result = s < 0 ? -static_cast<long long>(unsigned_result - 1u) - 1 : static_cast<long long>(unsigned_result); } #endif #ifdef BOOST_HAS_INT128 inline void eval_convert_to(uint128_type* result, const gmp_int& val) { if (mpz_sgn(val.data()) < 0) { BOOST_MP_THROW_EXCEPTION(std::range_error("Conversion from negative integer to an unsigned type results in undefined behaviour")); } *result = 0; gmp_int t(val); unsigned parts = sizeof(uint128_type) / sizeof(unsigned long); for (unsigned i = 0; i < parts; ++i) { uint128_type part = mpz_get_ui(t.data()); if (i) *result |= part << (i * sizeof(unsigned long) * CHAR_BIT); else *result = part; mpz_tdiv_q_2exp(t.data(), t.data(), sizeof(unsigned long) * CHAR_BIT); } } inline void eval_convert_to(int128_type* result, const gmp_int& val) { int s = mpz_sgn(val.data()); *result = 0; gmp_int t(val); unsigned parts = sizeof(uint128_type) / sizeof(unsigned long); uint128_type unsigned_result = 0; for (unsigned i = 0; i < parts; ++i) { uint128_type part = mpz_get_ui(t.data()); if (i) unsigned_result |= part << (i * sizeof(unsigned long) * CHAR_BIT); else unsigned_result = part; mpz_tdiv_q_2exp(t.data(), t.data(), sizeof(unsigned long) * CHAR_BIT); } // // Overflow check: // constexpr int128_type int128_max = static_cast<int128_type>((static_cast<uint128_type>(1u) << 127) - 1); constexpr int128_type int128_min = static_cast<int128_type>(static_cast<int128_type>(-int128_max) -1); bool overflow = false; if (mpz_sgn(t.data())) { overflow = true; } if ((s > 0) && (unsigned_result > static_cast<uint128_type>(int128_max))) overflow = true; if ((s < 0) && (unsigned_result > 1u - static_cast<uint128_type>(int128_min + 1))) overflow = true; if (overflow) *result = s < 0 ? int128_min : int128_max; else *result = s < 0 ? -static_cast<int128_type>(unsigned_result - 1u) - 1 : static_cast<int128_type>(unsigned_result); } template <unsigned digits10> inline void eval_convert_to(int128_type* result, const gmp_float<digits10>& val) { gmp_int i; mpz_set_f(i.data(), val.data()); eval_convert_to(result, i); } template <unsigned digits10> inline void eval_convert_to(uint128_type* result, const gmp_float<digits10>& val) { gmp_int i; mpz_set_f(i.data(), val.data()); eval_convert_to(result, i); } #endif #ifdef BOOST_HAS_FLOAT128 inline void eval_convert_to(float128_type* result, const gmp_int& val) { *result = float128_procs::strtoflt128(val.str(0, std::ios_base::fixed).c_str(), nullptr); } #endif inline void eval_abs(gmp_int& result, const gmp_int& val) { mpz_abs(result.data(), val.data()); } inline void eval_gcd(gmp_int& result, const gmp_int& a, const gmp_int& b) { mpz_gcd(result.data(), a.data(), b.data()); } inline void eval_lcm(gmp_int& result, const gmp_int& a, const gmp_int& b) { mpz_lcm(result.data(), a.data(), b.data()); } template <class I> inline typename std::enable_if<(boost::multiprecision::detail::is_unsigned<I>::value && (sizeof(I) <= sizeof(unsigned long)))>::type eval_gcd(gmp_int& result, const gmp_int& a, const I b) { mpz_gcd_ui(result.data(), a.data(), b); } template <class I> inline typename std::enable_if<(boost::multiprecision::detail::is_unsigned<I>::value && (sizeof(I) <= sizeof(unsigned long)))>::type eval_lcm(gmp_int& result, const gmp_int& a, const I b) { mpz_lcm_ui(result.data(), a.data(), b); } template <class I> inline typename std::enable_if<(boost::multiprecision::detail::is_signed<I>::value && boost::multiprecision::detail::is_integral<I>::value && (sizeof(I) <= sizeof(long)))>::type eval_gcd(gmp_int& result, const gmp_int& a, const I b) { mpz_gcd_ui(result.data(), a.data(), boost::multiprecision::detail::unsigned_abs(b)); } template <class I> inline typename std::enable_if<boost::multiprecision::detail::is_signed<I>::value && boost::multiprecision::detail::is_integral<I>::value && ((sizeof(I) <= sizeof(long)))>::type eval_lcm(gmp_int& result, const gmp_int& a, const I b) { mpz_lcm_ui(result.data(), a.data(), boost::multiprecision::detail::unsigned_abs(b)); } inline void eval_integer_sqrt(gmp_int& s, gmp_int& r, const gmp_int& x) { mpz_sqrtrem(s.data(), r.data(), x.data()); } inline std::size_t eval_lsb(const gmp_int& val) { int c = eval_get_sign(val); if (c == 0) { BOOST_MP_THROW_EXCEPTION(std::domain_error("No bits were set in the operand.")); } if (c < 0) { BOOST_MP_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined.")); } return static_cast<unsigned>(mpz_scan1(val.data(), 0)); } inline std::size_t eval_msb(const gmp_int& val) { int c = eval_get_sign(val); if (c == 0) { BOOST_MP_THROW_EXCEPTION(std::domain_error("No bits were set in the operand.")); } if (c < 0) { BOOST_MP_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined.")); } return static_cast<unsigned>(mpz_sizeinbase(val.data(), 2) - 1); } inline bool eval_bit_test(const gmp_int& val, std::size_t index) { return mpz_tstbit(val.data(), index) ? true : false; } inline void eval_bit_set(gmp_int& val, std::size_t index) { mpz_setbit(val.data(), index); } inline void eval_bit_unset(gmp_int& val, std::size_t index) { mpz_clrbit(val.data(), index); } inline void eval_bit_flip(gmp_int& val, std::size_t index) { mpz_combit(val.data(), index); } inline void eval_qr(const gmp_int& x, const gmp_int& y, gmp_int& q, gmp_int& r) { mpz_tdiv_qr(q.data(), r.data(), x.data(), y.data()); } template <class Integer> inline typename std::enable_if<boost::multiprecision::detail::is_unsigned<Integer>::value, Integer>::type eval_integer_modulus(const gmp_int& x, Integer val) { #if defined(__MPIR_VERSION) && (__MPIR_VERSION >= 3) if ((sizeof(Integer) <= sizeof(mpir_ui)) || (val <= (std::numeric_limits<mpir_ui>::max)())) #else if ((sizeof(Integer) <= sizeof(long)) || (val <= (std::numeric_limits<unsigned long>::max)())) #endif { return static_cast<Integer>(mpz_tdiv_ui(x.data(), val)); } else { return default_ops::eval_integer_modulus(x, val); } } template <class Integer> inline typename std::enable_if<boost::multiprecision::detail::is_signed<Integer>::value && boost::multiprecision::detail::is_integral<Integer>::value, Integer>::type eval_integer_modulus(const gmp_int& x, Integer val) { return eval_integer_modulus(x, boost::multiprecision::detail::unsigned_abs(val)); } inline void eval_powm(gmp_int& result, const gmp_int& base, const gmp_int& p, const gmp_int& m) { if (eval_get_sign(p) < 0) { BOOST_MP_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent.")); } mpz_powm(result.data(), base.data(), p.data(), m.data()); } template <class Integer> inline typename std::enable_if< boost::multiprecision::detail::is_unsigned<Integer>::value && (sizeof(Integer) <= sizeof(unsigned long))>::type eval_powm(gmp_int& result, const gmp_int& base, Integer p, const gmp_int& m) { mpz_powm_ui(result.data(), base.data(), p, m.data()); } template <class Integer> inline typename std::enable_if<boost::multiprecision::detail::is_signed<Integer>::value && boost::multiprecision::detail::is_integral<Integer>::value && (sizeof(Integer) <= sizeof(unsigned long))>::type eval_powm(gmp_int& result, const gmp_int& base, Integer p, const gmp_int& m) { if (p < 0) { BOOST_MP_THROW_EXCEPTION(std::runtime_error("powm requires a positive exponent.")); } mpz_powm_ui(result.data(), base.data(), p, m.data()); } inline std::size_t hash_value(const gmp_int& val) { // We should really use mpz_limbs_read here, but that's unsupported on older versions: std::size_t result = 0; for (int i = 0; i < std::abs(val.data()[0]._mp_size); ++i) boost::multiprecision::detail::hash_combine(result, val.data()[0]._mp_d[i]); boost::multiprecision::detail::hash_combine(result, val.data()[0]._mp_size); return result; } struct gmp_rational; void eval_add(gmp_rational& t, const gmp_rational& o); struct gmp_rational { #ifdef BOOST_HAS_LONG_LONG using signed_types = std::tuple<long, long long> ; using unsigned_types = std::tuple<unsigned long, unsigned long long>; #else using signed_types = std::tuple<long> ; using unsigned_types = std::tuple<unsigned long>; #endif using float_types = std::tuple<double, long double>; gmp_rational() { mpq_init(this->m_data); } gmp_rational(const gmp_rational& o) { mpq_init(m_data); if (o.m_data[0]._mp_num._mp_d) mpq_set(m_data, o.m_data); } gmp_rational(const gmp_int& o) { mpz_init_set(&m_data[0]._mp_num, o.data()); mpz_init_set_ui(&m_data[0]._mp_den, 1u); } gmp_rational(long i) { mpz_init_set_si(&m_data[0]._mp_num, i); mpz_init_set_ui(&m_data[0]._mp_den, 1u); } gmp_rational(unsigned long ui) { mpz_init_set_ui(&m_data[0]._mp_num, ui); mpz_init_set_ui(&m_data[0]._mp_den, 1u); } // 2-arg constructors: template <class T, class U> gmp_rational(const T& a, const U& b, typename std::enable_if<std::is_constructible<gmp_int, T>::value && std::is_constructible<gmp_int, U>::value>::type* = nullptr) { gmp_int i(a), j(b); if (eval_is_zero(j)) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); m_data[0]._mp_num = i.data()[0]; m_data[0]._mp_den = j.data()[0]; mpq_canonicalize(m_data); i.data()[0]._mp_d = nullptr; j.data()[0]._mp_d = nullptr; } template <class U> gmp_rational(const gmp_int& a, const U& b, typename std::enable_if<std::is_constructible<gmp_int, U>::value>::type* = nullptr) { gmp_int j(b); if (eval_is_zero(j)) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); mpz_init_set(&m_data[0]._mp_num, a.data()); m_data[0]._mp_den = j.data()[0]; if (boost::multiprecision::detail::unsigned_abs(b) > 1) mpq_canonicalize(m_data); j.data()[0]._mp_d = nullptr; } template <class U> gmp_rational(gmp_int&& a, const U& b, typename std::enable_if<std::is_constructible<gmp_int, U>::value>::type* = nullptr) { gmp_int j(b); if (eval_is_zero(j)) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); m_data[0]._mp_num = a.data()[0]; m_data[0]._mp_den = j.data()[0]; if (boost::multiprecision::detail::unsigned_abs(b) > 1) mpq_canonicalize(m_data); a.data()[0]._mp_d = nullptr; j.data()[0]._mp_d = nullptr; } template <class T> gmp_rational(const T& a, const gmp_int& b, typename std::enable_if<std::is_constructible<gmp_int, T>::value>::type* = nullptr) { if (eval_is_zero(b)) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); gmp_int i(a); m_data[0]._mp_num = i.data()[0]; mpz_init_set(&m_data[0]._mp_den, b.data()); if(boost::multiprecision::detail::unsigned_abs(a) > 1) mpq_canonicalize(m_data); i.data()[0]._mp_d = nullptr; } template <class T> gmp_rational(const T& a, gmp_int&& b, typename std::enable_if<std::is_constructible<gmp_int, T>::value>::type* = nullptr) { if (eval_is_zero(static_cast<gmp_int&&>(b))) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); gmp_int i(a); m_data[0]._mp_num = i.data()[0]; m_data[0]._mp_den = b.data()[0]; if(boost::multiprecision::detail::unsigned_abs(a) > 1) mpq_canonicalize(m_data); i.data()[0]._mp_d = nullptr; b.data()[0]._mp_d = nullptr; } gmp_rational(const gmp_int& a, const gmp_int& b) { if (eval_is_zero(b)) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); mpz_init_set(&m_data[0]._mp_num, a.data()); mpz_init_set(&m_data[0]._mp_den, b.data()); mpq_canonicalize(m_data); } gmp_rational(const gmp_int& a, gmp_int&& b) { if (eval_is_zero(static_cast<gmp_int&&>(b))) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); mpz_init_set(&m_data[0]._mp_num, a.data()); m_data[0]._mp_den = b.data()[0]; mpq_canonicalize(m_data); b.data()[0]._mp_d = nullptr; } gmp_rational(gmp_int&& a, const gmp_int& b) { if (eval_is_zero(b)) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); m_data[0]._mp_num = a.data()[0]; mpz_init_set(&m_data[0]._mp_den, b.data()); mpq_canonicalize(m_data); a.data()[0]._mp_d = nullptr; } gmp_rational(gmp_int&& a, gmp_int&& b) { if (eval_is_zero(static_cast<gmp_int&&>(b))) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); m_data[0]._mp_num = a.data()[0]; m_data[0]._mp_den = b.data()[0]; mpq_canonicalize(m_data); a.data()[0]._mp_d = nullptr; b.data()[0]._mp_d = nullptr; } // rvalue copy gmp_rational(gmp_rational&& o) noexcept { m_data[0] = o.m_data[0]; o.m_data[0]._mp_num._mp_d = nullptr; o.m_data[0]._mp_den._mp_d = nullptr; } gmp_rational(const mpq_t o) { mpq_init(m_data); mpq_set(m_data, o); } gmp_rational(const mpz_t o) { mpq_init(m_data); mpq_set_z(m_data, o); } gmp_rational& operator=(const gmp_rational& o) { if (m_data[0]._mp_den._mp_d == nullptr) mpq_init(m_data); mpq_set(m_data, o.m_data); return *this; } // rvalue assign gmp_rational& operator=(gmp_rational&& o) noexcept { mpq_swap(m_data, o.m_data); return *this; } #ifdef BOOST_HAS_LONG_LONG #if defined(ULLONG_MAX) && (ULLONG_MAX == ULONG_MAX) gmp_rational& operator=(unsigned long long i) { *this = static_cast<unsigned long>(i); return *this; } #else gmp_rational& operator=(unsigned long long i) { if (m_data[0]._mp_den._mp_d == nullptr) mpq_init(m_data); gmp_int zi; zi = i; mpq_set_z(m_data, zi.data()); return *this; } gmp_rational& operator=(long long i) { if (m_data[0]._mp_den._mp_d == nullptr) mpq_init(m_data); bool neg = i < 0; *this = boost::multiprecision::detail::unsigned_abs(i); if (neg) mpq_neg(m_data, m_data); return *this; } #endif #endif gmp_rational& operator=(unsigned long i) { if (m_data[0]._mp_den._mp_d == nullptr) mpq_init(m_data); mpq_set_ui(m_data, i, 1); return *this; } gmp_rational& operator=(long i) { if (m_data[0]._mp_den._mp_d == nullptr) mpq_init(m_data); mpq_set_si(m_data, i, 1); return *this; } gmp_rational& operator=(double d) { if (m_data[0]._mp_den._mp_d == nullptr) mpq_init(m_data); mpq_set_d(m_data, d); return *this; } template <class F> gmp_rational& assign_float(F a) { using default_ops::eval_add; using default_ops::eval_subtract; BOOST_MP_FLOAT128_USING using std::floor; using std::frexp; using std::ldexp; if (m_data[0]._mp_den._mp_d == nullptr) mpq_init(m_data); if (a == 0) { mpq_set_si(m_data, 0, 1); return *this; } if (a == 1) { mpq_set_si(m_data, 1, 1); return *this; } BOOST_MP_ASSERT(!BOOST_MP_ISINF(a)); BOOST_MP_ASSERT(!BOOST_MP_ISNAN(a)); int e; F f, term; mpq_set_ui(m_data, 0, 1); mpq_set_ui(m_data, 0u, 1); gmp_rational t; f = frexp(a, &e); constexpr int shift = std::numeric_limits<int>::digits - 1; while (f != static_cast<F>(0.0f)) { // extract int sized bits from f: f = ldexp(f, shift); term = floor(f); e -= shift; mpq_mul_2exp(m_data, m_data, shift); t = static_cast<long>(term); eval_add(*this, t); f -= term; } if (e > 0) mpq_mul_2exp(m_data, m_data, static_cast<mp_bitcnt_t>(e)); else if (e < 0) mpq_div_2exp(m_data, m_data, static_cast<mp_bitcnt_t>(-e)); return *this; } gmp_rational& operator=(long double a) { return assign_float(a); } #ifdef BOOST_HAS_FLOAT128 gmp_rational& operator=(float128_type a) { return assign_float(a); } #endif #ifdef BOOST_HAS_INT128 gmp_rational& operator=(uint128_type i) { gmp_int gi; gi = i; return *this = gi; } gmp_rational& operator=(int128_type i) { gmp_int gi; gi = i; return *this = gi; } #endif gmp_rational& operator=(const char* s) { if (m_data[0]._mp_den._mp_d == nullptr) mpq_init(m_data); if (0 != mpq_set_str(m_data, s, 10)) BOOST_MP_THROW_EXCEPTION(std::runtime_error(std::string("The string \"") + s + std::string("\"could not be interpreted as a valid rational number."))); return *this; } gmp_rational& operator=(const gmp_int& o) { if (m_data[0]._mp_den._mp_d == nullptr) mpq_init(m_data); mpq_set_z(m_data, o.data()); return *this; } gmp_rational& operator=(const mpq_t o) { if (m_data[0]._mp_den._mp_d == nullptr) mpq_init(m_data); mpq_set(m_data, o); return *this; } gmp_rational& operator=(const mpz_t o) { if (m_data[0]._mp_den._mp_d == nullptr) mpq_init(m_data); mpq_set_z(m_data, o); return *this; } void swap(gmp_rational& o) { mpq_swap(m_data, o.m_data); } std::string str(std::streamsize /*digits*/, std::ios_base::fmtflags /*f*/) const { BOOST_MP_ASSERT(m_data[0]._mp_num._mp_d); // TODO make a better job of this including handling of f!! void* (*alloc_func_ptr)(size_t); void* (*realloc_func_ptr)(void*, size_t, size_t); void (*free_func_ptr)(void*, size_t); const char* ps = mpq_get_str(nullptr, 10, m_data); std::string s = ps; mp_get_memory_functions(&alloc_func_ptr, &realloc_func_ptr, &free_func_ptr); (*free_func_ptr)((void*)ps, std::strlen(ps) + 1); return s; } ~gmp_rational() { if (m_data[0]._mp_num._mp_d || m_data[0]._mp_den._mp_d) mpq_clear(m_data); } void negate() { BOOST_MP_ASSERT(m_data[0]._mp_num._mp_d); mpq_neg(m_data, m_data); } int compare(const gmp_rational& o) const { BOOST_MP_ASSERT(m_data[0]._mp_num._mp_d && o.m_data[0]._mp_num._mp_d); return mpq_cmp(m_data, o.m_data); } template <class V> int compare(V v) const { gmp_rational d; d = v; return compare(d); } int compare(unsigned long v) const { BOOST_MP_ASSERT(m_data[0]._mp_num._mp_d); return mpq_cmp_ui(m_data, v, 1); } int compare(long v) const { BOOST_MP_ASSERT(m_data[0]._mp_num._mp_d); return mpq_cmp_si(m_data, v, 1); } mpq_t& data() { BOOST_MP_ASSERT(m_data[0]._mp_num._mp_d); return m_data; } const mpq_t& data() const { BOOST_MP_ASSERT(m_data[0]._mp_num._mp_d); return m_data; } protected: mpq_t m_data; }; inline bool eval_is_zero(const gmp_rational& val) { return mpq_sgn(val.data()) == 0; } template <class T> inline bool eval_eq(gmp_rational& a, const T& b) { return a.compare(b) == 0; } template <class T> inline bool eval_lt(gmp_rational& a, const T& b) { return a.compare(b) < 0; } template <class T> inline bool eval_gt(gmp_rational& a, const T& b) { return a.compare(b) > 0; } inline void eval_add(gmp_rational& t, const gmp_rational& o) { mpq_add(t.data(), t.data(), o.data()); } inline void eval_subtract(gmp_rational& t, const gmp_rational& o) { mpq_sub(t.data(), t.data(), o.data()); } inline void eval_multiply(gmp_rational& t, const gmp_rational& o) { mpq_mul(t.data(), t.data(), o.data()); } inline void eval_divide(gmp_rational& t, const gmp_rational& o) { if (eval_is_zero(o)) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); mpq_div(t.data(), t.data(), o.data()); } inline void eval_add(gmp_rational& t, const gmp_rational& p, const gmp_rational& o) { mpq_add(t.data(), p.data(), o.data()); } inline void eval_subtract(gmp_rational& t, const gmp_rational& p, const gmp_rational& o) { mpq_sub(t.data(), p.data(), o.data()); } inline void eval_multiply(gmp_rational& t, const gmp_rational& p, const gmp_rational& o) { mpq_mul(t.data(), p.data(), o.data()); } inline void eval_divide(gmp_rational& t, const gmp_rational& p, const gmp_rational& o) { if (eval_is_zero(o)) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); mpq_div(t.data(), p.data(), o.data()); } // // operator with scalars: // inline void eval_add(gmp_rational& result, gmp_rational const& a, gmp_int const& b) { // we allow result and a to be the same object here: if (&a != &result) { mpz_set(mpq_numref(result.data()), mpq_numref(a.data())); mpz_set(mpq_denref(result.data()), mpq_denref(a.data())); } mpz_addmul(mpq_numref(result.data()), mpq_denref(a.data()), b.data()); // no need to normalize, there can be no common divisor as long as a is already normalized. } inline void eval_add(gmp_rational& result, gmp_rational const& a, unsigned long b) { // we allow result and a to be the same object here: if (&a != &result) { mpz_set(mpq_numref(result.data()), mpq_numref(a.data())); mpz_set(mpq_denref(result.data()), mpq_denref(a.data())); } mpz_addmul_ui(mpq_numref(result.data()), mpq_denref(a.data()), b); // no need to normalize, there can be no common divisor as long as a is already normalized. } inline void eval_add(gmp_rational& result, gmp_rational const& a, long b) { // we allow result and a to be the same object here: if (&a != &result) { mpz_set(mpq_numref(result.data()), mpq_numref(a.data())); mpz_set(mpq_denref(result.data()), mpq_denref(a.data())); } using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; if(b > 0) mpz_addmul_ui(mpq_numref(result.data()), mpq_denref(a.data()), static_cast<local_uint_type>(b)); else mpz_submul_ui(mpq_numref(result.data()), mpq_denref(a.data()), static_cast<local_uint_type>(-b)); // no need to normalize, there can be no common divisor as long as a is already normalized. } template <class T> inline typename std::enable_if<boost::multiprecision::detail::is_integral<T>::value>::type eval_add(gmp_rational& result, gmp_rational const& a, const T& b) { gmp_int t; t = b; eval_add(result, a, t); } template <class T> inline typename std::enable_if<boost::multiprecision::detail::is_integral<T>::value>::type eval_add(gmp_rational& result, const T& b, gmp_rational const& a) { eval_add(result, a, b); } template <class T> inline typename std::enable_if<boost::multiprecision::detail::is_integral<T>::value>::type eval_add(gmp_rational& result, const T& b) { eval_add(result, result, b); } inline void eval_subtract(gmp_rational& result, gmp_rational const& a, gmp_int const& b) { // we allow result and a to be the same object here: if (&a != &result) { mpz_set(mpq_numref(result.data()), mpq_numref(a.data())); mpz_set(mpq_denref(result.data()), mpq_denref(a.data())); } mpz_submul(mpq_numref(result.data()), mpq_denref(a.data()), b.data()); // no need to normalize, there can be no common divisor as long as a is already normalized. } inline void eval_subtract(gmp_rational& result, gmp_rational const& a, unsigned long b) { // we allow result and a to be the same object here: if (&a != &result) { mpz_set(mpq_numref(result.data()), mpq_numref(a.data())); mpz_set(mpq_denref(result.data()), mpq_denref(a.data())); } mpz_submul_ui(mpq_numref(result.data()), mpq_denref(a.data()), b); // no need to normalize, there can be no common divisor as long as a is already normalized. } inline void eval_subtract(gmp_rational& result, gmp_rational const& a, long b) { // we allow result and a to be the same object here: if (&a != &result) { mpz_set(mpq_numref(result.data()), mpq_numref(a.data())); mpz_set(mpq_denref(result.data()), mpq_denref(a.data())); } using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; if(b > 0) mpz_submul_ui(mpq_numref(result.data()), mpq_denref(a.data()), static_cast<local_uint_type>(b)); else mpz_addmul_ui(mpq_numref(result.data()), mpq_denref(a.data()), static_cast<local_uint_type>(-b)); // no need to normalize, there can be no common divisor as long as a is already normalized. } template <class T> inline typename std::enable_if<boost::multiprecision::detail::is_integral<T>::value>::type eval_subtract(gmp_rational& result, gmp_rational const& a, const T& b) { gmp_int t; t = b; eval_subtract(result, a, t); } template <class T> inline typename std::enable_if<boost::multiprecision::detail::is_integral<T>::value>::type eval_subtract(gmp_rational& result, const T& b, gmp_rational const& a) { eval_subtract(result, a, b); result.negate(); } template <class T> inline typename std::enable_if<boost::multiprecision::detail::is_integral<T>::value>::type eval_subtract(gmp_rational& result, const T& b) { eval_subtract(result, result, b); } inline void eval_multiply(gmp_rational& result, gmp_rational const& a, gmp_int const& b) { gmp_int g, t; mpz_gcd(g.data(), mpq_denref(a.data()), b.data()); if (!mpz_fits_uint_p(g.data()) || (mpz_get_ui(g.data()) != 1)) { // We get here if the gcd is not unity, this is true if the number is // too large for an unsigned long, or if we get an unsigned long and check against 1. eval_divide(t, b, g); mpz_mul(mpq_numref(result.data()), t.data(), mpq_numref(a.data())); mpz_divexact(mpq_denref(result.data()), mpq_denref(a.data()), g.data()); } else { // gcd is 1. mpz_mul(mpq_numref(result.data()), mpq_numref(a.data()), b.data()); if (&result != &a) mpz_set(mpq_denref(result.data()), mpq_denref(a.data())); } } inline void eval_multiply(gmp_rational& result, gmp_rational const& a, unsigned long b) { if (b == 0) { mpq_set_ui(result.data(), b, 1); return; } if (mpz_sgn(mpq_numref(a.data())) == 0) { result = a; return; } unsigned long g = static_cast<unsigned long>(mpz_gcd_ui(nullptr, mpq_denref(a.data()), b)); if (g != 1) { BOOST_MP_ASSERT(g); b /= g; mpz_mul_ui(mpq_numref(result.data()), mpq_numref(a.data()), b); mpz_divexact_ui(mpq_denref(result.data()), mpq_denref(a.data()), g); } else { mpz_mul_ui(mpq_numref(result.data()), mpq_numref(a.data()), b); if (&result != &a) mpz_set(mpq_denref(result.data()), mpq_denref(a.data())); } } inline void eval_multiply(gmp_rational& result, gmp_rational const& a, long b) { eval_multiply(result, a, boost::multiprecision::detail::unsigned_abs(b)); if (b < 0) result.negate(); } template <class T> inline typename std::enable_if<boost::multiprecision::detail::is_integral<T>::value>::type eval_multiply(gmp_rational& result, gmp_rational const& a, const T& b) { gmp_int t; t = b; eval_multiply(result, a, t); } template <class T> inline typename std::enable_if<boost::multiprecision::detail::is_integral<T>::value>::type eval_multiply(gmp_rational& result, const T& b, gmp_rational const& a) { eval_multiply(result, a, b); } template <class T> inline typename std::enable_if<boost::multiprecision::detail::is_integral<T>::value>::type eval_multiply(gmp_rational& result, const T& b) { eval_multiply(result, result, b); } inline int eval_get_sign(const gmp_rational& val) { return mpq_sgn(val.data()); } template <class R> inline typename std::enable_if<number_category<R>::value == number_kind_floating_point>::type eval_convert_to(R* result, const gmp_rational& backend) { // // The generic conversion is as good as anything we can write here: // // This does not round correctly: // //*result = mpq_get_d(val.data()); // // This does: // ::boost::multiprecision::detail::generic_convert_rational_to_float(*result, backend); } #ifdef BOOST_HAS_FLOAT128 inline void eval_convert_to(float128_type* result, const gmp_rational& val) { using default_ops::eval_convert_to; gmp_int n, d; float128_type fn, fd; mpz_set(n.data(), mpq_numref(val.data())); mpz_set(d.data(), mpq_denref(val.data())); eval_convert_to(&fn, n); eval_convert_to(&fd, d); *result = fn / fd; } #endif template <class R> inline typename std::enable_if<number_category<R>::value == number_kind_integer>::type eval_convert_to(R* result, const gmp_rational& backend) { gmp_int n(mpq_numref(backend.data())); gmp_int d(mpq_denref(backend.data())); using default_ops::eval_divide; eval_divide(n, d); using default_ops::eval_convert_to; eval_convert_to(result, n); } inline void eval_abs(gmp_rational& result, const gmp_rational& val) { mpq_abs(result.data(), val.data()); } inline void assign_components(gmp_rational& result, unsigned long v1, unsigned long v2) { if (v2 == 0u) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); mpq_set_ui(result.data(), v1, v2); mpq_canonicalize(result.data()); } inline void assign_components(gmp_rational& result, long v1, long v2) { if (v2 == 0) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); using local_uint_type = typename boost::multiprecision::detail::make_unsigned<long>::type; if (v2 < 0) mpq_set_si(result.data(), -v1, static_cast<local_uint_type>(-v2)); else mpq_set_si(result.data(), v1, static_cast<local_uint_type>(v2)); mpq_canonicalize(result.data()); } inline void assign_components(gmp_rational& result, gmp_int const& v1, gmp_int const& v2) { if (eval_is_zero(v2)) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); mpz_set(mpq_numref(result.data()), v1.data()); mpz_set(mpq_denref(result.data()), v2.data()); mpq_canonicalize(result.data()); } template <class T, class U> void assign_components(gmp_rational& result, const T& a, const U& b) { gmp_int x, y; x = a; y = b; if (eval_is_zero(y)) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); std::swap(result.data()[0]._mp_num, x.data()[0]); std::swap(result.data()[0]._mp_den, y.data()[0]); mpq_canonicalize(result.data()); } template <class U> void assign_components(gmp_rational& result, const gmp_int& a, const U& b) { gmp_int y; y = b; if (eval_is_zero(y)) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); mpz_set(&result.data()[0]._mp_num, a.data()); std::swap(result.data()[0]._mp_den, y.data()[0]); mpq_canonicalize(result.data()); } template <class T> void assign_components(gmp_rational& result, const T& a, const gmp_int& b) { if (eval_is_zero(b)) BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero.")); gmp_int x; x = a; std::swap(result.data()[0]._mp_num, x.data()[0]); mpz_set(&result.data()[0]._mp_den, b.data()); mpq_canonicalize(result.data()); } inline std::size_t hash_value(const gmp_rational& val) { std::size_t result = 0; for (int i = 0; i < std::abs(val.data()[0]._mp_num._mp_size); ++i) boost::multiprecision::detail::hash_combine(result, val.data()[0]._mp_num._mp_d[i]); for (int i = 0; i < std::abs(val.data()[0]._mp_den._mp_size); ++i) boost::multiprecision::detail::hash_combine(result, val.data()[0]._mp_den._mp_d[i]); boost::multiprecision::detail::hash_combine(result, val.data()[0]._mp_num._mp_size); return result; } // // Some useful helpers: // inline std::size_t used_gmp_int_bits(const gmp_int& val) { return eval_msb(val) - eval_lsb(val) + 1; } inline std::size_t used_gmp_rational_bits(const gmp_rational& val) { unsigned d2_d = static_cast<unsigned>(mpz_sizeinbase(mpq_denref(val.data()), 2) - mpz_scan1(mpq_denref(val.data()), 0)); unsigned d2_n = static_cast<unsigned>(mpz_sizeinbase(mpq_numref(val.data()), 2) - mpz_scan1(mpq_numref(val.data()), 0)); return (std::max)(d2_d, d2_n); } // // Some member functions that are dependent upon previous code go here: // template <unsigned Digits10> template <unsigned D> inline gmp_float<Digits10>::gmp_float(const gmp_float<D>& o, typename std::enable_if<D <= Digits10>::type*) { mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(Digits10 ? Digits10 : (unsigned)this->get_default_precision())); mpf_set(this->m_data, o.data()); } template <unsigned Digits10> template <unsigned D> inline gmp_float<Digits10>::gmp_float(const gmp_float<D>& o, typename std::enable_if< !(D <= Digits10)>::type*) { mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(Digits10 ? Digits10 : (unsigned)this->get_default_precision())); mpf_set(this->m_data, o.data()); } template <unsigned Digits10> inline gmp_float<Digits10>::gmp_float(const gmp_int& o) { mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(Digits10 ? Digits10 : (unsigned)this->get_default_precision())); mpf_set_z(this->data(), o.data()); } template <unsigned Digits10> inline gmp_float<Digits10>::gmp_float(const gmp_rational& o) { mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(Digits10 ? Digits10 : (unsigned)this->get_default_precision())); mpf_set_q(this->data(), o.data()); } template <unsigned Digits10> template <unsigned D> inline gmp_float<Digits10>& gmp_float<Digits10>::operator=(const gmp_float<D>& o) { if (this->m_data[0]._mp_d == nullptr) mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(Digits10 ? Digits10 : (unsigned)this->get_default_precision())); mpf_set(this->m_data, o.data()); return *this; } template <unsigned Digits10> inline gmp_float<Digits10>& gmp_float<Digits10>::operator=(const gmp_int& o) { if (this->m_data[0]._mp_d == nullptr) mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(Digits10 ? Digits10 : (unsigned)this->get_default_precision())); mpf_set_z(this->data(), o.data()); return *this; } template <unsigned Digits10> inline gmp_float<Digits10>& gmp_float<Digits10>::operator=(const gmp_rational& o) { if (this->m_data[0]._mp_d == nullptr) mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(Digits10 ? Digits10 : (unsigned)this->get_default_precision())); mpf_set_q(this->data(), o.data()); return *this; } inline gmp_float<0>::gmp_float(const gmp_int& o) : requested_precision(get_default_precision()) { if (thread_default_variable_precision_options() >= variable_precision_options::preserve_all_precision) { std::size_t d2 = used_gmp_int_bits(o); std::size_t d10 = 1 + multiprecision::detail::digits2_2_10(d2); if (d10 > requested_precision) requested_precision = static_cast<unsigned>(d10); } mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(requested_precision)); mpf_set_z(this->data(), o.data()); } inline gmp_float<0>::gmp_float(const gmp_rational& o) : requested_precision(get_default_precision()) { if (thread_default_variable_precision_options() >= variable_precision_options::preserve_all_precision) { std::size_t d10 = 1 + multiprecision::detail::digits2_2_10(used_gmp_rational_bits(o)); if (d10 > requested_precision) requested_precision = static_cast<unsigned>(d10); } mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(requested_precision)); mpf_set_q(this->data(), o.data()); } inline gmp_float<0>& gmp_float<0>::operator=(const gmp_int& o) { if (this->m_data[0]._mp_d == nullptr) { requested_precision = this->get_default_precision(); if (thread_default_variable_precision_options() >= variable_precision_options::preserve_all_precision) { std::size_t d2 = used_gmp_int_bits(o); std::size_t d10 = 1 + multiprecision::detail::digits2_2_10(d2); if (d10 > requested_precision) requested_precision = static_cast<unsigned>(d10); } mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(requested_precision)); } else if (thread_default_variable_precision_options() >= variable_precision_options::preserve_all_precision) { std::size_t d2 = used_gmp_int_bits(o); std::size_t d10 = 1 + multiprecision::detail::digits2_2_10(d2); if (d10 > requested_precision) this->precision(static_cast<unsigned>(d10)); } mpf_set_z(this->data(), o.data()); return *this; } inline gmp_float<0>& gmp_float<0>::operator=(const gmp_rational& o) { if (this->m_data[0]._mp_d == nullptr) { requested_precision = this->get_default_precision(); if (thread_default_variable_precision_options() >= variable_precision_options::preserve_all_precision) { std::size_t d10 = 1 + multiprecision::detail::digits2_2_10(used_gmp_rational_bits(o)); if (d10 > requested_precision) requested_precision = static_cast<unsigned>(d10); } mpf_init2(this->m_data, multiprecision::detail::digits10_2_2(requested_precision)); } else if (thread_default_variable_precision_options() >= variable_precision_options::preserve_all_precision) { std::size_t d10 = 1 + multiprecision::detail::digits2_2_10(used_gmp_rational_bits(o)); if (d10 > requested_precision) this->precision(static_cast<unsigned>(d10)); } mpf_set_q(this->data(), o.data()); return *this; } inline gmp_int::gmp_int(const gmp_rational& o) { mpz_init(this->m_data); mpz_set_q(this->m_data, o.data()); } inline gmp_int& gmp_int::operator=(const gmp_rational& o) { if (this->m_data[0]._mp_d == nullptr) mpz_init(this->m_data); mpz_set_q(this->m_data, o.data()); return *this; } } //namespace backends template <expression_template_option ExpressionTemplates> struct component_type<number<gmp_rational, ExpressionTemplates> > { using type = number<gmp_int, ExpressionTemplates>; }; template <expression_template_option ET> inline number<gmp_int, ET> numerator(const number<gmp_rational, ET>& val) { number<gmp_int, ET> result; mpz_set(result.backend().data(), (mpq_numref(val.backend().data()))); return result; } template <expression_template_option ET> inline number<gmp_int, ET> denominator(const number<gmp_rational, ET>& val) { number<gmp_int, ET> result; mpz_set(result.backend().data(), (mpq_denref(val.backend().data()))); return result; } namespace detail { template <> struct digits2<number<gmp_float<0>, et_on> > { static long value() { return static_cast<long>(multiprecision::detail::digits10_2_2(gmp_float<0>::thread_default_precision())); } }; template <> struct digits2<number<gmp_float<0>, et_off> > { static long value() { return static_cast<long>(multiprecision::detail::digits10_2_2(gmp_float<0>::thread_default_precision())); } }; template <> struct digits2<number<debug_adaptor<gmp_float<0> >, et_on> > { static long value() { return static_cast<long>(multiprecision::detail::digits10_2_2(gmp_float<0>::thread_default_precision())); } }; template <> struct digits2<number<debug_adaptor<gmp_float<0> >, et_off> > { static long value() { return static_cast<long>(multiprecision::detail::digits10_2_2(gmp_float<0>::thread_default_precision())); } }; template <unsigned Digits10> struct transcendental_reduction_type<boost::multiprecision::backends::gmp_float<Digits10> > { // // The type used for trigonometric reduction needs 3 times the precision of the base type. // This is double the precision of the original type, plus the largest exponent supported. // As a practical measure the largest argument supported is 1/eps, as supporting larger // arguments requires the division of argument by PI/2 to also be done at higher precision, // otherwise the result (an integer) can not be represented exactly. // // See ARGUMENT REDUCTION FOR HUGE ARGUMENTS. K C Ng. // using type = boost::multiprecision::backends::gmp_float<Digits10 * 3>; }; } // namespace detail template <> struct number_category<detail::canonical<mpz_t, gmp_int>::type> : public std::integral_constant<int, number_kind_integer> {}; template <> struct number_category<detail::canonical<mpq_t, gmp_rational>::type> : public std::integral_constant<int, number_kind_rational> {}; template <> struct number_category<detail::canonical<mpf_t, gmp_float<0> >::type> : public std::integral_constant<int, number_kind_floating_point> {}; namespace detail { template <> struct is_variable_precision<backends::gmp_float<0> > : public std::integral_constant<bool, true> {}; } // namespace detail } // namespace multiprecision namespace math { namespace tools { #ifndef BOOST_MP_MATH_AVAILABLE template <typename T> inline int digits(); template <typename T> inline T max_value(); template <typename T> inline T min_value(); #endif // BOOST_MP_MATH_AVAILABLE inline void set_output_precision(const boost::multiprecision::mpf_float& val, std::ostream& os) { const int sz_prec = static_cast<int>(val.precision()); os << std::setprecision(sz_prec); } template <> inline int digits<boost::multiprecision::mpf_float>() #ifdef BOOST_MATH_NOEXCEPT noexcept #endif { return static_cast<int>(multiprecision::detail::digits10_2_2(boost::multiprecision::mpf_float::thread_default_precision())); } template <> inline int digits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, boost::multiprecision::et_off> >() #ifdef BOOST_MATH_NOEXCEPT noexcept #endif { return static_cast<int>(multiprecision::detail::digits10_2_2(boost::multiprecision::mpf_float::thread_default_precision())); } template <> inline boost::multiprecision::mpf_float max_value<boost::multiprecision::mpf_float>() { boost::multiprecision::mpf_float result(0.5); mpf_mul_2exp(result.backend().data(), result.backend().data(), (std::numeric_limits<mp_exp_t>::max)() / 64 + 1); return result; } template <> inline boost::multiprecision::mpf_float min_value<boost::multiprecision::mpf_float>() { boost::multiprecision::mpf_float result(0.5); mpf_div_2exp(result.backend().data(), result.backend().data(), (std::numeric_limits<mp_exp_t>::max)() / 64 + 1); return result; } template <> inline boost::multiprecision::number<boost::multiprecision::gmp_float<0>, boost::multiprecision::et_off> max_value<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, boost::multiprecision::et_off> >() { boost::multiprecision::number<boost::multiprecision::gmp_float<0>, boost::multiprecision::et_off> result(0.5); mpf_mul_2exp(result.backend().data(), result.backend().data(), (std::numeric_limits<mp_exp_t>::max)() / 64 + 1); return result; } template <> inline boost::multiprecision::number<boost::multiprecision::gmp_float<0>, boost::multiprecision::et_off> min_value<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, boost::multiprecision::et_off> >() { boost::multiprecision::number<boost::multiprecision::gmp_float<0>, boost::multiprecision::et_off> result(0.5); mpf_div_2exp(result.backend().data(), result.backend().data(), (std::numeric_limits<mp_exp_t>::max)() / 64 + 1); return result; } template <> inline int digits<boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::mpf_float::backend_type> > >() #ifdef BOOST_MATH_NOEXCEPT noexcept #endif { return static_cast<int>(multiprecision::detail::digits10_2_2(boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::mpf_float::backend_type> >::thread_default_precision())); } template <> inline int digits<boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::gmp_float<0> >, boost::multiprecision::et_off> >() #ifdef BOOST_MATH_NOEXCEPT noexcept #endif { return static_cast<int>(multiprecision::detail::digits10_2_2(boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::mpf_float::backend_type> >::thread_default_precision())); } template <> inline boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::mpf_float::backend_type> > max_value<boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::mpf_float::backend_type> > >() { return max_value<boost::multiprecision::mpf_float>().backend(); } template <> inline boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::mpf_float::backend_type> > min_value<boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::mpf_float::backend_type> > >() { return min_value<boost::multiprecision::mpf_float>().backend(); } template <> inline boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::gmp_float<0> >, boost::multiprecision::et_off> max_value<boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::gmp_float<0> >, boost::multiprecision::et_off> >() { return max_value<boost::multiprecision::mpf_float>().backend(); } template <> inline boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::gmp_float<0> >, boost::multiprecision::et_off> min_value<boost::multiprecision::number<boost::multiprecision::debug_adaptor<boost::multiprecision::gmp_float<0> >, boost::multiprecision::et_off> >() { return min_value<boost::multiprecision::mpf_float>().backend(); } }} // namespace math::tools } // namespace boost namespace std { // // numeric_limits [partial] specializations for the types declared in this header: // template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> class numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> > { using number_type = boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates>; // // min and max values chosen so as to not cause segfaults when calling // mpf_get_str on 64-bit Linux builds. Possibly we could use larger // exponent values elsewhere. // static number_type calc_min() { number_type result(1); mpf_div_2exp(result.backend().data(), result.backend().data(), (std::numeric_limits<mp_exp_t>::max)() / 64 + 1); return result; } static number_type calc_max() { number_type result(1); mpf_mul_2exp(result.backend().data(), result.backend().data(), (std::numeric_limits<mp_exp_t>::max)() / 64 + 1); return result; } static number_type calc_epsilon() { number_type result(1); mpf_div_2exp(result.backend().data(), result.backend().data(), std::numeric_limits<number_type>::digits - 1); return result; } public: static constexpr bool is_specialized = true; static number_type(min)() { // rely on C++11 thread safe initialization of statics: static const number_type value{calc_min()}; return value; } static number_type(max)() { static number_type value{calc_max()}; return value; } static constexpr number_type lowest() { return -(max)(); } static constexpr int digits = static_cast<int>((Digits10 * 1000L) / 301L + ((Digits10 * 1000L) % 301L ? 2 : 1)); static constexpr int digits10 = Digits10; // Have to allow for a possible extra limb inside the gmp data structure: static constexpr int max_digits10 = Digits10 + 3 + ((GMP_LIMB_BITS * 301L) / 1000L); static constexpr bool is_signed = true; static constexpr bool is_integer = false; static constexpr bool is_exact = false; static constexpr int radix = 2; static number_type epsilon() { static const number_type value{calc_epsilon()}; return value; } // What value should this be???? static number_type round_error() { return 1; } static constexpr long min_exponent = LONG_MIN; static constexpr long min_exponent10 = (LONG_MIN / 1000) * 301L; static constexpr long max_exponent = LONG_MAX; static constexpr long max_exponent10 = (LONG_MAX / 1000) * 301L; static constexpr bool has_infinity = false; static constexpr bool has_quiet_NaN = false; static constexpr bool has_signaling_NaN = false; static constexpr float_denorm_style has_denorm = denorm_absent; static constexpr bool has_denorm_loss = false; static constexpr number_type infinity() { return number_type(); } static constexpr number_type quiet_NaN() { return number_type(); } static constexpr number_type signaling_NaN() { return number_type(); } static constexpr number_type denorm_min() { return number_type(); } static constexpr bool is_iec559 = false; static constexpr bool is_bounded = true; static constexpr bool is_modulo = false; static constexpr bool traps = true; static constexpr bool tinyness_before = false; static constexpr float_round_style round_style = round_indeterminate; }; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::digits; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::digits10; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::max_digits10; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::is_signed; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::is_integer; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::is_exact; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::radix; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr long numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::min_exponent; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr long numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::min_exponent10; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr long numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::max_exponent; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr long numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::max_exponent10; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::has_infinity; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::has_quiet_NaN; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::has_signaling_NaN; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr float_denorm_style numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::has_denorm; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::has_denorm_loss; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::is_iec559; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::is_bounded; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::is_modulo; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::traps; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::tinyness_before; template <unsigned Digits10, boost::multiprecision::expression_template_option ExpressionTemplates> constexpr float_round_style numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<Digits10>, ExpressionTemplates> >::round_style; template <boost::multiprecision::expression_template_option ExpressionTemplates> class numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> > { using number_type = boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates>; public: static constexpr bool is_specialized = false; static number_type(min)() { return number_type(); } static number_type(max)() { return number_type(); } static number_type lowest() { return number_type(); } static constexpr int digits = 0; static constexpr int digits10 = 0; static constexpr int max_digits10 = 0; static constexpr bool is_signed = false; static constexpr bool is_integer = false; static constexpr bool is_exact = false; static constexpr int radix = 0; static number_type epsilon() { return number_type(); } static number_type round_error() { return number_type(); } static constexpr int min_exponent = 0; static constexpr int min_exponent10 = 0; static constexpr int max_exponent = 0; static constexpr int max_exponent10 = 0; static constexpr bool has_infinity = false; static constexpr bool has_quiet_NaN = false; static constexpr bool has_signaling_NaN = false; static constexpr float_denorm_style has_denorm = denorm_absent; static constexpr bool has_denorm_loss = false; static number_type infinity() { return number_type(); } static number_type quiet_NaN() { return number_type(); } static number_type signaling_NaN() { return number_type(); } static number_type denorm_min() { return number_type(); } static constexpr bool is_iec559 = false; static constexpr bool is_bounded = false; static constexpr bool is_modulo = false; static constexpr bool traps = false; static constexpr bool tinyness_before = false; static constexpr float_round_style round_style = round_indeterminate; }; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::digits; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::digits10; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::max_digits10; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::is_signed; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::is_integer; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::is_exact; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::radix; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::min_exponent; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::min_exponent10; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::max_exponent; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::max_exponent10; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::has_infinity; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::has_quiet_NaN; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::has_signaling_NaN; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr float_denorm_style numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::has_denorm; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::has_denorm_loss; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::is_iec559; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::is_bounded; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::is_modulo; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::traps; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::tinyness_before; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr float_round_style numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> >::round_style; template <boost::multiprecision::expression_template_option ExpressionTemplates> class numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> > { using number_type = boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates>; public: static constexpr bool is_specialized = true; // // Largest and smallest numbers are bounded only by available memory, set // to zero: // static number_type(min)() { return number_type(); } static number_type(max)() { return number_type(); } static number_type lowest() { return (min)(); } static constexpr int digits = INT_MAX; static constexpr int digits10 = (INT_MAX / 1000) * 301L; static constexpr int max_digits10 = digits10 + 3; static constexpr bool is_signed = true; static constexpr bool is_integer = true; static constexpr bool is_exact = true; static constexpr int radix = 2; static number_type epsilon() { return number_type(); } static number_type round_error() { return number_type(); } static constexpr int min_exponent = 0; static constexpr int min_exponent10 = 0; static constexpr int max_exponent = 0; static constexpr int max_exponent10 = 0; static constexpr bool has_infinity = false; static constexpr bool has_quiet_NaN = false; static constexpr bool has_signaling_NaN = false; static constexpr float_denorm_style has_denorm = denorm_absent; static constexpr bool has_denorm_loss = false; static number_type infinity() { return number_type(); } static number_type quiet_NaN() { return number_type(); } static number_type signaling_NaN() { return number_type(); } static number_type denorm_min() { return number_type(); } static constexpr bool is_iec559 = false; static constexpr bool is_bounded = false; static constexpr bool is_modulo = false; static constexpr bool traps = false; static constexpr bool tinyness_before = false; static constexpr float_round_style round_style = round_toward_zero; }; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::digits; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::digits10; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::max_digits10; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::is_signed; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::is_integer; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::is_exact; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::radix; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::min_exponent; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::min_exponent10; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::max_exponent; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::max_exponent10; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::has_infinity; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::has_quiet_NaN; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::has_signaling_NaN; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr float_denorm_style numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::has_denorm; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::has_denorm_loss; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::is_iec559; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::is_bounded; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::is_modulo; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::traps; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::tinyness_before; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr float_round_style numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_int, ExpressionTemplates> >::round_style; template <boost::multiprecision::expression_template_option ExpressionTemplates> class numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> > { using number_type = boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates>; public: static constexpr bool is_specialized = true; // // Largest and smallest numbers are bounded only by available memory, set // to zero: // static number_type(min)() { return number_type(); } static number_type(max)() { return number_type(); } static number_type lowest() { return (min)(); } // Digits are unbounded, use zero for now: static constexpr int digits = INT_MAX; static constexpr int digits10 = (INT_MAX / 1000) * 301L; static constexpr int max_digits10 = digits10 + 3; static constexpr bool is_signed = true; static constexpr bool is_integer = false; static constexpr bool is_exact = true; static constexpr int radix = 2; static number_type epsilon() { return number_type(); } static number_type round_error() { return number_type(); } static constexpr int min_exponent = 0; static constexpr int min_exponent10 = 0; static constexpr int max_exponent = 0; static constexpr int max_exponent10 = 0; static constexpr bool has_infinity = false; static constexpr bool has_quiet_NaN = false; static constexpr bool has_signaling_NaN = false; static constexpr float_denorm_style has_denorm = denorm_absent; static constexpr bool has_denorm_loss = false; static number_type infinity() { return number_type(); } static number_type quiet_NaN() { return number_type(); } static number_type signaling_NaN() { return number_type(); } static number_type denorm_min() { return number_type(); } static constexpr bool is_iec559 = false; static constexpr bool is_bounded = false; static constexpr bool is_modulo = false; static constexpr bool traps = false; static constexpr bool tinyness_before = false; static constexpr float_round_style round_style = round_toward_zero; }; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::digits; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::digits10; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::max_digits10; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::is_signed; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::is_integer; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::is_exact; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::radix; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::min_exponent; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::min_exponent10; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::max_exponent; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr int numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::max_exponent10; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::has_infinity; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::has_quiet_NaN; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::has_signaling_NaN; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr float_denorm_style numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::has_denorm; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::has_denorm_loss; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::is_iec559; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::is_bounded; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::is_modulo; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::traps; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr bool numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::tinyness_before; template <boost::multiprecision::expression_template_option ExpressionTemplates> constexpr float_round_style numeric_limits<boost::multiprecision::number<boost::multiprecision::gmp_rational, ExpressionTemplates> >::round_style; #ifdef BOOST_MSVC #pragma warning(pop) #endif } // namespace std namespace Eigen { template <class B1, class B2> struct NumTraitsImp; template <boost::multiprecision::expression_template_option ExpressionTemplates> struct NumTraitsImp<boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates>, boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates> > { using self_type = boost::multiprecision::number<boost::multiprecision::gmp_float<0>, ExpressionTemplates>; using Real = typename boost::multiprecision::scalar_result_from_possible_complex<self_type>::type; using NonInteger = self_type; // Not correct but we can't do much better?? using Literal = double; using Nested = self_type; enum { IsComplex = boost::multiprecision::number_category<self_type>::value == boost::multiprecision::number_kind_complex, IsInteger = boost::multiprecision::number_category<self_type>::value == boost::multiprecision::number_kind_integer, ReadCost = 1, AddCost = 4, MulCost = 8, IsSigned = std::numeric_limits<self_type>::is_specialized ? std::numeric_limits<self_type>::is_signed : true, RequireInitialization = 1, }; static Real highest() noexcept { return boost::math::tools::max_value<Real>(); } static Real lowest() noexcept { return boost::math::tools::min_value<Real>(); } static int digits() noexcept { return boost::math::tools::digits<Real>(); } static int digits10() { return Real::thread_default_precision(); } static Real epsilon() { return ldexp(Real(1), 1 - digits()); } static Real dummy_precision() { return 1000 * epsilon(); } static constexpr long min_exponent() noexcept { return LONG_MIN; } static constexpr long max_exponent() noexcept { return LONG_MAX; } static Real infinity() { return Real(); } static Real quiet_NaN() { return Real(); } }; } #endif