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