multiply.hpp 45 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848
  1. ///////////////////////////////////////////////////////////////
  2. // Copyright 2012-20 John Maddock.
  3. // Copyright 2019-20 Christopher Kormanyos.
  4. // Copyright 2019-20 Madhur Chauhan.
  5. // Distributed under the Boost Software License, Version 1.0.
  6. // (See accompanying file LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
  7. //
  8. // Comparison operators for cpp_int_backend:
  9. //
  10. #ifndef BOOST_MP_CPP_INT_MULTIPLY_HPP
  11. #define BOOST_MP_CPP_INT_MULTIPLY_HPP
  12. #include <limits>
  13. #include <boost/multiprecision/detail/standalone_config.hpp>
  14. #include <boost/multiprecision/detail/endian.hpp>
  15. #include <boost/multiprecision/detail/assert.hpp>
  16. #include <boost/multiprecision/integer.hpp>
  17. namespace boost { namespace multiprecision { namespace backends {
  18. #ifdef BOOST_MSVC
  19. #pragma warning(push)
  20. #pragma warning(disable : 4127) // conditional expression is constant
  21. #endif
  22. //
  23. // Multiplication by a single limb:
  24. //
  25. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, std::size_t MinBits2, std::size_t MaxBits2, cpp_integer_type SignType2, cpp_int_check_type Checked2, class Allocator2>
  26. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && !is_trivial_cpp_int<cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2> >::value>::type
  27. eval_multiply(
  28. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  29. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
  30. const limb_type& val) noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
  31. {
  32. if (!val)
  33. {
  34. result = static_cast<limb_type>(0);
  35. return;
  36. }
  37. if ((void*)&a != (void*)&result)
  38. result.resize(a.size(), a.size());
  39. double_limb_type carry = 0;
  40. typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_pointer p = result.limbs();
  41. typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_pointer pe = result.limbs() + result.size();
  42. typename cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>::const_limb_pointer pa = a.limbs();
  43. while (p != pe)
  44. {
  45. carry += static_cast<double_limb_type>(*pa) * static_cast<double_limb_type>(val);
  46. #ifdef __MSVC_RUNTIME_CHECKS
  47. *p = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  48. #else
  49. *p = static_cast<limb_type>(carry);
  50. #endif
  51. carry >>= cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
  52. ++p, ++pa;
  53. }
  54. if (carry)
  55. {
  56. std::size_t i = result.size();
  57. result.resize(i + 1, i + 1);
  58. if (result.size() > i)
  59. result.limbs()[i] = static_cast<limb_type>(carry);
  60. }
  61. result.sign(a.sign());
  62. if (is_fixed_precision<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value)
  63. result.normalize();
  64. }
  65. //
  66. // resize_for_carry forces a resize of the underlying buffer only if a previous request
  67. // for "required" elements could possibly have failed, *and* we have checking enabled.
  68. // This will cause an overflow error inside resize():
  69. //
  70. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  71. inline BOOST_MP_CXX14_CONSTEXPR void resize_for_carry(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& /*result*/, std::size_t /*required*/) {}
  72. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, class Allocator1>
  73. inline BOOST_MP_CXX14_CONSTEXPR void resize_for_carry(cpp_int_backend<MinBits1, MaxBits1, SignType1, checked, Allocator1>& result, std::size_t required)
  74. {
  75. if (result.size() < required)
  76. result.resize(required, required);
  77. }
  78. //
  79. // Minimum number of limbs required for Karatsuba to be worthwhile:
  80. //
  81. #ifdef BOOST_MP_KARATSUBA_CUTOFF
  82. const size_t karatsuba_cutoff = BOOST_MP_KARATSUBA_CUTOFF;
  83. #else
  84. const size_t karatsuba_cutoff = 40;
  85. #endif
  86. //
  87. // Core (recursive) Karatsuba multiplication, all the storage required is allocated upfront and
  88. // passed down the stack in this routine. Note that all the cpp_int_backend's must be the same type
  89. // and full variable precision. Karatsuba really doesn't play nice with fixed-size integers. If necessary
  90. // fixed precision integers will get aliased as variable-precision types before this is called.
  91. //
  92. template <std::size_t MinBits, std::size_t MaxBits, cpp_int_check_type Checked, class Allocator>
  93. inline void multiply_karatsuba(
  94. cpp_int_backend<MinBits, MaxBits, signed_magnitude, Checked, Allocator>& result,
  95. const cpp_int_backend<MinBits, MaxBits, signed_magnitude, Checked, Allocator>& a,
  96. const cpp_int_backend<MinBits, MaxBits, signed_magnitude, Checked, Allocator>& b,
  97. typename cpp_int_backend<MinBits, MaxBits, signed_magnitude, Checked, Allocator>::scoped_shared_storage& storage)
  98. {
  99. using cpp_int_type = cpp_int_backend<MinBits, MaxBits, signed_magnitude, Checked, Allocator>;
  100. std::size_t as = a.size();
  101. std::size_t bs = b.size();
  102. //
  103. // Termination condition: if either argument is smaller than karatsuba_cutoff
  104. // then schoolboy multiplication will be faster:
  105. //
  106. if ((as < karatsuba_cutoff) || (bs < karatsuba_cutoff))
  107. {
  108. eval_multiply(result, a, b);
  109. return;
  110. }
  111. //
  112. // Partitioning size: split the larger of a and b into 2 halves
  113. //
  114. std::size_t n = (as > bs ? as : bs) / 2 + 1;
  115. //
  116. // Partition a and b into high and low parts.
  117. // ie write a, b as a = a_h * 2^n + a_l, b = b_h * 2^n + b_l
  118. //
  119. // We could copy the high and low parts into new variables, but we'll
  120. // use aliasing to reference the internal limbs of a and b. There is one wart here:
  121. // if a and b are mismatched in size, then n may be larger than the smaller
  122. // of a and b. In that situation the high part is zero, and we have no limbs
  123. // to alias, so instead alias a local variable.
  124. // This raises 2 questions:
  125. // * Is this the best way to partition a and b?
  126. // * Since we have one high part zero, the arithmetic simplifies considerably,
  127. // so should we have a special routine for this?
  128. //
  129. std::size_t sz = (std::min)(as, n);
  130. const cpp_int_type a_l(a.limbs(), 0, sz);
  131. sz = (std::min)(bs, n);
  132. const cpp_int_type b_l(b.limbs(), 0, sz);
  133. limb_type zero = 0;
  134. const cpp_int_type a_h(as > n ? a.limbs() + n : &zero, 0, as > n ? as - n : 1);
  135. const cpp_int_type b_h(bs > n ? b.limbs() + n : &zero, 0, bs > n ? bs - n : 1);
  136. //
  137. // The basis for the Karatsuba algorithm is as follows:
  138. //
  139. // let x = a_h * b_ h
  140. // y = a_l * b_l
  141. // z = (a_h + a_l)*(b_h + b_l) - x - y
  142. // and therefore a * b = x * (2 ^ (2 * n))+ z * (2 ^ n) + y
  143. //
  144. // Begin by allocating our temporaries, these alias the memory already allocated in the shared storage:
  145. //
  146. cpp_int_type t1(storage, 2 * n + 2);
  147. cpp_int_type t2(storage, n + 1);
  148. cpp_int_type t3(storage, n + 1);
  149. //
  150. // Now we want:
  151. //
  152. // result = | a_h*b_h | a_l*b_l |
  153. // (bits) <-- 2*n -->
  154. //
  155. // We create aliases for the low and high parts of result, and multiply directly into them:
  156. //
  157. cpp_int_type result_low(result.limbs(), 0, 2 * n);
  158. cpp_int_type result_high(result.limbs(), 2 * n, result.size() - 2 * n);
  159. //
  160. // low part of result is a_l * b_l:
  161. //
  162. multiply_karatsuba(result_low, a_l, b_l, storage);
  163. //
  164. // We haven't zeroed out memory in result, so set to zero any unused limbs,
  165. // if a_l and b_l have mostly random bits then nothing happens here, but if
  166. // one is zero or nearly so, then a memset might be faster... it's not clear
  167. // that it's worth the extra logic though (and is darn hard to measure
  168. // what the "average" case is).
  169. //
  170. for (std::size_t i = result_low.size(); i < 2 * n; ++i)
  171. result.limbs()[i] = 0;
  172. //
  173. // Set the high part of result to a_h * b_h:
  174. //
  175. multiply_karatsuba(result_high, a_h, b_h, storage);
  176. for (std::size_t i = result_high.size() + 2 * n; i < result.size(); ++i)
  177. result.limbs()[i] = 0;
  178. //
  179. // Now calculate (a_h+a_l)*(b_h+b_l):
  180. //
  181. add_unsigned(t2, a_l, a_h);
  182. add_unsigned(t3, b_l, b_h);
  183. multiply_karatsuba(t1, t2, t3, storage); // t1 = (a_h+a_l)*(b_h+b_l)
  184. //
  185. // There is now a slight deviation from Karatsuba, we want to subtract
  186. // a_l*b_l + a_h*b_h from t1, but rather than use an addition and a subtraction
  187. // plus one temporary, we'll use 2 subtractions. On the minus side, a subtraction
  188. // is on average slightly slower than an addition, but we save a temporary (ie memory)
  189. // and also hammer the same piece of memory over and over rather than 2 disparate
  190. // memory regions. Overall it seems to be a slight win.
  191. //
  192. subtract_unsigned(t1, t1, result_high);
  193. subtract_unsigned(t1, t1, result_low);
  194. //
  195. // The final step is to left shift t1 by n bits and add to the result.
  196. // Rather than do an actual left shift, we can simply alias the result
  197. // and add to the alias:
  198. //
  199. cpp_int_type result_alias(result.limbs(), n, result.size() - n);
  200. add_unsigned(result_alias, result_alias, t1);
  201. //
  202. // Free up storage for use by sister branches to this one:
  203. //
  204. storage.deallocate(t1.capacity() + t2.capacity() + t3.capacity());
  205. result.normalize();
  206. }
  207. inline std::size_t karatsuba_storage_size(std::size_t s)
  208. {
  209. //
  210. // This estimates how much memory we will need based on
  211. // s-limb multiplication. In an ideal world the number of limbs
  212. // would halve with each recursion, and our storage requirements
  213. // would be 4s in the limit, and rather less in practice since
  214. // we bail out long before we reach one limb. In the real world
  215. // we don't quite halve s in each recursion, so this is an heuristic
  216. // which over-estimates how much we need. We could compute an exact
  217. // value, but it would be rather time consuming.
  218. //
  219. return 5 * s;
  220. }
  221. //
  222. // There are 2 entry point routines for Karatsuba multiplication:
  223. // one for variable precision types, and one for fixed precision types.
  224. // These are responsible for allocating all the storage required for the recursive
  225. // routines above, and are always at the outermost level.
  226. //
  227. // Normal variable precision case comes first:
  228. //
  229. template <std::size_t MinBits, std::size_t MaxBits, cpp_integer_type SignType, cpp_int_check_type Checked, class Allocator>
  230. inline typename std::enable_if<!is_fixed_precision<cpp_int_backend<MinBits, MaxBits, SignType, Checked, Allocator> >::value>::type
  231. setup_karatsuba(
  232. cpp_int_backend<MinBits, MaxBits, SignType, Checked, Allocator>& result,
  233. const cpp_int_backend<MinBits, MaxBits, SignType, Checked, Allocator>& a,
  234. const cpp_int_backend<MinBits, MaxBits, SignType, Checked, Allocator>& b)
  235. {
  236. std::size_t as = a.size();
  237. std::size_t bs = b.size();
  238. std::size_t s = as > bs ? as : bs;
  239. std::size_t storage_size = karatsuba_storage_size(s);
  240. if (storage_size < 300)
  241. {
  242. //
  243. // Special case: if we don't need too much memory, we can use stack based storage
  244. // and save a call to the allocator, this allows us to use Karatsuba multiply
  245. // at lower limb counts than would otherwise be possible:
  246. //
  247. limb_type limbs[300];
  248. typename cpp_int_backend<MinBits, MaxBits, SignType, Checked, Allocator>::scoped_shared_storage storage(limbs, storage_size);
  249. multiply_karatsuba(result, a, b, storage);
  250. }
  251. else
  252. {
  253. typename cpp_int_backend<MinBits, MaxBits, SignType, Checked, Allocator>::scoped_shared_storage storage(result.allocator(), storage_size);
  254. multiply_karatsuba(result, a, b, storage);
  255. }
  256. }
  257. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, std::size_t MinBits2, std::size_t MaxBits2, cpp_integer_type SignType2, cpp_int_check_type Checked2, class Allocator2, std::size_t MinBits3, std::size_t MaxBits3, cpp_integer_type SignType3, cpp_int_check_type Checked3, class Allocator3>
  258. inline typename std::enable_if<is_fixed_precision<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value || is_fixed_precision<cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2> >::value || is_fixed_precision<cpp_int_backend<MinBits3, MaxBits3, SignType3, Checked3, Allocator3> >::value>::type
  259. setup_karatsuba(
  260. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  261. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
  262. const cpp_int_backend<MinBits3, MaxBits3, SignType3, Checked3, Allocator3>& b)
  263. {
  264. //
  265. // Now comes the fixed precision case.
  266. // In fact Karatsuba doesn't really work with fixed precision since the logic
  267. // requires that we calculate all the bits of the result (especially in the
  268. // temporaries used internally). So... we'll convert all the arguments
  269. // to variable precision types by aliasing them, this also
  270. // reduce the number of template instantations:
  271. //
  272. using variable_precision_type = cpp_int_backend<0, 0, signed_magnitude, unchecked, std::allocator<limb_type> >;
  273. variable_precision_type a_t(a.limbs(), 0, a.size()), b_t(b.limbs(), 0, b.size());
  274. std::size_t as = a.size();
  275. std::size_t bs = b.size();
  276. std::size_t s = as > bs ? as : bs;
  277. std::size_t sz = as + bs;
  278. std::size_t storage_size = karatsuba_storage_size(s);
  279. if (!is_fixed_precision<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value || (sz * sizeof(limb_type) * CHAR_BIT <= MaxBits1))
  280. {
  281. // Result is large enough for all the bits of the result, so we can use aliasing:
  282. result.resize(sz, sz);
  283. variable_precision_type t(result.limbs(), 0, result.size());
  284. typename variable_precision_type::scoped_shared_storage storage(t.allocator(), storage_size);
  285. multiply_karatsuba(t, a_t, b_t, storage);
  286. result.resize(t.size(), t.size());
  287. }
  288. else
  289. {
  290. //
  291. // Not enough bit in result for the answer, so we must use a temporary
  292. // and then truncate (ie modular arithmetic):
  293. //
  294. typename variable_precision_type::scoped_shared_storage storage(variable_precision_type::allocator_type(), sz + storage_size);
  295. variable_precision_type t(storage, sz);
  296. multiply_karatsuba(t, a_t, b_t, storage);
  297. //
  298. // If there is truncation, and result is a checked type then this will throw:
  299. //
  300. result = t;
  301. }
  302. }
  303. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, std::size_t MinBits2, std::size_t MaxBits2, cpp_integer_type SignType2, cpp_int_check_type Checked2, class Allocator2, std::size_t MinBits3, std::size_t MaxBits3, cpp_integer_type SignType3, cpp_int_check_type Checked3, class Allocator3>
  304. inline typename std::enable_if<!is_fixed_precision<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && !is_fixed_precision<cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2> >::value && !is_fixed_precision<cpp_int_backend<MinBits3, MaxBits3, SignType3, Checked3, Allocator3> >::value>::type
  305. setup_karatsuba(
  306. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  307. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
  308. const cpp_int_backend<MinBits3, MaxBits3, SignType3, Checked3, Allocator3>& b)
  309. {
  310. //
  311. // Variable precision, mixed arguments, just alias and forward:
  312. //
  313. using variable_precision_type = cpp_int_backend<0, 0, signed_magnitude, unchecked, std::allocator<limb_type> >;
  314. variable_precision_type a_t(a.limbs(), 0, a.size()), b_t(b.limbs(), 0, b.size());
  315. std::size_t as = a.size();
  316. std::size_t bs = b.size();
  317. std::size_t s = as > bs ? as : bs;
  318. std::size_t sz = as + bs;
  319. std::size_t storage_size = karatsuba_storage_size(s);
  320. result.resize(sz, sz);
  321. variable_precision_type t(result.limbs(), 0, result.size());
  322. typename variable_precision_type::scoped_shared_storage storage(t.allocator(), storage_size);
  323. multiply_karatsuba(t, a_t, b_t, storage);
  324. }
  325. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, std::size_t MinBits2, std::size_t MaxBits2, cpp_integer_type SignType2, cpp_int_check_type Checked2, class Allocator2, std::size_t MinBits3, std::size_t MaxBits3, cpp_integer_type SignType3, cpp_int_check_type Checked3, class Allocator3>
  326. inline BOOST_MP_CXX14_CONSTEXPR void
  327. eval_multiply_comba(
  328. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  329. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
  330. const cpp_int_backend<MinBits3, MaxBits3, SignType3, Checked3, Allocator3>& b) noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
  331. {
  332. //
  333. // see PR #182
  334. // Comba Multiplier - based on Paul Comba's
  335. // Exponentiation cryptosystems on the IBM PC, 1990
  336. //
  337. std::ptrdiff_t as = a.size(),
  338. bs = b.size(),
  339. rs = result.size();
  340. typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_pointer pr = result.limbs();
  341. double_limb_type carry = 0,
  342. temp = 0;
  343. limb_type overflow = 0;
  344. const std::size_t limb_bits = sizeof(limb_type) * CHAR_BIT;
  345. const bool must_throw = rs < as + bs - 1;
  346. for (std::ptrdiff_t r = 0, lim = (std::min)(rs, as + bs - 1); r < lim; ++r, overflow = 0)
  347. {
  348. std::ptrdiff_t i = r >= as ? as - 1 : r,
  349. j = r - i,
  350. k = i < bs - j ? i + 1 : bs - j; // min(i+1, bs-j);
  351. typename cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>::const_limb_pointer pa = a.limbs() + i;
  352. typename cpp_int_backend<MinBits3, MaxBits3, SignType3, Checked3, Allocator3>::const_limb_pointer pb = b.limbs() + j;
  353. temp = carry;
  354. carry += static_cast<double_limb_type>(*(pa)) * (*(pb));
  355. overflow += carry < temp;
  356. for (--k; k; k--)
  357. {
  358. temp = carry;
  359. carry += static_cast<double_limb_type>(*(--pa)) * (*(++pb));
  360. overflow += carry < temp;
  361. }
  362. *(pr++) = static_cast<limb_type>(carry);
  363. carry = (static_cast<double_limb_type>(overflow) << limb_bits) | (carry >> limb_bits);
  364. }
  365. if (carry || must_throw)
  366. {
  367. resize_for_carry(result, as + bs);
  368. if (static_cast<int>(result.size()) >= as + bs)
  369. *pr = static_cast<limb_type>(carry);
  370. }
  371. }
  372. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, std::size_t MinBits2, std::size_t MaxBits2, cpp_integer_type SignType2, cpp_int_check_type Checked2, class Allocator2, std::size_t MinBits3, std::size_t MaxBits3, cpp_integer_type SignType3, cpp_int_check_type Checked3, class Allocator3>
  373. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && !is_trivial_cpp_int<cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2> >::value && !is_trivial_cpp_int<cpp_int_backend<MinBits3, MaxBits3, SignType3, Checked3, Allocator3> >::value>::type
  374. eval_multiply(
  375. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  376. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
  377. const cpp_int_backend<MinBits3, MaxBits3, SignType3, Checked3, Allocator3>& b)
  378. noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value
  379. && (karatsuba_cutoff * sizeof(limb_type) * CHAR_BIT > MaxBits1)
  380. && (karatsuba_cutoff * sizeof(limb_type)* CHAR_BIT > MaxBits2)
  381. && (karatsuba_cutoff * sizeof(limb_type)* CHAR_BIT > MaxBits3)))
  382. {
  383. // Uses simple (O(n^2)) multiplication when the limbs are less
  384. // otherwise switches to karatsuba algorithm based on experimental value (~40 limbs)
  385. //
  386. // Trivial cases first:
  387. //
  388. std::size_t as = a.size();
  389. std::size_t bs = b.size();
  390. typename cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>::const_limb_pointer pa = a.limbs();
  391. typename cpp_int_backend<MinBits3, MaxBits3, SignType3, Checked3, Allocator3>::const_limb_pointer pb = b.limbs();
  392. if (as == 1)
  393. {
  394. bool s = b.sign() != a.sign();
  395. if (bs == 1)
  396. {
  397. result = static_cast<double_limb_type>(*pa) * static_cast<double_limb_type>(*pb);
  398. }
  399. else
  400. {
  401. limb_type l = *pa;
  402. eval_multiply(result, b, l);
  403. }
  404. result.sign(s);
  405. return;
  406. }
  407. if (bs == 1)
  408. {
  409. bool s = b.sign() != a.sign();
  410. limb_type l = *pb;
  411. eval_multiply(result, a, l);
  412. result.sign(s);
  413. return;
  414. }
  415. if ((void*)&result == (void*)&a)
  416. {
  417. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t(a);
  418. eval_multiply(result, t, b);
  419. return;
  420. }
  421. if ((void*)&result == (void*)&b)
  422. {
  423. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t(b);
  424. eval_multiply(result, a, t);
  425. return;
  426. }
  427. constexpr double_limb_type limb_max = static_cast<double_limb_type>(~static_cast<limb_type>(0u));
  428. constexpr double_limb_type double_limb_max = static_cast<double_limb_type>(~static_cast<double_limb_type>(0u));
  429. result.resize(as + bs, as + bs - 1);
  430. #ifndef BOOST_MP_NO_CONSTEXPR_DETECTION
  431. if (!BOOST_MP_IS_CONST_EVALUATED(as) && (as >= karatsuba_cutoff && bs >= karatsuba_cutoff))
  432. #else
  433. if (as >= karatsuba_cutoff && bs >= karatsuba_cutoff)
  434. #endif
  435. {
  436. setup_karatsuba(result, a, b);
  437. //
  438. // Set the sign of the result:
  439. //
  440. result.sign(a.sign() != b.sign());
  441. return;
  442. }
  443. typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_pointer pr = result.limbs();
  444. static_assert(double_limb_max - 2 * limb_max >= limb_max * limb_max, "failed limb size sanity check");
  445. #ifndef BOOST_MP_NO_CONSTEXPR_DETECTION
  446. if (BOOST_MP_IS_CONST_EVALUATED(as))
  447. {
  448. for (std::size_t i = 0; i < result.size(); ++i)
  449. pr[i] = 0;
  450. }
  451. else
  452. #endif
  453. std::memset(pr, 0, result.size() * sizeof(limb_type));
  454. #if defined(BOOST_MP_COMBA)
  455. //
  456. // Comba Multiplier might not be efficient because of less efficient assembly
  457. // by the compiler as of 09/01/2020 (DD/MM/YY). See PR #182
  458. // Till then this will lay dormant :(
  459. //
  460. eval_multiply_comba(result, a, b);
  461. #else
  462. double_limb_type carry = 0;
  463. for (std::size_t i = 0; i < as; ++i)
  464. {
  465. BOOST_MP_ASSERT(result.size() > i);
  466. std::size_t inner_limit = !is_fixed_precision<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value ? bs : (std::min)(result.size() - i, bs);
  467. std::size_t j = 0;
  468. for (; j < inner_limit; ++j)
  469. {
  470. BOOST_MP_ASSERT(i + j < result.size());
  471. #if (!defined(__GLIBCXX__) && !defined(__GLIBCPP__)) || !BOOST_WORKAROUND(BOOST_GCC_VERSION, <= 50100)
  472. BOOST_MP_ASSERT(!std::numeric_limits<double_limb_type>::is_specialized || ((std::numeric_limits<double_limb_type>::max)() - carry >
  473. static_cast<double_limb_type>(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::max_limb_value) * static_cast<double_limb_type>(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::max_limb_value)));
  474. #endif
  475. carry += static_cast<double_limb_type>(pa[i]) * static_cast<double_limb_type>(pb[j]);
  476. BOOST_MP_ASSERT(!std::numeric_limits<double_limb_type>::is_specialized || ((std::numeric_limits<double_limb_type>::max)() - carry >= pr[i + j]));
  477. carry += pr[i + j];
  478. #ifdef __MSVC_RUNTIME_CHECKS
  479. pr[i + j] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  480. #else
  481. pr[i + j] = static_cast<limb_type>(carry);
  482. #endif
  483. carry >>= cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
  484. BOOST_MP_ASSERT(carry <= (cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::max_limb_value));
  485. }
  486. if (carry)
  487. {
  488. resize_for_carry(result, i + j + 1); // May throw if checking is enabled
  489. if (i + j < result.size())
  490. #ifdef __MSVC_RUNTIME_CHECKS
  491. pr[i + j] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  492. #else
  493. pr[i + j] = static_cast<limb_type>(carry);
  494. #endif
  495. }
  496. carry = 0;
  497. }
  498. #endif // ifdef(BOOST_MP_COMBA) ends
  499. result.normalize();
  500. //
  501. // Set the sign of the result:
  502. //
  503. result.sign(a.sign() != b.sign());
  504. }
  505. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, std::size_t MinBits2, std::size_t MaxBits2, cpp_integer_type SignType2, cpp_int_check_type Checked2, class Allocator2>
  506. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && !is_trivial_cpp_int<cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2> >::value>::type
  507. eval_multiply(
  508. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  509. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a)
  510. noexcept((noexcept(eval_multiply(std::declval<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&>(), std::declval<const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&>(), std::declval<const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>&>()))))
  511. {
  512. eval_multiply(result, result, a);
  513. }
  514. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  515. 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
  516. eval_multiply(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, const limb_type& val)
  517. noexcept((noexcept(eval_multiply(std::declval<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&>(), std::declval<const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&>(), std::declval<const limb_type&>()))))
  518. {
  519. eval_multiply(result, result, val);
  520. }
  521. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, std::size_t MinBits2, std::size_t MaxBits2, cpp_integer_type SignType2, cpp_int_check_type Checked2, class Allocator2>
  522. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && !is_trivial_cpp_int<cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2> >::value>::type
  523. eval_multiply(
  524. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  525. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
  526. const double_limb_type& val)
  527. noexcept(
  528. (noexcept(eval_multiply(std::declval<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&>(), std::declval<const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>&>(), std::declval<const limb_type&>())))
  529. && (noexcept(eval_multiply(std::declval<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&>(), std::declval<const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>&>(), std::declval<const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&>())))
  530. )
  531. {
  532. if (val <= (std::numeric_limits<limb_type>::max)())
  533. {
  534. eval_multiply(result, a, static_cast<limb_type>(val));
  535. }
  536. else
  537. {
  538. #if BOOST_MP_ENDIAN_LITTLE_BYTE && !defined(BOOST_MP_TEST_NO_LE)
  539. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t(val);
  540. #else
  541. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t;
  542. t = val;
  543. #endif
  544. eval_multiply(result, a, t);
  545. }
  546. }
  547. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  548. 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
  549. eval_multiply(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, const double_limb_type& val)
  550. noexcept((noexcept(eval_multiply(std::declval<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&>(), std::declval<const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&>(), std::declval<const double_limb_type&>()))))
  551. {
  552. eval_multiply(result, result, val);
  553. }
  554. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, std::size_t MinBits2, std::size_t MaxBits2, cpp_integer_type SignType2, cpp_int_check_type Checked2, class Allocator2>
  555. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && !is_trivial_cpp_int<cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2> >::value>::type
  556. eval_multiply(
  557. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  558. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
  559. const signed_limb_type& val)
  560. noexcept((noexcept(eval_multiply(std::declval<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&>(), std::declval<const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>&>(), std::declval<const limb_type&>()))))
  561. {
  562. if (val > 0)
  563. eval_multiply(result, a, static_cast<limb_type>(val));
  564. else
  565. {
  566. eval_multiply(result, a, static_cast<limb_type>(boost::multiprecision::detail::unsigned_abs(val)));
  567. result.negate();
  568. }
  569. }
  570. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  571. 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
  572. eval_multiply(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, const signed_limb_type& val)
  573. noexcept((noexcept(eval_multiply(std::declval<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&>(), std::declval<const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&>(), std::declval<const limb_type&>()))))
  574. {
  575. eval_multiply(result, result, val);
  576. }
  577. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, std::size_t MinBits2, std::size_t MaxBits2, cpp_integer_type SignType2, cpp_int_check_type Checked2, class Allocator2>
  578. inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && !is_trivial_cpp_int<cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2> >::value>::type
  579. eval_multiply(
  580. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  581. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
  582. const signed_double_limb_type& val)
  583. noexcept(
  584. (noexcept(eval_multiply(std::declval<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&>(), std::declval<const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>&>(), std::declval<const limb_type&>())))
  585. && (noexcept(eval_multiply(std::declval<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&>(), std::declval<const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>&>(), std::declval<const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&>())))
  586. )
  587. {
  588. if (val > 0)
  589. {
  590. if (val <= (std::numeric_limits<limb_type>::max)())
  591. {
  592. eval_multiply(result, a, static_cast<limb_type>(val));
  593. return;
  594. }
  595. }
  596. else if (val >= -static_cast<signed_double_limb_type>((std::numeric_limits<limb_type>::max)()))
  597. {
  598. eval_multiply(result, a, static_cast<limb_type>(boost::multiprecision::detail::unsigned_abs(val)));
  599. result.negate();
  600. return;
  601. }
  602. #if BOOST_MP_ENDIAN_LITTLE_BYTE && !defined(BOOST_MP_TEST_NO_LE)
  603. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t(val);
  604. #else
  605. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t;
  606. t = val;
  607. #endif
  608. eval_multiply(result, a, t);
  609. }
  610. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  611. 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
  612. eval_multiply(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, const signed_double_limb_type& val)
  613. noexcept(
  614. (noexcept(eval_multiply(std::declval<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&>(), std::declval<const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&>(), std::declval<const limb_type&>())))
  615. && (noexcept(eval_multiply(std::declval<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&>(), std::declval<const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&>(), std::declval<const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>&>())))
  616. )
  617. {
  618. eval_multiply(result, result, val);
  619. }
  620. //
  621. // Now over again for trivial cpp_int's:
  622. //
  623. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  624. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
  625. is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && 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 || is_signed_number<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value)>::type
  626. eval_multiply(
  627. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  628. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& o) noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
  629. {
  630. *result.limbs() = detail::checked_multiply(*result.limbs(), *o.limbs(), typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type());
  631. result.sign(result.sign() != o.sign());
  632. result.normalize();
  633. }
  634. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  635. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
  636. 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>::type
  637. eval_multiply(
  638. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  639. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& o) noexcept((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
  640. {
  641. *result.limbs() = detail::checked_multiply(*result.limbs(), *o.limbs(), typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type());
  642. result.normalize();
  643. }
  644. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  645. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
  646. is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && 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 || is_signed_number<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value)>::type
  647. eval_multiply(
  648. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  649. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
  650. 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))
  651. {
  652. *result.limbs() = detail::checked_multiply(*a.limbs(), *b.limbs(), typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type());
  653. result.sign(a.sign() != b.sign());
  654. result.normalize();
  655. }
  656. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  657. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
  658. 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>::type
  659. eval_multiply(
  660. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  661. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
  662. 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))
  663. {
  664. *result.limbs() = detail::checked_multiply(*a.limbs(), *b.limbs(), typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type());
  665. result.normalize();
  666. }
  667. //
  668. // Special routines for multiplying two integers to obtain a multiprecision result:
  669. //
  670. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  671. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
  672. !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
  673. eval_multiply(
  674. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  675. signed_double_limb_type a, signed_double_limb_type b)
  676. {
  677. constexpr signed_double_limb_type mask = static_cast<signed_double_limb_type>(~static_cast<limb_type>(0));
  678. constexpr std::size_t limb_bits = static_cast<std::size_t>(sizeof(limb_type) * CHAR_BIT);
  679. bool s = false;
  680. if (a < 0)
  681. {
  682. a = -a;
  683. s = true;
  684. }
  685. if (b < 0)
  686. {
  687. b = -b;
  688. s = !s;
  689. }
  690. double_limb_type w = a & mask;
  691. double_limb_type x = static_cast<double_limb_type>(a >> limb_bits);
  692. double_limb_type y = b & mask;
  693. double_limb_type z = static_cast<double_limb_type>(b >> limb_bits);
  694. result.resize(4, 4);
  695. limb_type* pr = result.limbs();
  696. double_limb_type carry = w * y;
  697. #ifdef __MSVC_RUNTIME_CHECKS
  698. pr[0] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  699. carry >>= limb_bits;
  700. carry += w * z + x * y;
  701. pr[1] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  702. carry >>= limb_bits;
  703. carry += x * z;
  704. pr[2] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  705. pr[3] = static_cast<limb_type>(carry >> limb_bits);
  706. #else
  707. pr[0] = static_cast<limb_type>(carry);
  708. carry >>= limb_bits;
  709. carry += w * z + x * y;
  710. pr[1] = static_cast<limb_type>(carry);
  711. carry >>= limb_bits;
  712. carry += x * z;
  713. pr[2] = static_cast<limb_type>(carry);
  714. pr[3] = static_cast<limb_type>(carry >> limb_bits);
  715. #endif
  716. result.sign(s);
  717. result.normalize();
  718. }
  719. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  720. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
  721. !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
  722. eval_multiply(
  723. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  724. double_limb_type a, double_limb_type b)
  725. {
  726. constexpr signed_double_limb_type mask = static_cast<signed_double_limb_type>(~static_cast<limb_type>(0));
  727. constexpr std::size_t limb_bits = static_cast<std::size_t>(sizeof(limb_type) * CHAR_BIT);
  728. double_limb_type w = a & mask;
  729. double_limb_type x = a >> limb_bits;
  730. double_limb_type y = b & mask;
  731. double_limb_type z = b >> limb_bits;
  732. result.resize(4, 4);
  733. limb_type* pr = result.limbs();
  734. double_limb_type carry = w * y;
  735. #ifdef __MSVC_RUNTIME_CHECKS
  736. pr[0] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  737. carry >>= limb_bits;
  738. carry += w * z;
  739. pr[1] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  740. carry >>= limb_bits;
  741. pr[2] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  742. carry = x * y + pr[1];
  743. pr[1] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  744. carry >>= limb_bits;
  745. carry += pr[2] + x * z;
  746. pr[2] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  747. pr[3] = static_cast<limb_type>(carry >> limb_bits);
  748. #else
  749. pr[0] = static_cast<limb_type>(carry);
  750. carry >>= limb_bits;
  751. carry += w * z;
  752. pr[1] = static_cast<limb_type>(carry);
  753. carry >>= limb_bits;
  754. pr[2] = static_cast<limb_type>(carry);
  755. carry = x * y + pr[1];
  756. pr[1] = static_cast<limb_type>(carry);
  757. carry >>= limb_bits;
  758. carry += pr[2] + x * z;
  759. pr[2] = static_cast<limb_type>(carry);
  760. pr[3] = static_cast<limb_type>(carry >> limb_bits);
  761. #endif
  762. result.sign(false);
  763. result.normalize();
  764. }
  765. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1,
  766. std::size_t MinBits2, std::size_t MaxBits2, cpp_integer_type SignType2, cpp_int_check_type Checked2, class Allocator2>
  767. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
  768. !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && is_trivial_cpp_int<cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2> >::value && is_trivial_cpp_int<cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2> >::value>::type
  769. eval_multiply(
  770. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  771. cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2> const& a,
  772. cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2> const& b)
  773. {
  774. using canonical_type = typename boost::multiprecision::detail::canonical<typename cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>::local_limb_type, cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::type;
  775. eval_multiply(result, static_cast<canonical_type>(*a.limbs()), static_cast<canonical_type>(*b.limbs()));
  776. result.sign(a.sign() != b.sign());
  777. }
  778. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class SI>
  779. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_signed<SI>::value && boost::multiprecision::detail::is_integral<SI>::value && (sizeof(SI) <= sizeof(signed_double_limb_type) / 2)>::type
  780. eval_multiply(
  781. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  782. SI a, SI b)
  783. {
  784. result = static_cast<signed_double_limb_type>(a) * static_cast<signed_double_limb_type>(b);
  785. }
  786. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class UI>
  787. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<boost::multiprecision::detail::is_unsigned<UI>::value && (sizeof(UI) <= sizeof(signed_double_limb_type) / 2)>::type
  788. eval_multiply(
  789. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  790. UI a, UI b)
  791. {
  792. result = static_cast<double_limb_type>(a) * static_cast<double_limb_type>(b);
  793. }
  794. #ifdef BOOST_MSVC
  795. #pragma warning(pop)
  796. #endif
  797. }}} // namespace boost::multiprecision::backends
  798. #endif