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
27namespace boost {
28namespace random {
29
30template<class IntType, IntType m>
31class const_mod
32{
33public:
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
103private:
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