annotate DEPENDENCIES/generic/include/boost/random/detail/large_arithmetic.hpp @ 133:4acb5d8d80b6 tip

Don't fail environmental check if README.md exists (but .txt and no-suffix don't)
author Chris Cannam
date Tue, 30 Jul 2019 12:25:44 +0100
parents c530137014c0
children
rev   line source
Chris@16 1 /* boost random/detail/large_arithmetic.hpp header file
Chris@16 2 *
Chris@16 3 * Copyright Steven Watanabe 2011
Chris@16 4 * Distributed under the Boost Software License, Version 1.0. (See
Chris@16 5 * accompanying file LICENSE_1_0.txt or copy at
Chris@16 6 * http://www.boost.org/LICENSE_1_0.txt)
Chris@16 7 *
Chris@16 8 * See http://www.boost.org for most recent version including documentation.
Chris@16 9 *
Chris@101 10 * $Id$
Chris@16 11 */
Chris@16 12
Chris@16 13 #ifndef BOOST_RANDOM_DETAIL_LARGE_ARITHMETIC_HPP
Chris@16 14 #define BOOST_RANDOM_DETAIL_LARGE_ARITHMETIC_HPP
Chris@16 15
Chris@16 16 #include <boost/cstdint.hpp>
Chris@16 17 #include <boost/integer.hpp>
Chris@16 18 #include <boost/limits.hpp>
Chris@16 19 #include <boost/random/detail/integer_log2.hpp>
Chris@16 20
Chris@16 21 #include <boost/random/detail/disable_warnings.hpp>
Chris@16 22
Chris@16 23 namespace boost {
Chris@16 24 namespace random {
Chris@16 25 namespace detail {
Chris@16 26
Chris@16 27 struct div_t {
Chris@16 28 boost::uintmax_t quotient;
Chris@16 29 boost::uintmax_t remainder;
Chris@16 30 };
Chris@16 31
Chris@16 32 inline div_t muldivmod(boost::uintmax_t a, boost::uintmax_t b, boost::uintmax_t m)
Chris@16 33 {
Chris@101 34 const int bits =
Chris@16 35 ::std::numeric_limits< ::boost::uintmax_t>::digits / 2;
Chris@101 36 const ::boost::uintmax_t mask = (::boost::uintmax_t(1) << bits) - 1;
Chris@16 37 typedef ::boost::uint_t<bits>::fast digit_t;
Chris@16 38
Chris@16 39 int shift = std::numeric_limits< ::boost::uintmax_t>::digits - 1
Chris@16 40 - detail::integer_log2(m);
Chris@16 41
Chris@16 42 a <<= shift;
Chris@16 43 m <<= shift;
Chris@16 44
Chris@16 45 digit_t product[4] = { 0, 0, 0, 0 };
Chris@16 46 digit_t a_[2] = { digit_t(a & mask), digit_t((a >> bits) & mask) };
Chris@16 47 digit_t b_[2] = { digit_t(b & mask), digit_t((b >> bits) & mask) };
Chris@16 48 digit_t m_[2] = { digit_t(m & mask), digit_t((m >> bits) & mask) };
Chris@16 49
Chris@16 50 // multiply a * b
Chris@16 51 for(int i = 0; i < 2; ++i) {
Chris@16 52 digit_t carry = 0;
Chris@16 53 for(int j = 0; j < 2; ++j) {
Chris@16 54 ::boost::uint64_t temp = ::boost::uintmax_t(a_[i]) * b_[j] +
Chris@16 55 carry + product[i + j];
Chris@16 56 product[i + j] = digit_t(temp & mask);
Chris@16 57 carry = digit_t(temp >> bits);
Chris@16 58 }
Chris@16 59 if(carry != 0) {
Chris@16 60 product[i + 2] += carry;
Chris@16 61 }
Chris@16 62 }
Chris@16 63
Chris@16 64 digit_t quotient[2];
Chris@16 65
Chris@16 66 if(m == 0) {
Chris@16 67 div_t result = {
Chris@16 68 ((::boost::uintmax_t(product[3]) << bits) | product[2]),
Chris@16 69 ((::boost::uintmax_t(product[1]) << bits) | product[0]) >> shift,
Chris@16 70 };
Chris@16 71 return result;
Chris@16 72 }
Chris@16 73
Chris@16 74 // divide product / m
Chris@16 75 for(int i = 3; i >= 2; --i) {
Chris@16 76 ::boost::uintmax_t temp =
Chris@16 77 ::boost::uintmax_t(product[i]) << bits | product[i - 1];
Chris@16 78
Chris@16 79 digit_t q = digit_t((product[i] == m_[1]) ? mask : temp / m_[1]);
Chris@16 80
Chris@16 81 ::boost::uintmax_t rem =
Chris@16 82 ((temp - ::boost::uintmax_t(q) * m_[1]) << bits) + product[i - 2];
Chris@16 83
Chris@16 84 ::boost::uintmax_t diff = m_[0] * ::boost::uintmax_t(q);
Chris@16 85
Chris@16 86 int error = 0;
Chris@16 87 if(diff > rem) {
Chris@16 88 if(diff - rem > m) {
Chris@16 89 error = 2;
Chris@16 90 } else {
Chris@16 91 error = 1;
Chris@16 92 }
Chris@16 93 }
Chris@16 94 q -= error;
Chris@16 95 rem = rem + error * m - diff;
Chris@16 96
Chris@16 97 quotient[i - 2] = q;
Chris@16 98 product[i] = 0;
Chris@101 99 product[i-1] = static_cast<digit_t>((rem >> bits) & mask);
Chris@101 100 product[i-2] = static_cast<digit_t>(rem & mask);
Chris@16 101 }
Chris@16 102
Chris@16 103 div_t result = {
Chris@16 104 ((::boost::uintmax_t(quotient[1]) << bits) | quotient[0]),
Chris@16 105 ((::boost::uintmax_t(product[1]) << bits) | product[0]) >> shift,
Chris@16 106 };
Chris@16 107 return result;
Chris@16 108 }
Chris@16 109
Chris@16 110 inline boost::uintmax_t muldiv(boost::uintmax_t a, boost::uintmax_t b, boost::uintmax_t m)
Chris@16 111 { return detail::muldivmod(a, b, m).quotient; }
Chris@16 112
Chris@16 113 inline boost::uintmax_t mulmod(boost::uintmax_t a, boost::uintmax_t b, boost::uintmax_t m)
Chris@16 114 { return detail::muldivmod(a, b, m).remainder; }
Chris@16 115
Chris@16 116 } // namespace detail
Chris@16 117 } // namespace random
Chris@16 118 } // namespace boost
Chris@16 119
Chris@16 120 #include <boost/random/detail/enable_warnings.hpp>
Chris@16 121
Chris@16 122 #endif // BOOST_RANDOM_DETAIL_LARGE_ARITHMETIC_HPP