Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-30 09:59:09

0001 /* boost random/detail/large_arithmetic.hpp header file
0002  *
0003  * Copyright Steven Watanabe 2011
0004  * Distributed under the Boost Software License, Version 1.0. (See
0005  * accompanying file LICENSE_1_0.txt or copy at
0006  * http://www.boost.org/LICENSE_1_0.txt)
0007  *
0008  * See http://www.boost.org for most recent version including documentation.
0009  *
0010  * $Id$
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     // multiply a * b
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     // divide product / m
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 } // namespace detail
0117 } // namespace random
0118 } // namespace boost
0119 
0120 #include <boost/random/detail/enable_warnings.hpp>
0121 
0122 #endif // BOOST_RANDOM_DETAIL_LARGE_ARITHMETIC_HPP