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
|