blitz Version 0.10
random/mt.h
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 /*
00004  * $Id: mt.h,v 1.9 2011/03/25 22:41:17 julianc Exp $
00005  *
00006  * A C-program for MT19937: Integer version (1998/4/6)            
00007  *  genrand() generates one pseudorandom unsigned integer (32bit) 
00008  * which is uniformly distributed among 0 to 2^32-1  for each    
00009  * call. sgenrand(seed) set initial values to the working area   
00010  * of 624 words. Before genrand(), sgenrand(seed) must be        
00011  * called once. (seed is any 32-bit integer except for 0).
00012  *   Coded by Takuji Nishimura, considering the suggestions by    
00013  * Topher Cooper and Marc Rieffel in July-Aug. 1997.             
00014  *
00015  * This library is free software; you can redistribute it and/or 
00016  * modify it under the terms of the GNU Library General Public   
00017  * License as published by the Free Software Foundation; either    
00018  * version 2 of the License, or (at your option) any later         
00019  * version.                                                        
00020  * This library is distributed in the hope that it will be useful, 
00021  * but WITHOUT ANY WARRANTY; without even the implied warranty of  
00022  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.            
00023  * See the GNU Library General Public License for more details.    
00024  * You should have received a copy of the GNU Library General      
00025  * Public License along with this library; if not, write to the    
00026  * Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA    
00027  * 02111-1307  USA                                                 
00028  *
00029  * Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.       
00030  * When you use this, send an email to: matumoto@math.keio.ac.jp   
00031  * with an appropriate reference to your work.                     
00032  *
00033  * REFERENCE                                                       
00034  * M. Matsumoto and T. Nishimura,                                  
00035  * "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform
00036  * Pseudo-Random Number Generator",                                
00037  * ACM Transactions on Modeling and Computer Simulation,           
00038  * Vol. 8, No. 1, January 1998, pp 3--30.                          
00039  *
00040  * See 
00041  *     http://www.math.keio.ac.jp/~matumoto/emt.html
00042  * and
00043  *     http://www.acm.org/pubs/citations/journals/tomacs/1998-8-1/p3-matsumoto/
00044  *
00045  */
00046 
00047 #ifndef BZ_RAND_MT
00048 #define BZ_RAND_MT
00049 
00050 #include <blitz/blitz.h>
00051 
00052 #include <vector>
00053 #include <string>
00054 #include <sstream> 
00055 #include <iostream>
00056 #include <iterator>
00057 
00058 #ifndef UINT_MAX
00059   #include <limits.h>
00060 #endif
00061 
00062 BZ_NAMESPACE(ranlib)
00063 
00064 #if UINT_MAX < 4294967295U
00065   typedef unsigned long twist_int;  // must be at least 32 bits
00066 #else
00067   typedef unsigned int twist_int;
00068 #endif
00069 
00070 class MersenneTwister
00071 {
00072 public:
00073 
00074 private:
00075 
00076 #if defined(BZ_HAVE_NAMESPACES) && defined(BZ_HAVE_STD)
00077   typedef std::vector<twist_int> State;
00078 #else
00079   typedef vector<twist_int> State;
00080 #endif
00081   typedef State::size_type SizeType;
00082   typedef State::iterator Iter;
00083 
00084   // Implements Step 2 and half of Step 3 in the MN98 description
00085   struct BitMixer {
00086     BitMixer() : s0(0), K(0x9908b0df) {};
00087     BitMixer(twist_int k) : s0(0), K(k) {};
00088 
00089     // Return 0 if lsb of s1=0, a if lsb of s1=1.
00090     inline twist_int low_mask (twist_int s1) const {
00091       // This does not actually result in a branch because it's
00092       // equivalent to ( -(s1 & 1u) ) & K
00093       return (s1&1u) ? K : 0u;
00094     }
00095     
00096     // Return y>>1 in MN98 (step 2 + part of 3).
00097     // y = (x[i] AND u) OR (x[i+1 mod n] AND ll), where s0=x[i] and
00098     // s1=x[i+1 mod n]
00099     inline twist_int high_mask (twist_int s1) const {
00100       return ( (s0&0x80000000) | (s1&0x7fffffff) ) >>1;
00101     }
00102 
00103     // Calculate (half of) step 3 in MN98.
00104     inline twist_int operator() (twist_int s1) {
00105       // (y>>1) XOR (0 if lsb of s1=0, a if lsb of s1=1)
00106       // x[i+m] is XORed in reload
00107       // (Note that it is OK to call low_mask with s1 (x[i+1]) and not
00108       // with y, like MN98 does, because s1&1 == y&1 by construction.
00109       twist_int r = high_mask(s1) ^ low_mask(s1);
00110       s0 = s1;
00111       return r;
00112     }
00113     twist_int s0; // this is "x[i]" in the MN98 description
00114     const twist_int K; // MN98 "a" vector
00115   };
00116 
00117 enum { N = 624, 
00118        PF = 397, // MN98 "m"
00119        reference_seed = 4357 }; 
00120  
00121   void initialize()
00122   {
00123     S.resize(N);
00124     I = S.end();
00125   }
00126  
00127 public: 
00128   MersenneTwister() : b_(0x9D2C5680), c_(0xEFC60000)
00129   {
00130     initialize();
00131     seed();
00132 
00133     // There is a problem: static initialization + templates do not
00134     // mix very well in C++.  If you have a static member
00135     // of a class template, there is no guarantee on its order iin
00136     // static initialization.  This MersenneTwister class is used
00137     // elsewhere as a static member of a template class, and it is
00138     // possible (in fact, I've done so) to create a static initializer
00139     // that will invoke the seed() method of this object before its
00140     // ctor has been called (result: crash). 
00141     // ANSI C++ is stranger than fiction.
00142     // Currently the documentation forbids using RNGs from
00143     // static initializers.  There doesn't seem to be a good
00144     // fix.
00145   }
00146 
00147   MersenneTwister(twist_int aa, twist_int bb, twist_int cc) : 
00148     twist_(aa), b_(bb), c_(cc)
00149   {
00150     initialize();
00151     seed();
00152   }
00153 
00154   MersenneTwister(twist_int initial_seed) : b_(0x9D2C5680), c_(0xEFC60000)
00155   {
00156     initialize();
00157     seed(initial_seed);
00158   }
00159 
00160   MersenneTwister(twist_int aa, twist_int bb, twist_int cc, 
00161       twist_int initial_seed) : twist_(aa), b_(bb), c_(cc)
00162   {
00163     initialize();
00164     seed(initial_seed);
00165   }
00166 
00167   // Seed. Uses updated seed algorithm from mt19937ar. The old
00168   // algorithm would yield sequences that were close from seeds that
00169   // were close.
00170   void seed (twist_int seed = reference_seed)
00171   {
00172     // seed cannot equal 0
00173     if (seed == 0)
00174       seed = reference_seed;
00175 
00176     S[0] = seed & 0xFFFFFFFF;
00177     for (SizeType mti=1; mti<S.size(); ++mti) {
00178       S[mti] = (1812433253U * (S[mti-1] ^ (S[mti-1] >> 30)) + mti); 
00179       S[mti] &= 0xffffffffU;
00180     }
00181 
00182     reload();
00183   }
00184 
00185   // Seed by array, swiped directly from mt19937ar. Gives a larger
00186   // initial seed space.
00187   void seed (State seed_vector)
00188   {
00189     SizeType i, j, k;
00190     seed(19650218U);
00191     i=1; j=0;
00192     const SizeType N=S.size();
00193   const SizeType n=seed_vector.size();
00194     k = (N>n ? N : n);
00195     for (; k; k--) {
00196       S[i] = (S[i] ^ ((S[i-1] ^ (S[i-1] >> 30)) * 1664525U))
00197   + seed_vector[j] + j; /* non linear */
00198       S[i] &= 0xffffffffU; /* for WORDSIZE > 32 machines */
00199       i++; j++;
00200       if (i>=N) { S[0] = S[N-1]; i=1; }
00201       if (j>=n) j=0;
00202     }
00203     for (k=N-1; k; k--) {
00204       S[i] = (S[i] ^ ((S[i-1] ^ (S[i-1] >> 30)) * 1566083941UL))
00205   - i; /* non linear */
00206       S[i] &= 0xffffffffU; /* for WORDSIZE > 32 machines */
00207       i++;
00208       if (i>=N) { S[0] = S[N-1]; i=1; }
00209     }
00210 
00211     S[0] = 0x80000000U; /* MSB is 1; assuring non-zero initial array */ 
00212 
00213     reload();
00214   }
00215 
00216   // generate a full new x array
00217   void reload (void)
00218   {
00219     // This check is required because it is possible to call random()
00220     // before the constructor.  See the note above about static
00221     // initialization.
00222 
00223     Iter p0 = S.begin();
00224     Iter pM = p0 + PF;
00225     twist_ (S[0]); // set x[i]=x[0] in the twister (prime the pump)
00226     for (Iter pf_end = S.begin()+(N-PF); p0 != pf_end; ++p0, ++pM)
00227       // mt[kk] = mt[kk+m] XOR ((y>>1)XOR(mag01), as calc by BitMixer)
00228       *p0 = *pM ^ twist_ (p0[1]);
00229 
00230     // This is the "modulo part" where kk+m rolls over
00231     pM = S.begin();
00232     for (Iter s_end = S.begin()+(N-1); p0 != s_end; ++p0, ++pM)
00233       *p0 = *pM ^ twist_ (p0[1]);
00234     // and final element where kk+1 rolls over
00235     *p0 = *pM ^ twist_ (S[0]);
00236 
00237     I = S.begin();
00238   }
00239 
00240   inline twist_int random (void)
00241   {
00242     if (I >= S.end()) reload();
00243     // get next word from array
00244     twist_int y = *I++;
00245     // Step 4+5 in MN98, multiply by tempering matrix
00246     y ^= (y >> 11);
00247     y ^= (y <<  7) & b_; 
00248     y ^= (y << 15) & c_; 
00249     y ^= (y >> 18);
00250     return y;
00251   }
00252 
00253   // functions for getting/setting state
00254   class mt_state {
00255     friend class MersenneTwister;
00256   private:
00257     State S;
00258   State::difference_type I;
00259   public: 
00260     mt_state() { }
00261   mt_state(State s, State::difference_type i) : S(s), I(i) { }
00262     mt_state(const std::string& s) {
00263       std::istringstream is(s);
00264       is >> I;
00265       S = State(std::istream_iterator<twist_int>(is),
00266     std::istream_iterator<twist_int>());
00267       assert(!S.empty());
00268     }
00269     operator bool() const { return !S.empty(); }
00270     std::string str() const {
00271       if (S.empty())
00272   return std::string();
00273       std::ostringstream os;
00274       os << I << " ";
00275       std::copy(S.begin(), S.end(),
00276     std::ostream_iterator<twist_int>(os," "));
00277       return os.str();
00278     }
00279   };
00280   
00281   typedef mt_state T_state;
00282   T_state getState() const { return T_state(S, I-S.begin()); }
00283   std::string getStateString() const {
00284     T_state tmp(S, I-S.begin());
00285     return tmp.str();
00286   }
00287   void setState(const T_state& s) {
00288     if (!s) {
00289       std::cerr << "Error: state is empty" << std::endl;
00290       return;
00291     } 
00292     S = s.S;
00293     I = S.begin() + s.I;
00294   }
00295   void setState(const std::string& s) {
00296     T_state tmp(s);
00297     setState(tmp);
00298   }
00299   
00300 private:
00301   BitMixer twist_;
00302   const twist_int b_,c_;
00303 
00304   State   S;
00305   Iter    I;
00306 };
00307 
00308 
00311 class MersenneTwisterCreator
00312 {
00313 public:
00314   static MersenneTwister create(unsigned int i) {
00315     // We only have n different parameter sets
00316     i = i % n;
00317     return MersenneTwister(a_[i], b_[i], c_[i]);
00318   };
00319 
00320 private:
00321   static const unsigned int n=48;
00322   static const twist_int a_[n];
00323   static const twist_int b_[n];
00324   static const twist_int c_[n];
00325 };
00326 
00327 BZ_NAMESPACE_END
00328 
00329 #endif // BZ_RAND_MT
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Friends Defines