PageRenderTime 50ms CodeModel.GetById 13ms app.highlight 32ms RepoModel.GetById 1ms app.codeStats 0ms

/Src/Dependencies/Boost/boost/random/detail/const_mod.hpp

http://hadesmem.googlecode.com/
C++ Header | 216 lines | 174 code | 24 blank | 18 comment | 64 complexity | e04a4003ba3c42eeecca7a08bbae9181 MD5 | raw file
  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 71018 2011-04-05 21:27:52Z steven_watanabe $
 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 supress_warnings = (m == 0);
 40      BOOST_ASSERT(supress_warnings == 0);
 41      return x % (m + supress_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 supress_warnings = (m == 0);
 81      BOOST_ASSERT(supress_warnings == 0);
 82      return (a*x+c) % (m + supress_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 supress_warnings = (m == 0);
112    BOOST_ASSERT(supress_warnings == 0);
113    return a*x % (m + supress_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