multiply.hpp 43 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833
  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 enable_if_c<!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) BOOST_MP_NOEXCEPT_IF((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. typedef cpp_int_backend<MinBits, MaxBits, signed_magnitude, Checked, Allocator> cpp_int_type;
  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 enable_if_c<!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 enable_if_c<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. typedef cpp_int_backend<0, 0, signed_magnitude, unchecked, std::allocator<limb_type> > variable_precision_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 (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 BOOST_MP_CXX14_CONSTEXPR void
  300. eval_multiply_comba(
  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) BOOST_MP_NOEXCEPT_IF((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
  304. {
  305. //
  306. // see PR #182
  307. // Comba Multiplier - based on Paul Comba's
  308. // Exponentiation cryptosystems on the IBM PC, 1990
  309. //
  310. int as = a.size(),
  311. bs = b.size(),
  312. rs = result.size();
  313. typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_pointer pr = result.limbs();
  314. double_limb_type carry = 0,
  315. temp = 0;
  316. limb_type overflow = 0;
  317. const unsigned limb_bits = sizeof(limb_type) * CHAR_BIT;
  318. const bool must_throw = rs < as + bs - 1;
  319. for (int r = 0, lim = (std::min)(rs, as + bs - 1); r < lim; ++r, overflow = 0)
  320. {
  321. int i = r >= as ? as - 1 : r,
  322. j = r - i,
  323. k = i < bs - j ? i + 1 : bs - j; // min(i+1, bs-j);
  324. typename cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>::const_limb_pointer pa = a.limbs() + i;
  325. typename cpp_int_backend<MinBits3, MaxBits3, SignType3, Checked3, Allocator3>::const_limb_pointer pb = b.limbs() + j;
  326. temp = carry;
  327. carry += static_cast<double_limb_type>(*(pa)) * (*(pb));
  328. overflow += carry < temp;
  329. for (--k; k; k--)
  330. {
  331. temp = carry;
  332. carry += static_cast<double_limb_type>(*(--pa)) * (*(++pb));
  333. overflow += carry < temp;
  334. }
  335. *(pr++) = static_cast<limb_type>(carry);
  336. carry = (static_cast<double_limb_type>(overflow) << limb_bits) | (carry >> limb_bits);
  337. }
  338. if (carry || must_throw)
  339. {
  340. resize_for_carry(result, as + bs);
  341. if ((int)result.size() >= as + bs)
  342. *pr = static_cast<limb_type>(carry);
  343. }
  344. }
  345. 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>
  346. inline BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<!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
  347. eval_multiply(
  348. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  349. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
  350. const cpp_int_backend<MinBits3, MaxBits3, SignType3, Checked3, Allocator3>& b)
  351. BOOST_MP_NOEXCEPT_IF((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value
  352. && (karatsuba_cutoff * sizeof(limb_type) * CHAR_BIT > MaxBits1)
  353. && (karatsuba_cutoff * sizeof(limb_type)* CHAR_BIT > MaxBits2)
  354. && (karatsuba_cutoff * sizeof(limb_type)* CHAR_BIT > MaxBits3)))
  355. {
  356. // Uses simple (O(n^2)) multiplication when the limbs are less
  357. // otherwise switches to karatsuba algorithm based on experimental value (~40 limbs)
  358. //
  359. // Trivial cases first:
  360. //
  361. unsigned as = a.size();
  362. unsigned bs = b.size();
  363. typename cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>::const_limb_pointer pa = a.limbs();
  364. typename cpp_int_backend<MinBits3, MaxBits3, SignType3, Checked3, Allocator3>::const_limb_pointer pb = b.limbs();
  365. if (as == 1)
  366. {
  367. bool s = b.sign() != a.sign();
  368. if (bs == 1)
  369. {
  370. result = static_cast<double_limb_type>(*pa) * static_cast<double_limb_type>(*pb);
  371. }
  372. else
  373. {
  374. limb_type l = *pa;
  375. eval_multiply(result, b, l);
  376. }
  377. result.sign(s);
  378. return;
  379. }
  380. if (bs == 1)
  381. {
  382. bool s = b.sign() != a.sign();
  383. limb_type l = *pb;
  384. eval_multiply(result, a, l);
  385. result.sign(s);
  386. return;
  387. }
  388. if ((void*)&result == (void*)&a)
  389. {
  390. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t(a);
  391. eval_multiply(result, t, b);
  392. return;
  393. }
  394. if ((void*)&result == (void*)&b)
  395. {
  396. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t(b);
  397. eval_multiply(result, a, t);
  398. return;
  399. }
  400. #ifdef BOOST_NO_CXX14_CONSTEXPR
  401. static const double_limb_type limb_max = ~static_cast<limb_type>(0u);
  402. static const double_limb_type double_limb_max = ~static_cast<double_limb_type>(0u);
  403. #else
  404. constexpr const double_limb_type limb_max = ~static_cast<limb_type>(0u);
  405. constexpr const double_limb_type double_limb_max = ~static_cast<double_limb_type>(0u);
  406. #endif
  407. result.resize(as + bs, as + bs - 1);
  408. #ifndef BOOST_MP_NO_CONSTEXPR_DETECTION
  409. if (!BOOST_MP_IS_CONST_EVALUATED(as) && (as >= karatsuba_cutoff && bs >= karatsuba_cutoff))
  410. #else
  411. if (as >= karatsuba_cutoff && bs >= karatsuba_cutoff)
  412. #endif
  413. {
  414. setup_karatsuba(result, a, b);
  415. //
  416. // Set the sign of the result:
  417. //
  418. result.sign(a.sign() != b.sign());
  419. return;
  420. }
  421. typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_pointer pr = result.limbs();
  422. BOOST_STATIC_ASSERT(double_limb_max - 2 * limb_max >= limb_max * limb_max);
  423. #ifndef BOOST_MP_NO_CONSTEXPR_DETECTION
  424. if (BOOST_MP_IS_CONST_EVALUATED(as))
  425. {
  426. for (unsigned i = 0; i < result.size(); ++i)
  427. pr[i] = 0;
  428. }
  429. else
  430. #endif
  431. std::memset(pr, 0, result.size() * sizeof(limb_type));
  432. #if defined(BOOST_MP_COMBA)
  433. //
  434. // Comba Multiplier might not be efficient because of less efficient assembly
  435. // by the compiler as of 09/01/2020 (DD/MM/YY). See PR #182
  436. // Till then this will lay dormant :(
  437. //
  438. eval_multiply_comba(result, a, b);
  439. #else
  440. double_limb_type carry = 0;
  441. for (unsigned i = 0; i < as; ++i)
  442. {
  443. BOOST_ASSERT(result.size() > i);
  444. unsigned inner_limit = !is_fixed_precision<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value ? bs : (std::min)(result.size() - i, bs);
  445. unsigned j = 0;
  446. for (; j < inner_limit; ++j)
  447. {
  448. BOOST_ASSERT(i + j < result.size());
  449. #if (!defined(__GLIBCXX__) && !defined(__GLIBCPP__)) || !BOOST_WORKAROUND(BOOST_GCC_VERSION, <= 50100)
  450. BOOST_ASSERT(!std::numeric_limits<double_limb_type>::is_specialized || ((std::numeric_limits<double_limb_type>::max)() - carry >
  451. 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)));
  452. #endif
  453. carry += static_cast<double_limb_type>(pa[i]) * static_cast<double_limb_type>(pb[j]);
  454. BOOST_ASSERT(!std::numeric_limits<double_limb_type>::is_specialized || ((std::numeric_limits<double_limb_type>::max)() - carry >= pr[i + j]));
  455. carry += pr[i + j];
  456. #ifdef __MSVC_RUNTIME_CHECKS
  457. pr[i + j] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  458. #else
  459. pr[i + j] = static_cast<limb_type>(carry);
  460. #endif
  461. carry >>= cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_bits;
  462. BOOST_ASSERT(carry <= (cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::max_limb_value));
  463. }
  464. if (carry)
  465. {
  466. resize_for_carry(result, i + j + 1); // May throw if checking is enabled
  467. if (i + j < result.size())
  468. #ifdef __MSVC_RUNTIME_CHECKS
  469. pr[i + j] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  470. #else
  471. pr[i + j] = static_cast<limb_type>(carry);
  472. #endif
  473. }
  474. carry = 0;
  475. }
  476. #endif // ifdef(BOOST_MP_COMBA) ends
  477. result.normalize();
  478. //
  479. // Set the sign of the result:
  480. //
  481. result.sign(a.sign() != b.sign());
  482. }
  483. 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>
  484. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<!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
  485. eval_multiply(
  486. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  487. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a)
  488. BOOST_MP_NOEXCEPT_IF((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>&>()))))
  489. {
  490. eval_multiply(result, result, a);
  491. }
  492. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  493. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
  494. eval_multiply(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, const limb_type& val)
  495. BOOST_MP_NOEXCEPT_IF((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&>()))))
  496. {
  497. eval_multiply(result, result, val);
  498. }
  499. 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>
  500. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<!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
  501. eval_multiply(
  502. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  503. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
  504. const double_limb_type& val)
  505. BOOST_MP_NOEXCEPT_IF(
  506. (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&>())))
  507. && (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>&>())))
  508. )
  509. {
  510. if (val <= (std::numeric_limits<limb_type>::max)())
  511. {
  512. eval_multiply(result, a, static_cast<limb_type>(val));
  513. }
  514. else
  515. {
  516. #if BOOST_ENDIAN_LITTLE_BYTE && !defined(BOOST_MP_TEST_NO_LE)
  517. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t(val);
  518. #else
  519. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t;
  520. t = val;
  521. #endif
  522. eval_multiply(result, a, t);
  523. }
  524. }
  525. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  526. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
  527. eval_multiply(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, const double_limb_type& val)
  528. BOOST_MP_NOEXCEPT_IF((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&>()))))
  529. {
  530. eval_multiply(result, result, val);
  531. }
  532. 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>
  533. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<!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
  534. eval_multiply(
  535. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  536. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
  537. const signed_limb_type& val)
  538. BOOST_MP_NOEXCEPT_IF((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&>()))))
  539. {
  540. if (val > 0)
  541. eval_multiply(result, a, static_cast<limb_type>(val));
  542. else
  543. {
  544. eval_multiply(result, a, static_cast<limb_type>(boost::multiprecision::detail::unsigned_abs(val)));
  545. result.negate();
  546. }
  547. }
  548. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  549. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
  550. eval_multiply(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, const signed_limb_type& val)
  551. BOOST_MP_NOEXCEPT_IF((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&>()))))
  552. {
  553. eval_multiply(result, result, val);
  554. }
  555. 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>
  556. inline BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<!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
  557. eval_multiply(
  558. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  559. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
  560. const signed_double_limb_type& val)
  561. BOOST_MP_NOEXCEPT_IF(
  562. (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&>())))
  563. && (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>&>())))
  564. )
  565. {
  566. if (val > 0)
  567. {
  568. if (val <= (std::numeric_limits<limb_type>::max)())
  569. {
  570. eval_multiply(result, a, static_cast<limb_type>(val));
  571. return;
  572. }
  573. }
  574. else if (val >= -static_cast<signed_double_limb_type>((std::numeric_limits<limb_type>::max)()))
  575. {
  576. eval_multiply(result, a, static_cast<limb_type>(boost::multiprecision::detail::unsigned_abs(val)));
  577. result.negate();
  578. return;
  579. }
  580. #if BOOST_ENDIAN_LITTLE_BYTE && !defined(BOOST_MP_TEST_NO_LE)
  581. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t(val);
  582. #else
  583. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t;
  584. t = val;
  585. #endif
  586. eval_multiply(result, a, t);
  587. }
  588. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  589. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
  590. eval_multiply(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, const signed_double_limb_type& val)
  591. BOOST_MP_NOEXCEPT_IF(
  592. (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&>())))
  593. && (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>&>())))
  594. )
  595. {
  596. eval_multiply(result, result, val);
  597. }
  598. //
  599. // Now over again for trivial cpp_int's:
  600. //
  601. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  602. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<
  603. 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
  604. eval_multiply(
  605. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  606. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& o) BOOST_MP_NOEXCEPT_IF((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
  607. {
  608. *result.limbs() = detail::checked_multiply(*result.limbs(), *o.limbs(), typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type());
  609. result.sign(result.sign() != o.sign());
  610. result.normalize();
  611. }
  612. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  613. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<
  614. 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
  615. eval_multiply(
  616. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  617. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& o) BOOST_MP_NOEXCEPT_IF((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
  618. {
  619. *result.limbs() = detail::checked_multiply(*result.limbs(), *o.limbs(), typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type());
  620. result.normalize();
  621. }
  622. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  623. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<
  624. 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
  625. eval_multiply(
  626. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  627. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a,
  628. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& b) BOOST_MP_NOEXCEPT_IF((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
  629. {
  630. *result.limbs() = detail::checked_multiply(*a.limbs(), *b.limbs(), typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type());
  631. result.sign(a.sign() != b.sign());
  632. result.normalize();
  633. }
  634. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  635. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<
  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>& a,
  640. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& b) BOOST_MP_NOEXCEPT_IF((is_non_throwing_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value))
  641. {
  642. *result.limbs() = detail::checked_multiply(*a.limbs(), *b.limbs(), typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::checked_type());
  643. result.normalize();
  644. }
  645. //
  646. // Special routines for multiplying two integers to obtain a multiprecision result:
  647. //
  648. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  649. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<
  650. !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
  651. eval_multiply(
  652. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  653. signed_double_limb_type a, signed_double_limb_type b)
  654. {
  655. #ifdef BOOST_NO_CXX14_CONSTEXPR
  656. static const signed_double_limb_type mask = ~static_cast<limb_type>(0);
  657. static const unsigned limb_bits = sizeof(limb_type) * CHAR_BIT;
  658. #else
  659. constexpr const signed_double_limb_type mask = ~static_cast<limb_type>(0);
  660. constexpr const unsigned limb_bits = sizeof(limb_type) * CHAR_BIT;
  661. #endif
  662. bool s = false;
  663. if (a < 0)
  664. {
  665. a = -a;
  666. s = true;
  667. }
  668. if (b < 0)
  669. {
  670. b = -b;
  671. s = !s;
  672. }
  673. double_limb_type w = a & mask;
  674. double_limb_type x = a >> limb_bits;
  675. double_limb_type y = b & mask;
  676. double_limb_type z = b >> limb_bits;
  677. result.resize(4, 4);
  678. limb_type* pr = result.limbs();
  679. double_limb_type carry = w * y;
  680. #ifdef __MSVC_RUNTIME_CHECKS
  681. pr[0] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  682. carry >>= limb_bits;
  683. carry += w * z + x * y;
  684. pr[1] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  685. carry >>= limb_bits;
  686. carry += x * z;
  687. pr[2] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  688. pr[3] = static_cast<limb_type>(carry >> limb_bits);
  689. #else
  690. pr[0] = static_cast<limb_type>(carry);
  691. carry >>= limb_bits;
  692. carry += w * z + x * y;
  693. pr[1] = static_cast<limb_type>(carry);
  694. carry >>= limb_bits;
  695. carry += x * z;
  696. pr[2] = static_cast<limb_type>(carry);
  697. pr[3] = static_cast<limb_type>(carry >> limb_bits);
  698. #endif
  699. result.sign(s);
  700. result.normalize();
  701. }
  702. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  703. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<
  704. !is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
  705. eval_multiply(
  706. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  707. double_limb_type a, double_limb_type b)
  708. {
  709. #ifdef BOOST_NO_CXX14_CONSTEXPR
  710. static const signed_double_limb_type mask = ~static_cast<limb_type>(0);
  711. static const unsigned limb_bits = sizeof(limb_type) * CHAR_BIT;
  712. #else
  713. constexpr const signed_double_limb_type mask = ~static_cast<limb_type>(0);
  714. constexpr const unsigned limb_bits = sizeof(limb_type) * CHAR_BIT;
  715. #endif
  716. double_limb_type w = a & mask;
  717. double_limb_type x = a >> limb_bits;
  718. double_limb_type y = b & mask;
  719. double_limb_type z = b >> limb_bits;
  720. result.resize(4, 4);
  721. limb_type* pr = result.limbs();
  722. double_limb_type carry = w * y;
  723. #ifdef __MSVC_RUNTIME_CHECKS
  724. pr[0] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  725. carry >>= limb_bits;
  726. carry += w * z;
  727. pr[1] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  728. carry >>= limb_bits;
  729. pr[2] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  730. carry = x * y + pr[1];
  731. pr[1] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  732. carry >>= limb_bits;
  733. carry += pr[2] + x * z;
  734. pr[2] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  735. pr[3] = static_cast<limb_type>(carry >> limb_bits);
  736. #else
  737. pr[0] = static_cast<limb_type>(carry);
  738. carry >>= limb_bits;
  739. carry += w * z;
  740. pr[1] = static_cast<limb_type>(carry);
  741. carry >>= limb_bits;
  742. pr[2] = static_cast<limb_type>(carry);
  743. carry = x * y + pr[1];
  744. pr[1] = static_cast<limb_type>(carry);
  745. carry >>= limb_bits;
  746. carry += pr[2] + x * z;
  747. pr[2] = static_cast<limb_type>(carry);
  748. pr[3] = static_cast<limb_type>(carry >> limb_bits);
  749. #endif
  750. result.sign(false);
  751. result.normalize();
  752. }
  753. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1,
  754. unsigned MinBits2, unsigned MaxBits2, cpp_integer_type SignType2, cpp_int_check_type Checked2, class Allocator2>
  755. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<
  756. !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
  757. eval_multiply(
  758. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  759. cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2> const& a,
  760. cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2> const& b)
  761. {
  762. typedef 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 canonical_type;
  763. eval_multiply(result, static_cast<canonical_type>(*a.limbs()), static_cast<canonical_type>(*b.limbs()));
  764. result.sign(a.sign() != b.sign());
  765. }
  766. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class SI>
  767. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<is_signed<SI>::value && (sizeof(SI) <= sizeof(signed_double_limb_type) / 2)>::type
  768. eval_multiply(
  769. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  770. SI a, SI b)
  771. {
  772. result = static_cast<signed_double_limb_type>(a) * static_cast<signed_double_limb_type>(b);
  773. }
  774. template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, class UI>
  775. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<is_unsigned<UI>::value && (sizeof(UI) <= sizeof(signed_double_limb_type) / 2)>::type
  776. eval_multiply(
  777. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  778. UI a, UI b)
  779. {
  780. result = static_cast<double_limb_type>(a) * static_cast<double_limb_type>(b);
  781. }
  782. #ifdef _MSC_VER
  783. #pragma warning(pop)
  784. #endif
  785. }}} // namespace boost::multiprecision::backends
  786. #endif