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