misc.hpp 55 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293
  1. ///////////////////////////////////////////////////////////////
  2. // Copyright 2012-2020 John Maddock.
  3. // Copyright 2020 Madhur Chauhan.
  4. // Distributed under the Boost Software License, Version 1.0.
  5. // (See accompanying file LICENSE_1_0.txt or copy at
  6. // https://www.boost.org/LICENSE_1_0.txt)
  7. //
  8. // Comparison operators for cpp_int_backend:
  9. //
  10. #ifndef BOOST_MP_CPP_INT_MISC_HPP
  11. #define BOOST_MP_CPP_INT_MISC_HPP
  12. #include <boost/multiprecision/detail/constexpr.hpp>
  13. #include <boost/multiprecision/detail/bitscan.hpp> // lsb etc
  14. #include <boost/integer/common_factor_rt.hpp> // gcd/lcm
  15. #include <boost/functional/hash_fwd.hpp>
  16. #include <numeric> // std::gcd
  17. #ifdef BOOST_MSVC
  18. #pragma warning(push)
  19. #pragma warning(disable : 4702)
  20. #pragma warning(disable : 4127) // conditional expression is constant
  21. #pragma warning(disable : 4146) // unary minus operator applied to unsigned type, result still unsigned
  22. #endif
  23. namespace boost { namespace multiprecision { namespace backends {
  24. template <class T, bool has_limits = std::numeric_limits<T>::is_specialized>
  25. struct numeric_limits_workaround : public std::numeric_limits<T>
  26. {
  27. };
  28. template <class R>
  29. struct numeric_limits_workaround<R, false>
  30. {
  31. static constexpr unsigned digits = ~static_cast<R>(0) < 0 ? sizeof(R) * CHAR_BIT - 1 : sizeof(R) * CHAR_BIT;
  32. static constexpr R (min)(){ return (static_cast<R>(-1) < 0) ? static_cast<R>(1) << digits : 0; }
  33. static constexpr R (max)() { return (static_cast<R>(-1) < 0) ? ~(static_cast<R>(1) << digits) : ~static_cast<R>(0); }
  34. };
  35. template <class R, class CppInt>
  36. BOOST_MP_CXX14_CONSTEXPR void check_in_range(const CppInt& val, const std::integral_constant<int, checked>&)
  37. {
  38. using cast_type = typename boost::multiprecision::detail::canonical<R, CppInt>::type;
  39. if (val.sign())
  40. {
  41. BOOST_IF_CONSTEXPR (boost::multiprecision::detail::is_signed<R>::value == false)
  42. BOOST_THROW_EXCEPTION(std::range_error("Attempt to assign a negative value to an unsigned type."));
  43. if (val.compare(static_cast<cast_type>((numeric_limits_workaround<R>::min)())) < 0)
  44. BOOST_THROW_EXCEPTION(std::overflow_error("Could not convert to the target type - -value is out of range."));
  45. }
  46. else
  47. {
  48. if (val.compare(static_cast<cast_type>((numeric_limits_workaround<R>::max)())) > 0)
  49. BOOST_THROW_EXCEPTION(std::overflow_error("Could not convert to the target type - -value is out of range."));
  50. }
  51. }
  52. template <class R, class CppInt>
  53. inline BOOST_MP_CXX14_CONSTEXPR void check_in_range(const CppInt& /*val*/, const std::integral_constant<int, unchecked>&) noexcept {}
  54. inline BOOST_MP_CXX14_CONSTEXPR void check_is_negative(const std::integral_constant<bool, true>&) noexcept {}
  55. inline void check_is_negative(const std::integral_constant<bool, false>&)
  56. {
  57. BOOST_THROW_EXCEPTION(std::range_error("Attempt to assign a negative value to an unsigned type."));
  58. }
  59. template <class Integer>
  60. inline BOOST_MP_CXX14_CONSTEXPR Integer negate_integer(Integer i, const std::integral_constant<bool, true>&) noexcept
  61. {
  62. return -i;
  63. }
  64. template <class Integer>
  65. inline BOOST_MP_CXX14_CONSTEXPR Integer negate_integer(Integer i, const std::integral_constant<bool, false>&) noexcept
  66. {
  67. return ~(i - 1);
  68. }
  69. template <class R, unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  70. 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
  71. eval_convert_to(R* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& backend)
  72. {
  73. using checked_type = std::integral_constant<int, Checked1>;
  74. check_in_range<R>(backend, checked_type());
  75. BOOST_IF_CONSTEXPR(numeric_limits_workaround<R>::digits < cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits)
  76. {
  77. 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]))
  78. {
  79. *result = (numeric_limits_workaround<R>::min)();
  80. return;
  81. }
  82. 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])
  83. {
  84. *result = (numeric_limits_workaround<R>::max)();
  85. return;
  86. }
  87. else
  88. *result = static_cast<R>(backend.limbs()[0]);
  89. }
  90. else
  91. *result = static_cast<R>(backend.limbs()[0]);
  92. unsigned shift = cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
  93. unsigned i = 1;
  94. BOOST_IF_CONSTEXPR(numeric_limits_workaround<R>::digits > cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits)
  95. {
  96. while ((i < backend.size()) && (shift < static_cast<unsigned>(numeric_limits_workaround<R>::digits - cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits)))
  97. {
  98. *result += static_cast<R>(backend.limbs()[i]) << shift;
  99. shift += cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
  100. ++i;
  101. }
  102. //
  103. // We have one more limb to extract, but may not need all the bits, so treat this as a special case:
  104. //
  105. if (i < backend.size())
  106. {
  107. 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;
  108. *result += (static_cast<R>(backend.limbs()[i]) & mask) << shift;
  109. if ((static_cast<R>(backend.limbs()[i]) & static_cast<limb_type>(~mask)) || (i + 1 < backend.size()))
  110. {
  111. // Overflow:
  112. if (backend.sign())
  113. {
  114. check_is_negative(boost::multiprecision::detail::is_signed<R>());
  115. *result = (numeric_limits_workaround<R>::min)();
  116. }
  117. else if (boost::multiprecision::detail::is_signed<R>::value)
  118. *result = (numeric_limits_workaround<R>::max)();
  119. return;
  120. }
  121. }
  122. }
  123. else if (backend.size() > 1)
  124. {
  125. // Overflow:
  126. if (backend.sign())
  127. {
  128. check_is_negative(boost::multiprecision::detail::is_signed<R>());
  129. *result = (numeric_limits_workaround<R>::min)();
  130. }
  131. else if (boost::multiprecision::detail::is_signed<R>::value)
  132. *result = (numeric_limits_workaround<R>::max)();
  133. return;
  134. }
  135. if (backend.sign())
  136. {
  137. check_is_negative(std::integral_constant<bool, boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value>());
  138. *result = negate_integer(*result, std::integral_constant<bool, boost::multiprecision::detail::is_signed<R>::value && boost::multiprecision::detail::is_integral<R>::value>());
  139. }
  140. }
  141. template <class R, unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  142. 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
  143. eval_convert_to(R* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& backend) noexcept(boost::multiprecision::detail::is_arithmetic<R>::value)
  144. {
  145. typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::const_limb_pointer p = backend.limbs();
  146. unsigned shift = cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
  147. *result = static_cast<R>(*p);
  148. for (unsigned i = 1; i < backend.size(); ++i)
  149. {
  150. *result += static_cast<R>(std::ldexp(static_cast<long double>(p[i]), shift));
  151. shift += cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
  152. }
  153. if (backend.sign())
  154. *result = -*result;
  155. }
  156. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  157. 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
  158. eval_is_zero(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) noexcept
  159. {
  160. return (val.size() == 1) && (val.limbs()[0] == 0);
  161. }
  162. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  163. 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
  164. eval_get_sign(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) noexcept
  165. {
  166. return eval_is_zero(val) ? 0 : val.sign() ? -1 : 1;
  167. }
  168. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  169. 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
  170. 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))
  171. {
  172. result = val;
  173. result.sign(false);
  174. }
  175. //
  176. // Get the location of the least-significant-bit:
  177. //
  178. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  179. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, unsigned>::type
  180. eval_lsb(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
  181. {
  182. using default_ops::eval_get_sign;
  183. if (eval_get_sign(a) == 0)
  184. {
  185. BOOST_THROW_EXCEPTION(std::domain_error("No bits were set in the operand."));
  186. }
  187. if (a.sign())
  188. {
  189. BOOST_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined."));
  190. }
  191. //
  192. // Find the index of the least significant limb that is non-zero:
  193. //
  194. unsigned index = 0;
  195. while (!a.limbs()[index] && (index < a.size()))
  196. ++index;
  197. //
  198. // Find the index of the least significant bit within that limb:
  199. //
  200. unsigned result = boost::multiprecision::detail::find_lsb(a.limbs()[index]);
  201. return result + index * cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
  202. }
  203. //
  204. // Get the location of the most-significant-bit:
  205. //
  206. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  207. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, unsigned>::type
  208. eval_msb_imp(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
  209. {
  210. //
  211. // Find the index of the most significant bit that is non-zero:
  212. //
  213. return (a.size() - 1) * cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits + boost::multiprecision::detail::find_msb(a.limbs()[a.size() - 1]);
  214. }
  215. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  216. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, unsigned>::type
  217. eval_msb(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
  218. {
  219. using default_ops::eval_get_sign;
  220. if (eval_get_sign(a) == 0)
  221. {
  222. BOOST_THROW_EXCEPTION(std::domain_error("No bits were set in the operand."));
  223. }
  224. if (a.sign())
  225. {
  226. BOOST_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined."));
  227. }
  228. return eval_msb_imp(a);
  229. }
  230. #ifdef BOOST_GCC
  231. //
  232. // We really shouldn't need to be disabling this warning, but it really does appear to be
  233. // spurious. The warning appears only when in release mode, and asserts are on.
  234. //
  235. #pragma GCC diagnostic push
  236. #pragma GCC diagnostic ignored "-Warray-bounds"
  237. #endif
  238. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  239. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, bool>::type
  240. eval_bit_test(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val, unsigned index) noexcept
  241. {
  242. unsigned offset = index / cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
  243. unsigned shift = index % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
  244. limb_type mask = shift ? limb_type(1u) << shift : limb_type(1u);
  245. if (offset >= val.size())
  246. return false;
  247. return val.limbs()[offset] & mask ? true : false;
  248. }
  249. #ifdef BOOST_GCC
  250. #pragma GCC diagnostic pop
  251. #endif
  252. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  253. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
  254. eval_bit_set(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val, unsigned index)
  255. {
  256. unsigned offset = index / cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
  257. unsigned shift = index % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
  258. limb_type mask = shift ? limb_type(1u) << shift : limb_type(1u);
  259. if (offset >= val.size())
  260. {
  261. unsigned os = val.size();
  262. val.resize(offset + 1, offset + 1);
  263. if (offset >= val.size())
  264. return; // fixed precision overflow
  265. for (unsigned i = os; i <= offset; ++i)
  266. val.limbs()[i] = 0;
  267. }
  268. val.limbs()[offset] |= mask;
  269. }
  270. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  271. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
  272. eval_bit_unset(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val, unsigned index) noexcept
  273. {
  274. unsigned offset = index / cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
  275. unsigned shift = index % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
  276. limb_type mask = shift ? limb_type(1u) << shift : limb_type(1u);
  277. if (offset >= val.size())
  278. return;
  279. val.limbs()[offset] &= ~mask;
  280. val.normalize();
  281. }
  282. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  283. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
  284. eval_bit_flip(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val, unsigned index)
  285. {
  286. unsigned offset = index / cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
  287. unsigned shift = index % cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
  288. limb_type mask = shift ? limb_type(1u) << shift : limb_type(1u);
  289. if (offset >= val.size())
  290. {
  291. unsigned os = val.size();
  292. val.resize(offset + 1, offset + 1);
  293. if (offset >= val.size())
  294. return; // fixed precision overflow
  295. for (unsigned i = os; i <= offset; ++i)
  296. val.limbs()[i] = 0;
  297. }
  298. val.limbs()[offset] ^= mask;
  299. val.normalize();
  300. }
  301. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  302. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
  303. eval_qr(
  304. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x,
  305. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& y,
  306. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& q,
  307. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& r) noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
  308. {
  309. divide_unsigned_helper(&q, x, y, r);
  310. q.sign(x.sign() != y.sign());
  311. r.sign(x.sign());
  312. }
  313. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  314. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
  315. eval_qr(
  316. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x,
  317. limb_type y,
  318. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& q,
  319. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& r) noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
  320. {
  321. divide_unsigned_helper(&q, x, y, r);
  322. q.sign(x.sign());
  323. r.sign(x.sign());
  324. }
  325. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class U>
  326. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_integral<U>::value>::type eval_qr(
  327. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x,
  328. U y,
  329. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& q,
  330. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& r) noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
  331. {
  332. using default_ops::eval_qr;
  333. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t;
  334. t = y;
  335. eval_qr(x, t, q, r);
  336. }
  337. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Integer>
  338. 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
  339. eval_integer_modulus(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a, Integer mod)
  340. {
  341. BOOST_IF_CONSTEXPR (sizeof(Integer) <= sizeof(limb_type))
  342. {
  343. if (mod <= (std::numeric_limits<limb_type>::max)())
  344. {
  345. const int n = a.size();
  346. const double_limb_type two_n_mod = static_cast<limb_type>(1u) + (~static_cast<limb_type>(0u) - mod) % mod;
  347. limb_type res = a.limbs()[n - 1] % mod;
  348. for (int i = n - 2; i >= 0; --i)
  349. res = static_cast<limb_type>((res * two_n_mod + a.limbs()[i]) % mod);
  350. return res;
  351. }
  352. else
  353. return default_ops::eval_integer_modulus(a, mod);
  354. }
  355. else
  356. {
  357. return default_ops::eval_integer_modulus(a, mod);
  358. }
  359. }
  360. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Integer>
  361. 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
  362. eval_integer_modulus(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& x, Integer val)
  363. {
  364. return eval_integer_modulus(x, boost::multiprecision::detail::unsigned_abs(val));
  365. }
  366. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR limb_type eval_gcd(limb_type u, limb_type v)
  367. {
  368. // boundary cases
  369. if (!u || !v)
  370. return u | v;
  371. #if __cpp_lib_gcd_lcm >= 201606L
  372. return std::gcd(u, v);
  373. #else
  374. unsigned shift = boost::multiprecision::detail::find_lsb(u | v);
  375. u >>= boost::multiprecision::detail::find_lsb(u);
  376. do
  377. {
  378. v >>= boost::multiprecision::detail::find_lsb(v);
  379. if (u > v)
  380. std_constexpr::swap(u, v);
  381. v -= u;
  382. } while (v);
  383. return u << shift;
  384. #endif
  385. }
  386. inline BOOST_MP_CXX14_CONSTEXPR double_limb_type eval_gcd(double_limb_type u, double_limb_type v)
  387. {
  388. #if (__cpp_lib_gcd_lcm >= 201606L) && (!defined(BOOST_HAS_INT128) || !defined(__STRICT_ANSI__))
  389. return std::gcd(u, v);
  390. #else
  391. unsigned shift = boost::multiprecision::detail::find_lsb(u | v);
  392. u >>= boost::multiprecision::detail::find_lsb(u);
  393. do
  394. {
  395. v >>= boost::multiprecision::detail::find_lsb(v);
  396. if (u > v)
  397. std_constexpr::swap(u, v);
  398. v -= u;
  399. } while (v);
  400. return u << shift;
  401. #endif
  402. }
  403. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  404. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
  405. eval_gcd(
  406. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  407. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
  408. limb_type b)
  409. {
  410. int s = eval_get_sign(a);
  411. if (!b || !s)
  412. {
  413. result = a;
  414. *result.limbs() |= b;
  415. }
  416. else
  417. {
  418. eval_modulus(result, a, b);
  419. limb_type& res = *result.limbs();
  420. res = eval_gcd(res, b);
  421. }
  422. result.sign(false);
  423. }
  424. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  425. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
  426. eval_gcd(
  427. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  428. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
  429. double_limb_type b)
  430. {
  431. int s = eval_get_sign(a);
  432. if (!b || !s)
  433. {
  434. if (!s)
  435. result = b;
  436. else
  437. result = a;
  438. return;
  439. }
  440. double_limb_type res = 0;
  441. if(a.sign() == 0)
  442. res = eval_integer_modulus(a, b);
  443. else
  444. {
  445. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t(a);
  446. t.negate();
  447. res = eval_integer_modulus(t, b);
  448. }
  449. res = eval_gcd(res, b);
  450. result = res;
  451. result.sign(false);
  452. }
  453. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  454. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
  455. eval_gcd(
  456. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  457. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
  458. signed_double_limb_type v)
  459. {
  460. eval_gcd(result, a, static_cast<double_limb_type>(v < 0 ? -v : v));
  461. }
  462. //
  463. // These 2 overloads take care of gcd against an (unsigned) short etc:
  464. //
  465. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Integer>
  466. 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
  467. eval_gcd(
  468. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  469. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
  470. const Integer& v)
  471. {
  472. eval_gcd(result, a, static_cast<limb_type>(v));
  473. }
  474. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Integer>
  475. 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
  476. eval_gcd(
  477. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  478. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
  479. const Integer& v)
  480. {
  481. eval_gcd(result, a, static_cast<limb_type>(v < 0 ? -v : v));
  482. }
  483. //
  484. // What follows is Lehmer's GCD algorithm:
  485. // Essentially this uses the leading digit(s) of U and V
  486. // only to run a "simulated" Euclid algorithm. It stops
  487. // when the calculated quotient differs from what would have been
  488. // the true quotient. At that point the cosequences are used to
  489. // calculate the new U and V. A nice lucid description appears
  490. // in "An Analysis of Lehmer's Euclidean GCD Algorithm",
  491. // by Jonathan Sorenson. https://www.researchgate.net/publication/2424634_An_Analysis_of_Lehmer%27s_Euclidean_GCD_Algorithm
  492. // DOI: 10.1145/220346.220378.
  493. //
  494. // There are two versions of this algorithm here, and both are "double digit"
  495. // variations: which is to say if there are k bits per limb, then they extract
  496. // 2k bits into a double_limb_type and then run the algorithm on that. The first
  497. // version is a straightforward version of the algorithm, and is designed for
  498. // situations where double_limb_type is a native integer (for example where
  499. // limb_type is a 32-bit integer on a 64-bit machine). For 32-bit limbs it
  500. // reduces the size of U by about 30 bits per call. The second is a more complex
  501. // version for situations where double_limb_type is a synthetic type: for example
  502. // __int128. For 64 bit limbs it reduces the size of U by about 62 bits per call.
  503. //
  504. // The complexity of the algorithm given by Sorenson is roughly O(ln^2(N)) for
  505. // two N bit numbers.
  506. //
  507. // The original double-digit version of the algorithm is described in:
  508. //
  509. // "A Double Digit Lehmer-Euclid Algorithm for Finding the GCD of Long Integers",
  510. // Tudor Jebelean, J Symbolic Computation, 1995 (19), 145.
  511. //
  512. #ifndef BOOST_HAS_INT128
  513. //
  514. // When double_limb_type is a native integer type then we should just use it and not worry about the consequences.
  515. // This can eliminate approximately a full limb with each call.
  516. //
  517. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Storage>
  518. 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)
  519. {
  520. //
  521. // Extract the leading 2 * bits_per_limb bits from U and V:
  522. //
  523. unsigned h = lu % bits_per_limb;
  524. double_limb_type u = (static_cast<double_limb_type>((U.limbs()[U.size() - 1])) << bits_per_limb) | U.limbs()[U.size() - 2];
  525. 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];
  526. if (h)
  527. {
  528. u <<= bits_per_limb - h;
  529. u |= U.limbs()[U.size() - 3] >> h;
  530. v <<= bits_per_limb - h;
  531. v |= V.limbs()[U.size() - 3] >> h;
  532. }
  533. //
  534. // Co-sequences x an y: we need only the last 3 values of these,
  535. // the first 2 values are known correct, the third gets checked
  536. // in each loop operation, and we terminate when they go wrong.
  537. //
  538. // x[i+0] is positive for even i.
  539. // y[i+0] is positive for odd i.
  540. //
  541. // However we track only absolute values here:
  542. //
  543. double_limb_type x[3] = {1, 0};
  544. double_limb_type y[3] = {0, 1};
  545. unsigned i = 0;
  546. #ifdef BOOST_MP_GCD_DEBUG
  547. cpp_int UU, VV;
  548. UU = U;
  549. VV = V;
  550. #endif
  551. while (true)
  552. {
  553. double_limb_type q = u / v;
  554. x[2] = x[0] + q * x[1];
  555. y[2] = y[0] + q * y[1];
  556. double_limb_type tu = u;
  557. u = v;
  558. v = tu - q * v;
  559. ++i;
  560. //
  561. // We must make sure that y[2] occupies a single limb otherwise
  562. // the multiprecision multiplications below would be much more expensive.
  563. // This can sometimes lose us one iteration, but is worth it for improved
  564. // calculation efficiency.
  565. //
  566. if (y[2] >> bits_per_limb)
  567. break;
  568. //
  569. // These are Jebelean's exact termination conditions:
  570. //
  571. if ((i & 1u) == 0)
  572. {
  573. BOOST_ASSERT(u > v);
  574. if ((v < x[2]) || ((u - v) < (y[2] + y[1])))
  575. break;
  576. }
  577. else
  578. {
  579. BOOST_ASSERT(u > v);
  580. if ((v < y[2]) || ((u - v) < (x[2] + x[1])))
  581. break;
  582. }
  583. #ifdef BOOST_MP_GCD_DEBUG
  584. BOOST_ASSERT(q == UU / VV);
  585. UU %= VV;
  586. UU.swap(VV);
  587. #endif
  588. x[0] = x[1];
  589. x[1] = x[2];
  590. y[0] = y[1];
  591. y[1] = y[2];
  592. }
  593. if (i == 1)
  594. {
  595. // No change to U and V we've stalled!
  596. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t;
  597. eval_modulus(t, U, V);
  598. U.swap(V);
  599. V.swap(t);
  600. return;
  601. }
  602. //
  603. // Update U and V.
  604. // We have:
  605. //
  606. // U = x[0]U + y[0]V and
  607. // V = x[1]U + y[1]V.
  608. //
  609. // But since we track only absolute values of x and y
  610. // we have to take account of the implied signs and perform
  611. // the appropriate subtraction depending on the whether i is
  612. // even or odd:
  613. //
  614. unsigned ts = U.size() + 1;
  615. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t1(storage, ts), t2(storage, ts), t3(storage, ts);
  616. eval_multiply(t1, U, static_cast<limb_type>(x[0]));
  617. eval_multiply(t2, V, static_cast<limb_type>(y[0]));
  618. eval_multiply(t3, U, static_cast<limb_type>(x[1]));
  619. if ((i & 1u) == 0)
  620. {
  621. if (x[0] == 0)
  622. U = t2;
  623. else
  624. {
  625. BOOST_ASSERT(t2.compare(t1) >= 0);
  626. eval_subtract(U, t2, t1);
  627. BOOST_ASSERT(U.sign() == false);
  628. }
  629. }
  630. else
  631. {
  632. BOOST_ASSERT(t1.compare(t2) >= 0);
  633. eval_subtract(U, t1, t2);
  634. BOOST_ASSERT(U.sign() == false);
  635. }
  636. eval_multiply(t2, V, static_cast<limb_type>(y[1]));
  637. if (i & 1u)
  638. {
  639. if (x[1] == 0)
  640. V = t2;
  641. else
  642. {
  643. BOOST_ASSERT(t2.compare(t3) >= 0);
  644. eval_subtract(V, t2, t3);
  645. BOOST_ASSERT(V.sign() == false);
  646. }
  647. }
  648. else
  649. {
  650. BOOST_ASSERT(t3.compare(t2) >= 0);
  651. eval_subtract(V, t3, t2);
  652. BOOST_ASSERT(V.sign() == false);
  653. }
  654. BOOST_ASSERT(U.compare(V) >= 0);
  655. BOOST_ASSERT(lu > eval_msb(U));
  656. #ifdef BOOST_MP_GCD_DEBUG
  657. BOOST_ASSERT(UU == U);
  658. BOOST_ASSERT(VV == V);
  659. extern unsigned total_lehmer_gcd_calls;
  660. extern unsigned total_lehmer_gcd_bits_saved;
  661. extern unsigned total_lehmer_gcd_cycles;
  662. ++total_lehmer_gcd_calls;
  663. total_lehmer_gcd_bits_saved += lu - eval_msb(U);
  664. total_lehmer_gcd_cycles += i;
  665. #endif
  666. if (lu < 2048)
  667. {
  668. //
  669. // Since we have stripped all common powers of 2 from U and V at the start
  670. // if either are even at this point, we can remove stray powers of 2 now.
  671. // Note that it is not possible for *both* U and V to be even at this point.
  672. //
  673. // This has an adverse effect on performance for high bit counts, but has
  674. // a significant positive effect for smaller counts.
  675. //
  676. if ((U.limbs()[0] & 1u) == 0)
  677. {
  678. eval_right_shift(U, eval_lsb(U));
  679. if (U.compare(V) < 0)
  680. U.swap(V);
  681. }
  682. else if ((V.limbs()[0] & 1u) == 0)
  683. {
  684. eval_right_shift(V, eval_lsb(V));
  685. }
  686. }
  687. storage.deallocate(ts * 3);
  688. }
  689. #else
  690. //
  691. // This branch is taken when double_limb_type is a synthetic type with no native hardware support.
  692. // For example __int128. The assumption is that add/subtract/multiply of double_limb_type are efficient,
  693. // but that division is very slow.
  694. //
  695. // We begin with a specialized routine for division.
  696. // We know that u > v > ~limb_type(0), and therefore
  697. // that the result will fit into a single limb_type.
  698. // We also know that most of the time this is called the result will be 1.
  699. // For small limb counts, this almost doubles the performance of Lehmer's routine!
  700. //
  701. BOOST_FORCEINLINE void divide_subtract(limb_type& q, double_limb_type& u, const double_limb_type& v)
  702. {
  703. BOOST_ASSERT(q == 1); // precondition on entry.
  704. u -= v;
  705. while (u >= v)
  706. {
  707. u -= v;
  708. if (++q > 30)
  709. {
  710. limb_type t = u / v;
  711. u -= t * v;
  712. q += t;
  713. }
  714. }
  715. }
  716. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class Storage>
  717. 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)
  718. {
  719. //
  720. // Extract the leading 2*bits_per_limb bits from U and V:
  721. //
  722. unsigned h = lu % bits_per_limb;
  723. double_limb_type u, v;
  724. if (h)
  725. {
  726. u = (static_cast<double_limb_type>((U.limbs()[U.size() - 1])) << bits_per_limb) | U.limbs()[U.size() - 2];
  727. v = (static_cast<double_limb_type>((V.size() < U.size() ? 0 : V.limbs()[V.size() - 1])) << bits_per_limb) | V.limbs()[U.size() - 2];
  728. u <<= bits_per_limb - h;
  729. u |= U.limbs()[U.size() - 3] >> h;
  730. v <<= bits_per_limb - h;
  731. v |= V.limbs()[U.size() - 3] >> h;
  732. }
  733. else
  734. {
  735. u = (static_cast<double_limb_type>(U.limbs()[U.size() - 1]) << bits_per_limb) | U.limbs()[U.size() - 2];
  736. v = (static_cast<double_limb_type>(V.limbs()[U.size() - 1]) << bits_per_limb) | V.limbs()[U.size() - 2];
  737. }
  738. //
  739. // Cosequences are stored as limb_types, we take care not to overflow these:
  740. //
  741. // x[i+0] is positive for even i.
  742. // y[i+0] is positive for odd i.
  743. //
  744. // However we track only absolute values here:
  745. //
  746. limb_type x[3] = { 1, 0 };
  747. limb_type y[3] = { 0, 1 };
  748. unsigned i = 0;
  749. #ifdef BOOST_MP_GCD_DEBUG
  750. cpp_int UU, VV;
  751. UU = U;
  752. VV = V;
  753. #endif
  754. //
  755. // We begine by running a single digit version of Lehmer's algorithm, we still have
  756. // to track u and v at double precision, but this adds only a tiny performance penalty.
  757. // What we gain is fast division, and fast termination testing.
  758. // When you see static_cast<limb_type>(u >> bits_per_limb) here, this is really just
  759. // a direct access to the upper bits_per_limb of the double limb type. For __int128
  760. // this is simple a load of the upper 64 bits and the "shift" is optimised away.
  761. //
  762. double_limb_type old_u, old_v;
  763. while (true)
  764. {
  765. limb_type q = static_cast<limb_type>(u >> bits_per_limb) / static_cast<limb_type>(v >> bits_per_limb);
  766. x[2] = x[0] + q * x[1];
  767. y[2] = y[0] + q * y[1];
  768. double_limb_type tu = u;
  769. old_u = u;
  770. old_v = v;
  771. u = v;
  772. double_limb_type t = q * v;
  773. if (tu < t)
  774. {
  775. ++i;
  776. break;
  777. }
  778. v = tu - t;
  779. ++i;
  780. if ((i & 1u) == 0)
  781. {
  782. BOOST_ASSERT(u > v);
  783. 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])))
  784. break;
  785. }
  786. else
  787. {
  788. BOOST_ASSERT(u > v);
  789. 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])))
  790. break;
  791. }
  792. #ifdef BOOST_MP_GCD_DEBUG
  793. BOOST_ASSERT(q == UU / VV);
  794. UU %= VV;
  795. UU.swap(VV);
  796. #endif
  797. x[0] = x[1];
  798. x[1] = x[2];
  799. y[0] = y[1];
  800. y[1] = y[2];
  801. }
  802. //
  803. // We get here when the single digit algorithm has gone wrong, back up i, u and v:
  804. //
  805. --i;
  806. u = old_u;
  807. v = old_v;
  808. //
  809. // Now run the full double-digit algorithm:
  810. //
  811. while (true)
  812. {
  813. limb_type q = 1;
  814. double_limb_type tt = u;
  815. divide_subtract(q, u, v);
  816. std::swap(u, v);
  817. tt = y[0] + q * static_cast<double_limb_type>(y[1]);
  818. //
  819. // If calculation of y[2] would overflow a single limb, then we *must* terminate.
  820. // Note that x[2] < y[2] so there is no need to check that as well:
  821. //
  822. if (tt >> bits_per_limb)
  823. {
  824. ++i;
  825. break;
  826. }
  827. x[2] = x[0] + q * x[1];
  828. y[2] = tt;
  829. ++i;
  830. if ((i & 1u) == 0)
  831. {
  832. BOOST_ASSERT(u > v);
  833. if ((v < x[2]) || ((u - v) < (static_cast<double_limb_type>(y[2]) + y[1])))
  834. break;
  835. }
  836. else
  837. {
  838. BOOST_ASSERT(u > v);
  839. if ((v < y[2]) || ((u - v) < (static_cast<double_limb_type>(x[2]) + x[1])))
  840. break;
  841. }
  842. #ifdef BOOST_MP_GCD_DEBUG
  843. BOOST_ASSERT(q == UU / VV);
  844. UU %= VV;
  845. UU.swap(VV);
  846. #endif
  847. x[0] = x[1];
  848. x[1] = x[2];
  849. y[0] = y[1];
  850. y[1] = y[2];
  851. }
  852. if (i == 1)
  853. {
  854. // No change to U and V we've stalled!
  855. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t;
  856. eval_modulus(t, U, V);
  857. U.swap(V);
  858. V.swap(t);
  859. return;
  860. }
  861. //
  862. // Update U and V.
  863. // We have:
  864. //
  865. // U = x[0]U + y[0]V and
  866. // V = x[1]U + y[1]V.
  867. //
  868. // But since we track only absolute values of x and y
  869. // we have to take account of the implied signs and perform
  870. // the appropriate subtraction depending on the whether i is
  871. // even or odd:
  872. //
  873. unsigned ts = U.size() + 1;
  874. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t1(storage, ts), t2(storage, ts), t3(storage, ts);
  875. eval_multiply(t1, U, x[0]);
  876. eval_multiply(t2, V, y[0]);
  877. eval_multiply(t3, U, x[1]);
  878. if ((i & 1u) == 0)
  879. {
  880. if (x[0] == 0)
  881. U = t2;
  882. else
  883. {
  884. BOOST_ASSERT(t2.compare(t1) >= 0);
  885. eval_subtract(U, t2, t1);
  886. BOOST_ASSERT(U.sign() == false);
  887. }
  888. }
  889. else
  890. {
  891. BOOST_ASSERT(t1.compare(t2) >= 0);
  892. eval_subtract(U, t1, t2);
  893. BOOST_ASSERT(U.sign() == false);
  894. }
  895. eval_multiply(t2, V, y[1]);
  896. if (i & 1u)
  897. {
  898. if (x[1] == 0)
  899. V = t2;
  900. else
  901. {
  902. BOOST_ASSERT(t2.compare(t3) >= 0);
  903. eval_subtract(V, t2, t3);
  904. BOOST_ASSERT(V.sign() == false);
  905. }
  906. }
  907. else
  908. {
  909. BOOST_ASSERT(t3.compare(t2) >= 0);
  910. eval_subtract(V, t3, t2);
  911. BOOST_ASSERT(V.sign() == false);
  912. }
  913. BOOST_ASSERT(U.compare(V) >= 0);
  914. BOOST_ASSERT(lu > eval_msb(U));
  915. #ifdef BOOST_MP_GCD_DEBUG
  916. BOOST_ASSERT(UU == U);
  917. BOOST_ASSERT(VV == V);
  918. extern unsigned total_lehmer_gcd_calls;
  919. extern unsigned total_lehmer_gcd_bits_saved;
  920. extern unsigned total_lehmer_gcd_cycles;
  921. ++total_lehmer_gcd_calls;
  922. total_lehmer_gcd_bits_saved += lu - eval_msb(U);
  923. total_lehmer_gcd_cycles += i;
  924. #endif
  925. if (lu < 2048)
  926. {
  927. //
  928. // Since we have stripped all common powers of 2 from U and V at the start
  929. // if either are even at this point, we can remove stray powers of 2 now.
  930. // Note that it is not possible for *both* U and V to be even at this point.
  931. //
  932. // This has an adverse effect on performance for high bit counts, but has
  933. // a significant positive effect for smaller counts.
  934. //
  935. if ((U.limbs()[0] & 1u) == 0)
  936. {
  937. eval_right_shift(U, eval_lsb(U));
  938. if (U.compare(V) < 0)
  939. U.swap(V);
  940. }
  941. else if ((V.limbs()[0] & 1u) == 0)
  942. {
  943. eval_right_shift(V, eval_lsb(V));
  944. }
  945. }
  946. storage.deallocate(ts * 3);
  947. }
  948. #endif
  949. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  950. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
  951. eval_gcd(
  952. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  953. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
  954. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& b)
  955. {
  956. using default_ops::eval_get_sign;
  957. using default_ops::eval_is_zero;
  958. using default_ops::eval_lsb;
  959. if (a.size() == 1)
  960. {
  961. eval_gcd(result, b, *a.limbs());
  962. return;
  963. }
  964. if (b.size() == 1)
  965. {
  966. eval_gcd(result, a, *b.limbs());
  967. return;
  968. }
  969. unsigned temp_size = (std::max)(a.size(), b.size()) + 1;
  970. typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::scoped_shared_storage storage(a, temp_size * 6);
  971. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> U(storage, temp_size);
  972. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> V(storage, temp_size);
  973. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t(storage, temp_size);
  974. U = a;
  975. V = b;
  976. int s = eval_get_sign(U);
  977. /* GCD(0,x) := x */
  978. if (s < 0)
  979. {
  980. U.negate();
  981. }
  982. else if (s == 0)
  983. {
  984. result = V;
  985. return;
  986. }
  987. s = eval_get_sign(V);
  988. if (s < 0)
  989. {
  990. V.negate();
  991. }
  992. else if (s == 0)
  993. {
  994. result = U;
  995. return;
  996. }
  997. //
  998. // Remove common factors of 2:
  999. //
  1000. unsigned us = eval_lsb(U);
  1001. unsigned vs = eval_lsb(V);
  1002. int shift = (std::min)(us, vs);
  1003. if (us)
  1004. eval_right_shift(U, us);
  1005. if (vs)
  1006. eval_right_shift(V, vs);
  1007. if (U.compare(V) < 0)
  1008. U.swap(V);
  1009. while (!eval_is_zero(V))
  1010. {
  1011. if (U.size() <= 2)
  1012. {
  1013. //
  1014. // Special case: if V has no more than 2 limbs
  1015. // then we can reduce U and V to a pair of integers and perform
  1016. // direct integer gcd:
  1017. //
  1018. if (U.size() == 1)
  1019. U = eval_gcd(*V.limbs(), *U.limbs());
  1020. else
  1021. {
  1022. double_limb_type i = U.limbs()[0] | (static_cast<double_limb_type>(U.limbs()[1]) << sizeof(limb_type) * CHAR_BIT);
  1023. 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);
  1024. U = eval_gcd(i, j);
  1025. }
  1026. break;
  1027. }
  1028. unsigned lu = eval_msb(U) + 1;
  1029. unsigned lv = eval_msb(V) + 1;
  1030. #ifndef BOOST_MP_NO_CONSTEXPR_DETECTION
  1031. if (!BOOST_MP_IS_CONST_EVALUATED(lu) && (lu - lv <= bits_per_limb / 2))
  1032. #else
  1033. if (lu - lv <= bits_per_limb / 2)
  1034. #endif
  1035. {
  1036. eval_gcd_lehmer(U, V, lu, storage);
  1037. }
  1038. else
  1039. {
  1040. eval_modulus(t, U, V);
  1041. U.swap(V);
  1042. V.swap(t);
  1043. }
  1044. }
  1045. result = U;
  1046. if (shift)
  1047. eval_left_shift(result, shift);
  1048. }
  1049. //
  1050. // Now again for trivial backends:
  1051. //
  1052. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  1053. 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
  1054. 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
  1055. {
  1056. *result.limbs() = boost::integer::gcd(*a.limbs(), *b.limbs());
  1057. }
  1058. // This one is only enabled for unchecked cpp_int's, for checked int's we need the checking in the default version:
  1059. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  1060. 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
  1061. 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))
  1062. {
  1063. *result.limbs() = boost::integer::lcm(*a.limbs(), *b.limbs());
  1064. result.normalize(); // result may overflow the specified number of bits
  1065. }
  1066. inline void conversion_overflow(const std::integral_constant<int, checked>&)
  1067. {
  1068. BOOST_THROW_EXCEPTION(std::overflow_error("Overflow in conversion to narrower type"));
  1069. }
  1070. inline BOOST_MP_CXX14_CONSTEXPR void conversion_overflow(const std::integral_constant<int, unchecked>&) {}
  1071. #if defined(__clang__) && defined(__MINGW32__)
  1072. //
  1073. // clang-11 on Mingw segfaults on conversion of __int128 -> float.
  1074. // See: https://bugs.llvm.org/show_bug.cgi?id=48941
  1075. // These workarounds pass everything through an intermediate uint64_t.
  1076. //
  1077. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  1078. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
  1079. 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
  1080. eval_convert_to(float* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val)
  1081. {
  1082. float f = static_cast<std::uint64_t>((*val.limbs()) >> 64);
  1083. *result = std::ldexp(f, 64);
  1084. *result += static_cast<std::uint64_t>((*val.limbs()));
  1085. if(val.sign())
  1086. *result = -*result;
  1087. }
  1088. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  1089. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
  1090. 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
  1091. eval_convert_to(double* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val)
  1092. {
  1093. float f = static_cast<std::uint64_t>((*val.limbs()) >> 64);
  1094. *result = std::ldexp(f, 64);
  1095. *result += static_cast<std::uint64_t>((*val.limbs()));
  1096. if(val.sign())
  1097. *result = -*result;
  1098. }
  1099. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  1100. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
  1101. 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
  1102. eval_convert_to(long double* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val)
  1103. {
  1104. float f = static_cast<std::uint64_t>((*val.limbs()) >> 64);
  1105. *result = std::ldexp(f, 64);
  1106. *result += static_cast<std::uint64_t>((*val.limbs()));
  1107. if(val.sign())
  1108. *result = -*result;
  1109. }
  1110. #endif
  1111. template <class R, unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  1112. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
  1113. 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
  1114. eval_convert_to(R* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val)
  1115. {
  1116. using common_type = typename std::common_type<R, typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::local_limb_type>::type;
  1117. BOOST_IF_CONSTEXPR(std::numeric_limits<R>::is_specialized)
  1118. {
  1119. if (static_cast<common_type>(*val.limbs()) > static_cast<common_type>((std::numeric_limits<R>::max)()))
  1120. {
  1121. if (val.isneg())
  1122. {
  1123. 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) > ());
  1124. if (static_cast<common_type>(*val.limbs()) > -static_cast<common_type>((std::numeric_limits<R>::min)()))
  1125. conversion_overflow(typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type());
  1126. *result = (std::numeric_limits<R>::min)();
  1127. }
  1128. else
  1129. {
  1130. conversion_overflow(typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type());
  1131. *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());
  1132. }
  1133. }
  1134. else
  1135. {
  1136. *result = static_cast<R>(*val.limbs());
  1137. if (val.isneg())
  1138. {
  1139. 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) > ());
  1140. *result = negate_integer(*result, std::integral_constant < bool, is_signed_number<R>::value || (number_category<R>::value == number_kind_floating_point) > ());
  1141. }
  1142. }
  1143. }
  1144. else
  1145. {
  1146. *result = static_cast<R>(*val.limbs());
  1147. if (val.isneg())
  1148. {
  1149. 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) > ());
  1150. *result = negate_integer(*result, std::integral_constant<bool, is_signed_number<R>::value || (number_category<R>::value == number_kind_floating_point) > ());
  1151. }
  1152. }
  1153. }
  1154. template <class R, unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  1155. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
  1156. 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
  1157. eval_convert_to(R* result, const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val)
  1158. {
  1159. using common_type = typename std::common_type<R, typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::local_limb_type>::type;
  1160. BOOST_IF_CONSTEXPR(std::numeric_limits<R>::is_specialized)
  1161. {
  1162. if(static_cast<common_type>(*val.limbs()) > static_cast<common_type>((std::numeric_limits<R>::max)()))
  1163. {
  1164. conversion_overflow(typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type());
  1165. *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());
  1166. }
  1167. else
  1168. *result = static_cast<R>(*val.limbs());
  1169. }
  1170. else
  1171. *result = static_cast<R>(*val.limbs());
  1172. }
  1173. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  1174. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, unsigned>::type
  1175. eval_lsb(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
  1176. {
  1177. using default_ops::eval_get_sign;
  1178. if (eval_get_sign(a) == 0)
  1179. {
  1180. BOOST_THROW_EXCEPTION(std::domain_error("No bits were set in the operand."));
  1181. }
  1182. if (a.sign())
  1183. {
  1184. BOOST_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined."));
  1185. }
  1186. //
  1187. // Find the index of the least significant bit within that limb:
  1188. //
  1189. return boost::multiprecision::detail::find_lsb(*a.limbs());
  1190. }
  1191. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  1192. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, unsigned>::type
  1193. eval_msb_imp(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
  1194. {
  1195. //
  1196. // Find the index of the least significant bit within that limb:
  1197. //
  1198. return boost::multiprecision::detail::find_msb(*a.limbs());
  1199. }
  1200. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  1201. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value, unsigned>::type
  1202. eval_msb(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a)
  1203. {
  1204. using default_ops::eval_get_sign;
  1205. if (eval_get_sign(a) == 0)
  1206. {
  1207. BOOST_THROW_EXCEPTION(std::domain_error("No bits were set in the operand."));
  1208. }
  1209. if (a.sign())
  1210. {
  1211. BOOST_THROW_EXCEPTION(std::domain_error("Testing individual bits in negative values is not supported - results are undefined."));
  1212. }
  1213. return eval_msb_imp(a);
  1214. }
  1215. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  1216. inline BOOST_MP_CXX14_CONSTEXPR std::size_t hash_value(const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& val) noexcept
  1217. {
  1218. std::size_t result = 0;
  1219. for (unsigned i = 0; i < val.size(); ++i)
  1220. {
  1221. boost::hash_combine(result, val.limbs()[i]);
  1222. }
  1223. boost::hash_combine(result, val.sign());
  1224. return result;
  1225. }
  1226. #ifdef BOOST_MSVC
  1227. #pragma warning(pop)
  1228. #endif
  1229. }}} // namespace boost::multiprecision::backends
  1230. #endif