123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293 |
- ///////////////////////////////////////////////////////////////
- // Copyright 2012-2020 John Maddock.
- // Copyright 2020 Madhur Chauhan.
- // Distributed under the Boost Software License, Version 1.0.
- // (See accompanying file LICENSE_1_0.txt or copy at
- // https://www.boost.org/LICENSE_1_0.txt)
- //
- // Comparison operators for cpp_int_backend:
- //
- #ifndef BOOST_MP_CPP_INT_MISC_HPP
- #define BOOST_MP_CPP_INT_MISC_HPP
- #include <boost/multiprecision/detail/constexpr.hpp>
- #include <boost/multiprecision/detail/bitscan.hpp> // lsb etc
- #include <boost/integer/common_factor_rt.hpp> // gcd/lcm
- #include <boost/functional/hash_fwd.hpp>
- #include <numeric> // std::gcd
- #ifdef BOOST_MSVC
- #pragma warning(push)
- #pragma warning(disable : 4702)
- #pragma warning(disable : 4127) // conditional expression is constant
- #pragma warning(disable : 4146) // unary minus operator applied to unsigned type, result still unsigned
- #endif
- namespace boost { namespace multiprecision { namespace backends {
- template <class T, bool has_limits = std::numeric_limits<T>::is_specialized>
- struct numeric_limits_workaround : public std::numeric_limits<T>
- {
- };
- template <class R>
- struct numeric_limits_workaround<R, false>
- {
- static constexpr unsigned digits = ~static_cast<R>(0) < 0 ? sizeof(R) * CHAR_BIT - 1 : sizeof(R) * CHAR_BIT;
- static constexpr R (min)(){ return (static_cast<R>(-1) < 0) ? static_cast<R>(1) << digits : 0; }
- static constexpr R (max)() { return (static_cast<R>(-1) < 0) ? ~(static_cast<R>(1) << digits) : ~static_cast<R>(0); }
- };
- template <class R, class CppInt>
- BOOST_MP_CXX14_CONSTEXPR void check_in_range(const CppInt& val, const std::integral_constant<int, checked>&)
- {
- using cast_type = typename boost::multiprecision::detail::canonical<R, CppInt>::type;
- if (val.sign())
- {
- BOOST_IF_CONSTEXPR (boost::multiprecision::detail::is_signed<R>::value == false)
- BOOST_THROW_EXCEPTION(std::range_error("Attempt to assign a negative value to an unsigned type."));
- if (val.compare(static_cast<cast_type>((numeric_limits_workaround<R>::min)())) < 0)
- BOOST_THROW_EXCEPTION(std::overflow_error("Could not convert to the target type - -value is out of range."));
- }
- else
- {
- if (val.compare(static_cast<cast_type>((numeric_limits_workaround<R>::max)())) > 0)
- BOOST_THROW_EXCEPTION(std::overflow_error("Could not convert to the target type - -value is out of range."));
- }
- }
- template <class R, class CppInt>
- inline BOOST_MP_CXX14_CONSTEXPR void check_in_range(const CppInt& /*val*/, const std::integral_constant<int, unchecked>&) noexcept {}
- inline BOOST_MP_CXX14_CONSTEXPR void check_is_negative(const std::integral_constant<bool, true>&) noexcept {}
- inline void check_is_negative(const std::integral_constant<bool, false>&)
- {
- BOOST_THROW_EXCEPTION(std::range_error("Attempt to assign a negative value to an unsigned type."));
- }
- template <class Integer>
- inline BOOST_MP_CXX14_CONSTEXPR Integer negate_integer(Integer i, const std::integral_constant<bool, true>&) noexcept
- {
- return -i;
- }
- template <class Integer>
- inline BOOST_MP_CXX14_CONSTEXPR Integer negate_integer(Integer i, const std::integral_constant<bool, false>&) noexcept
- {
- return ~(i - 1);
- }
- template <class R, unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- 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
- eval_convert_to(R* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& backend)
- {
- using checked_type = std::integral_constant<int, Checked1>;
- check_in_range<R>(backend, checked_type());
- BOOST_IF_CONSTEXPR(numeric_limits_workaround<R>::digits < cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits)
- {
- 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]))
- {
- *result = (numeric_limits_workaround<R>::min)();
- return;
- }
- 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])
- {
- *result = (numeric_limits_workaround<R>::max)();
- return;
- }
- else
- *result = static_cast<R>(backend.limbs()[0]);
- }
- else
- *result = static_cast<R>(backend.limbs()[0]);
- unsigned shift = cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
- unsigned i = 1;
- BOOST_IF_CONSTEXPR(numeric_limits_workaround<R>::digits > cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits)
- {
- while ((i < backend.size()) && (shift < static_cast<unsigned>(numeric_limits_workaround<R>::digits - cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits)))
- {
- *result += static_cast<R>(backend.limbs()[i]) << shift;
- shift += cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
- ++i;
- }
- //
- // We have one more limb to extract, but may not need all the bits, so treat this as a special case:
- //
- if (i < backend.size())
- {
- 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>(1u) << (numeric_limits_workaround<R>::digits - shift)) - 1;
- *result += (static_cast<R>(backend.limbs()[i]) & mask) << shift;
- if ((static_cast<R>(backend.limbs()[i]) & static_cast<limb_type>(~mask)) || (i + 1 < backend.size()))
- {
- // Overflow:
- if (backend.sign())
- {
- check_is_negative(boost::multiprecision::detail::is_signed<R>());
- *result = (numeric_limits_workaround<R>::min)();
- }
- else if (boost::multiprecision::detail::is_signed<R>::value)
- *result = (numeric_limits_workaround<R>::max)();
- return;
- }
- }
- }
- else if (backend.size() > 1)
- {
- // Overflow:
- if (backend.sign())
- {
- check_is_negative(boost::multiprecision::detail::is_signed<R>());
- *result = (numeric_limits_workaround<R>::min)();
- }
- else if (boost::multiprecision::detail::is_signed<R>::value)
- *result = (numeric_limits_workaround<R>::max)();
- return;
- }
- if (backend.sign())
- {
- check_is_negative(std::integral_constant<bool, boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value>());
- *result = negate_integer(*result, std::integral_constant<bool, boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value>());
- }
- }
- template <class R, unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- 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
- eval_convert_to(R* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& backend) noexcept(boost::multiprecision::detail::is_arithmetic<R>::value)
- {
- typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::const_limb_pointer p = backend.limbs();
- unsigned shift = cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
- *result = static_cast<R>(*p);
- for (unsigned i = 1; i < backend.size(); ++i)
- {
- *result += static_cast<R>(std::ldexp(static_cast<long double>(p[i]), shift));
- shift += cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
- }
- if (backend.sign())
- *result = -*result;
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- 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
- eval_is_zero(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) noexcept
- {
- return (val.size() == 1) && (val.limbs()[0] == 0);
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- 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
- eval_get_sign(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) noexcept
- {
- return eval_is_zero(val) ? 0 : val.sign() ? -1 : 1;
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- 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
- 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))
- {
- result = val;
- result.sign(false);
- }
- //
- // Get the location of the least-significant-bit:
- //
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, unsigned>::type
- eval_lsb(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
- {
- using default_ops::eval_get_sign;
- if (eval_get_sign(a) == 0)
- {
- BOOST_THROW_EXCEPTION(std::domain_error("No bits were set in the operand."));
- }
- if (a.sign())
- {
- BOOST_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined."));
- }
- //
- // Find the index of the least significant limb that is non-zero:
- //
- unsigned index = 0;
- while (!a.limbs()[index] && (index < a.size()))
- ++index;
- //
- // Find the index of the least significant bit within that limb:
- //
- unsigned result = boost::multiprecision::detail::find_lsb(a.limbs()[index]);
- return result + index * cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
- }
- //
- // Get the location of the most-significant-bit:
- //
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, unsigned>::type
- eval_msb_imp(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
- {
- //
- // Find the index of the most significant bit that is non-zero:
- //
- return (a.size() - 1) * cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits + boost::multiprecision::detail::find_msb(a.limbs()[a.size() - 1]);
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, unsigned>::type
- eval_msb(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
- {
- using default_ops::eval_get_sign;
- if (eval_get_sign(a) == 0)
- {
- BOOST_THROW_EXCEPTION(std::domain_error("No bits were set in the operand."));
- }
- if (a.sign())
- {
- BOOST_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined."));
- }
- return eval_msb_imp(a);
- }
- #ifdef BOOST_GCC
- //
- // We really shouldn't need to be disabling this warning, but it really does appear to be
- // spurious. The warning appears only when in release mode, and asserts are on.
- //
- #pragma GCC diagnostic push
- #pragma GCC diagnostic ignored "-Warray-bounds"
- #endif
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, bool>::type
- eval_bit_test(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val, unsigned index) noexcept
- {
- unsigned offset = index / cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
- unsigned shift = index % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
- limb_type mask = shift ? limb_type(1u) << shift : limb_type(1u);
- if (offset >= val.size())
- return false;
- return val.limbs()[offset] & mask ? true : false;
- }
- #ifdef BOOST_GCC
- #pragma GCC diagnostic pop
- #endif
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
- eval_bit_set(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val, unsigned index)
- {
- unsigned offset = index / cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
- unsigned shift = index % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
- limb_type mask = shift ? limb_type(1u) << shift : limb_type(1u);
- if (offset >= val.size())
- {
- unsigned os = val.size();
- val.resize(offset + 1, offset + 1);
- if (offset >= val.size())
- return; // fixed precision overflow
- for (unsigned i = os; i <= offset; ++i)
- val.limbs()[i] = 0;
- }
- val.limbs()[offset] |= mask;
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
- eval_bit_unset(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val, unsigned index) noexcept
- {
- unsigned offset = index / cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
- unsigned shift = index % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
- limb_type mask = shift ? limb_type(1u) << shift : limb_type(1u);
- if (offset >= val.size())
- return;
- val.limbs()[offset] &= ~mask;
- val.normalize();
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
- eval_bit_flip(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val, unsigned index)
- {
- unsigned offset = index / cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
- unsigned shift = index % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
- limb_type mask = shift ? limb_type(1u) << shift : limb_type(1u);
- if (offset >= val.size())
- {
- unsigned os = val.size();
- val.resize(offset + 1, offset + 1);
- if (offset >= val.size())
- return; // fixed precision overflow
- for (unsigned i = os; i <= offset; ++i)
- val.limbs()[i] = 0;
- }
- val.limbs()[offset] ^= mask;
- val.normalize();
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
- eval_qr(
- const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x,
- const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& y,
- cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& q,
- cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& r) noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
- {
- divide_unsigned_helper(&q, x, y, r);
- q.sign(x.sign() != y.sign());
- r.sign(x.sign());
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
- eval_qr(
- const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x,
- limb_type y,
- cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& q,
- cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& r) noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
- {
- divide_unsigned_helper(&q, x, y, r);
- q.sign(x.sign());
- r.sign(x.sign());
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class U>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<U>::value>::type eval_qr(
- const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x,
- U y,
- cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& q,
- cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& r) noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
- {
- using default_ops::eval_qr;
- cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t;
- t = y;
- eval_qr(x, t, q, r);
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Integer>
- 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
- eval_integer_modulus(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a, Integer mod)
- {
- BOOST_IF_CONSTEXPR (sizeof(Integer) <= sizeof(limb_type))
- {
- if (mod <= (std::numeric_limits<limb_type>::max)())
- {
- const int n = a.size();
- const double_limb_type two_n_mod = static_cast<limb_type>(1u) + (~static_cast<limb_type>(0u) - mod) % mod;
- limb_type res = a.limbs()[n - 1] % mod;
- for (int i = n - 2; i >= 0; --i)
- res = static_cast<limb_type>((res * two_n_mod + a.limbs()[i]) % mod);
- return res;
- }
- else
- return default_ops::eval_integer_modulus(a, mod);
- }
- else
- {
- return default_ops::eval_integer_modulus(a, mod);
- }
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Integer>
- 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
- eval_integer_modulus(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x, Integer val)
- {
- return eval_integer_modulus(x, boost::multiprecision::detail::unsigned_abs(val));
- }
- BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR limb_type eval_gcd(limb_type u, limb_type v)
- {
- // boundary cases
- if (!u || !v)
- return u | v;
- #if __cpp_lib_gcd_lcm >= 201606L
- return std::gcd(u, v);
- #else
- unsigned shift = boost::multiprecision::detail::find_lsb(u | v);
- u >>= boost::multiprecision::detail::find_lsb(u);
- do
- {
- v >>= boost::multiprecision::detail::find_lsb(v);
- if (u > v)
- std_constexpr::swap(u, v);
- v -= u;
- } while (v);
- return u << shift;
- #endif
- }
- inline BOOST_MP_CXX14_CONSTEXPR double_limb_type eval_gcd(double_limb_type u, double_limb_type v)
- {
- #if (__cpp_lib_gcd_lcm >= 201606L) && (!defined(BOOST_HAS_INT128) || !defined(__STRICT_ANSI__))
- return std::gcd(u, v);
- #else
- unsigned shift = boost::multiprecision::detail::find_lsb(u | v);
- u >>= boost::multiprecision::detail::find_lsb(u);
- do
- {
- v >>= boost::multiprecision::detail::find_lsb(v);
- if (u > v)
- std_constexpr::swap(u, v);
- v -= u;
- } while (v);
- return u << shift;
- #endif
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
- eval_gcd(
- cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
- const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
- limb_type b)
- {
- int s = eval_get_sign(a);
- if (!b || !s)
- {
- result = a;
- *result.limbs() |= b;
- }
- else
- {
- eval_modulus(result, a, b);
- limb_type& res = *result.limbs();
- res = eval_gcd(res, b);
- }
- result.sign(false);
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
- eval_gcd(
- cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
- const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
- double_limb_type b)
- {
- int s = eval_get_sign(a);
- if (!b || !s)
- {
- if (!s)
- result = b;
- else
- result = a;
- return;
- }
- double_limb_type res = 0;
- if(a.sign() == 0)
- res = eval_integer_modulus(a, b);
- else
- {
- cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t(a);
- t.negate();
- res = eval_integer_modulus(t, b);
- }
- res = eval_gcd(res, b);
- result = res;
- result.sign(false);
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
- eval_gcd(
- cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
- const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
- signed_double_limb_type v)
- {
- eval_gcd(result, a, static_cast<double_limb_type>(v < 0 ? -v : v));
- }
- //
- // These 2 overloads take care of gcd against an (unsigned) short etc:
- //
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Integer>
- 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
- eval_gcd(
- cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
- const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
- const Integer& v)
- {
- eval_gcd(result, a, static_cast<limb_type>(v));
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Integer>
- 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
- eval_gcd(
- cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
- const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
- const Integer& v)
- {
- eval_gcd(result, a, static_cast<limb_type>(v < 0 ? -v : v));
- }
- //
- // What follows is Lehmer's GCD algorithm:
- // Essentially this uses the leading digit(s) of U and V
- // only to run a "simulated" Euclid algorithm. It stops
- // when the calculated quotient differs from what would have been
- // the true quotient. At that point the cosequences are used to
- // calculate the new U and V. A nice lucid description appears
- // in "An Analysis of Lehmer's Euclidean GCD Algorithm",
- // by Jonathan Sorenson. https://www.researchgate.net/publication/2424634_An_Analysis_of_Lehmer%27s_Euclidean_GCD_Algorithm
- // DOI: 10.1145/220346.220378.
- //
- // There are two versions of this algorithm here, and both are "double digit"
- // variations: which is to say if there are k bits per limb, then they extract
- // 2k bits into a double_limb_type and then run the algorithm on that. The first
- // version is a straightforward version of the algorithm, and is designed for
- // situations where double_limb_type is a native integer (for example where
- // limb_type is a 32-bit integer on a 64-bit machine). For 32-bit limbs it
- // reduces the size of U by about 30 bits per call. The second is a more complex
- // version for situations where double_limb_type is a synthetic type: for example
- // __int128. For 64 bit limbs it reduces the size of U by about 62 bits per call.
- //
- // The complexity of the algorithm given by Sorenson is roughly O(ln^2(N)) for
- // two N bit numbers.
- //
- // The original double-digit version of the algorithm is described in:
- //
- // "A Double Digit Lehmer-Euclid Algorithm for Finding the GCD of Long Integers",
- // Tudor Jebelean, J Symbolic Computation, 1995 (19), 145.
- //
- #ifndef BOOST_HAS_INT128
- //
- // When double_limb_type is a native integer type then we should just use it and not worry about the consequences.
- // This can eliminate approximately a full limb with each call.
- //
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Storage>
- void eval_gcd_lehmer(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& U, cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& V, unsigned lu, Storage& storage)
- {
- //
- // Extract the leading 2 * bits_per_limb bits from U and V:
- //
- unsigned h = lu % bits_per_limb;
- double_limb_type u = (static_cast<double_limb_type>((U.limbs()[U.size() - 1])) << bits_per_limb) | U.limbs()[U.size() - 2];
- 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];
- if (h)
- {
- u <<= bits_per_limb - h;
- u |= U.limbs()[U.size() - 3] >> h;
- v <<= bits_per_limb - h;
- v |= V.limbs()[U.size() - 3] >> h;
- }
- //
- // Co-sequences x an y: we need only the last 3 values of these,
- // the first 2 values are known correct, the third gets checked
- // in each loop operation, and we terminate when they go wrong.
- //
- // x[i+0] is positive for even i.
- // y[i+0] is positive for odd i.
- //
- // However we track only absolute values here:
- //
- double_limb_type x[3] = {1, 0};
- double_limb_type y[3] = {0, 1};
- unsigned i = 0;
- #ifdef BOOST_MP_GCD_DEBUG
- cpp_int UU, VV;
- UU = U;
- VV = V;
- #endif
- while (true)
- {
- double_limb_type q = u / v;
- x[2] = x[0] + q * x[1];
- y[2] = y[0] + q * y[1];
- double_limb_type tu = u;
- u = v;
- v = tu - q * v;
- ++i;
- //
- // We must make sure that y[2] occupies a single limb otherwise
- // the multiprecision multiplications below would be much more expensive.
- // This can sometimes lose us one iteration, but is worth it for improved
- // calculation efficiency.
- //
- if (y[2] >> bits_per_limb)
- break;
- //
- // These are Jebelean's exact termination conditions:
- //
- if ((i & 1u) == 0)
- {
- BOOST_ASSERT(u > v);
- if ((v < x[2]) || ((u - v) < (y[2] + y[1])))
- break;
- }
- else
- {
- BOOST_ASSERT(u > v);
- if ((v < y[2]) || ((u - v) < (x[2] + x[1])))
- break;
- }
- #ifdef BOOST_MP_GCD_DEBUG
- BOOST_ASSERT(q == UU / VV);
- UU %= VV;
- UU.swap(VV);
- #endif
- x[0] = x[1];
- x[1] = x[2];
- y[0] = y[1];
- y[1] = y[2];
- }
- if (i == 1)
- {
- // No change to U and V we've stalled!
- cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t;
- eval_modulus(t, U, V);
- U.swap(V);
- V.swap(t);
- return;
- }
- //
- // Update U and V.
- // We have:
- //
- // U = x[0]U + y[0]V and
- // V = x[1]U + y[1]V.
- //
- // But since we track only absolute values of x and y
- // we have to take account of the implied signs and perform
- // the appropriate subtraction depending on the whether i is
- // even or odd:
- //
- unsigned ts = U.size() + 1;
- cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t1(storage, ts), t2(storage, ts), t3(storage, ts);
- eval_multiply(t1, U, static_cast<limb_type>(x[0]));
- eval_multiply(t2, V, static_cast<limb_type>(y[0]));
- eval_multiply(t3, U, static_cast<limb_type>(x[1]));
- if ((i & 1u) == 0)
- {
- if (x[0] == 0)
- U = t2;
- else
- {
- BOOST_ASSERT(t2.compare(t1) >= 0);
- eval_subtract(U, t2, t1);
- BOOST_ASSERT(U.sign() == false);
- }
- }
- else
- {
- BOOST_ASSERT(t1.compare(t2) >= 0);
- eval_subtract(U, t1, t2);
- BOOST_ASSERT(U.sign() == false);
- }
- eval_multiply(t2, V, static_cast<limb_type>(y[1]));
- if (i & 1u)
- {
- if (x[1] == 0)
- V = t2;
- else
- {
- BOOST_ASSERT(t2.compare(t3) >= 0);
- eval_subtract(V, t2, t3);
- BOOST_ASSERT(V.sign() == false);
- }
- }
- else
- {
- BOOST_ASSERT(t3.compare(t2) >= 0);
- eval_subtract(V, t3, t2);
- BOOST_ASSERT(V.sign() == false);
- }
- BOOST_ASSERT(U.compare(V) >= 0);
- BOOST_ASSERT(lu > eval_msb(U));
- #ifdef BOOST_MP_GCD_DEBUG
- BOOST_ASSERT(UU == U);
- BOOST_ASSERT(VV == V);
- extern unsigned total_lehmer_gcd_calls;
- extern unsigned total_lehmer_gcd_bits_saved;
- extern unsigned total_lehmer_gcd_cycles;
- ++total_lehmer_gcd_calls;
- total_lehmer_gcd_bits_saved += lu - eval_msb(U);
- total_lehmer_gcd_cycles += i;
- #endif
- if (lu < 2048)
- {
- //
- // Since we have stripped all common powers of 2 from U and V at the start
- // if either are even at this point, we can remove stray powers of 2 now.
- // Note that it is not possible for *both* U and V to be even at this point.
- //
- // This has an adverse effect on performance for high bit counts, but has
- // a significant positive effect for smaller counts.
- //
- if ((U.limbs()[0] & 1u) == 0)
- {
- eval_right_shift(U, eval_lsb(U));
- if (U.compare(V) < 0)
- U.swap(V);
- }
- else if ((V.limbs()[0] & 1u) == 0)
- {
- eval_right_shift(V, eval_lsb(V));
- }
- }
- storage.deallocate(ts * 3);
- }
- #else
- //
- // This branch is taken when double_limb_type is a synthetic type with no native hardware support.
- // For example __int128. The assumption is that add/subtract/multiply of double_limb_type are efficient,
- // but that division is very slow.
- //
- // We begin with a specialized routine for division.
- // We know that u > v > ~limb_type(0), and therefore
- // that the result will fit into a single limb_type.
- // We also know that most of the time this is called the result will be 1.
- // For small limb counts, this almost doubles the performance of Lehmer's routine!
- //
- BOOST_FORCEINLINE void divide_subtract(limb_type& q, double_limb_type& u, const double_limb_type& v)
- {
- BOOST_ASSERT(q == 1); // precondition on entry.
- u -= v;
- while (u >= v)
- {
- u -= v;
- if (++q > 30)
- {
- limb_type t = u / v;
- u -= t * v;
- q += t;
- }
- }
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Storage>
- void eval_gcd_lehmer(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& U, cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& V, unsigned lu, Storage& storage)
- {
- //
- // Extract the leading 2*bits_per_limb bits from U and V:
- //
- unsigned h = lu % bits_per_limb;
- double_limb_type u, v;
- if (h)
- {
- u = (static_cast<double_limb_type>((U.limbs()[U.size() - 1])) << bits_per_limb) | U.limbs()[U.size() - 2];
- v = (static_cast<double_limb_type>((V.size() < U.size() ? 0 : V.limbs()[V.size() - 1])) << bits_per_limb) | V.limbs()[U.size() - 2];
- u <<= bits_per_limb - h;
- u |= U.limbs()[U.size() - 3] >> h;
- v <<= bits_per_limb - h;
- v |= V.limbs()[U.size() - 3] >> h;
- }
- else
- {
- u = (static_cast<double_limb_type>(U.limbs()[U.size() - 1]) << bits_per_limb) | U.limbs()[U.size() - 2];
- v = (static_cast<double_limb_type>(V.limbs()[U.size() - 1]) << bits_per_limb) | V.limbs()[U.size() - 2];
- }
- //
- // Cosequences are stored as limb_types, we take care not to overflow these:
- //
- // x[i+0] is positive for even i.
- // y[i+0] is positive for odd i.
- //
- // However we track only absolute values here:
- //
- limb_type x[3] = { 1, 0 };
- limb_type y[3] = { 0, 1 };
- unsigned i = 0;
- #ifdef BOOST_MP_GCD_DEBUG
- cpp_int UU, VV;
- UU = U;
- VV = V;
- #endif
- //
- // We begine by running a single digit version of Lehmer's algorithm, we still have
- // to track u and v at double precision, but this adds only a tiny performance penalty.
- // What we gain is fast division, and fast termination testing.
- // When you see static_cast<limb_type>(u >> bits_per_limb) here, this is really just
- // a direct access to the upper bits_per_limb of the double limb type. For __int128
- // this is simple a load of the upper 64 bits and the "shift" is optimised away.
- //
- double_limb_type old_u, old_v;
- while (true)
- {
- limb_type q = static_cast<limb_type>(u >> bits_per_limb) / static_cast<limb_type>(v >> bits_per_limb);
- x[2] = x[0] + q * x[1];
- y[2] = y[0] + q * y[1];
- double_limb_type tu = u;
- old_u = u;
- old_v = v;
- u = v;
- double_limb_type t = q * v;
- if (tu < t)
- {
- ++i;
- break;
- }
- v = tu - t;
- ++i;
- if ((i & 1u) == 0)
- {
- BOOST_ASSERT(u > v);
- 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])))
- break;
- }
- else
- {
- BOOST_ASSERT(u > v);
- 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])))
- break;
- }
- #ifdef BOOST_MP_GCD_DEBUG
- BOOST_ASSERT(q == UU / VV);
- UU %= VV;
- UU.swap(VV);
- #endif
- x[0] = x[1];
- x[1] = x[2];
- y[0] = y[1];
- y[1] = y[2];
- }
- //
- // We get here when the single digit algorithm has gone wrong, back up i, u and v:
- //
- --i;
- u = old_u;
- v = old_v;
- //
- // Now run the full double-digit algorithm:
- //
- while (true)
- {
- limb_type q = 1;
- double_limb_type tt = u;
- divide_subtract(q, u, v);
- std::swap(u, v);
- tt = y[0] + q * static_cast<double_limb_type>(y[1]);
- //
- // If calculation of y[2] would overflow a single limb, then we *must* terminate.
- // Note that x[2] < y[2] so there is no need to check that as well:
- //
- if (tt >> bits_per_limb)
- {
- ++i;
- break;
- }
- x[2] = x[0] + q * x[1];
- y[2] = tt;
- ++i;
- if ((i & 1u) == 0)
- {
- BOOST_ASSERT(u > v);
- if ((v < x[2]) || ((u - v) < (static_cast<double_limb_type>(y[2]) + y[1])))
- break;
- }
- else
- {
- BOOST_ASSERT(u > v);
- if ((v < y[2]) || ((u - v) < (static_cast<double_limb_type>(x[2]) + x[1])))
- break;
- }
- #ifdef BOOST_MP_GCD_DEBUG
- BOOST_ASSERT(q == UU / VV);
- UU %= VV;
- UU.swap(VV);
- #endif
- x[0] = x[1];
- x[1] = x[2];
- y[0] = y[1];
- y[1] = y[2];
- }
- if (i == 1)
- {
- // No change to U and V we've stalled!
- cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t;
- eval_modulus(t, U, V);
- U.swap(V);
- V.swap(t);
- return;
- }
- //
- // Update U and V.
- // We have:
- //
- // U = x[0]U + y[0]V and
- // V = x[1]U + y[1]V.
- //
- // But since we track only absolute values of x and y
- // we have to take account of the implied signs and perform
- // the appropriate subtraction depending on the whether i is
- // even or odd:
- //
- unsigned ts = U.size() + 1;
- cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t1(storage, ts), t2(storage, ts), t3(storage, ts);
- eval_multiply(t1, U, x[0]);
- eval_multiply(t2, V, y[0]);
- eval_multiply(t3, U, x[1]);
- if ((i & 1u) == 0)
- {
- if (x[0] == 0)
- U = t2;
- else
- {
- BOOST_ASSERT(t2.compare(t1) >= 0);
- eval_subtract(U, t2, t1);
- BOOST_ASSERT(U.sign() == false);
- }
- }
- else
- {
- BOOST_ASSERT(t1.compare(t2) >= 0);
- eval_subtract(U, t1, t2);
- BOOST_ASSERT(U.sign() == false);
- }
- eval_multiply(t2, V, y[1]);
- if (i & 1u)
- {
- if (x[1] == 0)
- V = t2;
- else
- {
- BOOST_ASSERT(t2.compare(t3) >= 0);
- eval_subtract(V, t2, t3);
- BOOST_ASSERT(V.sign() == false);
- }
- }
- else
- {
- BOOST_ASSERT(t3.compare(t2) >= 0);
- eval_subtract(V, t3, t2);
- BOOST_ASSERT(V.sign() == false);
- }
- BOOST_ASSERT(U.compare(V) >= 0);
- BOOST_ASSERT(lu > eval_msb(U));
- #ifdef BOOST_MP_GCD_DEBUG
- BOOST_ASSERT(UU == U);
- BOOST_ASSERT(VV == V);
- extern unsigned total_lehmer_gcd_calls;
- extern unsigned total_lehmer_gcd_bits_saved;
- extern unsigned total_lehmer_gcd_cycles;
- ++total_lehmer_gcd_calls;
- total_lehmer_gcd_bits_saved += lu - eval_msb(U);
- total_lehmer_gcd_cycles += i;
- #endif
- if (lu < 2048)
- {
- //
- // Since we have stripped all common powers of 2 from U and V at the start
- // if either are even at this point, we can remove stray powers of 2 now.
- // Note that it is not possible for *both* U and V to be even at this point.
- //
- // This has an adverse effect on performance for high bit counts, but has
- // a significant positive effect for smaller counts.
- //
- if ((U.limbs()[0] & 1u) == 0)
- {
- eval_right_shift(U, eval_lsb(U));
- if (U.compare(V) < 0)
- U.swap(V);
- }
- else if ((V.limbs()[0] & 1u) == 0)
- {
- eval_right_shift(V, eval_lsb(V));
- }
- }
- storage.deallocate(ts * 3);
- }
- #endif
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
- 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)
- {
- using default_ops::eval_get_sign;
- using default_ops::eval_is_zero;
- using default_ops::eval_lsb;
- if (a.size() == 1)
- {
- eval_gcd(result, b, *a.limbs());
- return;
- }
- if (b.size() == 1)
- {
- eval_gcd(result, a, *b.limbs());
- return;
- }
- unsigned temp_size = (std::max)(a.size(), b.size()) + 1;
- typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::scoped_shared_storage storage(a, temp_size * 6);
- cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> U(storage, temp_size);
- cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> V(storage, temp_size);
- cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t(storage, temp_size);
- U = a;
- V = b;
- int s = eval_get_sign(U);
- /* GCD(0,x) := x */
- if (s < 0)
- {
- U.negate();
- }
- else if (s == 0)
- {
- result = V;
- return;
- }
- s = eval_get_sign(V);
- if (s < 0)
- {
- V.negate();
- }
- else if (s == 0)
- {
- result = U;
- return;
- }
- //
- // Remove common factors of 2:
- //
- unsigned us = eval_lsb(U);
- unsigned vs = eval_lsb(V);
- int shift = (std::min)(us, vs);
- if (us)
- eval_right_shift(U, us);
- if (vs)
- eval_right_shift(V, vs);
- if (U.compare(V) < 0)
- U.swap(V);
- while (!eval_is_zero(V))
- {
- if (U.size() <= 2)
- {
- //
- // Special case: if V has no more than 2 limbs
- // then we can reduce U and V to a pair of integers and perform
- // direct integer gcd:
- //
- if (U.size() == 1)
- U = eval_gcd(*V.limbs(), *U.limbs());
- else
- {
- double_limb_type i = U.limbs()[0] | (static_cast<double_limb_type>(U.limbs()[1]) << sizeof(limb_type) * CHAR_BIT);
- 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);
- U = eval_gcd(i, j);
- }
- break;
- }
- unsigned lu = eval_msb(U) + 1;
- unsigned lv = eval_msb(V) + 1;
- #ifndef BOOST_MP_NO_CONSTEXPR_DETECTION
- if (!BOOST_MP_IS_CONST_EVALUATED(lu) && (lu - lv <= bits_per_limb / 2))
- #else
- if (lu - lv <= bits_per_limb / 2)
- #endif
- {
- eval_gcd_lehmer(U, V, lu, storage);
- }
- else
- {
- eval_modulus(t, U, V);
- U.swap(V);
- V.swap(t);
- }
- }
- result = U;
- if (shift)
- eval_left_shift(result, shift);
- }
- //
- // Now again for trivial backends:
- //
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- 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
- 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
- {
- *result.limbs() = boost::integer::gcd(*a.limbs(), *b.limbs());
- }
- // This one is only enabled for unchecked cpp_int's, for checked int's we need the checking in the default version:
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- 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
- 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))
- {
- *result.limbs() = boost::integer::lcm(*a.limbs(), *b.limbs());
- result.normalize(); // result may overflow the specified number of bits
- }
- inline void conversion_overflow(const std::integral_constant<int, checked>&)
- {
- BOOST_THROW_EXCEPTION(std::overflow_error("Overflow in conversion to narrower type"));
- }
- inline BOOST_MP_CXX14_CONSTEXPR void conversion_overflow(const std::integral_constant<int, unchecked>&) {}
- #if defined(__clang__) && defined(__MINGW32__)
- //
- // clang-11 on Mingw segfaults on conversion of __int128 -> float.
- // See: https://bugs.llvm.org/show_bug.cgi?id=48941
- // These workarounds pass everything through an intermediate uint64_t.
- //
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
- 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
- eval_convert_to(float* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val)
- {
- float f = static_cast<std::uint64_t>((*val.limbs()) >> 64);
- *result = std::ldexp(f, 64);
- *result += static_cast<std::uint64_t>((*val.limbs()));
- if(val.sign())
- *result = -*result;
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
- 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
- eval_convert_to(double* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val)
- {
- float f = static_cast<std::uint64_t>((*val.limbs()) >> 64);
- *result = std::ldexp(f, 64);
- *result += static_cast<std::uint64_t>((*val.limbs()));
- if(val.sign())
- *result = -*result;
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
- 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
- eval_convert_to(long double* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val)
- {
- float f = static_cast<std::uint64_t>((*val.limbs()) >> 64);
- *result = std::ldexp(f, 64);
- *result += static_cast<std::uint64_t>((*val.limbs()));
- if(val.sign())
- *result = -*result;
- }
- #endif
- template <class R, unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
- 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
- eval_convert_to(R* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val)
- {
- using common_type = typename std::common_type<R, typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::local_limb_type>::type;
- BOOST_IF_CONSTEXPR(std::numeric_limits<R>::is_specialized)
- {
- if (static_cast<common_type>(*val.limbs()) > static_cast<common_type>((std::numeric_limits<R>::max)()))
- {
- if (val.isneg())
- {
- 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) > ());
- if (static_cast<common_type>(*val.limbs()) > -static_cast<common_type>((std::numeric_limits<R>::min)()))
- conversion_overflow(typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type());
- *result = (std::numeric_limits<R>::min)();
- }
- else
- {
- conversion_overflow(typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type());
- *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());
- }
- }
- else
- {
- *result = static_cast<R>(*val.limbs());
- if (val.isneg())
- {
- 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) > ());
- *result = negate_integer(*result, std::integral_constant < bool, is_signed_number<R>::value || (number_category<R>::value == number_kind_floating_point) > ());
- }
- }
- }
- else
- {
- *result = static_cast<R>(*val.limbs());
- if (val.isneg())
- {
- 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) > ());
- *result = negate_integer(*result, std::integral_constant<bool, is_signed_number<R>::value || (number_category<R>::value == number_kind_floating_point) > ());
- }
- }
- }
- template <class R, unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
- 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
- eval_convert_to(R* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val)
- {
- using common_type = typename std::common_type<R, typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::local_limb_type>::type;
- BOOST_IF_CONSTEXPR(std::numeric_limits<R>::is_specialized)
- {
- if(static_cast<common_type>(*val.limbs()) > static_cast<common_type>((std::numeric_limits<R>::max)()))
- {
- conversion_overflow(typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type());
- *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());
- }
- else
- *result = static_cast<R>(*val.limbs());
- }
- else
- *result = static_cast<R>(*val.limbs());
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, unsigned>::type
- eval_lsb(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
- {
- using default_ops::eval_get_sign;
- if (eval_get_sign(a) == 0)
- {
- BOOST_THROW_EXCEPTION(std::domain_error("No bits were set in the operand."));
- }
- if (a.sign())
- {
- BOOST_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined."));
- }
- //
- // Find the index of the least significant bit within that limb:
- //
- return boost::multiprecision::detail::find_lsb(*a.limbs());
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, unsigned>::type
- eval_msb_imp(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
- {
- //
- // Find the index of the least significant bit within that limb:
- //
- return boost::multiprecision::detail::find_msb(*a.limbs());
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, unsigned>::type
- eval_msb(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
- {
- using default_ops::eval_get_sign;
- if (eval_get_sign(a) == 0)
- {
- BOOST_THROW_EXCEPTION(std::domain_error("No bits were set in the operand."));
- }
- if (a.sign())
- {
- BOOST_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined."));
- }
- return eval_msb_imp(a);
- }
- template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
- inline BOOST_MP_CXX14_CONSTEXPR std::size_t hash_value(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) noexcept
- {
- std::size_t result = 0;
- for (unsigned i = 0; i < val.size(); ++i)
- {
- boost::hash_combine(result, val.limbs()[i]);
- }
- boost::hash_combine(result, val.sign());
- return result;
- }
- #ifdef BOOST_MSVC
- #pragma warning(pop)
- #endif
- }}} // namespace boost::multiprecision::backends
- #endif
|