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