X-Git-Url: https://git.donarmstrong.com/?p=rsem.git;a=blobdiff_plain;f=boost%2Frandom%2Fdetail%2Flarge_arithmetic.hpp;fp=boost%2Frandom%2Fdetail%2Flarge_arithmetic.hpp;h=24177dcdda43379bcbd0e30f9ce8b03f78cb59e4;hp=0000000000000000000000000000000000000000;hb=2d71eb92104693ca9baa5a2e1c23eeca776d8fd3;hpb=da57529b92adbb7ae74a89861cb39fb35ac7c62d diff --git a/boost/random/detail/large_arithmetic.hpp b/boost/random/detail/large_arithmetic.hpp new file mode 100644 index 0000000..24177dc --- /dev/null +++ b/boost/random/detail/large_arithmetic.hpp @@ -0,0 +1,122 @@ +/* boost random/detail/large_arithmetic.hpp header file + * + * Copyright Steven Watanabe 2011 + * 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_0.txt) + * + * See http://www.boost.org for most recent version including documentation. + * + * $Id: large_arithmetic.hpp 71018 2011-04-05 21:27:52Z steven_watanabe $ + */ + +#ifndef BOOST_RANDOM_DETAIL_LARGE_ARITHMETIC_HPP +#define BOOST_RANDOM_DETAIL_LARGE_ARITHMETIC_HPP + +#include +#include +#include +#include + +#include + +namespace boost { +namespace random { +namespace detail { + +struct div_t { + boost::uintmax_t quotient; + boost::uintmax_t remainder; +}; + +inline div_t muldivmod(boost::uintmax_t a, boost::uintmax_t b, boost::uintmax_t m) +{ + static const int bits = + ::std::numeric_limits< ::boost::uintmax_t>::digits / 2; + static const ::boost::uintmax_t mask = (::boost::uintmax_t(1) << bits) - 1; + typedef ::boost::uint_t::fast digit_t; + + int shift = std::numeric_limits< ::boost::uintmax_t>::digits - 1 + - detail::integer_log2(m); + + a <<= shift; + m <<= shift; + + digit_t product[4] = { 0, 0, 0, 0 }; + digit_t a_[2] = { digit_t(a & mask), digit_t((a >> bits) & mask) }; + digit_t b_[2] = { digit_t(b & mask), digit_t((b >> bits) & mask) }; + digit_t m_[2] = { digit_t(m & mask), digit_t((m >> bits) & mask) }; + + // multiply a * b + for(int i = 0; i < 2; ++i) { + digit_t carry = 0; + for(int j = 0; j < 2; ++j) { + ::boost::uint64_t temp = ::boost::uintmax_t(a_[i]) * b_[j] + + carry + product[i + j]; + product[i + j] = digit_t(temp & mask); + carry = digit_t(temp >> bits); + } + if(carry != 0) { + product[i + 2] += carry; + } + } + + digit_t quotient[2]; + + if(m == 0) { + div_t result = { + ((::boost::uintmax_t(product[3]) << bits) | product[2]), + ((::boost::uintmax_t(product[1]) << bits) | product[0]) >> shift, + }; + return result; + } + + // divide product / m + for(int i = 3; i >= 2; --i) { + ::boost::uintmax_t temp = + ::boost::uintmax_t(product[i]) << bits | product[i - 1]; + + digit_t q = digit_t((product[i] == m_[1]) ? mask : temp / m_[1]); + + ::boost::uintmax_t rem = + ((temp - ::boost::uintmax_t(q) * m_[1]) << bits) + product[i - 2]; + + ::boost::uintmax_t diff = m_[0] * ::boost::uintmax_t(q); + + int error = 0; + if(diff > rem) { + if(diff - rem > m) { + error = 2; + } else { + error = 1; + } + } + q -= error; + rem = rem + error * m - diff; + + quotient[i - 2] = q; + product[i] = 0; + product[i-1] = (rem >> bits) & mask; + product[i-2] = rem & mask; + } + + div_t result = { + ((::boost::uintmax_t(quotient[1]) << bits) | quotient[0]), + ((::boost::uintmax_t(product[1]) << bits) | product[0]) >> shift, + }; + return result; +} + +inline boost::uintmax_t muldiv(boost::uintmax_t a, boost::uintmax_t b, boost::uintmax_t m) +{ return detail::muldivmod(a, b, m).quotient; } + +inline boost::uintmax_t mulmod(boost::uintmax_t a, boost::uintmax_t b, boost::uintmax_t m) +{ return detail::muldivmod(a, b, m).remainder; } + +} // namespace detail +} // namespace random +} // namespace boost + +#include + +#endif // BOOST_RANDOM_DETAIL_LARGE_ARITHMETIC_HPP