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 http://www.boost.org/LICENSE_1_
6 // Comparison operators for cpp_int_backend:
8 #ifndef BOOST_MP_CPP_INT_DIV_HPP
9 #define BOOST_MP_CPP_INT_DIV_HPP
11 namespace boost{ namespace multiprecision{ namespace backends{
13 template <class CppInt1, class CppInt2, class CppInt3>
14 void divide_unsigned_helper(
20 if(((void*)result == (void*)&x) || ((void*)&r == (void*)&x))
23 divide_unsigned_helper(result, t, y, r);
26 if(((void*)result == (void*)&y) || ((void*)&r == (void*)&y))
29 divide_unsigned_helper(result, x, t, r);
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.
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.
52 using default_ops::eval_subtract;
57 divide_unsigned_helper(result, x, y, rem);
63 // Find the most significant words of numerator and denominator.
65 limb_type y_order = y.size() - 1;
70 // Only a single non-zero limb in the denominator, in this case
71 // we can use a specialized divide-by-single-limb routine which is
72 // much faster. This also handles division by zero:
74 divide_unsigned_helper(result, x, y.limbs()[y_order], r);
78 typename CppInt2::const_limb_pointer px = x.limbs();
79 typename CppInt3::const_limb_pointer py = y.limbs();
81 limb_type r_order = x.size() - 1;
82 if((r_order == 0) && (*px == 0))
84 // x is zero, so is the result:
94 *result = static_cast<limb_type>(0u);
96 // Check if the remainder is already less than the divisor, if so
97 // we already have the result. Note we try and avoid a full compare
100 if(r_order <= y_order)
102 if((r_order < y_order) || (r.compare_unsigned(y) < 0))
112 // See if we can short-circuit long division, and use basic arithmetic instead:
118 *result = px[0] / py[0];
123 else if(r_order == 1)
125 double_limb_type a, b;
126 a = (static_cast<double_limb_type>(px[1]) << CppInt1::limb_bits) | px[0];
128 (static_cast<double_limb_type>(py[1]) << CppInt1::limb_bits) | py[0]
141 result->resize(1 + r_order - y_order, 1 + r_order - y_order);
142 typename CppInt1::const_limb_pointer prem = r.limbs();
143 // This is initialised just to keep the compiler from emitting useless warnings later on:
144 typename CppInt1::limb_pointer pr
145 = typename CppInt1::limb_pointer();
148 pr = result->limbs();
149 for(unsigned i = 1; i < 1 + r_order - y_order; ++i)
152 bool first_pass = true;
157 // Calculate our best guess for how many times y divides into r:
160 if((prem[r_order] <= py[y_order]) && (r_order > 0))
162 double_limb_type a, b, v;
163 a = (static_cast<double_limb_type>(prem[r_order]) << CppInt1::limb_bits) | prem[r_order - 1];
166 if(v > CppInt1::max_limb_value)
170 guess = static_cast<limb_type>(v);
174 else if(r_order == 0)
176 guess = prem[0] / py[y_order];
180 double_limb_type a, b, v;
181 a = (static_cast<double_limb_type>(prem[r_order]) << CppInt1::limb_bits) | prem[r_order - 1];
182 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);
184 guess = static_cast<limb_type>(v);
186 BOOST_ASSERT(guess); // If the guess ever gets to zero we go on forever....
190 limb_type shift = r_order - y_order;
195 if(pr[shift] > guess)
199 t.resize(shift + 1, shift + 1);
200 t.limbs()[shift] = guess;
201 for(unsigned i = 0; i < shift; ++i)
203 eval_subtract(*result, t);
206 else if(CppInt1::max_limb_value - pr[shift] > guess)
210 t.resize(shift + 1, shift + 1);
211 t.limbs()[shift] = guess;
212 for(unsigned i = 0; i < shift; ++i)
214 eval_add(*result, t);
218 // Calculate guess * y, we use a fused mutiply-shift O(N) for this
219 // rather than a full O(N^2) multiply:
221 double_limb_type carry = 0;
222 t.resize(y.size() + shift + 1, y.size() + shift);
223 bool truncated_t = (t.size() != y.size() + shift + 1);
224 typename CppInt1::limb_pointer pt = t.limbs();
225 for(unsigned i = 0; i < shift; ++i)
227 for(unsigned i = 0; i < y.size(); ++i)
229 carry += static_cast<double_limb_type>(py[i]) * static_cast<double_limb_type>(guess);
230 #ifdef __MSVC_RUNTIME_CHECKS
231 pt[i + shift] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
233 pt[i + shift] = static_cast<limb_type>(carry);
235 carry >>= CppInt1::limb_bits;
237 if(carry && !truncated_t)
239 #ifdef __MSVC_RUNTIME_CHECKS
240 pt[t.size() - 1] = static_cast<limb_type>(carry & ~static_cast<limb_type>(0));
242 pt[t.size() - 1] = static_cast<limb_type>(carry);
245 else if(!truncated_t)
247 t.resize(t.size() - 1, t.size() - 1);
250 // Update r in a way that won't actually produce a negative result
251 // in case the argument types are unsigned:
253 if(truncated_t && carry)
255 // We need to calculate 2^n + t - r
256 // where n is the number of bits in this type.
257 // Simplest way is to get 2^n - r by complementing
258 // r, then add t to it. Note that we can't call eval_complement
259 // in case this is a signed checked type:
260 for(unsigned i = 0; i <= r_order; ++i)
261 r.limbs()[i] = ~prem[i];
267 else if(r.compare(t) > 0)
279 // First time through we need to strip any leading zero, otherwise
280 // the termination condition goes belly-up:
282 if(result && first_pass)
285 while(pr[result->size() - 1] == 0)
286 result->resize(result->size() - 1, result->size() - 1);
291 r_order = r.size() - 1;
292 if(r_order < y_order)
295 // Termination condition is really just a check that r > y, but with a common
296 // short-circuit case handled first:
297 while((r_order > y_order) || (r.compare_unsigned(y) >= 0));
300 // We now just have to normalise the result:
302 if(r_neg && eval_get_sign(r))
304 // We have one too many in the result:
306 eval_decrement(*result);
313 eval_subtract(r, y, r);
316 BOOST_ASSERT(r.compare_unsigned(y) < 0); // remainder must be less than the divisor or our code has failed
319 template <class CppInt1, class CppInt2>
320 void divide_unsigned_helper(
326 if(((void*)result == (void*)&x) || ((void*)&r == (void*)&x))
329 divide_unsigned_helper(result, t, y, r);
336 divide_unsigned_helper(result, x, y, rem);
341 // As above, but simplified for integer divisor:
343 using default_ops::eval_subtract;
347 BOOST_THROW_EXCEPTION(std::overflow_error("Integer Division by zero."));
350 // Find the most significant word of numerator.
352 limb_type r_order = x.size() - 1;
355 // Set remainder and result to their initial values:
359 typename CppInt1::limb_pointer pr = r.limbs();
362 // check for x < y, try to do this without actually having to
363 // do a full comparison:
365 if((r_order == 0) && (*pr < y))
368 *result = static_cast<limb_type>(0u);
373 // See if we can short-circuit long division, and use basic arithmetic instead:
380 result->sign(x.sign());
386 else if(r_order == 1)
389 a = (static_cast<double_limb_type>(pr[r_order]) << CppInt1::limb_bits) | pr[0];
393 result->sign(x.sign());
400 // This is initialised just to keep the compiler from emitting useless warnings later on:
401 typename CppInt1::limb_pointer pres = typename CppInt1::limb_pointer();
404 result->resize(r_order + 1, r_order + 1);
405 pres = result->limbs();
406 if(result->size() > r_order)
407 pres[r_order] = 0; // just in case we don't set the most significant limb below.
413 // Calculate our best guess for how many times y divides into r:
415 if((pr[r_order] < y) && r_order)
417 double_limb_type a, b;
418 a = (static_cast<double_limb_type>(pr[r_order]) << CppInt1::limb_bits) | pr[r_order - 1];
420 r.resize(r.size() - 1, r.size() - 1);
422 pr[r_order] = static_cast<limb_type>(b);
424 pres[r_order] = static_cast<limb_type>(a / y);
425 if(r_order && pr[r_order] == 0)
427 --r_order; // No remainder, division was exact.
428 r.resize(r.size() - 1, r.size() - 1);
430 pres[r_order] = static_cast<limb_type>(0u);
436 pres[r_order] = pr[r_order] / y;
438 if(r_order && pr[r_order] == 0)
440 --r_order; // No remainder, division was exact.
441 r.resize(r.size() - 1, r.size() - 1);
443 pres[r_order] = static_cast<limb_type>(0u);
447 // Termination condition is really just a check that r >= y, but with two common
448 // short-circuit cases handled first:
449 while(r_order || (pr[r_order] >= y));
454 result->sign(x.sign());
459 BOOST_ASSERT(r.compare(y) < 0); // remainder must be less than the divisor or our code has failed
462 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>
463 BOOST_MP_FORCEINLINE 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
465 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
466 const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
467 const cpp_int_backend<MinBits3, MaxBits3, SignType3, Checked3, Allocator3>& b)
469 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> r;
470 bool s = a.sign() != b.sign();
471 divide_unsigned_helper(&result, a, b, r);
475 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>
476 BOOST_MP_FORCEINLINE 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
478 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
479 const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
482 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> r;
484 divide_unsigned_helper(&result, a, b, r);
488 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>
489 BOOST_MP_FORCEINLINE 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
491 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
492 const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
495 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> r;
496 bool s = a.sign() != (b < 0);
497 divide_unsigned_helper(&result, a, static_cast<limb_type>(boost::multiprecision::detail::unsigned_abs(b)), r);
501 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>
502 BOOST_MP_FORCEINLINE 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
504 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
505 const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& b)
507 // There is no in place divide:
508 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> a(result);
509 eval_divide(result, a, b);
512 template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
513 BOOST_MP_FORCEINLINE typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
515 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
518 // There is no in place divide:
519 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> a(result);
520 eval_divide(result, a, b);
523 template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
524 BOOST_MP_FORCEINLINE typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
526 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
529 // There is no in place divide:
530 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> a(result);
531 eval_divide(result, a, b);
534 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>
535 BOOST_MP_FORCEINLINE 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
537 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
538 const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
539 const cpp_int_backend<MinBits3, MaxBits3, SignType3, Checked3, Allocator3>& b)
542 divide_unsigned_helper(static_cast<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>* >(0), a, b, result);
546 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>
547 BOOST_MP_FORCEINLINE 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
549 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
550 const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a, limb_type b)
553 divide_unsigned_helper(static_cast<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>* >(0), a, b, result);
557 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>
558 BOOST_MP_FORCEINLINE 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
560 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
561 const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& a,
565 divide_unsigned_helper(static_cast<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>* >(0), a, static_cast<limb_type>(boost::multiprecision::detail::unsigned_abs(b)), result);
569 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>
570 BOOST_MP_FORCEINLINE 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
572 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
573 const cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2>& b)
575 // There is no in place divide:
576 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> a(result);
577 eval_modulus(result, a, b);
580 template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
581 BOOST_MP_FORCEINLINE typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
583 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
586 // There is no in place divide:
587 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> a(result);
588 eval_modulus(result, a, b);
591 template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
592 BOOST_MP_FORCEINLINE typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
594 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
597 // There is no in place divide:
598 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> a(result);
599 eval_modulus(result, a, b);
603 // Over again for trivial cpp_int's:
605 template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
606 BOOST_MP_FORCEINLINE typename enable_if_c<
607 is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value
608 && is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value
609 && (is_signed_number<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value
610 || is_signed_number<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value)
613 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
614 const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& o)
617 BOOST_THROW_EXCEPTION(std::overflow_error("Division by zero."));
618 *result.limbs() /= *o.limbs();
619 result.sign(result.sign() != o.sign());
622 template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
623 BOOST_MP_FORCEINLINE typename enable_if_c<
624 is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value
625 && is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value
626 && is_unsigned_number<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value
627 && is_unsigned_number<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value
630 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
631 const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& o)
634 BOOST_THROW_EXCEPTION(std::overflow_error("Division by zero."));
635 *result.limbs() /= *o.limbs();
638 template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
639 BOOST_MP_FORCEINLINE typename enable_if_c<
640 is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value
641 && is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value
644 cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result,
645 const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& o)
648 BOOST_THROW_EXCEPTION(std::overflow_error("Division by zero."));
649 *result.limbs() %= *o.limbs();
650 result.sign(result.sign());