multiply.hpp 44 KB

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