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