1 /* boost random/mersenne_twister.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: mersenne_twister.hpp,v 1.20 2005/07/21 22:04:31 jmaurer Exp $

    11  *

    12  * Revision history

    13  *  2001-02-18  moved to individual header files

    14  */

    15 

    16 #ifndef BOOST_RANDOM_MERSENNE_TWISTER_HPP

    17 #define BOOST_RANDOM_MERSENNE_TWISTER_HPP

    18 

    19 #include <iostream>

    20 #include <algorithm>     // std::copy

    21 #include <stdexcept>

    22 #include <boost/config.hpp>

    23 #include <boost/limits.hpp>

    24 #include <boost/static_assert.hpp>

    25 #include <boost/integer_traits.hpp>

    26 #include <boost/cstdint.hpp>

    27 #include <boost/random/linear_congruential.hpp>

    28 #include <boost/detail/workaround.hpp>

    29 #include <boost/random/detail/ptr_helper.hpp>

    30 

    31 namespace boost {

    32 namespace random {

    33 

    34 // http://www.math.keio.ac.jp/matumoto/emt.html

    35 template<class UIntType, int w, int n, int m, int r, UIntType a, int u,

    36   int s, UIntType b, int t, UIntType c, int l, UIntType val>

    37 class mersenne_twister

    38 {

    39 public:

    40   typedef UIntType result_type;

    41   BOOST_STATIC_CONSTANT(int, word_size = w);

    42   BOOST_STATIC_CONSTANT(int, state_size = n);

    43   BOOST_STATIC_CONSTANT(int, shift_size = m);

    44   BOOST_STATIC_CONSTANT(int, mask_bits = r);

    45   BOOST_STATIC_CONSTANT(UIntType, parameter_a = a);

    46   BOOST_STATIC_CONSTANT(int, output_u = u);

    47   BOOST_STATIC_CONSTANT(int, output_s = s);

    48   BOOST_STATIC_CONSTANT(UIntType, output_b = b);

    49   BOOST_STATIC_CONSTANT(int, output_t = t);

    50   BOOST_STATIC_CONSTANT(UIntType, output_c = c);

    51   BOOST_STATIC_CONSTANT(int, output_l = l);

    52 

    53   BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);

    54   

    55   mersenne_twister() { seed(); }

    56 

    57 #if defined(__SUNPRO_CC) && (__SUNPRO_CC <= 0x520)

    58   // Work around overload resolution problem (Gennadiy E. Rozental)

    59   explicit mersenne_twister(const UIntType& value)

    60 #else

    61   explicit mersenne_twister(UIntType value)

    62 #endif

    63   { seed(value); }

    64   template<class It> mersenne_twister(It& first, It last) { seed(first,last); }

    65 

    66   template<class Generator>

    67   explicit mersenne_twister(Generator & gen) { seed(gen); }

    68 

    69   // compiler-generated copy ctor and assignment operator are fine

    70 

    71   void seed() { seed(UIntType(5489)); }

    72 

    73 #if defined(__SUNPRO_CC) && (__SUNPRO_CC <= 0x520)

    74   // Work around overload resolution problem (Gennadiy E. Rozental)

    75   void seed(const UIntType& value)

    76 #else

    77   void seed(UIntType value)

    78 #endif

    79   {

    80     // New seeding algorithm from 

    81     // http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/MT2002/emt19937ar.html

    82     // In the previous versions, MSBs of the seed affected only MSBs of the

    83     // state x[].

    84     const UIntType mask = ~0u;

    85     x[0] = value & mask;

    86     for (i = 1; i < n; i++) {

    87       // See Knuth "The Art of Computer Programming" Vol. 2, 3rd ed., page 106

    88       x[i] = (1812433253UL * (x[i-1] ^ (x[i-1] >> (w-2))) + i) & mask;

    89     }

    90   }

    91 

    92   // For GCC, moving this function out-of-line prevents inlining, which may

    93   // reduce overall object code size.  However, MSVC does not grok

    94   // out-of-line definitions of member function templates.

    95   template<class Generator>

    96   void seed(Generator & gen)

    97   {

    98 #ifndef BOOST_NO_LIMITS_COMPILE_TIME_CONSTANTS

    99     BOOST_STATIC_ASSERT(!std::numeric_limits<result_type>::is_signed);

   100 #endif

   101     // I could have used std::generate_n, but it takes "gen" by value

   102     for(int j = 0; j < n; j++)

   103       x[j] = gen();

   104     i = n;

   105   }

   106 

   107   template<class It>

   108   void seed(It& first, It last)

   109   {

   110     int j;

   111     for(j = 0; j < n && first != last; ++j, ++first)

   112       x[j] = *first;

   113     i = n;

   114     if(first == last && j < n)

   115       throw std::invalid_argument("mersenne_twister::seed");

   116   }

   117   

   118   result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 0; }

   119   result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const

   120   {

   121     // avoid "left shift count >= with of type" warning

   122     result_type res = 0;

   123     for(int i = 0; i < w; ++i)

   124       res |= (1u << i);

   125     return res;

   126   }

   127 

   128   result_type operator()();

   129   static bool validation(result_type v) { return val == v; }

   130 

   131 #ifndef BOOST_NO_OPERATORS_IN_NAMESPACE

   132 

   133 #ifndef BOOST_NO_MEMBER_TEMPLATE_FRIENDS

   134   template<class CharT, class Traits>

   135   friend std::basic_ostream<CharT,Traits>&

   136   operator<<(std::basic_ostream<CharT,Traits>& os, const mersenne_twister& mt)

   137   {

   138     for(int j = 0; j < mt.state_size; ++j)

   139       os << mt.compute(j) << " ";

   140     return os;

   141   }

   142 

   143   template<class CharT, class Traits>

   144   friend std::basic_istream<CharT,Traits>&

   145   operator>>(std::basic_istream<CharT,Traits>& is, mersenne_twister& mt)

   146   {

   147     for(int j = 0; j < mt.state_size; ++j)

   148       is >> mt.x[j] >> std::ws;

   149     // MSVC (up to 7.1) and Borland (up to 5.64) don't handle the template

   150     // value parameter "n" available from the class template scope, so use

   151     // the static constant with the same value

   152     mt.i = mt.state_size;

   153     return is;

   154   }

   155 #endif

   156 

   157   friend bool operator==(const mersenne_twister& x, const mersenne_twister& y)

   158   {

   159     for(int j = 0; j < state_size; ++j)

   160       if(x.compute(j) != y.compute(j))

   161         return false;

   162     return true;

   163   }

   164 

   165   friend bool operator!=(const mersenne_twister& x, const mersenne_twister& y)

   166   { return !(x == y); }

   167 #else

   168   // Use a member function; Streamable concept not supported.

   169   bool operator==(const mersenne_twister& rhs) const

   170   {

   171     for(int j = 0; j < state_size; ++j)

   172       if(compute(j) != rhs.compute(j))

   173         return false;

   174     return true;

   175   }

   176 

   177   bool operator!=(const mersenne_twister& rhs) const

   178   { return !(*this == rhs); }

   179 #endif

   180 

   181 private:

   182   // returns x(i-n+index), where index is in 0..n-1

   183   UIntType compute(unsigned int index) const

   184   {

   185     // equivalent to (i-n+index) % 2n, but doesn't produce negative numbers

   186     return x[ (i + n + index) % (2*n) ];

   187   }

   188   void twist(int block);

   189 

   190   // state representation: next output is o(x(i))

   191   //   x[0]  ... x[k] x[k+1] ... x[n-1]     x[n]     ... x[2*n-1]   represents

   192   //  x(i-k) ... x(i) x(i+1) ... x(i-k+n-1) x(i-k-n) ... x[i(i-k-1)]

   193   // The goal is to always have x(i-n) ... x(i-1) available for

   194   // operator== and save/restore.

   195 

   196   UIntType x[2*n]; 

   197   int i;

   198 };

   199 

   200 #ifndef BOOST_NO_INCLASS_MEMBER_INITIALIZATION

   201 //  A definition is required even for integral static constants

   202 template<class UIntType, int w, int n, int m, int r, UIntType a, int u,

   203   int s, UIntType b, int t, UIntType c, int l, UIntType val>

   204 const bool mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::has_fixed_range;

   205 template<class UIntType, int w, int n, int m, int r, UIntType a, int u,

   206   int s, UIntType b, int t, UIntType c, int l, UIntType val>

   207 const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::state_size;

   208 template<class UIntType, int w, int n, int m, int r, UIntType a, int u,

   209   int s, UIntType b, int t, UIntType c, int l, UIntType val>

   210 const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::shift_size;

   211 template<class UIntType, int w, int n, int m, int r, UIntType a, int u,

   212   int s, UIntType b, int t, UIntType c, int l, UIntType val>

   213 const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::mask_bits;

   214 template<class UIntType, int w, int n, int m, int r, UIntType a, int u,

   215   int s, UIntType b, int t, UIntType c, int l, UIntType val>

   216 const UIntType mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::parameter_a;

   217 template<class UIntType, int w, int n, int m, int r, UIntType a, int u,

   218   int s, UIntType b, int t, UIntType c, int l, UIntType val>

   219 const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_u;

   220 template<class UIntType, int w, int n, int m, int r, UIntType a, int u,

   221   int s, UIntType b, int t, UIntType c, int l, UIntType val>

   222 const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_s;

   223 template<class UIntType, int w, int n, int m, int r, UIntType a, int u,

   224   int s, UIntType b, int t, UIntType c, int l, UIntType val>

   225 const UIntType mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_b;

   226 template<class UIntType, int w, int n, int m, int r, UIntType a, int u,

   227   int s, UIntType b, int t, UIntType c, int l, UIntType val>

   228 const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_t;

   229 template<class UIntType, int w, int n, int m, int r, UIntType a, int u,

   230   int s, UIntType b, int t, UIntType c, int l, UIntType val>

   231 const UIntType mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_c;

   232 template<class UIntType, int w, int n, int m, int r, UIntType a, int u,

   233   int s, UIntType b, int t, UIntType c, int l, UIntType val>

   234 const int mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::output_l;

   235 #endif

   236 

   237 template<class UIntType, int w, int n, int m, int r, UIntType a, int u,

   238   int s, UIntType b, int t, UIntType c, int l, UIntType val>

   239 void mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::twist(int block)

   240 {

   241   const UIntType upper_mask = (~0u) << r;

   242   const UIntType lower_mask = ~upper_mask;

   243 

   244   if(block == 0) {

   245     for(int j = n; j < 2*n; j++) {

   246       UIntType y = (x[j-n] & upper_mask) | (x[j-(n-1)] & lower_mask);

   247       x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0);

   248     }

   249   } else if (block == 1) {

   250     // split loop to avoid costly modulo operations

   251     {  // extra scope for MSVC brokenness w.r.t. for scope

   252       for(int j = 0; j < n-m; j++) {

   253         UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask);

   254         x[j] = x[j+n+m] ^ (y >> 1) ^ (y&1 ? a : 0);

   255       }

   256     }

   257     

   258     for(int j = n-m; j < n-1; j++) {

   259       UIntType y = (x[j+n] & upper_mask) | (x[j+n+1] & lower_mask);

   260       x[j] = x[j-(n-m)] ^ (y >> 1) ^ (y&1 ? a : 0);

   261     }

   262     // last iteration

   263     UIntType y = (x[2*n-1] & upper_mask) | (x[0] & lower_mask);

   264     x[n-1] = x[m-1] ^ (y >> 1) ^ (y&1 ? a : 0);

   265     i = 0;

   266   }

   267 }

   268 

   269 template<class UIntType, int w, int n, int m, int r, UIntType a, int u,

   270   int s, UIntType b, int t, UIntType c, int l, UIntType val>

   271 inline typename mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::result_type

   272 mersenne_twister<UIntType,w,n,m,r,a,u,s,b,t,c,l,val>::operator()()

   273 {

   274   if(i == n)

   275     twist(0);

   276   else if(i >= 2*n)

   277     twist(1);

   278   // Step 4

   279   UIntType z = x[i];

   280   ++i;

   281   z ^= (z >> u);

   282   z ^= ((z << s) & b);

   283   z ^= ((z << t) & c);

   284   z ^= (z >> l);

   285   return z;

   286 }

   287 

   288 } // namespace random

   289 

   290 

   291 typedef random::mersenne_twister<uint32_t,32,351,175,19,0xccab8ee7,11,

   292   7,0x31b6ab00,15,0xffe50000,17, 0xa37d3c92> mt11213b;

   293 

   294 // validation by experiment from mt19937.c

   295 typedef random::mersenne_twister<uint32_t,32,624,397,31,0x9908b0df,11,

   296   7,0x9d2c5680,15,0xefc60000,18, 3346425566U> mt19937;

   297 

   298 } // namespace boost

   299 

   300 BOOST_RANDOM_PTR_HELPER_SPEC(boost::mt19937)

   301 

   302 #endif // BOOST_RANDOM_MERSENNE_TWISTER_HPP