]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 |