williamr@2
|
1 |
/* boost random/detail/const_mod.hpp header file
|
williamr@2
|
2 |
*
|
williamr@2
|
3 |
* Copyright Jens Maurer 2000-2001
|
williamr@2
|
4 |
* Distributed under the Boost Software License, Version 1.0. (See
|
williamr@2
|
5 |
* accompanying file LICENSE_1_0.txt or copy at
|
williamr@2
|
6 |
* http://www.boost.org/LICENSE_1_0.txt)
|
williamr@2
|
7 |
*
|
williamr@2
|
8 |
* See http://www.boost.org for most recent version including documentation.
|
williamr@2
|
9 |
*
|
williamr@2
|
10 |
* $Id: const_mod.hpp,v 1.8 2004/07/27 03:43:32 dgregor Exp $
|
williamr@2
|
11 |
*
|
williamr@2
|
12 |
* Revision history
|
williamr@2
|
13 |
* 2001-02-18 moved to individual header files
|
williamr@2
|
14 |
*/
|
williamr@2
|
15 |
|
williamr@2
|
16 |
#ifndef BOOST_RANDOM_CONST_MOD_HPP
|
williamr@2
|
17 |
#define BOOST_RANDOM_CONST_MOD_HPP
|
williamr@2
|
18 |
|
williamr@2
|
19 |
#include <cassert>
|
williamr@2
|
20 |
#include <boost/static_assert.hpp>
|
williamr@2
|
21 |
#include <boost/cstdint.hpp>
|
williamr@2
|
22 |
#include <boost/integer_traits.hpp>
|
williamr@2
|
23 |
#include <boost/detail/workaround.hpp>
|
williamr@2
|
24 |
|
williamr@2
|
25 |
namespace boost {
|
williamr@2
|
26 |
namespace random {
|
williamr@2
|
27 |
|
williamr@2
|
28 |
/*
|
williamr@2
|
29 |
* Some random number generators require modular arithmetic. Put
|
williamr@2
|
30 |
* everything we need here.
|
williamr@2
|
31 |
* IntType must be an integral type.
|
williamr@2
|
32 |
*/
|
williamr@2
|
33 |
|
williamr@2
|
34 |
namespace detail {
|
williamr@2
|
35 |
|
williamr@2
|
36 |
template<bool is_signed>
|
williamr@2
|
37 |
struct do_add
|
williamr@2
|
38 |
{ };
|
williamr@2
|
39 |
|
williamr@2
|
40 |
template<>
|
williamr@2
|
41 |
struct do_add<true>
|
williamr@2
|
42 |
{
|
williamr@2
|
43 |
template<class IntType>
|
williamr@2
|
44 |
static IntType add(IntType m, IntType x, IntType c)
|
williamr@2
|
45 |
{
|
williamr@2
|
46 |
x += (c-m);
|
williamr@2
|
47 |
if(x < 0)
|
williamr@2
|
48 |
x += m;
|
williamr@2
|
49 |
return x;
|
williamr@2
|
50 |
}
|
williamr@2
|
51 |
};
|
williamr@2
|
52 |
|
williamr@2
|
53 |
template<>
|
williamr@2
|
54 |
struct do_add<false>
|
williamr@2
|
55 |
{
|
williamr@2
|
56 |
template<class IntType>
|
williamr@2
|
57 |
static IntType add(IntType, IntType, IntType)
|
williamr@2
|
58 |
{
|
williamr@2
|
59 |
// difficult
|
williamr@2
|
60 |
assert(!"const_mod::add with c too large");
|
williamr@2
|
61 |
return 0;
|
williamr@2
|
62 |
}
|
williamr@2
|
63 |
};
|
williamr@2
|
64 |
} // namespace detail
|
williamr@2
|
65 |
|
williamr@2
|
66 |
#if !(defined(__BORLANDC__) && (__BORLANDC__ == 0x560))
|
williamr@2
|
67 |
|
williamr@2
|
68 |
template<class IntType, IntType m>
|
williamr@2
|
69 |
class const_mod
|
williamr@2
|
70 |
{
|
williamr@2
|
71 |
public:
|
williamr@2
|
72 |
static IntType add(IntType x, IntType c)
|
williamr@2
|
73 |
{
|
williamr@2
|
74 |
if(c == 0)
|
williamr@2
|
75 |
return x;
|
williamr@2
|
76 |
else if(c <= traits::const_max - m) // i.e. m+c < max
|
williamr@2
|
77 |
return add_small(x, c);
|
williamr@2
|
78 |
else
|
williamr@2
|
79 |
return detail::do_add<traits::is_signed>::add(m, x, c);
|
williamr@2
|
80 |
}
|
williamr@2
|
81 |
|
williamr@2
|
82 |
static IntType mult(IntType a, IntType x)
|
williamr@2
|
83 |
{
|
williamr@2
|
84 |
if(a == 1)
|
williamr@2
|
85 |
return x;
|
williamr@2
|
86 |
else if(m <= traits::const_max/a) // i.e. a*m <= max
|
williamr@2
|
87 |
return mult_small(a, x);
|
williamr@2
|
88 |
else if(traits::is_signed && (m%a < m/a))
|
williamr@2
|
89 |
return mult_schrage(a, x);
|
williamr@2
|
90 |
else {
|
williamr@2
|
91 |
// difficult
|
williamr@2
|
92 |
assert(!"const_mod::mult with a too large");
|
williamr@2
|
93 |
return 0;
|
williamr@2
|
94 |
}
|
williamr@2
|
95 |
}
|
williamr@2
|
96 |
|
williamr@2
|
97 |
static IntType mult_add(IntType a, IntType x, IntType c)
|
williamr@2
|
98 |
{
|
williamr@2
|
99 |
if(m <= (traits::const_max-c)/a) // i.e. a*m+c <= max
|
williamr@2
|
100 |
return (a*x+c) % m;
|
williamr@2
|
101 |
else
|
williamr@2
|
102 |
return add(mult(a, x), c);
|
williamr@2
|
103 |
}
|
williamr@2
|
104 |
|
williamr@2
|
105 |
static IntType invert(IntType x)
|
williamr@2
|
106 |
{ return x == 0 ? 0 : invert_euclidian(x); }
|
williamr@2
|
107 |
|
williamr@2
|
108 |
private:
|
williamr@2
|
109 |
typedef integer_traits<IntType> traits;
|
williamr@2
|
110 |
|
williamr@2
|
111 |
const_mod(); // don't instantiate
|
williamr@2
|
112 |
|
williamr@2
|
113 |
static IntType add_small(IntType x, IntType c)
|
williamr@2
|
114 |
{
|
williamr@2
|
115 |
x += c;
|
williamr@2
|
116 |
if(x >= m)
|
williamr@2
|
117 |
x -= m;
|
williamr@2
|
118 |
return x;
|
williamr@2
|
119 |
}
|
williamr@2
|
120 |
|
williamr@2
|
121 |
static IntType mult_small(IntType a, IntType x)
|
williamr@2
|
122 |
{
|
williamr@2
|
123 |
return a*x % m;
|
williamr@2
|
124 |
}
|
williamr@2
|
125 |
|
williamr@2
|
126 |
static IntType mult_schrage(IntType a, IntType value)
|
williamr@2
|
127 |
{
|
williamr@2
|
128 |
const IntType q = m / a;
|
williamr@2
|
129 |
const IntType r = m % a;
|
williamr@2
|
130 |
|
williamr@2
|
131 |
assert(r < q); // check that overflow cannot happen
|
williamr@2
|
132 |
|
williamr@2
|
133 |
value = a*(value%q) - r*(value/q);
|
williamr@2
|
134 |
// An optimizer bug in the SGI MIPSpro 7.3.1.x compiler requires this
|
williamr@2
|
135 |
// convoluted formulation of the loop (Synge Todo)
|
williamr@2
|
136 |
for(;;) {
|
williamr@2
|
137 |
if (value > 0)
|
williamr@2
|
138 |
break;
|
williamr@2
|
139 |
value += m;
|
williamr@2
|
140 |
}
|
williamr@2
|
141 |
return value;
|
williamr@2
|
142 |
}
|
williamr@2
|
143 |
|
williamr@2
|
144 |
// invert c in the finite field (mod m) (m must be prime)
|
williamr@2
|
145 |
static IntType invert_euclidian(IntType c)
|
williamr@2
|
146 |
{
|
williamr@2
|
147 |
// we are interested in the gcd factor for c, because this is our inverse
|
williamr@2
|
148 |
BOOST_STATIC_ASSERT(m > 0);
|
williamr@2
|
149 |
#if BOOST_WORKAROUND(__MWERKS__, BOOST_TESTED_AT(0x3003))
|
williamr@2
|
150 |
assert(boost::integer_traits<IntType>::is_signed);
|
williamr@2
|
151 |
#elif !defined(BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS)
|
williamr@2
|
152 |
BOOST_STATIC_ASSERT(boost::integer_traits<IntType>::is_signed);
|
williamr@2
|
153 |
#endif
|
williamr@2
|
154 |
assert(c > 0);
|
williamr@2
|
155 |
IntType l1 = 0;
|
williamr@2
|
156 |
IntType l2 = 1;
|
williamr@2
|
157 |
IntType n = c;
|
williamr@2
|
158 |
IntType p = m;
|
williamr@2
|
159 |
for(;;) {
|
williamr@2
|
160 |
IntType q = p / n;
|
williamr@2
|
161 |
l1 -= q * l2; // this requires a signed IntType!
|
williamr@2
|
162 |
p -= q * n;
|
williamr@2
|
163 |
if(p == 0)
|
williamr@2
|
164 |
return (l2 < 1 ? l2 + m : l2);
|
williamr@2
|
165 |
IntType q2 = n / p;
|
williamr@2
|
166 |
l2 -= q2 * l1;
|
williamr@2
|
167 |
n -= q2 * p;
|
williamr@2
|
168 |
if(n == 0)
|
williamr@2
|
169 |
return (l1 < 1 ? l1 + m : l1);
|
williamr@2
|
170 |
}
|
williamr@2
|
171 |
}
|
williamr@2
|
172 |
};
|
williamr@2
|
173 |
|
williamr@2
|
174 |
// The modulus is exactly the word size: rely on machine overflow handling.
|
williamr@2
|
175 |
// Due to a GCC bug, we cannot partially specialize in the presence of
|
williamr@2
|
176 |
// template value parameters.
|
williamr@2
|
177 |
template<>
|
williamr@2
|
178 |
class const_mod<unsigned int, 0>
|
williamr@2
|
179 |
{
|
williamr@2
|
180 |
typedef unsigned int IntType;
|
williamr@2
|
181 |
public:
|
williamr@2
|
182 |
static IntType add(IntType x, IntType c) { return x+c; }
|
williamr@2
|
183 |
static IntType mult(IntType a, IntType x) { return a*x; }
|
williamr@2
|
184 |
static IntType mult_add(IntType a, IntType x, IntType c) { return a*x+c; }
|
williamr@2
|
185 |
|
williamr@2
|
186 |
// m is not prime, thus invert is not useful
|
williamr@2
|
187 |
private: // don't instantiate
|
williamr@2
|
188 |
const_mod();
|
williamr@2
|
189 |
};
|
williamr@2
|
190 |
|
williamr@2
|
191 |
template<>
|
williamr@2
|
192 |
class const_mod<unsigned long, 0>
|
williamr@2
|
193 |
{
|
williamr@2
|
194 |
typedef unsigned long IntType;
|
williamr@2
|
195 |
public:
|
williamr@2
|
196 |
static IntType add(IntType x, IntType c) { return x+c; }
|
williamr@2
|
197 |
static IntType mult(IntType a, IntType x) { return a*x; }
|
williamr@2
|
198 |
static IntType mult_add(IntType a, IntType x, IntType c) { return a*x+c; }
|
williamr@2
|
199 |
|
williamr@2
|
200 |
// m is not prime, thus invert is not useful
|
williamr@2
|
201 |
private: // don't instantiate
|
williamr@2
|
202 |
const_mod();
|
williamr@2
|
203 |
};
|
williamr@2
|
204 |
|
williamr@2
|
205 |
// the modulus is some power of 2: rely partly on machine overflow handling
|
williamr@2
|
206 |
// we only specialize for rand48 at the moment
|
williamr@2
|
207 |
#ifndef BOOST_NO_INT64_T
|
williamr@2
|
208 |
template<>
|
williamr@2
|
209 |
class const_mod<uint64_t, uint64_t(1) << 48>
|
williamr@2
|
210 |
{
|
williamr@2
|
211 |
typedef uint64_t IntType;
|
williamr@2
|
212 |
public:
|
williamr@2
|
213 |
static IntType add(IntType x, IntType c) { return c == 0 ? x : mod(x+c); }
|
williamr@2
|
214 |
static IntType mult(IntType a, IntType x) { return mod(a*x); }
|
williamr@2
|
215 |
static IntType mult_add(IntType a, IntType x, IntType c)
|
williamr@2
|
216 |
{ return mod(a*x+c); }
|
williamr@2
|
217 |
static IntType mod(IntType x) { return x &= ((uint64_t(1) << 48)-1); }
|
williamr@2
|
218 |
|
williamr@2
|
219 |
// m is not prime, thus invert is not useful
|
williamr@2
|
220 |
private: // don't instantiate
|
williamr@2
|
221 |
const_mod();
|
williamr@2
|
222 |
};
|
williamr@2
|
223 |
#endif /* !BOOST_NO_INT64_T */
|
williamr@2
|
224 |
|
williamr@2
|
225 |
#else
|
williamr@2
|
226 |
|
williamr@2
|
227 |
//
|
williamr@2
|
228 |
// for some reason Borland C++ Builder 6 has problems with
|
williamr@2
|
229 |
// the full specialisations of const_mod, define a generic version
|
williamr@2
|
230 |
// instead, the compiler will optimise away the const-if statements:
|
williamr@2
|
231 |
//
|
williamr@2
|
232 |
|
williamr@2
|
233 |
template<class IntType, IntType m>
|
williamr@2
|
234 |
class const_mod
|
williamr@2
|
235 |
{
|
williamr@2
|
236 |
public:
|
williamr@2
|
237 |
static IntType add(IntType x, IntType c)
|
williamr@2
|
238 |
{
|
williamr@2
|
239 |
if(0 == m)
|
williamr@2
|
240 |
{
|
williamr@2
|
241 |
return x+c;
|
williamr@2
|
242 |
}
|
williamr@2
|
243 |
else
|
williamr@2
|
244 |
{
|
williamr@2
|
245 |
if(c == 0)
|
williamr@2
|
246 |
return x;
|
williamr@2
|
247 |
else if(c <= traits::const_max - m) // i.e. m+c < max
|
williamr@2
|
248 |
return add_small(x, c);
|
williamr@2
|
249 |
else
|
williamr@2
|
250 |
return detail::do_add<traits::is_signed>::add(m, x, c);
|
williamr@2
|
251 |
}
|
williamr@2
|
252 |
}
|
williamr@2
|
253 |
|
williamr@2
|
254 |
static IntType mult(IntType a, IntType x)
|
williamr@2
|
255 |
{
|
williamr@2
|
256 |
if(x == 0)
|
williamr@2
|
257 |
{
|
williamr@2
|
258 |
return a*x;
|
williamr@2
|
259 |
}
|
williamr@2
|
260 |
else
|
williamr@2
|
261 |
{
|
williamr@2
|
262 |
if(a == 1)
|
williamr@2
|
263 |
return x;
|
williamr@2
|
264 |
else if(m <= traits::const_max/a) // i.e. a*m <= max
|
williamr@2
|
265 |
return mult_small(a, x);
|
williamr@2
|
266 |
else if(traits::is_signed && (m%a < m/a))
|
williamr@2
|
267 |
return mult_schrage(a, x);
|
williamr@2
|
268 |
else {
|
williamr@2
|
269 |
// difficult
|
williamr@2
|
270 |
assert(!"const_mod::mult with a too large");
|
williamr@2
|
271 |
return 0;
|
williamr@2
|
272 |
}
|
williamr@2
|
273 |
}
|
williamr@2
|
274 |
}
|
williamr@2
|
275 |
|
williamr@2
|
276 |
static IntType mult_add(IntType a, IntType x, IntType c)
|
williamr@2
|
277 |
{
|
williamr@2
|
278 |
if(m == 0)
|
williamr@2
|
279 |
{
|
williamr@2
|
280 |
return a*x+c;
|
williamr@2
|
281 |
}
|
williamr@2
|
282 |
else
|
williamr@2
|
283 |
{
|
williamr@2
|
284 |
if(m <= (traits::const_max-c)/a) // i.e. a*m+c <= max
|
williamr@2
|
285 |
return (a*x+c) % m;
|
williamr@2
|
286 |
else
|
williamr@2
|
287 |
return add(mult(a, x), c);
|
williamr@2
|
288 |
}
|
williamr@2
|
289 |
}
|
williamr@2
|
290 |
|
williamr@2
|
291 |
static IntType invert(IntType x)
|
williamr@2
|
292 |
{ return x == 0 ? 0 : invert_euclidian(x); }
|
williamr@2
|
293 |
|
williamr@2
|
294 |
private:
|
williamr@2
|
295 |
typedef integer_traits<IntType> traits;
|
williamr@2
|
296 |
|
williamr@2
|
297 |
const_mod(); // don't instantiate
|
williamr@2
|
298 |
|
williamr@2
|
299 |
static IntType add_small(IntType x, IntType c)
|
williamr@2
|
300 |
{
|
williamr@2
|
301 |
x += c;
|
williamr@2
|
302 |
if(x >= m)
|
williamr@2
|
303 |
x -= m;
|
williamr@2
|
304 |
return x;
|
williamr@2
|
305 |
}
|
williamr@2
|
306 |
|
williamr@2
|
307 |
static IntType mult_small(IntType a, IntType x)
|
williamr@2
|
308 |
{
|
williamr@2
|
309 |
return a*x % m;
|
williamr@2
|
310 |
}
|
williamr@2
|
311 |
|
williamr@2
|
312 |
static IntType mult_schrage(IntType a, IntType value)
|
williamr@2
|
313 |
{
|
williamr@2
|
314 |
const IntType q = m / a;
|
williamr@2
|
315 |
const IntType r = m % a;
|
williamr@2
|
316 |
|
williamr@2
|
317 |
assert(r < q); // check that overflow cannot happen
|
williamr@2
|
318 |
|
williamr@2
|
319 |
value = a*(value%q) - r*(value/q);
|
williamr@2
|
320 |
while(value <= 0)
|
williamr@2
|
321 |
value += m;
|
williamr@2
|
322 |
return value;
|
williamr@2
|
323 |
}
|
williamr@2
|
324 |
|
williamr@2
|
325 |
// invert c in the finite field (mod m) (m must be prime)
|
williamr@2
|
326 |
static IntType invert_euclidian(IntType c)
|
williamr@2
|
327 |
{
|
williamr@2
|
328 |
// we are interested in the gcd factor for c, because this is our inverse
|
williamr@2
|
329 |
BOOST_STATIC_ASSERT(m > 0);
|
williamr@2
|
330 |
#ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS
|
williamr@2
|
331 |
BOOST_STATIC_ASSERT(boost::integer_traits<IntType>::is_signed);
|
williamr@2
|
332 |
#endif
|
williamr@2
|
333 |
assert(c > 0);
|
williamr@2
|
334 |
IntType l1 = 0;
|
williamr@2
|
335 |
IntType l2 = 1;
|
williamr@2
|
336 |
IntType n = c;
|
williamr@2
|
337 |
IntType p = m;
|
williamr@2
|
338 |
for(;;) {
|
williamr@2
|
339 |
IntType q = p / n;
|
williamr@2
|
340 |
l1 -= q * l2; // this requires a signed IntType!
|
williamr@2
|
341 |
p -= q * n;
|
williamr@2
|
342 |
if(p == 0)
|
williamr@2
|
343 |
return (l2 < 1 ? l2 + m : l2);
|
williamr@2
|
344 |
IntType q2 = n / p;
|
williamr@2
|
345 |
l2 -= q2 * l1;
|
williamr@2
|
346 |
n -= q2 * p;
|
williamr@2
|
347 |
if(n == 0)
|
williamr@2
|
348 |
return (l1 < 1 ? l1 + m : l1);
|
williamr@2
|
349 |
}
|
williamr@2
|
350 |
}
|
williamr@2
|
351 |
};
|
williamr@2
|
352 |
|
williamr@2
|
353 |
|
williamr@2
|
354 |
#endif
|
williamr@2
|
355 |
|
williamr@2
|
356 |
} // namespace random
|
williamr@2
|
357 |
} // namespace boost
|
williamr@2
|
358 |
|
williamr@2
|
359 |
#endif // BOOST_RANDOM_CONST_MOD_HPP
|