/////////////////////////////////////////////////////////////// // Copyright 2012 John Maddock. Distributed under the Boost // Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_ // // Comparison operators for cpp_int_backend: // #ifndef BOOST_MP_CPP_INT_DIV_HPP #define BOOST_MP_CPP_INT_DIV_HPP namespace boost{ namespace multiprecision{ namespace backends{ template void divide_unsigned_helper( CppInt1* result, const CppInt2& x, const CppInt3& y, CppInt1& r) { if(((void*)result == (void*)&x) || ((void*)&r == (void*)&x)) { CppInt2 t(x); divide_unsigned_helper(result, t, y, r); return; } if(((void*)result == (void*)&y) || ((void*)&r == (void*)&y)) { CppInt3 t(y); divide_unsigned_helper(result, x, t, r); return; } /* Very simple, fairly braindead long division. Start by setting the remainder equal to x, and the result equal to 0. Then in each loop we calculate our "best guess" for how many times y divides into r, add our guess to the result, and subtract guess*y from the remainder r. One wrinkle is that the remainder may go negative, in which case we subtract the current guess from the result rather than adding. The value of the guess is determined by dividing the most-significant-limb of the current remainder by the most-significant-limb of y. Note that there are more efficient algorithms than this available, in particular see Knuth Vol 2. However for small numbers of limbs this generally outperforms the alternatives and avoids the normalisation step which would require extra storage. */ using default_ops::eval_subtract; if(result == &r) { CppInt1 rem; divide_unsigned_helper(result, x, y, rem); r = rem; return; } // // Find the most significant words of numerator and denominator. // limb_type y_order = y.size() - 1; if(y_order == 0) { // // Only a single non-zero limb in the denominator, in this case // we can use a specialized divide-by-single-limb routine which is // much faster. This also handles division by zero: // divide_unsigned_helper(result, x, y.limbs()[y_order], r); return; } typename CppInt2::const_limb_pointer px = x.limbs(); typename CppInt3::const_limb_pointer py = y.limbs(); limb_type r_order = x.size() - 1; if((r_order == 0) && (*px == 0)) { // x is zero, so is the result: r = x; if(result) *result = x; return; } r = x; r.sign(false); if(result) *result = static_cast(0u); // // Check if the remainder is already less than the divisor, if so // we already have the result. Note we try and avoid a full compare // if we can: // if(r_order <= y_order) { if((r_order < y_order) || (r.compare_unsigned(y) < 0)) { return; } } CppInt1 t; bool r_neg = false; // // See if we can short-circuit long division, and use basic arithmetic instead: // if(r_order == 0) { if(result) { *result = px[0] / py[0]; } r = px[0] % py[0]; return; } else if(r_order == 1) { double_limb_type a, b; a = (static_cast(px[1]) << CppInt1::limb_bits) | px[0]; b = y_order ? (static_cast(py[1]) << CppInt1::limb_bits) | py[0] : py[0]; if(result) { *result = a / b; } r = a % b; return; } // // prepare result: // if(result) result->resize(1 + r_order - y_order, 1 + r_order - y_order); typename CppInt1::const_limb_pointer prem = r.limbs(); // This is initialised just to keep the compiler from emitting useless warnings later on: typename CppInt1::limb_pointer pr = typename CppInt1::limb_pointer(); if(result) { pr = result->limbs(); for(unsigned i = 1; i < 1 + r_order - y_order; ++i) pr[i] = 0; } bool first_pass = true; do { // // Calculate our best guess for how many times y divides into r: // limb_type guess; if((prem[r_order] <= py[y_order]) && (r_order > 0)) { double_limb_type a, b, v; a = (static_cast(prem[r_order]) << CppInt1::limb_bits) | prem[r_order - 1]; b = py[y_order]; v = a / b; if(v > CppInt1::max_limb_value) guess = 1; else { guess = static_cast(v); --r_order; } } else if(r_order == 0) { guess = prem[0] / py[y_order]; } else { double_limb_type a, b, v; a = (static_cast(prem[r_order]) << CppInt1::limb_bits) | prem[r_order - 1]; b = (y_order > 0) ? (static_cast(py[y_order]) << CppInt1::limb_bits) | py[y_order - 1] : (static_cast(py[y_order]) << CppInt1::limb_bits); v = a / b; guess = static_cast(v); } BOOST_ASSERT(guess); // If the guess ever gets to zero we go on forever.... // // Update result: // limb_type shift = r_order - y_order; if(result) { if(r_neg) { if(pr[shift] > guess) pr[shift] -= guess; else { t.resize(shift + 1, shift + 1); t.limbs()[shift] = guess; for(unsigned i = 0; i < shift; ++i) t.limbs()[i] = 0; eval_subtract(*result, t); } } else if(CppInt1::max_limb_value - pr[shift] > guess) pr[shift] += guess; else { t.resize(shift + 1, shift + 1); t.limbs()[shift] = guess; for(unsigned i = 0; i < shift; ++i) t.limbs()[i] = 0; eval_add(*result, t); } } // // Calculate guess * y, we use a fused mutiply-shift O(N) for this // rather than a full O(N^2) multiply: // double_limb_type carry = 0; t.resize(y.size() + shift + 1, y.size() + shift); bool truncated_t = !CppInt1::variable && (t.size() != y.size() + shift + 1); typename CppInt1::limb_pointer pt = t.limbs(); for(unsigned i = 0; i < shift; ++i) pt[i] = 0; for(unsigned i = 0; i < y.size(); ++i) { carry += static_cast(py[i]) * static_cast(guess); #ifdef __MSVC_RUNTIME_CHECKS pt[i + shift] = static_cast(carry & ~static_cast(0)); #else pt[i + shift] = static_cast(carry); #endif carry >>= CppInt1::limb_bits; } if(carry && !truncated_t) { #ifdef __MSVC_RUNTIME_CHECKS pt[t.size() - 1] = static_cast(carry & ~static_cast(0)); #else pt[t.size() - 1] = static_cast(carry); #endif } else if(!truncated_t) { t.resize(t.size() - 1, t.size() - 1); } // // Update r in a way that won't actually produce a negative result // in case the argument types are unsigned: // if(truncated_t && carry) { // We need to calculate 2^n + t - r // where n is the number of bits in this type. // Simplest way is to get 2^n - r by complementing // r, then add t to it. Note that we can't call eval_complement // in case this is a signed checked type: for(unsigned i = 0; i <= r_order; ++i) r.limbs()[i] = ~prem[i]; r.normalize(); eval_increment(r); eval_add(r, t); r_neg = !r_neg; } else if(r.compare(t) > 0) { eval_subtract(r, t); } else { r.swap(t); eval_subtract(r, t); prem = r.limbs(); r_neg = !r_neg; } // // First time through we need to strip any leading zero, otherwise // the termination condition goes belly-up: // if(result && first_pass) { first_pass = false; while(pr[result->size() - 1] == 0) result->resize(result->size() - 1, result->size() - 1); } // // Update r_order: // r_order = r.size() - 1; if(r_order < y_order) break; } // Termination condition is really just a check that r > y, but with a common // short-circuit case handled first: while((r_order > y_order) || (r.compare_unsigned(y) >= 0)); // // We now just have to normalise the result: // if(r_neg && eval_get_sign(r)) { // We have one too many in the result: if(result) eval_decrement(*result); if(y.sign()) { r.negate(); eval_subtract(r, y); } else eval_subtract(r, y, r); } BOOST_ASSERT(r.compare_unsigned(y) < 0); // remainder must be less than the divisor or our code has failed } template void divide_unsigned_helper( CppInt1* result, const CppInt2& x, limb_type y, CppInt1& r) { if(((void*)result == (void*)&x) || ((void*)&r == (void*)&x)) { CppInt2 t(x); divide_unsigned_helper(result, t, y, r); return; } if(result == &r) { CppInt1 rem; divide_unsigned_helper(result, x, y, rem); r = rem; return; } // As above, but simplified for integer divisor: using default_ops::eval_subtract; if(y == 0) { BOOST_THROW_EXCEPTION(std::overflow_error("Integer Division by zero.")); } // // Find the most significant word of numerator. // limb_type r_order = x.size() - 1; // // Set remainder and result to their initial values: // r = x; r.sign(false); typename CppInt1::limb_pointer pr = r.limbs(); // // check for x < y, try to do this without actually having to // do a full comparison: // if((r_order == 0) && (*pr < y)) { if(result) *result = static_cast(0u); return; } // // See if we can short-circuit long division, and use basic arithmetic instead: // if(r_order == 0) { if(result) { *result = *pr / y; result->sign(x.sign()); } *pr %= y; r.sign(x.sign()); return; } else if(r_order == 1) { double_limb_type a; a = (static_cast(pr[r_order]) << CppInt1::limb_bits) | pr[0]; if(result) { *result = a / y; result->sign(x.sign()); } r = a % y; r.sign(x.sign()); return; } // This is initialised just to keep the compiler from emitting useless warnings later on: typename CppInt1::limb_pointer pres = typename CppInt1::limb_pointer(); if(result) { result->resize(r_order + 1, r_order + 1); pres = result->limbs(); if(result->size() > r_order) pres[r_order] = 0; // just in case we don't set the most significant limb below. } do { // // Calculate our best guess for how many times y divides into r: // if((pr[r_order] < y) && r_order) { double_limb_type a, b; a = (static_cast(pr[r_order]) << CppInt1::limb_bits) | pr[r_order - 1]; b = a % y; r.resize(r.size() - 1, r.size() - 1); --r_order; pr[r_order] = static_cast(b); if(result) pres[r_order] = static_cast(a / y); if(r_order && pr[r_order] == 0) { --r_order; // No remainder, division was exact. r.resize(r.size() - 1, r.size() - 1); if(result) pres[r_order] = static_cast(0u); } } else { if(result) pres[r_order] = pr[r_order] / y; pr[r_order] %= y; if(r_order && pr[r_order] == 0) { --r_order; // No remainder, division was exact. r.resize(r.size() - 1, r.size() - 1); if(result) pres[r_order] = static_cast(0u); } } } // Termination condition is really just a check that r >= y, but with two common // short-circuit cases handled first: while(r_order || (pr[r_order] >= y)); if(result) { result->normalize(); result->sign(x.sign()); } r.normalize(); r.sign(x.sign()); BOOST_ASSERT(r.compare(y) < 0); // remainder must be less than the divisor or our code has failed } template BOOST_MP_FORCEINLINE typename enable_if_c >::value && !is_trivial_cpp_int >::value && !is_trivial_cpp_int >::value >::type eval_divide( cpp_int_backend& result, const cpp_int_backend& a, const cpp_int_backend& b) { cpp_int_backend r; bool s = a.sign() != b.sign(); divide_unsigned_helper(&result, a, b, r); result.sign(s); } template BOOST_MP_FORCEINLINE typename enable_if_c >::value && !is_trivial_cpp_int >::value >::type eval_divide( cpp_int_backend& result, const cpp_int_backend& a, limb_type& b) { cpp_int_backend r; bool s = a.sign(); divide_unsigned_helper(&result, a, b, r); result.sign(s); } template BOOST_MP_FORCEINLINE typename enable_if_c >::value && !is_trivial_cpp_int >::value >::type eval_divide( cpp_int_backend& result, const cpp_int_backend& a, signed_limb_type& b) { cpp_int_backend r; bool s = a.sign() != (b < 0); divide_unsigned_helper(&result, a, static_cast(boost::multiprecision::detail::unsigned_abs(b)), r); result.sign(s); } template BOOST_MP_FORCEINLINE typename enable_if_c >::value && !is_trivial_cpp_int >::value >::type eval_divide( cpp_int_backend& result, const cpp_int_backend& b) { // There is no in place divide: cpp_int_backend a(result); eval_divide(result, a, b); } template BOOST_MP_FORCEINLINE typename enable_if_c >::value>::type eval_divide( cpp_int_backend& result, limb_type b) { // There is no in place divide: cpp_int_backend a(result); eval_divide(result, a, b); } template BOOST_MP_FORCEINLINE typename enable_if_c >::value>::type eval_divide( cpp_int_backend& result, signed_limb_type b) { // There is no in place divide: cpp_int_backend a(result); eval_divide(result, a, b); } template BOOST_MP_FORCEINLINE typename enable_if_c >::value && !is_trivial_cpp_int >::value && !is_trivial_cpp_int >::value >::type eval_modulus( cpp_int_backend& result, const cpp_int_backend& a, const cpp_int_backend& b) { bool s = a.sign(); divide_unsigned_helper(static_cast* >(0), a, b, result); result.sign(s); } template BOOST_MP_FORCEINLINE typename enable_if_c >::value && !is_trivial_cpp_int >::value >::type eval_modulus( cpp_int_backend& result, const cpp_int_backend& a, limb_type b) { bool s = a.sign(); divide_unsigned_helper(static_cast* >(0), a, b, result); result.sign(s); } template BOOST_MP_FORCEINLINE typename enable_if_c >::value && !is_trivial_cpp_int >::value >::type eval_modulus( cpp_int_backend& result, const cpp_int_backend& a, signed_limb_type b) { bool s = a.sign(); divide_unsigned_helper(static_cast* >(0), a, static_cast(boost::multiprecision::detail::unsigned_abs(b)), result); result.sign(s); } template BOOST_MP_FORCEINLINE typename enable_if_c >::value && !is_trivial_cpp_int >::value >::type eval_modulus( cpp_int_backend& result, const cpp_int_backend& b) { // There is no in place divide: cpp_int_backend a(result); eval_modulus(result, a, b); } template BOOST_MP_FORCEINLINE typename enable_if_c >::value>::type eval_modulus( cpp_int_backend& result, limb_type b) { // There is no in place divide: cpp_int_backend a(result); eval_modulus(result, a, b); } template BOOST_MP_FORCEINLINE typename enable_if_c >::value>::type eval_modulus( cpp_int_backend& result, signed_limb_type b) { // There is no in place divide: cpp_int_backend a(result); eval_modulus(result, a, b); } // // Over again for trivial cpp_int's: // template BOOST_MP_FORCEINLINE typename enable_if_c< is_trivial_cpp_int >::value && is_trivial_cpp_int >::value && (is_signed_number >::value || is_signed_number >::value) >::type eval_divide( cpp_int_backend& result, const cpp_int_backend& o) { if(!*o.limbs()) BOOST_THROW_EXCEPTION(std::overflow_error("Division by zero.")); *result.limbs() /= *o.limbs(); result.sign(result.sign() != o.sign()); } template BOOST_MP_FORCEINLINE typename enable_if_c< is_trivial_cpp_int >::value && is_trivial_cpp_int >::value && is_unsigned_number >::value && is_unsigned_number >::value >::type eval_divide( cpp_int_backend& result, const cpp_int_backend& o) { if(!*o.limbs()) BOOST_THROW_EXCEPTION(std::overflow_error("Division by zero.")); *result.limbs() /= *o.limbs(); } template BOOST_MP_FORCEINLINE typename enable_if_c< is_trivial_cpp_int >::value && is_trivial_cpp_int >::value >::type eval_modulus( cpp_int_backend& result, const cpp_int_backend& o) { if(!*o.limbs()) BOOST_THROW_EXCEPTION(std::overflow_error("Division by zero.")); *result.limbs() %= *o.limbs(); result.sign(result.sign()); } }}} // namespaces #endif