multiprecision/cpp_int/misc.hpp

0001 ///////////////////////////////////////////////////////////////
0002 //  Copyright 2012-2020 John Maddock.
0003 //  Copyright 2020 Madhur Chauhan.
0004 //  Copyright 2021 Matt Borland.
0005 //  Distributed under the Boost Software License, Version 1.0.
0006 //  (See accompanying file LICENSE_1_0.txt or copy at
0007 //   https://www.boost.org/LICENSE_1_0.txt)
0008 //
0009 // Comparison operators for cpp_int_backend:
0010 //
0011 #ifndef BOOST_MP_CPP_INT_MISC_HPP
0012 #define BOOST_MP_CPP_INT_MISC_HPP
0013
0014 #include <boost/multiprecision/detail/standalone_config.hpp>
0015 #include <boost/multiprecision/detail/number_base.hpp>
0016 #include <boost/multiprecision/cpp_int/cpp_int_config.hpp>
0017 #include <boost/multiprecision/detail/float128_functions.hpp>
0018 #include <boost/multiprecision/detail/assert.hpp>
0019 #include <boost/multiprecision/detail/constexpr.hpp>
0020 #include <boost/multiprecision/detail/bitscan.hpp> // lsb etc
0021 #include <boost/multiprecision/detail/hash.hpp>
0022 #include <boost/multiprecision/detail/no_exceptions_support.hpp>
0023 #include <numeric> // std::gcd
0024 #include <type_traits>
0025 #include <stdexcept>
0026 #include <cmath>
0027
0028 #ifndef BOOST_MP_STANDALONE
0029 #include <boost/integer/common_factor_rt.hpp>
0030 #endif
0031
0032 #ifdef BOOST_MP_MATH_AVAILABLE
0033 #include <boost/math/special_functions/next.hpp>
0034 #endif
0035
0036 #ifdef BOOST_MSVC
0037 #pragma warning(push)
0038 #pragma warning(disable : 4702)
0039 #pragma warning(disable : 4127) // conditional expression is constant
0040 #pragma warning(disable : 4146) // unary minus operator applied to unsigned type, result still unsigned
0041 #endif
0042
0043 // Forward decleration of gcd and lcm functions
0044 namespace boost { namespace multiprecision { namespace detail {
0045
0046 template <typename T>
0047 inline BOOST_CXX14_CONSTEXPR T constexpr_gcd(T a, T b) noexcept;
0048
0049 template <typename T>
0050 inline BOOST_CXX14_CONSTEXPR T constexpr_lcm(T a, T b) noexcept;
0051
0052 }}} // namespace boost::multiprecision::detail
0053
0054 namespace boost { namespace multiprecision { namespace backends {
0055
0056 template <class T, bool has_limits = std::numeric_limits<T>::is_specialized>
0057 struct numeric_limits_workaround : public std::numeric_limits<T>
0058 {
0059 };
0060 template <class R>
0061 struct numeric_limits_workaround<R, false>
0062 {
0063    static constexpr unsigned digits = ~static_cast<R>(0) < 0 ? sizeof(R) * CHAR_BIT - 1 : sizeof(R) * CHAR_BIT;
0064    static constexpr R (min)(){ return (static_cast<R>(-1) < 0) ? static_cast<R>(1) << digits : 0; }
0065    static constexpr R (max)() { return (static_cast<R>(-1) < 0) ? ~(static_cast<R>(1) << digits) : ~static_cast<R>(0); }
0066 };
0067
0068 template <class R, class CppInt>
0069 BOOST_MP_CXX14_CONSTEXPR void check_in_range(const CppInt& val, const std::integral_constant<int, checked>&)
0070 {
0071    using cast_type = typename boost::multiprecision::detail::canonical<R, CppInt>::type;
0072
0073    if (val.sign())
0074    {
0075       BOOST_IF_CONSTEXPR (boost::multiprecision::detail::is_signed<R>::value == false)
0076          BOOST_MP_THROW_EXCEPTION(std::range_error("Attempt to assign a negative value to an unsigned type."));
0077       if (val.compare(static_cast<cast_type>((numeric_limits_workaround<R>::min)())) < 0)
0078          BOOST_MP_THROW_EXCEPTION(std::overflow_error("Could not convert to the target type - -value is out of range."));
0079    }
0080    else
0081    {
0082       if (val.compare(static_cast<cast_type>((numeric_limits_workaround<R>::max)())) > 0)
0083          BOOST_MP_THROW_EXCEPTION(std::overflow_error("Could not convert to the target type - -value is out of range."));
0084    }
0085 }
0086 template <class R, class CppInt>
0087 inline BOOST_MP_CXX14_CONSTEXPR void check_in_range(const CppInt& /*val*/, const std::integral_constant<int, unchecked>&) noexcept {}
0088
0089 inline BOOST_MP_CXX14_CONSTEXPR void check_is_negative(const std::integral_constant<bool, true>&) noexcept {}
0090 inline void                          check_is_negative(const std::integral_constant<bool, false>&)
0091 {
0092    BOOST_MP_THROW_EXCEPTION(std::range_error("Attempt to assign a negative value to an unsigned type."));
0093 }
0094
0095 template <class Integer>
0096 inline BOOST_MP_CXX14_CONSTEXPR Integer negate_integer(Integer i, const std::integral_constant<bool, true>&) noexcept
0097 {
0098    return -i;
0099 }
0100 template <class Integer>
0101 inline BOOST_MP_CXX14_CONSTEXPR Integer negate_integer(Integer i, const std::integral_constant<bool, false>&) noexcept
0102 {
0103    return ~(i - 1);
0104 }
0105
0106 template <class R, std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
0107 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<R>::value && !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, void>::type
0108 eval_convert_to(R* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& backend)
0109 {
0110    using checked_type = std::integral_constant<int, Checked1>;
0111    check_in_range<R>(backend, checked_type());
0112
0113    BOOST_IF_CONSTEXPR(numeric_limits_workaround<R>::digits < cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits)
0114    {
0115       if ((backend.sign() && boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value) && (1 + static_cast<boost::multiprecision::limb_type>((std::numeric_limits<R>::max)()) <= backend.limbs()[0]))
0116       {
0117          *result = (numeric_limits_workaround<R>::min)();
0118          return;
0119       }
0120       else if (boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value && !backend.sign() && static_cast<boost::multiprecision::limb_type>((std::numeric_limits<R>::max)()) <= backend.limbs()[0])
0121       {
0122          *result = (numeric_limits_workaround<R>::max)();
0123          return;
0124       }
0125       else
0126          *result = static_cast<R>(backend.limbs()[0]);
0127    }
0128    else
0129       *result = static_cast<R>(backend.limbs()[0]);
0130
0131    BOOST_IF_CONSTEXPR(numeric_limits_workaround<R>::digits > cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits)
0132    {
0133      std::size_t shift = cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
0134      std::size_t i     = 1u;
0135
0136       while ((i < backend.size()) && (shift < static_cast<unsigned>(numeric_limits_workaround<R>::digits - cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits)))
0137       {
0138          *result += static_cast<R>(backend.limbs()[i]) << shift;
0139          shift += cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
0140          ++i;
0141       }
0142       //
0143       // We have one more limb to extract, but may not need all the bits, so treat this as a special case:
0144       //
0145       if (i < backend.size())
0146       {
0147          const limb_type mask                 = ((numeric_limits_workaround<R>::digits - shift) == cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits) ? ~static_cast<limb_type>(0) : static_cast<limb_type>(static_cast<limb_type>(1u) << (numeric_limits_workaround<R>::digits - shift)) - 1u;
0148          const limb_type limb_at_index_masked = static_cast<limb_type>(backend.limbs()[i] & mask);
0149
0150          *result = static_cast<R>(*result + static_cast<R>(static_cast<R>(limb_at_index_masked) << shift));
0151
0152          if ((backend.limbs()[i] & static_cast<limb_type>(~mask)) || (i + 1 < backend.size()))
0153          {
0154             // Overflow:
0155             if (backend.sign())
0156             {
0157                check_is_negative(boost::multiprecision::detail::is_signed<R>());
0158                *result = (numeric_limits_workaround<R>::min)();
0159             }
0160             else if (boost::multiprecision::detail::is_signed<R>::value)
0161                *result = (numeric_limits_workaround<R>::max)();
0162             return;
0163          }
0164       }
0165    }
0166    else if (backend.size() > 1)
0167    {
0168       // Overflow:
0169       if (backend.sign())
0170       {
0171          check_is_negative(boost::multiprecision::detail::is_signed<R>());
0172          *result = (numeric_limits_workaround<R>::min)();
0173       }
0174       else if (boost::multiprecision::detail::is_signed<R>::value)
0175          *result = (numeric_limits_workaround<R>::max)();
0176       return;
0177    }
0178    if (backend.sign())
0179    {
0180       check_is_negative(std::integral_constant<bool, boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value>());
0181       *result = negate_integer(*result, std::integral_constant<bool, boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value>());
0182    }
0183 }
0184
0185 template <class R, std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
0186 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<std::is_floating_point<R>::value && !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, void>::type
0187 eval_convert_to(R* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& backend) noexcept(boost::multiprecision::detail::is_arithmetic<R>::value && std::numeric_limits<R>::has_infinity)
0188 {
0189    BOOST_MP_FLOAT128_USING using std::ldexp;
0190    if (eval_is_zero(backend))
0191    {
0192       *result = 0.0f;
0193       return;
0194    }
0195
0196 #ifdef BOOST_HAS_FLOAT128
0197    std::ptrdiff_t bits_to_keep = static_cast<std::ptrdiff_t>(std::is_same<R, float128_type>::value ? 113 : std::numeric_limits<R>::digits);
0198 #else
0199    std::ptrdiff_t bits_to_keep = static_cast<std::ptrdiff_t>(std::numeric_limits<R>::digits);
0200 #endif
0201    std::ptrdiff_t bits = static_cast<std::ptrdiff_t>(eval_msb_imp(backend) + 1);
0202
0203    if (bits > bits_to_keep)
0204    {
0205       // Extract the bits we need, and then manually round the result:
0206       *result = 0.0f;
0207       typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::const_limb_pointer p = backend.limbs();
0208       limb_type mask = ~static_cast<limb_type>(0u);
0209       std::size_t index = backend.size() - 1;
0210       std::size_t shift = cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits * index;
0211       while (bits_to_keep > 0)
0212       {
0213          if (bits_to_keep < (std::ptrdiff_t)cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits)
0214          {
0215             if(index != backend.size() - 1)
0216             {
0217                const std::ptrdiff_t left_shift_amount = static_cast<std::ptrdiff_t>(static_cast<std::ptrdiff_t>(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits) - bits_to_keep);
0218
0219                mask <<= left_shift_amount;
0220             }
0221             else
0222             {
0223                std::ptrdiff_t bits_in_first_limb = static_cast<std::ptrdiff_t>(bits % static_cast<std::ptrdiff_t>(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits));
0224                if (bits_in_first_limb == 0)
0225                   bits_in_first_limb = cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
0226                if (bits_in_first_limb > bits_to_keep)
0227                   mask <<= bits_in_first_limb - bits_to_keep;
0228             }
0229          }
0230          *result += ldexp(static_cast<R>(p[index] & mask), static_cast<int>(shift));
0231          shift -= cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
0232
0233          const bool bits_has_non_zero_remainder = (bits % static_cast<std::ptrdiff_t>(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits) != 0);
0234
0235          bits_to_keep -= ((index == backend.size() - 1) && bits_has_non_zero_remainder)
0236             ? bits % static_cast<std::ptrdiff_t>(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits)
0237             :        static_cast<std::ptrdiff_t>(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits);
0238          --index;
0239       }
0240       // Perform rounding:
0241       bits -= 1 + std::numeric_limits<R>::digits;
0242       if (eval_bit_test(backend, static_cast<unsigned>(bits)))
0243       {
0244          if ((eval_lsb_imp(backend) < static_cast<std::size_t>(bits)) || eval_bit_test(backend, static_cast<std::size_t>(bits + 1)))
0245          {
0246             #ifdef BOOST_MP_MATH_AVAILABLE
0247             BOOST_IF_CONSTEXPR(std::numeric_limits<R>::has_infinity)
0248             {
0249                // Must NOT throw:
0250                *result = boost::math::float_next(*result, boost::math::policies::make_policy(boost::math::policies::overflow_error<boost::math::policies::ignore_error>()));
0251             }
0252             else
0253             {
0254                *result = boost::math::float_next(*result);
0255             }
0256             #else
0257             using std::nextafter; BOOST_MP_FLOAT128_USING
0258             *result = nextafter(*result, *result * 2);
0259             #endif
0260          }
0261       }
0262    }
0263    else
0264    {
0265       typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::const_limb_pointer p = backend.limbs();
0266       std::size_t shift = cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
0267       *result = static_cast<R>(*p);
0268       for (std::size_t i = 1; i < backend.size(); ++i)
0269       {
0270          *result += static_cast<R>(ldexp(static_cast<long double>(p[i]), static_cast<int>(shift)));
0271          shift += cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
0272       }
0273    }
0274    if (backend.sign())
0275       *result = -*result;
0276 }
0277
0278 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
0279 BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, bool>::type
0280 eval_is_zero(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) noexcept
0281 {
0282    return (val.size() == 1) && (val.limbs()[0] == 0);
0283 }
0284 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
0285 BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, int>::type
0286 eval_get_sign(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) noexcept
0287 {
0288    return eval_is_zero(val) ? 0 : val.sign() ? -1 : 1;
0289 }
0290 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
0291 BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
0292 eval_abs(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
0293 {
0294    result = val;
0295    result.sign(false);
0296 }
0297
0298 //
0299 // Get the location of the least-significant-bit:
0300 //
0301 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
0302 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, std::size_t>::type
0303 eval_lsb_imp(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
0304 {
0305    //
0306    // Find the index of the least significant limb that is non-zero:
0307    //
0308    std::size_t index = 0;
0309    while (!a.limbs()[index] && (index < a.size()))
0310       ++index;
0311    //
0312    // Find the index of the least significant bit within that limb:
0313    //
0314    std::size_t result = boost::multiprecision::detail::find_lsb(a.limbs()[index]);
0315
0316    return result + index * cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
0317 }
0318
0319 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
0320 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, std::size_t>::type
0321 eval_lsb(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
0322 {
0323    using default_ops::eval_get_sign;
0324    if (eval_get_sign(a) == 0)
0325    {
0326       BOOST_MP_THROW_EXCEPTION(std::domain_error("No bits were set in the operand."));
0327    }
0328    if (a.sign())
0329    {
0330       BOOST_MP_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined."));
0331    }
0332    return eval_lsb_imp(a);
0333 }
0334
0335 //
0336 // Get the location of the most-significant-bit:
0337 //
0338 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
0339 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, std::size_t>::type
0340 eval_msb_imp(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
0341 {
0342    //
0343    // Find the index of the most significant bit that is non-zero:
0344    //
0345    return (a.size() - 1) * cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits + boost::multiprecision::detail::find_msb(a.limbs()[a.size() - 1]);
0346 }
0347
0348 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
0349 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, std::size_t>::type
0350 eval_msb(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
0351 {
0352    using default_ops::eval_get_sign;
0353    if (eval_get_sign(a) == 0)
0354    {
0355       BOOST_MP_THROW_EXCEPTION(std::domain_error("No bits were set in the operand."));
0356    }
0357    if (a.sign())
0358    {
0359       BOOST_MP_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined."));
0360    }
0361    return eval_msb_imp(a);
0362 }
0363
0364 #ifdef BOOST_GCC
0365 //
0366 // We really shouldn't need to be disabling this warning, but it really does appear to be
0367 // spurious.  The warning appears only when in release mode, and asserts are on.
0368 //
0369 #pragma GCC diagnostic push
0370 #pragma GCC diagnostic ignored "-Warray-bounds"
0371 #endif
0372
0373 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
0374 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, bool>::type
0375 eval_bit_test(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val, std::size_t index) noexcept
0376 {
0377    std::size_t  offset = index / cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
0378    std::size_t  shift  = index % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
0379    limb_type mask   = shift ? limb_type(1u) << shift : limb_type(1u);
0380    if (offset >= val.size())
0381       return false;
0382    return val.limbs()[offset] & mask ? true : false;
0383 }
0384
0385 #ifdef BOOST_GCC
0386 #pragma GCC diagnostic pop
0387 #endif
0388
0389 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
0390 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
0391 eval_bit_set(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val, std::size_t index)
0392 {
0393    std::size_t  offset = index / cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
0394    std::size_t  shift  = index % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
0395    limb_type mask   = shift ? limb_type(1u) << shift : limb_type(1u);
0396    if (offset >= val.size())
0397    {
0398       std::size_t os = val.size();
0399       val.resize(offset + 1, offset + 1);
0400       if (offset >= val.size())
0401          return; // fixed precision overflow
0402       for (std::size_t i = os; i <= offset; ++i)
0403          val.limbs()[i] = 0;
0404    }
0405    val.limbs()[offset] |= mask;
0406 }
0407
0408 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
0409 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
0410 eval_bit_unset(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val, std::size_t index) noexcept
0411 {
0412    std::size_t  offset = index / cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
0413    std::size_t  shift  = index % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
0414    limb_type mask   = shift ? limb_type(1u) << shift : limb_type(1u);
0415    if (offset >= val.size())
0416       return;
0417    val.limbs()[offset] &= ~mask;
0418    val.normalize();
0419 }
0420
0421 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
0422 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
0423 eval_bit_flip(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val, std::size_t index)
0424 {
0425    std::size_t  offset = index / cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
0426    std::size_t  shift  = index % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
0427    limb_type mask   = shift ? limb_type(1u) << shift : limb_type(1u);
0428    if (offset >= val.size())
0429    {
0430       std::size_t os = val.size();
0431       val.resize(offset + 1, offset + 1);
0432       if (offset >= val.size())
0433          return; // fixed precision overflow
0434       for (std::size_t i = os; i <= offset; ++i)
0435          val.limbs()[i] = 0;
0436    }
0437    val.limbs()[offset] ^= mask;
0438    val.normalize();
0439 }
0440
0441 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
0442 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
0443 eval_qr(
0444     const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x,
0445     const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& y,
0446     cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&       q,
0447     cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&       r) noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
0448 {
0449    divide_unsigned_helper(&q, x, y, r);
0450    q.sign(x.sign() != y.sign());
0451    r.sign(x.sign());
0452 }
0453
0454 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
0455 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
0456 eval_qr(
0457     const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x,
0458     limb_type                                                                   y,
0459     cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&       q,
0460     cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&       r) noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
0461 {
0462    divide_unsigned_helper(&q, x, y, r);
0463    q.sign(x.sign());
0464    r.sign(x.sign());
0465 }
0466
0467 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class U>
0468 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<U>::value>::type eval_qr(
0469     const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x,
0470     U                                                                           y,
0471     cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&       q,
0472     cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&       r) noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
0473 {
0474    using default_ops::eval_qr;
0475    cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t;
0476    t = y;
0477    eval_qr(x, t, q, r);
0478 }
0479
0480 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Integer>
0481 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_unsigned<Integer>::value && !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, Integer>::type
0482 eval_integer_modulus(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a, Integer mod)
0483 {
0484    BOOST_IF_CONSTEXPR (sizeof(Integer) <= sizeof(limb_type))
0485    {
0486       if (mod <= (std::numeric_limits<limb_type>::max)())
0487       {
0488          const std::ptrdiff_t n = a.size();
0489          const double_limb_type two_n_mod = static_cast<limb_type>(1u) + (~static_cast<limb_type>(0u) - mod) % mod;
0490          limb_type              res = a.limbs()[n - 1] % mod;
0491
0492          for (std::ptrdiff_t i = n - 2; i >= 0; --i)
0493             res = static_cast<limb_type>((res * two_n_mod + a.limbs()[i]) % mod);
0494          return res;
0495       }
0496       else
0497          return default_ops::eval_integer_modulus(a, mod);
0498    }
0499    else
0500    {
0501       return default_ops::eval_integer_modulus(a, mod);
0502    }
0503 }
0504
0505 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Integer>
0506 BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_signed<Integer>::value && boost::multiprecision::detail::is_integral<Integer>::value && !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, Integer>::type
0507 eval_integer_modulus(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x, Integer val)
0508 {
0509    return eval_integer_modulus(x, boost::multiprecision::detail::unsigned_abs(val));
0510 }
0511
0512 BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR limb_type eval_gcd(limb_type u, limb_type v)
0513 {
0514    // boundary cases
0515    if (!u || !v)
0516       return u | v;
0517 #if (defined(__cpp_lib_gcd_lcm) && (__cpp_lib_gcd_lcm >= 201606L))
0518    return std::gcd(u, v);
0519 #else
0520    std::size_t shift = boost::multiprecision::detail::find_lsb(u | v);
0521    u >>= boost::multiprecision::detail::find_lsb(u);
0522    do
0523    {
0524       v >>= boost::multiprecision::detail::find_lsb(v);
0525       if (u > v)
0526          std_constexpr::swap(u, v);
0527       v -= u;
0528    } while (v);
0529    return u << shift;
0530 #endif
0531 }
0532
0533 inline BOOST_MP_CXX14_CONSTEXPR double_limb_type eval_gcd(double_limb_type u, double_limb_type v)
0534 {
0535 #if (defined(__cpp_lib_gcd_lcm) && (__cpp_lib_gcd_lcm >= 201606L)) && (!defined(BOOST_HAS_INT128) || !defined(__STRICT_ANSI__))
0536    return std::gcd(u, v);
0537 #else
0538    if (u == 0)
0539       return v;
0540
0541    std::size_t shift = boost::multiprecision::detail::find_lsb(u | v);
0542    u >>= boost::multiprecision::detail::find_lsb(u);
0543    do
0544    {
0545       v >>= boost::multiprecision::detail::find_lsb(v);
0546       if (u > v)
0547          std_constexpr::swap(u, v);
0548       v -= u;
0549    } while (v);
0550    return u << shift;
0551 #endif
0552 }
0553
0554 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
0555 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
0556 eval_gcd(
0557     cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&       result,
0558     const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
0559     limb_type                                                                   b)
0560 {
0561    int s = eval_get_sign(a);
0562    if (!b || !s)
0563    {
0564       result = a;
0565       *result.limbs() |= b;
0566    }
0567    else
0568    {
0569       eval_modulus(result, a, b);
0570       limb_type& res = *result.limbs();
0571       res = eval_gcd(res, b);
0572    }
0573    result.sign(false);
0574 }
0575
0576 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
0577 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
0578 eval_gcd(
0579     cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&       result,
0580     const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
0581     double_limb_type                                                            b)
0582 {
0583    int s = eval_get_sign(a);
0584    if (!b || !s)
0585    {
0586       if (!s)
0587          result = b;
0588       else
0589          result = a;
0590       return;
0591    }
0592    double_limb_type res = 0;
0593    if(a.sign() == 0)
0594       res = eval_integer_modulus(a, b);
0595    else
0596    {
0597       cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t(a);
0598       t.negate();
0599       res = eval_integer_modulus(t, b);
0600    }
0601    res            = eval_gcd(res, b);
0602    result = res;
0603    result.sign(false);
0604 }
0605 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
0606 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
0607 eval_gcd(
0608    cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
0609    const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
0610    signed_double_limb_type                                                     v)
0611 {
0612    eval_gcd(result, a, static_cast<double_limb_type>(v < 0 ? -v : v));
0613 }
0614 //
0615 // These 2 overloads take care of gcd against an (unsigned) short etc:
0616 //
0617 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Integer>
0618 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_unsigned<Integer>::value && (sizeof(Integer) <= sizeof(limb_type)) && !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
0619 eval_gcd(
0620     cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&       result,
0621     const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
0622     const Integer&                                                              v)
0623 {
0624    eval_gcd(result, a, static_cast<limb_type>(v));
0625 }
0626 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Integer>
0627 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_signed<Integer>::value && boost::multiprecision::detail::is_integral<Integer>::value && (sizeof(Integer) <= sizeof(limb_type)) && !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
0628 eval_gcd(
0629     cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&       result,
0630     const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
0631     const Integer&                                                              v)
0632 {
0633    eval_gcd(result, a, static_cast<limb_type>(v < 0 ? -v : v));
0634 }
0635 //
0636 // What follows is Lehmer's GCD algorithm:
0637 // Essentially this uses the leading digit(s) of U and V
0638 // only to run a "simulated" Euclid algorithm.  It stops
0639 // when the calculated quotient differs from what would have been
0640 // the true quotient.  At that point the cosequences are used to
0641 // calculate the new U and V.  A nice lucid description appears
0642 // in "An Analysis of Lehmer's Euclidean GCD Algorithm",
0643 // by Jonathan Sorenson.  https://www.researchgate.net/publication/2424634_An_Analysis_of_Lehmer%27s_Euclidean_GCD_Algorithm
0644 // DOI: 10.1145/220346.220378.
0645 //
0646 // There are two versions of this algorithm here, and both are "double digit"
0647 // variations: which is to say if there are k bits per limb, then they extract
0648 // 2k bits into a double_limb_type and then run the algorithm on that.  The first
0649 // version is a straightforward version of the algorithm, and is designed for
0650 // situations where double_limb_type is a native integer (for example where
0651 // limb_type is a 32-bit integer on a 64-bit machine).  For 32-bit limbs it
0652 // reduces the size of U by about 30 bits per call.  The second is a more complex
0653 // version for situations where double_limb_type is a synthetic type: for example
0654 // __int128.  For 64 bit limbs it reduces the size of U by about 62 bits per call.
0655 //
0656 // The complexity of the algorithm given by Sorenson is roughly O(ln^2(N)) for
0657 // two N bit numbers.
0658 //
0659 // The original double-digit version of the algorithm is described in:
0660 //
0661 // "A Double Digit Lehmer-Euclid Algorithm for Finding the GCD of Long Integers",
0662 // Tudor Jebelean, J Symbolic Computation, 1995 (19), 145.
0663 //
0664 #ifndef BOOST_HAS_INT128
0665 //
0666 // When double_limb_type is a native integer type then we should just use it and not worry about the consequences.
0667 // This can eliminate approximately a full limb with each call.
0668 //
0669 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Storage>
0670 void eval_gcd_lehmer(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& U, cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& V, std::size_t lu, Storage& storage)
0671 {
0672    //
0673    // Extract the leading 2 * bits_per_limb bits from U and V:
0674    //
0675    std::size_t         h = lu % bits_per_limb;
0676    double_limb_type u = (static_cast<double_limb_type>((U.limbs()[U.size() - 1])) << bits_per_limb) | U.limbs()[U.size() - 2];
0677    double_limb_type v = (static_cast<double_limb_type>((V.size() < U.size() ? 0 : V.limbs()[V.size() - 1])) << bits_per_limb) | V.limbs()[U.size() - 2];
0678    if (h)
0679    {
0680       u <<= bits_per_limb - h;
0681       u |= U.limbs()[U.size() - 3] >> h;
0682       v <<= bits_per_limb - h;
0683       v |= V.limbs()[U.size() - 3] >> h;
0684    }
0685    //
0686    // Co-sequences x an y: we need only the last 3 values of these,
0687    // the first 2 values are known correct, the third gets checked
0688    // in each loop operation, and we terminate when they go wrong.
0689    //
0690    // x[i+0] is positive for even i.
0691    // y[i+0] is positive for odd i.
0692    //
0693    // However we track only absolute values here:
0694    //
0695    double_limb_type x[3] = {1, 0};
0696    double_limb_type y[3] = {0, 1};
0697    std::size_t         i    = 0;
0698
0699 #ifdef BOOST_MP_GCD_DEBUG
0700    cpp_int UU, VV;
0701    UU = U;
0702    VV = V;
0703 #endif
0704
0705    while (true)
0706    {
0707       double_limb_type q  = u / v;
0708       x[2]                = x[0] + q * x[1];
0709       y[2]                = y[0] + q * y[1];
0710       double_limb_type tu = u;
0711       u                   = v;
0712       v                   = tu - q * v;
0713       ++i;
0714       //
0715       // We must make sure that y[2] occupies a single limb otherwise
0716       // the multiprecision multiplications below would be much more expensive.
0717       // This can sometimes lose us one iteration, but is worth it for improved
0718       // calculation efficiency.
0719       //
0720       if (y[2] >> bits_per_limb)
0721          break;
0722       //
0723       // These are Jebelean's exact termination conditions:
0724       //
0725       if ((i & 1u) == 0)
0726       {
0727          BOOST_MP_ASSERT(u > v);
0728          if ((v < x[2]) || ((u - v) < (y[2] + y[1])))
0729             break;
0730       }
0731       else
0732       {
0733          BOOST_MP_ASSERT(u > v);
0734          if ((v < y[2]) || ((u - v) < (x[2] + x[1])))
0735             break;
0736       }
0737 #ifdef BOOST_MP_GCD_DEBUG
0738       BOOST_MP_ASSERT(q == UU / VV);
0739       UU %= VV;
0740       UU.swap(VV);
0741 #endif
0742       x[0] = x[1];
0743       x[1] = x[2];
0744       y[0] = y[1];
0745       y[1] = y[2];
0746    }
0747    if (i == 1)
0748    {
0749       // No change to U and V we've stalled!
0750       cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t;
0751       eval_modulus(t, U, V);
0752       U.swap(V);
0753       V.swap(t);
0754       return;
0755    }
0756    //
0757    // Update U and V.
0758    // We have:
0759    //
0760    // U = x[0]U + y[0]V and
0761    // V = x[1]U + y[1]V.
0762    //
0763    // But since we track only absolute values of x and y
0764    // we have to take account of the implied signs and perform
0765    // the appropriate subtraction depending on the whether i is
0766    // even or odd:
0767    //
0768    std::size_t                                                             ts = U.size() + 1;
0769    cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t1(storage, ts), t2(storage, ts), t3(storage, ts);
0770    eval_multiply(t1, U, static_cast<limb_type>(x[0]));
0771    eval_multiply(t2, V, static_cast<limb_type>(y[0]));
0772    eval_multiply(t3, U, static_cast<limb_type>(x[1]));
0773    if ((i & 1u) == 0)
0774    {
0775       if (x[0] == 0)
0776          U = t2;
0777       else
0778       {
0779          BOOST_MP_ASSERT(t2.compare(t1) >= 0);
0780          eval_subtract(U, t2, t1);
0781          BOOST_MP_ASSERT(U.sign() == false);
0782       }
0783    }
0784    else
0785    {
0786       BOOST_MP_ASSERT(t1.compare(t2) >= 0);
0787       eval_subtract(U, t1, t2);
0788       BOOST_MP_ASSERT(U.sign() == false);
0789    }
0790    eval_multiply(t2, V, static_cast<limb_type>(y[1]));
0791    if (i & 1u)
0792    {
0793       if (x[1] == 0)
0794          V = t2;
0795       else
0796       {
0797          BOOST_MP_ASSERT(t2.compare(t3) >= 0);
0798          eval_subtract(V, t2, t3);
0799          BOOST_MP_ASSERT(V.sign() == false);
0800       }
0801    }
0802    else
0803    {
0804       BOOST_MP_ASSERT(t3.compare(t2) >= 0);
0805       eval_subtract(V, t3, t2);
0806       BOOST_MP_ASSERT(V.sign() == false);
0807    }
0808    BOOST_MP_ASSERT(U.compare(V) >= 0);
0809    BOOST_MP_ASSERT(lu > eval_msb(U));
0810 #ifdef BOOST_MP_GCD_DEBUG
0811
0812    BOOST_MP_ASSERT(UU == U);
0813    BOOST_MP_ASSERT(VV == V);
0814
0815    extern std::size_t total_lehmer_gcd_calls;
0816    extern std::size_t total_lehmer_gcd_bits_saved;
0817    extern std::size_t total_lehmer_gcd_cycles;
0818
0819    ++total_lehmer_gcd_calls;
0820    total_lehmer_gcd_bits_saved += lu - eval_msb(U);
0821    total_lehmer_gcd_cycles += i;
0822 #endif
0823    if (lu < 2048)
0824    {
0825       //
0826       // Since we have stripped all common powers of 2 from U and V at the start
0827       // if either are even at this point, we can remove stray powers of 2 now.
0828       // Note that it is not possible for *both* U and V to be even at this point.
0829       //
0830       // This has an adverse effect on performance for high bit counts, but has
0831       // a significant positive effect for smaller counts.
0832       //
0833       if ((U.limbs()[0] & 1u) == 0)
0834       {
0835          eval_right_shift(U, eval_lsb(U));
0836          if (U.compare(V) < 0)
0837             U.swap(V);
0838       }
0839       else if ((V.limbs()[0] & 1u) == 0)
0840       {
0841          eval_right_shift(V, eval_lsb(V));
0842       }
0843    }
0844    storage.deallocate(ts * 3);
0845 }
0846
0847 #else
0848 //
0849 // This branch is taken when double_limb_type is a synthetic type with no native hardware support.
0850 // For example __int128.  The assumption is that add/subtract/multiply of double_limb_type are efficient,
0851 // but that division is very slow.
0852 //
0853 // We begin with a specialized routine for division.
0854 // We know that most of the time this is called the result will be 1.
0855 // For small limb counts, this almost doubles the performance of Lehmer's routine!
0856 //
0857 BOOST_FORCEINLINE void divide_subtract(double_limb_type& q, double_limb_type& u, const double_limb_type& v)
0858 {
0859    BOOST_MP_ASSERT(q == 1); // precondition on entry.
0860    u -= v;
0861    while (u >= v)
0862    {
0863       u -= v;
0864       if (++q > 30)
0865       {
0866          double_limb_type t = u / v;
0867          u -= t * v;
0868          q += t;
0869       }
0870    }
0871 }
0872
0873 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Storage>
0874 void eval_gcd_lehmer(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& U, cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& V, std::size_t lu, Storage& storage)
0875 {
0876    //
0877    // Extract the leading 2*bits_per_limb bits from U and V:
0878    //
0879    std::size_t  h = lu % bits_per_limb;
0880    double_limb_type u, v;
0881    if (h)
0882    {
0883       u = (static_cast<double_limb_type>((U.limbs()[U.size() - 1])) << bits_per_limb) | U.limbs()[U.size() - 2];
0884       v = (static_cast<double_limb_type>((V.size() < U.size() ? 0 : V.limbs()[V.size() - 1])) << bits_per_limb) | V.limbs()[U.size() - 2];
0885       u <<= bits_per_limb - h;
0886       u |= U.limbs()[U.size() - 3] >> h;
0887       v <<= bits_per_limb - h;
0888       v |= V.limbs()[U.size() - 3] >> h;
0889    }
0890    else
0891    {
0892       u = (static_cast<double_limb_type>(U.limbs()[U.size() - 1]) << bits_per_limb) | U.limbs()[U.size() - 2];
0893       v = (static_cast<double_limb_type>(V.limbs()[U.size() - 1]) << bits_per_limb) | V.limbs()[U.size() - 2];
0894    }
0895    //
0896    // Cosequences are stored as limb_types, we take care not to overflow these:
0897    //
0898    // x[i+0] is positive for even i.
0899    // y[i+0] is positive for odd i.
0900    //
0901    // However we track only absolute values here:
0902    //
0903    limb_type x[3] = { 1, 0 };
0904    limb_type y[3] = { 0, 1 };
0905    std::size_t  i = 0;
0906
0907 #ifdef BOOST_MP_GCD_DEBUG
0908    cpp_int UU, VV;
0909    UU = U;
0910    VV = V;
0911 #endif
0912    //
0913    // We begine by running a single digit version of Lehmer's algorithm, we still have
0914    // to track u and v at double precision, but this adds only a tiny performance penalty.
0915    // What we gain is fast division, and fast termination testing.
0916    // When you see static_cast<limb_type>(u >> bits_per_limb) here, this is really just
0917    // a direct access to the upper bits_per_limb of the double limb type.  For __int128
0918    // this is simple a load of the upper 64 bits and the "shift" is optimised away.
0919    //
0920    double_limb_type old_u, old_v;
0921    while (true)
0922    {
0923       limb_type q = static_cast<limb_type>(u >> bits_per_limb) / static_cast<limb_type>(v >> bits_per_limb);
0924       x[2] = x[0] + q * x[1];
0925       y[2] = y[0] + q * y[1];
0926       double_limb_type tu = u;
0927       old_u = u;
0928       old_v = v;
0929       u = v;
0930       double_limb_type t = q * v;
0931       if (tu < t)
0932       {
0933          ++i;
0934          break;
0935       }
0936       v = tu - t;
0937       ++i;
0938       BOOST_MP_ASSERT((u <= v) || (t / q == old_v));
0939       if (u <= v)
0940       {
0941          // We've gone terribly wrong, probably numeric overflow:
0942          break;
0943       }
0944       if ((i & 1u) == 0)
0945       {
0946          if ((static_cast<limb_type>(v >> bits_per_limb) < x[2]) || ((static_cast<limb_type>(u >> bits_per_limb) - static_cast<limb_type>(v >> bits_per_limb)) < (y[2] + y[1])))
0947             break;
0948       }
0949       else
0950       {
0951          if ((static_cast<limb_type>(v >> bits_per_limb) < y[2]) || ((static_cast<limb_type>(u >> bits_per_limb) - static_cast<limb_type>(v >> bits_per_limb)) < (x[2] + x[1])))
0952             break;
0953       }
0954 #ifdef BOOST_MP_GCD_DEBUG
0955       BOOST_MP_ASSERT(q == UU / VV);
0956       UU %= VV;
0957       UU.swap(VV);
0958 #endif
0959       x[0] = x[1];
0960       x[1] = x[2];
0961       y[0] = y[1];
0962       y[1] = y[2];
0963    }
0964    //
0965    // We get here when the single digit algorithm has gone wrong, back up i, u and v:
0966    //
0967    --i;
0968    u = old_u;
0969    v = old_v;
0970    //
0971    // Now run the full double-digit algorithm:
0972    //
0973    while (true)
0974    {
0975       double_limb_type q = 1u;
0976       double_limb_type tt = u;
0977       divide_subtract(q, u, v);
0978       std::swap(u, v);
0979       tt = y[0] + q * static_cast<double_limb_type>(y[1]);
0980       //
0981       // If calculation of y[2] would overflow a single limb, then we *must* terminate.
0982       // Note that x[2] < y[2] so there is no need to check that as well:
0983       //
0984       if (tt >> bits_per_limb)
0985       {
0986          ++i;
0987          break;
0988       }
0989       x[2] = static_cast<limb_type>(x[0] + static_cast<double_limb_type>(q * x[1]));
0990       y[2] = static_cast<limb_type>(tt);
0991       ++i;
0992       if ((i & 1u) == 0)
0993       {
0994          BOOST_MP_ASSERT(u > v);
0995          if ((v < x[2]) || ((u - v) < (static_cast<double_limb_type>(y[2]) + y[1])))
0996             break;
0997       }
0998       else
0999       {
1000          BOOST_MP_ASSERT(u > v);
1001          if ((v < y[2]) || ((u - v) < (static_cast<double_limb_type>(x[2]) + x[1])))
1002             break;
1003       }
1004 #ifdef BOOST_MP_GCD_DEBUG
1005       BOOST_MP_ASSERT(q == UU / VV);
1006       UU %= VV;
1007       UU.swap(VV);
1008 #endif
1009       x[0] = x[1];
1010       x[1] = x[2];
1011       y[0] = y[1];
1012       y[1] = y[2];
1013    }
1014    if (i == 1)
1015    {
1016       // No change to U and V we've stalled!
1017       cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t;
1018       eval_modulus(t, U, V);
1019       U.swap(V);
1020       V.swap(t);
1021       return;
1022    }
1023    //
1024    // Update U and V.
1025    // We have:
1026    //
1027    // U = x[0]U + y[0]V and
1028    // V = x[1]U + y[1]V.
1029    //
1030    // But since we track only absolute values of x and y
1031    // we have to take account of the implied signs and perform
1032    // the appropriate subtraction depending on the whether i is
1033    // even or odd:
1034    //
1035    std::size_t ts = U.size() + 1;
1036    cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t1(storage, ts), t2(storage, ts), t3(storage, ts);
1037    eval_multiply(t1, U, x[0]);
1038    eval_multiply(t2, V, y[0]);
1039    eval_multiply(t3, U, x[1]);
1040    if ((i & 1u) == 0)
1041    {
1042       if (x[0] == 0)
1043          U = t2;
1044       else
1045       {
1046          BOOST_MP_ASSERT(t2.compare(t1) >= 0);
1047          eval_subtract(U, t2, t1);
1048          BOOST_MP_ASSERT(U.sign() == false);
1049       }
1050    }
1051    else
1052    {
1053       BOOST_MP_ASSERT(t1.compare(t2) >= 0);
1054       eval_subtract(U, t1, t2);
1055       BOOST_MP_ASSERT(U.sign() == false);
1056    }
1057    eval_multiply(t2, V, y[1]);
1058    if (i & 1u)
1059    {
1060       if (x[1] == 0)
1061          V = t2;
1062       else
1063       {
1064          BOOST_MP_ASSERT(t2.compare(t3) >= 0);
1065          eval_subtract(V, t2, t3);
1066          BOOST_MP_ASSERT(V.sign() == false);
1067       }
1068    }
1069    else
1070    {
1071       BOOST_MP_ASSERT(t3.compare(t2) >= 0);
1072       eval_subtract(V, t3, t2);
1073       BOOST_MP_ASSERT(V.sign() == false);
1074    }
1075    BOOST_MP_ASSERT(U.compare(V) >= 0);
1076    BOOST_MP_ASSERT(lu > eval_msb(U));
1077 #ifdef BOOST_MP_GCD_DEBUG
1078
1079    BOOST_MP_ASSERT(UU == U);
1080    BOOST_MP_ASSERT(VV == V);
1081
1082    extern std::size_t total_lehmer_gcd_calls;
1083    extern std::size_t total_lehmer_gcd_bits_saved;
1084    extern std::size_t total_lehmer_gcd_cycles;
1085
1086    ++total_lehmer_gcd_calls;
1087    total_lehmer_gcd_bits_saved += lu - eval_msb(U);
1088    total_lehmer_gcd_cycles += i;
1089 #endif
1090    if (lu < 2048)
1091    {
1092       //
1093       // Since we have stripped all common powers of 2 from U and V at the start
1094       // if either are even at this point, we can remove stray powers of 2 now.
1095       // Note that it is not possible for *both* U and V to be even at this point.
1096       //
1097       // This has an adverse effect on performance for high bit counts, but has
1098       // a significant positive effect for smaller counts.
1099       //
1100       if ((U.limbs()[0] & 1u) == 0)
1101       {
1102          eval_right_shift(U, eval_lsb(U));
1103          if (U.compare(V) < 0)
1104             U.swap(V);
1105       }
1106       else if ((V.limbs()[0] & 1u) == 0)
1107       {
1108          eval_right_shift(V, eval_lsb(V));
1109       }
1110    }
1111    storage.deallocate(ts * 3);
1112 }
1113
1114 #endif
1115
1116 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
1117 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
1118 eval_gcd(
1119    cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
1120    const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
1121    const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& b)
1122 {
1123    using default_ops::eval_get_sign;
1124    using default_ops::eval_is_zero;
1125    using default_ops::eval_lsb;
1126
1127    if (a.size() == 1)
1128    {
1129       eval_gcd(result, b, *a.limbs());
1130       return;
1131    }
1132    if (b.size() == 1)
1133    {
1134       eval_gcd(result, a, *b.limbs());
1135       return;
1136    }
1137    std::size_t temp_size = (std::max)(a.size(), b.size()) + 1;
1138    typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::scoped_shared_storage storage(a, temp_size * 6);
1139
1140    cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> U(storage, temp_size);
1141    cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> V(storage, temp_size);
1142    cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t(storage, temp_size);
1143    U = a;
1144    V = b;
1145
1146    int s = eval_get_sign(U);
1147
1148    /* GCD(0,x) := x */
1149    if (s < 0)
1150    {
1151       U.negate();
1152    }
1153    else if (s == 0)
1154    {
1155       result = V;
1156       return;
1157    }
1158    s = eval_get_sign(V);
1159    if (s < 0)
1160    {
1161       V.negate();
1162    }
1163    else if (s == 0)
1164    {
1165       result = U;
1166       return;
1167    }
1168    //
1169    // Remove common factors of 2:
1170    //
1171    std::size_t us = eval_lsb(U);
1172    std::size_t vs = eval_lsb(V);
1173    std::size_t shift = (std::min)(us, vs);
1174    if (us)
1175       eval_right_shift(U, us);
1176    if (vs)
1177       eval_right_shift(V, vs);
1178
1179    if (U.compare(V) < 0)
1180       U.swap(V);
1181
1182    while (!eval_is_zero(V))
1183    {
1184       if (U.size() <= 2)
1185       {
1186          //
1187          // Special case: if V has no more than 2 limbs
1188          // then we can reduce U and V to a pair of integers and perform
1189          // direct integer gcd:
1190          //
1191          if (U.size() == 1)
1192             U = eval_gcd(*V.limbs(), *U.limbs());
1193          else
1194          {
1195             double_limb_type i = U.limbs()[0] | (static_cast<double_limb_type>(U.limbs()[1]) << sizeof(limb_type) * CHAR_BIT);
1196             double_limb_type j = (V.size() == 1) ? *V.limbs() : V.limbs()[0] | (static_cast<double_limb_type>(V.limbs()[1]) << sizeof(limb_type) * CHAR_BIT);
1197             U = eval_gcd(i, j);
1198          }
1199          break;
1200       }
1201       std::size_t lu = eval_msb(U) + 1;
1202       std::size_t lv = eval_msb(V) + 1;
1203 #ifndef BOOST_MP_NO_CONSTEXPR_DETECTION
1204       if (!BOOST_MP_IS_CONST_EVALUATED(lu) && (lu - lv <= bits_per_limb / 2))
1205 #else
1206       if (lu - lv <= bits_per_limb / 2)
1207 #endif
1208       {
1209          eval_gcd_lehmer(U, V, lu, storage);
1210       }
1211       else
1212       {
1213          eval_modulus(t, U, V);
1214          U.swap(V);
1215          V.swap(t);
1216       }
1217    }
1218    result = U;
1219    if (shift)
1220       eval_left_shift(result, shift);
1221 }
1222 //
1223 // Now again for trivial backends:
1224 //
1225 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
1226 BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
1227 eval_gcd(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& b) noexcept
1228 {
1229    *result.limbs() = boost::multiprecision::detail::constexpr_gcd(*a.limbs(), *b.limbs());
1230    result.sign(false);
1231 }
1232 // This one is only enabled for unchecked cpp_int's, for checked int's we need the checking in the default version:
1233 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
1234 BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && (Checked1 == unchecked)>::type
1235 eval_lcm(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& b) noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
1236 {
1237    *result.limbs() = boost::multiprecision::detail::constexpr_lcm(*a.limbs(), *b.limbs());
1238    result.normalize(); // result may overflow the specified number of bits
1239    result.sign(false);
1240 }
1241
1242 inline void conversion_overflow(const std::integral_constant<int, checked>&)
1243 {
1244    BOOST_MP_THROW_EXCEPTION(std::overflow_error("Overflow in conversion to narrower type"));
1245 }
1246 inline BOOST_MP_CXX14_CONSTEXPR void conversion_overflow(const std::integral_constant<int, unchecked>&) {}
1247
1248 #if defined(__clang__) && defined(__MINGW32__)
1249 //
1250 // clang-11 on Mingw segfaults on conversion of __int128 -> float.
1251 // See: https://bugs.llvm.org/show_bug.cgi?id=48941
1252 // These workarounds pass everything through an intermediate uint64_t.
1253 //
1254 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
1255 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
1256     is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && is_signed_number<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && std::is_same<typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::local_limb_type, double_limb_type>::value>::type
1257 eval_convert_to(float* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val)
1258 {
1259    float f = static_cast<std::uint64_t>((*val.limbs()) >> 64);
1260    *result = std::ldexp(f, 64);
1261    *result += static_cast<std::uint64_t>((*val.limbs()));
1262    if(val.sign())
1263       *result = -*result;
1264 }
1265 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
1266 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
1267     is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && is_signed_number<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && std::is_same<typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::local_limb_type, double_limb_type>::value>::type
1268 eval_convert_to(double* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val)
1269 {
1270    float f = static_cast<std::uint64_t>((*val.limbs()) >> 64);
1271    *result = std::ldexp(f, 64);
1272    *result += static_cast<std::uint64_t>((*val.limbs()));
1273    if(val.sign())
1274       *result = -*result;
1275 }
1276 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
1277 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
1278     is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && is_signed_number<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && std::is_same<typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::local_limb_type, double_limb_type>::value>::type
1279 eval_convert_to(long double* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val)
1280 {
1281    float f = static_cast<std::uint64_t>((*val.limbs()) >> 64);
1282    *result = std::ldexp(f, 64);
1283    *result += static_cast<std::uint64_t>((*val.limbs()));
1284    if(val.sign())
1285       *result = -*result;
1286 }
1287 #endif
1288
1289 template <class R, std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
1290 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
1291     is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && is_signed_number<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && std::is_convertible<typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::local_limb_type, R>::value>::type
1292 eval_convert_to(R* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val)
1293 {
1294    BOOST_IF_CONSTEXPR(std::numeric_limits<R>::is_specialized)
1295    {
1296       using common_type = typename std::common_type<R, typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::local_limb_type>::type;
1297
1298       if (static_cast<common_type>(*val.limbs()) > static_cast<common_type>((std::numeric_limits<R>::max)()))
1299       {
1300          if (val.isneg())
1301          {
1302             check_is_negative(std::integral_constant < bool, (boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value) || (number_category<R>::value == number_kind_floating_point) > ());
1303             if (static_cast<common_type>(*val.limbs()) > -static_cast<common_type>((std::numeric_limits<R>::min)()))
1304                conversion_overflow(typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type());
1305             *result = (std::numeric_limits<R>::min)();
1306          }
1307          else
1308          {
1309             conversion_overflow(typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type());
1310             *result = boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value ? (std::numeric_limits<R>::max)() : static_cast<R>(*val.limbs());
1311          }
1312       }
1313       else
1314       {
1315          *result = static_cast<R>(*val.limbs());
1316          if (val.isneg())
1317          {
1318             check_is_negative(std::integral_constant < bool, (boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value) || (number_category<R>::value == number_kind_floating_point) > ());
1319             *result = negate_integer(*result, std::integral_constant < bool, is_signed_number<R>::value || (number_category<R>::value == number_kind_floating_point) > ());
1320          }
1321       }
1322    }
1323    else
1324    {
1325       *result = static_cast<R>(*val.limbs());
1326       if (val.isneg())
1327       {
1328          check_is_negative(std::integral_constant<bool, (boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value) || (number_category<R>::value == number_kind_floating_point) > ());
1329          *result = negate_integer(*result, std::integral_constant<bool, is_signed_number<R>::value || (number_category<R>::value == number_kind_floating_point) > ());
1330       }
1331    }
1332 }
1333
1334 template <class R, std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
1335 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
1336     is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && is_unsigned_number<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && std::is_convertible<typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::local_limb_type, R>::value>::type
1337 eval_convert_to(R* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val)
1338 {
1339    BOOST_IF_CONSTEXPR(std::numeric_limits<R>::is_specialized)
1340    {
1341       using common_type = typename std::common_type<R, typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::local_limb_type>::type;
1342
1343       if(static_cast<common_type>(*val.limbs()) > static_cast<common_type>((std::numeric_limits<R>::max)()))
1344       {
1345          conversion_overflow(typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type());
1346          *result = boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value ? (std::numeric_limits<R>::max)() : static_cast<R>(*val.limbs());
1347       }
1348       else
1349          *result = static_cast<R>(*val.limbs());
1350    }
1351    else
1352       *result = static_cast<R>(*val.limbs());
1353 }
1354
1355 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
1356 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, std::size_t>::type
1357 eval_lsb(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
1358 {
1359    using default_ops::eval_get_sign;
1360    if (eval_get_sign(a) == 0)
1361    {
1362       BOOST_MP_THROW_EXCEPTION(std::domain_error("No bits were set in the operand."));
1363    }
1364    if (a.sign())
1365    {
1366       BOOST_MP_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined."));
1367    }
1368    //
1369    // Find the index of the least significant bit within that limb:
1370    //
1371    return boost::multiprecision::detail::find_lsb(*a.limbs());
1372 }
1373
1374 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
1375 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, std::size_t>::type
1376 eval_msb_imp(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
1377 {
1378    //
1379    // Find the index of the least significant bit within that limb:
1380    //
1381    return boost::multiprecision::detail::find_msb(*a.limbs());
1382 }
1383
1384 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
1385 inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, std::size_t>::type
1386 eval_msb(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
1387 {
1388    using default_ops::eval_get_sign;
1389    if (eval_get_sign(a) == 0)
1390    {
1391       BOOST_MP_THROW_EXCEPTION(std::domain_error("No bits were set in the operand."));
1392    }
1393    if (a.sign())
1394    {
1395       BOOST_MP_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined."));
1396    }
1397    return eval_msb_imp(a);
1398 }
1399
1400 template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
1401 inline BOOST_MP_CXX14_CONSTEXPR std::size_t hash_value(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) noexcept
1402 {
1403    std::size_t result = 0;
1404    for (std::size_t i = 0; i < val.size(); ++i)
1405    {
1406       boost::multiprecision::detail::hash_combine(result, val.limbs()[i]);
1407    }
1408    boost::multiprecision::detail::hash_combine(result, val.sign());
1409    return result;
1410 }
1411
1412 #ifdef BOOST_MSVC
1413 #pragma warning(pop)
1414 #endif
1415
1416 } // Namespace backends
1417
1418 namespace detail {
1419
1420 #ifndef BOOST_MP_STANDALONE
1421 template <typename T>
1422 inline BOOST_CXX14_CONSTEXPR T constexpr_gcd(T a, T b) noexcept
1423 {
1424    return boost::integer::gcd(a, b);
1425 }
1426
1427 template <typename T>
1428 inline BOOST_CXX14_CONSTEXPR T constexpr_lcm(T a, T b) noexcept
1429 {
1430    return boost::integer::lcm(a, b);
1431 }
1432
1433 #else
1434
1435 template <typename T>
1436 inline BOOST_CXX14_CONSTEXPR T constexpr_gcd(T a, T b) noexcept
1437 {
1438    return boost::multiprecision::backends::eval_gcd(a, b);
1439 }
1440
1441 template <typename T>
1442 inline BOOST_CXX14_CONSTEXPR T constexpr_lcm(T a, T b) noexcept
1443 {
1444    const T ab_gcd = boost::multiprecision::detail::constexpr_gcd(a, b);
1445    return (a * b) / ab_gcd;
1446 }
1447
1448 #endif // BOOST_MP_STANDALONE
1449
1450 }
1451
1452 }} // Namespace boost::multiprecision
1453
1454 #endif