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