source: NonGTP/Boost/boost/random/detail/const_mod.hpp @ 857

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