EIC code displayed by LXR master/​include/​boost/​multiprecision/​gmp.hpp	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2025-01-18 09:42:31

 ///////////////////////////////////////////////////////////////////////////////
 //  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 (min)(); }
    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 (min)(); }
    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