divide.hpp 27 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661
  1. ///////////////////////////////////////////////////////////////
  2. // Copyright 2012 John Maddock. Distributed under the Boost
  3. // Software License, Version 1.0. (See accompanying file
  4. // LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
  5. //
  6. // Comparison operators for cpp_int_backend:
  7. //
  8. #ifndef BOOST_MP_CPP_INT_DIVIDE_HPP
  9. #define BOOST_MP_CPP_INT_DIVIDE_HPP
  10. #include <boost/multiprecision/detail/no_exceptions_support.hpp>
  11. #include <boost/multiprecision/detail/assert.hpp>
  12. namespace boost { namespace multiprecision { namespace backends {
  13. template <class CppInt1, class CppInt2, class CppInt3>
  14. BOOST_MP_CXX14_CONSTEXPR void divide_unsigned_helper(
  15. CppInt1* result,
  16. const CppInt2& x,
  17. const CppInt3& y,
  18. CppInt1& r)
  19. {
  20. if (((void*)result == (void*)&x) || ((void*)&r == (void*)&x))
  21. {
  22. CppInt2 t(x);
  23. divide_unsigned_helper(result, t, y, r);
  24. return;
  25. }
  26. if (((void*)result == (void*)&y) || ((void*)&r == (void*)&y))
  27. {
  28. CppInt3 t(y);
  29. divide_unsigned_helper(result, x, t, r);
  30. return;
  31. }
  32. /*
  33. Very simple, fairly braindead long division.
  34. Start by setting the remainder equal to x, and the
  35. result equal to 0. Then in each loop we calculate our
  36. "best guess" for how many times y divides into r,
  37. add our guess to the result, and subtract guess*y
  38. from the remainder r. One wrinkle is that the remainder
  39. may go negative, in which case we subtract the current guess
  40. from the result rather than adding. The value of the guess
  41. is determined by dividing the most-significant-limb of the
  42. current remainder by the most-significant-limb of y.
  43. Note that there are more efficient algorithms than this
  44. available, in particular see Knuth Vol 2. However for small
  45. numbers of limbs this generally outperforms the alternatives
  46. and avoids the normalisation step which would require extra storage.
  47. */
  48. using default_ops::eval_subtract;
  49. if (result == &r)
  50. {
  51. CppInt1 rem;
  52. divide_unsigned_helper(result, x, y, rem);
  53. r = rem;
  54. return;
  55. }
  56. //
  57. // Find the most significant words of numerator and denominator.
  58. //
  59. std::size_t y_order = y.size() - 1;
  60. if (y_order == 0)
  61. {
  62. //
  63. // Only a single non-zero limb in the denominator, in this case
  64. // we can use a specialized divide-by-single-limb routine which is
  65. // much faster. This also handles division by zero:
  66. //
  67. divide_unsigned_helper(result, x, y.limbs()[y_order], r);
  68. return;
  69. }
  70. typename CppInt2::const_limb_pointer px = x.limbs();
  71. typename CppInt3::const_limb_pointer py = y.limbs();
  72. std::size_t r_order = x.size() - 1;
  73. if ((r_order == 0) && (*px == 0))
  74. {
  75. // x is zero, so is the result:
  76. r = x;
  77. if (result)
  78. *result = x;
  79. return;
  80. }
  81. r = x;
  82. r.sign(false);
  83. if (result)
  84. *result = static_cast<limb_type>(0u);
  85. //
  86. // Check if the remainder is already less than the divisor, if so
  87. // we already have the result. Note we try and avoid a full compare
  88. // if we can:
  89. //
  90. if (r_order <= y_order)
  91. {
  92. if ((r_order < y_order) || (r.compare_unsigned(y) < 0))
  93. {
  94. return;
  95. }
  96. }
  97. CppInt1 t;
  98. bool r_neg = false;
  99. //
  100. // See if we can short-circuit long division, and use basic arithmetic instead:
  101. //
  102. if (r_order == 0)
  103. {
  104. if (result)
  105. {
  106. *result = px[0] / py[0];
  107. }
  108. r = px[0] % py[0];
  109. return;
  110. }
  111. else if (r_order == 1)
  112. {
  113. double_limb_type a = (static_cast<double_limb_type>(px[1]) << CppInt1::limb_bits) | px[0];
  114. double_limb_type b = y_order ? (static_cast<double_limb_type>(py[1]) << CppInt1::limb_bits) | py[0]
  115. : py[0];
  116. if (result)
  117. {
  118. *result = a / b;
  119. }
  120. r = a % b;
  121. return;
  122. }
  123. //
  124. // prepare result:
  125. //
  126. if (result)
  127. result->resize(1 + r_order - y_order, 1 + r_order - y_order);
  128. typename CppInt1::const_limb_pointer prem = r.limbs();
  129. // This is initialised just to keep the compiler from emitting useless warnings later on:
  130. typename CppInt1::limb_pointer pr = typename CppInt1::limb_pointer();
  131. if (result)
  132. {
  133. pr = result->limbs();
  134. for (std::size_t i = 1; i < 1 + r_order - y_order; ++i)
  135. pr[i] = 0;
  136. }
  137. bool first_pass = true;
  138. do
  139. {
  140. //
  141. // Calculate our best guess for how many times y divides into r:
  142. //
  143. limb_type guess = 1;
  144. if ((prem[r_order] <= py[y_order]) && (r_order > 0))
  145. {
  146. double_limb_type a = (static_cast<double_limb_type>(prem[r_order]) << CppInt1::limb_bits) | prem[r_order - 1];
  147. double_limb_type b = py[y_order];
  148. double_limb_type v = a / b;
  149. if (v <= CppInt1::max_limb_value)
  150. {
  151. guess = static_cast<limb_type>(v);
  152. --r_order;
  153. }
  154. }
  155. else if (r_order == 0)
  156. {
  157. guess = prem[0] / py[y_order];
  158. }
  159. else
  160. {
  161. double_limb_type a = (static_cast<double_limb_type>(prem[r_order]) << CppInt1::limb_bits) | prem[r_order - 1];
  162. double_limb_type b = (y_order > 0) ? (static_cast<double_limb_type>(py[y_order]) << CppInt1::limb_bits) | py[y_order - 1] : (static_cast<double_limb_type>(py[y_order]) << CppInt1::limb_bits);
  163. BOOST_MP_ASSERT(b);
  164. double_limb_type v = a / b;
  165. guess = static_cast<limb_type>(v);
  166. }
  167. BOOST_MP_ASSERT(guess); // If the guess ever gets to zero we go on forever....
  168. //
  169. // Update result:
  170. //
  171. std::size_t shift = r_order - y_order;
  172. if (result)
  173. {
  174. if (r_neg)
  175. {
  176. if (pr[shift] > guess)
  177. pr[shift] -= guess;
  178. else
  179. {
  180. t.resize(shift + 1, shift + 1);
  181. t.limbs()[shift] = guess;
  182. for (std::size_t i = 0; i < shift; ++i)
  183. t.limbs()[i] = 0;
  184. eval_subtract(*result, t);
  185. }
  186. }
  187. else if (CppInt1::max_limb_value - pr[shift] > guess)
  188. pr[shift] += guess;
  189. else
  190. {
  191. t.resize(shift + 1, shift + 1);
  192. t.limbs()[shift] = guess;
  193. for (std::size_t i = 0; i < shift; ++i)
  194. t.limbs()[i] = 0;
  195. eval_add(*result, t);
  196. }
  197. }
  198. //
  199. // Calculate guess * y, we use a fused mutiply-shift O(N) for this
  200. // rather than a full O(N^2) multiply:
  201. //
  202. double_limb_type carry = 0;
  203. t.resize(y.size() + shift + 1, y.size() + shift);
  204. bool truncated_t = (t.size() != y.size() + shift + 1);
  205. typename CppInt1::limb_pointer pt = t.limbs();
  206. for (std::size_t i = 0; i < shift; ++i)
  207. pt[i] = 0;
  208. for (std::size_t i = 0; i < y.size(); ++i)
  209. {
  210. carry += static_cast<double_limb_type>(py[i]) * static_cast<double_limb_type>(guess);
  211. #ifdef __MSVC_RUNTIME_CHECKS
  212. pt[i + shift] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  213. #else
  214. pt[i + shift] = static_cast<limb_type>(carry);
  215. #endif
  216. carry >>= CppInt1::limb_bits;
  217. }
  218. if (carry && !truncated_t)
  219. {
  220. #ifdef __MSVC_RUNTIME_CHECKS
  221. pt[t.size() - 1] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
  222. #else
  223. pt[t.size() - 1] = static_cast<limb_type>(carry);
  224. #endif
  225. }
  226. else if (!truncated_t)
  227. {
  228. t.resize(t.size() - 1, t.size() - 1);
  229. }
  230. //
  231. // Update r in a way that won't actually produce a negative result
  232. // in case the argument types are unsigned:
  233. //
  234. if (truncated_t && carry)
  235. {
  236. // We need to calculate 2^n + t - r
  237. // where n is the number of bits in this type.
  238. // Simplest way is to get 2^n - r by complementing
  239. // r, then add t to it. Note that we can't call eval_complement
  240. // in case this is a signed checked type:
  241. for (std::size_t i = 0; i <= r_order; ++i)
  242. r.limbs()[i] = ~prem[i];
  243. r.normalize();
  244. eval_increment(r);
  245. eval_add(r, t);
  246. r_neg = !r_neg;
  247. }
  248. else if (r.compare(t) > 0)
  249. {
  250. eval_subtract(r, t);
  251. }
  252. else
  253. {
  254. r.swap(t);
  255. eval_subtract(r, t);
  256. prem = r.limbs();
  257. r_neg = !r_neg;
  258. }
  259. //
  260. // First time through we need to strip any leading zero, otherwise
  261. // the termination condition goes belly-up:
  262. //
  263. if (result && first_pass)
  264. {
  265. first_pass = false;
  266. while (pr[result->size() - 1] == 0)
  267. result->resize(result->size() - 1, result->size() - 1);
  268. }
  269. //
  270. // Update r_order:
  271. //
  272. r_order = r.size() - 1;
  273. if (r_order < y_order)
  274. break;
  275. }
  276. // Termination condition is really just a check that r > y, but with a common
  277. // short-circuit case handled first:
  278. while ((r_order > y_order) || (r.compare_unsigned(y) >= 0));
  279. //
  280. // We now just have to normalise the result:
  281. //
  282. if (r_neg && eval_get_sign(r))
  283. {
  284. // We have one too many in the result:
  285. if (result)
  286. eval_decrement(*result);
  287. if (y.sign())
  288. {
  289. r.negate();
  290. eval_subtract(r, y);
  291. }
  292. else
  293. eval_subtract(r, y, r);
  294. }
  295. BOOST_MP_ASSERT(r.compare_unsigned(y) < 0); // remainder must be less than the divisor or our code has failed
  296. }
  297. template <class CppInt1, class CppInt2>
  298. BOOST_MP_CXX14_CONSTEXPR void divide_unsigned_helper(
  299. CppInt1* result,
  300. const CppInt2& x,
  301. limb_type y,
  302. CppInt1& r)
  303. {
  304. if (((void*)result == (void*)&x) || ((void*)&r == (void*)&x))
  305. {
  306. CppInt2 t(x);
  307. divide_unsigned_helper(result, t, y, r);
  308. return;
  309. }
  310. if (result == &r)
  311. {
  312. CppInt1 rem;
  313. divide_unsigned_helper(result, x, y, rem);
  314. r = rem;
  315. return;
  316. }
  317. // As above, but simplified for integer divisor:
  318. using default_ops::eval_subtract;
  319. if (y == 0)
  320. {
  321. BOOST_MP_THROW_EXCEPTION(std::overflow_error("Integer Division by zero."));
  322. }
  323. //
  324. // Find the most significant word of numerator.
  325. //
  326. std::size_t r_order = x.size() - 1;
  327. //
  328. // Set remainder and result to their initial values:
  329. //
  330. r = x;
  331. r.sign(false);
  332. typename CppInt1::limb_pointer pr = r.limbs();
  333. //
  334. // check for x < y, try to do this without actually having to
  335. // do a full comparison:
  336. //
  337. if ((r_order == 0) && (*pr < y))
  338. {
  339. if (result)
  340. *result = static_cast<limb_type>(0u);
  341. return;
  342. }
  343. //
  344. // See if we can short-circuit long division, and use basic arithmetic instead:
  345. //
  346. if (r_order == 0)
  347. {
  348. if (result)
  349. {
  350. *result = *pr / y;
  351. result->sign(x.sign());
  352. }
  353. *pr %= y;
  354. r.sign(x.sign());
  355. return;
  356. }
  357. else if (r_order == 1)
  358. {
  359. double_limb_type a = (static_cast<double_limb_type>(pr[r_order]) << CppInt1::limb_bits) | pr[0];
  360. if (result)
  361. {
  362. *result = a / y;
  363. result->sign(x.sign());
  364. }
  365. r = a % y;
  366. r.sign(x.sign());
  367. return;
  368. }
  369. // This is initialised just to keep the compiler from emitting useless warnings later on:
  370. typename CppInt1::limb_pointer pres = typename CppInt1::limb_pointer();
  371. if (result)
  372. {
  373. result->resize(r_order + 1, r_order + 1);
  374. pres = result->limbs();
  375. if (result->size() > r_order)
  376. pres[r_order] = 0; // just in case we don't set the most significant limb below.
  377. }
  378. do
  379. {
  380. //
  381. // Calculate our best guess for how many times y divides into r:
  382. //
  383. if ((pr[r_order] < y) && r_order)
  384. {
  385. double_limb_type a = (static_cast<double_limb_type>(pr[r_order]) << CppInt1::limb_bits) | pr[r_order - 1];
  386. double_limb_type b = a % y;
  387. r.resize(r.size() - 1, r.size() - 1);
  388. --r_order;
  389. pr[r_order] = static_cast<limb_type>(b);
  390. if (result)
  391. pres[r_order] = static_cast<limb_type>(a / y);
  392. if (r_order && pr[r_order] == 0)
  393. {
  394. --r_order; // No remainder, division was exact.
  395. r.resize(r.size() - 1, r.size() - 1);
  396. if (result)
  397. pres[r_order] = static_cast<limb_type>(0u);
  398. }
  399. }
  400. else
  401. {
  402. if (result)
  403. pres[r_order] = pr[r_order] / y;
  404. pr[r_order] %= y;
  405. if (r_order && pr[r_order] == 0)
  406. {
  407. --r_order; // No remainder, division was exact.
  408. r.resize(r.size() - 1, r.size() - 1);
  409. if (result)
  410. pres[r_order] = static_cast<limb_type>(0u);
  411. }
  412. }
  413. }
  414. // Termination condition is really just a check that r >= y, but with two common
  415. // short-circuit cases handled first:
  416. while (r_order || (pr[r_order] >= y));
  417. if (result)
  418. {
  419. result->normalize();
  420. result->sign(x.sign());
  421. }
  422. r.normalize();
  423. r.sign(x.sign());
  424. BOOST_MP_ASSERT(r.compare(y) < 0); // remainder must be less than the divisor or our code has failed
  425. }
  426. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, std::size_t MinBits2, std::size_t MaxBits2, cpp_integer_type SignType2, cpp_int_check_type Checked2, class Allocator2, std::size_t MinBits3, std::size_t MaxBits3, cpp_integer_type SignType3, cpp_int_check_type Checked3, class Allocator3>
  427. 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 && !is_trivial_cpp_int<cpp_int_backend<MinBits3, MaxBits3, SignType3, Checked3, Allocator3> >::value>::type
  428. eval_divide(
  429. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  430. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
  431. const cpp_int_backend<MinBits3, MaxBits3, SignType3, Checked3, Allocator3>& b)
  432. {
  433. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> r;
  434. bool s = a.sign() != b.sign();
  435. divide_unsigned_helper(&result, a, b, r);
  436. result.sign(s);
  437. }
  438. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, std::size_t MinBits2, std::size_t MaxBits2, cpp_integer_type SignType2, cpp_int_check_type Checked2, class Allocator2>
  439. 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
  440. eval_divide(
  441. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  442. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
  443. limb_type& b)
  444. {
  445. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> r;
  446. bool s = a.sign();
  447. divide_unsigned_helper(&result, a, b, r);
  448. result.sign(s);
  449. }
  450. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, std::size_t MinBits2, std::size_t MaxBits2, cpp_integer_type SignType2, cpp_int_check_type Checked2, class Allocator2>
  451. 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
  452. eval_divide(
  453. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  454. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
  455. signed_limb_type& b)
  456. {
  457. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> r;
  458. bool s = a.sign() != (b < 0);
  459. divide_unsigned_helper(&result, a, static_cast<limb_type>(boost::multiprecision::detail::unsigned_abs(b)), r);
  460. result.sign(s);
  461. }
  462. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, std::size_t MinBits2, std::size_t MaxBits2, cpp_integer_type SignType2, cpp_int_check_type Checked2, class Allocator2>
  463. 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
  464. eval_divide(
  465. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  466. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& b)
  467. {
  468. // There is no in place divide:
  469. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> a(result);
  470. eval_divide(result, a, b);
  471. }
  472. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  473. 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
  474. eval_divide(
  475. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  476. limb_type b)
  477. {
  478. // There is no in place divide:
  479. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> a(result);
  480. eval_divide(result, a, b);
  481. }
  482. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  483. 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
  484. eval_divide(
  485. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  486. signed_limb_type b)
  487. {
  488. // There is no in place divide:
  489. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> a(result);
  490. eval_divide(result, a, b);
  491. }
  492. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, std::size_t MinBits2, std::size_t MaxBits2, cpp_integer_type SignType2, cpp_int_check_type Checked2, class Allocator2, std::size_t MinBits3, std::size_t MaxBits3, cpp_integer_type SignType3, cpp_int_check_type Checked3, class Allocator3>
  493. 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 && !is_trivial_cpp_int<cpp_int_backend<MinBits3, MaxBits3, SignType3, Checked3, Allocator3> >::value>::type
  494. eval_modulus(
  495. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  496. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
  497. const cpp_int_backend<MinBits3, MaxBits3, SignType3, Checked3, Allocator3>& b)
  498. {
  499. bool s = a.sign();
  500. if (b.size() == 1)
  501. eval_modulus(result, a, *b.limbs());
  502. else
  503. {
  504. using cpp_int_backend1_type = cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>;
  505. divide_unsigned_helper(static_cast<cpp_int_backend1_type*>(nullptr), a, b, result);
  506. }
  507. result.sign(s);
  508. }
  509. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, std::size_t MinBits2, std::size_t MaxBits2, cpp_integer_type SignType2, cpp_int_check_type Checked2, class Allocator2>
  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 && !is_trivial_cpp_int<cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2> >::value>::type
  511. eval_modulus(
  512. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  513. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
  514. const limb_type mod)
  515. {
  516. const std::ptrdiff_t n = static_cast<std::ptrdiff_t>(a.size());
  517. const double_limb_type two_n_mod =
  518. static_cast<double_limb_type>
  519. (
  520. static_cast<double_limb_type>(1u) + static_cast<limb_type>(static_cast<limb_type>(~static_cast<limb_type>(0u) - mod) % mod)
  521. );
  522. limb_type res = a.limbs()[n - 1] % mod;
  523. for (std::ptrdiff_t i = n - 2; i >= 0; --i)
  524. res = static_cast<limb_type>(static_cast<double_limb_type>(static_cast<double_limb_type>(res * two_n_mod) + a.limbs()[i]) % mod);
  525. //
  526. // We must not modify result until here in case
  527. // result and a are the same object:
  528. //
  529. result.resize(1, 1);
  530. *result.limbs() = res;
  531. result.sign(a.sign());
  532. }
  533. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, std::size_t MinBits2, std::size_t MaxBits2, cpp_integer_type SignType2, cpp_int_check_type Checked2, class Allocator2>
  534. 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
  535. eval_modulus(
  536. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  537. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
  538. signed_limb_type b)
  539. {
  540. const limb_type t = b < 0 ? static_cast<limb_type>(-b) : static_cast<limb_type>(b);
  541. eval_modulus(result, a, t);
  542. result.sign(a.sign());
  543. }
  544. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, std::size_t MinBits2, std::size_t MaxBits2, cpp_integer_type SignType2, cpp_int_check_type Checked2, class Allocator2>
  545. 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
  546. eval_modulus(
  547. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  548. const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& b)
  549. {
  550. // There is no in place divide:
  551. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> a(result);
  552. eval_modulus(result, a, b);
  553. }
  554. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  555. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
  556. eval_modulus(
  557. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  558. limb_type b)
  559. {
  560. // Single limb modulus is in place:
  561. eval_modulus(result, result, b);
  562. }
  563. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  564. 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
  565. eval_modulus(
  566. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  567. signed_limb_type b)
  568. {
  569. // Single limb modulus is in place:
  570. eval_modulus(result, result, b);
  571. }
  572. //
  573. // Over again for trivial cpp_int's:
  574. //
  575. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  576. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
  577. 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
  578. eval_divide(
  579. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  580. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& o)
  581. {
  582. if (!*o.limbs())
  583. BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero."));
  584. *result.limbs() /= *o.limbs();
  585. result.sign(result.sign() != o.sign());
  586. }
  587. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  588. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
  589. 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_unsigned_number<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && is_unsigned_number<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
  590. eval_divide(
  591. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  592. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& o)
  593. {
  594. if (!*o.limbs())
  595. BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero."));
  596. *result.limbs() /= *o.limbs();
  597. }
  598. template <std::size_t MinBits1, std::size_t MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
  599. BOOST_MP_FORCEINLINE BOOST_MP_CXX14_CONSTEXPR typename std::enable_if<
  600. 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>::type
  601. eval_modulus(
  602. cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
  603. const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& o)
  604. {
  605. if (!*o.limbs())
  606. BOOST_MP_THROW_EXCEPTION(std::overflow_error("Division by zero."));
  607. *result.limbs() %= *o.limbs();
  608. result.sign(result.sign());
  609. }
  610. }}} // namespace boost::multiprecision::backends
  611. #endif