File indexing completed on 2025-01-30 09:59:09
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013 #ifndef BOOST_RANDOM_DETAIL_LARGE_ARITHMETIC_HPP
0014 #define BOOST_RANDOM_DETAIL_LARGE_ARITHMETIC_HPP
0015
0016 #include <boost/cstdint.hpp>
0017 #include <boost/integer.hpp>
0018 #include <boost/limits.hpp>
0019 #include <boost/random/detail/integer_log2.hpp>
0020
0021 #include <boost/random/detail/disable_warnings.hpp>
0022
0023 namespace boost {
0024 namespace random {
0025 namespace detail {
0026
0027 struct div_t {
0028 boost::uintmax_t quotient;
0029 boost::uintmax_t remainder;
0030 };
0031
0032 inline div_t muldivmod(boost::uintmax_t a, boost::uintmax_t b, boost::uintmax_t m)
0033 {
0034 const int bits =
0035 ::std::numeric_limits< ::boost::uintmax_t>::digits / 2;
0036 const ::boost::uintmax_t mask = (::boost::uintmax_t(1) << bits) - 1;
0037 typedef ::boost::uint_t<bits>::fast digit_t;
0038
0039 int shift = std::numeric_limits< ::boost::uintmax_t>::digits - 1
0040 - detail::integer_log2(m);
0041
0042 a <<= shift;
0043 m <<= shift;
0044
0045 digit_t product[4] = { 0, 0, 0, 0 };
0046 digit_t a_[2] = { digit_t(a & mask), digit_t((a >> bits) & mask) };
0047 digit_t b_[2] = { digit_t(b & mask), digit_t((b >> bits) & mask) };
0048 digit_t m_[2] = { digit_t(m & mask), digit_t((m >> bits) & mask) };
0049
0050
0051 for(int i = 0; i < 2; ++i) {
0052 digit_t carry = 0;
0053 for(int j = 0; j < 2; ++j) {
0054 ::boost::uint64_t temp = ::boost::uintmax_t(a_[i]) * b_[j] +
0055 carry + product[i + j];
0056 product[i + j] = digit_t(temp & mask);
0057 carry = digit_t(temp >> bits);
0058 }
0059 if(carry != 0) {
0060 product[i + 2] += carry;
0061 }
0062 }
0063
0064 digit_t quotient[2];
0065
0066 if(m == 0) {
0067 div_t result = {
0068 ((::boost::uintmax_t(product[3]) << bits) | product[2]),
0069 ((::boost::uintmax_t(product[1]) << bits) | product[0]) >> shift,
0070 };
0071 return result;
0072 }
0073
0074
0075 for(int i = 3; i >= 2; --i) {
0076 ::boost::uintmax_t temp =
0077 ::boost::uintmax_t(product[i]) << bits | product[i - 1];
0078
0079 digit_t q = digit_t((product[i] == m_[1]) ? mask : temp / m_[1]);
0080
0081 ::boost::uintmax_t rem =
0082 ((temp - ::boost::uintmax_t(q) * m_[1]) << bits) + product[i - 2];
0083
0084 ::boost::uintmax_t diff = m_[0] * ::boost::uintmax_t(q);
0085
0086 int error = 0;
0087 if(diff > rem) {
0088 if(diff - rem > m) {
0089 error = 2;
0090 } else {
0091 error = 1;
0092 }
0093 }
0094 q -= error;
0095 rem = rem + error * m - diff;
0096
0097 quotient[i - 2] = q;
0098 product[i] = 0;
0099 product[i-1] = static_cast<digit_t>((rem >> bits) & mask);
0100 product[i-2] = static_cast<digit_t>(rem & mask);
0101 }
0102
0103 div_t result = {
0104 ((::boost::uintmax_t(quotient[1]) << bits) | quotient[0]),
0105 ((::boost::uintmax_t(product[1]) << bits) | product[0]) >> shift,
0106 };
0107 return result;
0108 }
0109
0110 inline boost::uintmax_t muldiv(boost::uintmax_t a, boost::uintmax_t b, boost::uintmax_t m)
0111 { return detail::muldivmod(a, b, m).quotient; }
0112
0113 inline boost::uintmax_t mulmod(boost::uintmax_t a, boost::uintmax_t b, boost::uintmax_t m)
0114 { return detail::muldivmod(a, b, m).remainder; }
0115
0116 }
0117 }
0118 }
0119
0120 #include <boost/random/detail/enable_warnings.hpp>
0121
0122 #endif