1 | /* boost random/detail/const_mod.hpp header file |
2 | * |
3 | * Copyright Jens Maurer 2000-2001 |
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 | * Revision history |
13 | * 2001-02-18 moved to individual header files |
14 | */ |
15 | |
16 | #ifndef BOOST_RANDOM_CONST_MOD_HPP |
17 | #define BOOST_RANDOM_CONST_MOD_HPP |
18 | |
19 | #include <boost/assert.hpp> |
20 | #include <boost/static_assert.hpp> |
21 | #include <boost/integer_traits.hpp> |
22 | #include <boost/type_traits/make_unsigned.hpp> |
23 | #include <boost/random/detail/large_arithmetic.hpp> |
24 | |
25 | #include <boost/random/detail/disable_warnings.hpp> |
26 | |
27 | namespace boost { |
28 | namespace random { |
29 | |
30 | template<class IntType, IntType m> |
31 | class const_mod |
32 | { |
33 | public: |
34 | static IntType apply(IntType x) |
35 | { |
36 | if(((unsigned_m() - 1) & unsigned_m()) == 0) |
37 | return (unsigned_type(x)) & (unsigned_m() - 1); |
38 | else { |
39 | IntType suppress_warnings = (m == 0); |
40 | BOOST_ASSERT(suppress_warnings == 0); |
41 | return x % (m + suppress_warnings); |
42 | } |
43 | } |
44 | |
45 | static IntType add(IntType x, IntType c) |
46 | { |
47 | if(((unsigned_m() - 1) & unsigned_m()) == 0) |
48 | return (unsigned_type(x) + unsigned_type(c)) & (unsigned_m() - 1); |
49 | else if(c == 0) |
50 | return x; |
51 | else if(x < m - c) |
52 | return x + c; |
53 | else |
54 | return x - (m - c); |
55 | } |
56 | |
57 | static IntType mult(IntType a, IntType x) |
58 | { |
59 | if(((unsigned_m() - 1) & unsigned_m()) == 0) |
60 | return unsigned_type(a) * unsigned_type(x) & (unsigned_m() - 1); |
61 | else if(a == 0) |
62 | return 0; |
63 | else if(a == 1) |
64 | return x; |
65 | else if(m <= traits::const_max/a) // i.e. a*m <= max |
66 | return mult_small(a, x); |
67 | else if(traits::is_signed && (m%a < m/a)) |
68 | return mult_schrage(a, x); |
69 | else |
70 | return mult_general(a, x); |
71 | } |
72 | |
73 | static IntType mult_add(IntType a, IntType x, IntType c) |
74 | { |
75 | if(((unsigned_m() - 1) & unsigned_m()) == 0) |
76 | return (unsigned_type(a) * unsigned_type(x) + unsigned_type(c)) & (unsigned_m() - 1); |
77 | else if(a == 0) |
78 | return c; |
79 | else if(m <= (traits::const_max-c)/a) { // i.e. a*m+c <= max |
80 | IntType suppress_warnings = (m == 0); |
81 | BOOST_ASSERT(suppress_warnings == 0); |
82 | return (a*x+c) % (m + suppress_warnings); |
83 | } else |
84 | return add(mult(a, x), c); |
85 | } |
86 | |
87 | static IntType pow(IntType a, boost::uintmax_t exponent) |
88 | { |
89 | IntType result = 1; |
90 | while(exponent != 0) { |
91 | if(exponent % 2 == 1) { |
92 | result = mult(result, a); |
93 | } |
94 | a = mult(a, a); |
95 | exponent /= 2; |
96 | } |
97 | return result; |
98 | } |
99 | |
100 | static IntType invert(IntType x) |
101 | { return x == 0 ? 0 : (m == 0? invert_euclidian0(x) : invert_euclidian(x)); } |
102 | |
103 | private: |
104 | typedef integer_traits<IntType> traits; |
105 | typedef typename make_unsigned<IntType>::type unsigned_type; |
106 | |
107 | const_mod(); // don't instantiate |
108 | |
109 | static IntType mult_small(IntType a, IntType x) |
110 | { |
111 | IntType suppress_warnings = (m == 0); |
112 | BOOST_ASSERT(suppress_warnings == 0); |
113 | return a*x % (m + suppress_warnings); |
114 | } |
115 | |
116 | static IntType mult_schrage(IntType a, IntType value) |
117 | { |
118 | const IntType q = m / a; |
119 | const IntType r = m % a; |
120 | |
121 | BOOST_ASSERT(r < q); // check that overflow cannot happen |
122 | |
123 | return sub(a*(value%q), r*(value/q)); |
124 | } |
125 | |
126 | static IntType mult_general(IntType a, IntType b) |
127 | { |
128 | IntType suppress_warnings = (m == 0); |
129 | BOOST_ASSERT(suppress_warnings == 0); |
130 | IntType modulus = m + suppress_warnings; |
131 | BOOST_ASSERT(modulus == m); |
132 | if(::boost::uintmax_t(modulus) <= |
133 | (::std::numeric_limits< ::boost::uintmax_t>::max)() / modulus) |
134 | { |
135 | return static_cast<IntType>(boost::uintmax_t(a) * b % modulus); |
136 | } else { |
137 | return static_cast<IntType>(detail::mulmod(a, b, modulus)); |
138 | } |
139 | } |
140 | |
141 | static IntType sub(IntType a, IntType b) |
142 | { |
143 | if(a < b) |
144 | return m - (b - a); |
145 | else |
146 | return a - b; |
147 | } |
148 | |
149 | static unsigned_type unsigned_m() |
150 | { |
151 | if(m == 0) { |
152 | return unsigned_type((std::numeric_limits<IntType>::max)()) + 1; |
153 | } else { |
154 | return unsigned_type(m); |
155 | } |
156 | } |
157 | |
158 | // invert c in the finite field (mod m) (m must be prime) |
159 | static IntType invert_euclidian(IntType c) |
160 | { |
161 | // we are interested in the gcd factor for c, because this is our inverse |
162 | BOOST_ASSERT(c > 0); |
163 | IntType l1 = 0; |
164 | IntType l2 = 1; |
165 | IntType n = c; |
166 | IntType p = m; |
167 | for(;;) { |
168 | IntType q = p / n; |
169 | l1 += q * l2; |
170 | p -= q * n; |
171 | if(p == 0) |
172 | return l2; |
173 | IntType q2 = n / p; |
174 | l2 += q2 * l1; |
175 | n -= q2 * p; |
176 | if(n == 0) |
177 | return m - l1; |
178 | } |
179 | } |
180 | |
181 | // invert c in the finite field (mod m) (c must be relatively prime to m) |
182 | static IntType invert_euclidian0(IntType c) |
183 | { |
184 | // we are interested in the gcd factor for c, because this is our inverse |
185 | BOOST_ASSERT(c > 0); |
186 | if(c == 1) return 1; |
187 | IntType l1 = 0; |
188 | IntType l2 = 1; |
189 | IntType n = c; |
190 | IntType p = m; |
191 | IntType max = (std::numeric_limits<IntType>::max)(); |
192 | IntType q = max / n; |
193 | BOOST_ASSERT(max % n != n - 1 && "c must be relatively prime to m." ); |
194 | l1 += q * l2; |
195 | p = max - q * n + 1; |
196 | for(;;) { |
197 | if(p == 0) |
198 | return l2; |
199 | IntType q2 = n / p; |
200 | l2 += q2 * l1; |
201 | n -= q2 * p; |
202 | if(n == 0) |
203 | return m - l1; |
204 | q = p / n; |
205 | l1 += q * l2; |
206 | p -= q * n; |
207 | } |
208 | } |
209 | }; |
210 | |
211 | } // namespace random |
212 | } // namespace boost |
213 | |
214 | #include <boost/random/detail/enable_warnings.hpp> |
215 | |
216 | #endif // BOOST_RANDOM_CONST_MOD_HPP |
217 | |