random/mt.h

Go to the documentation of this file.
00001 
00002 /*
00003  * $Id: mt.h,v 1.6 2005/10/13 20:08:53 julianc Exp $
00004  *
00005  * A C-program for MT19937: Integer version (1998/4/6)            
00006  *  genrand() generates one pseudorandom unsigned integer (32bit) 
00007  * which is uniformly distributed among 0 to 2^32-1  for each    
00008  * call. sgenrand(seed) set initial values to the working area   
00009  * of 624 words. Before genrand(), sgenrand(seed) must be        
00010  * called once. (seed is any 32-bit integer except for 0).
00011  *   Coded by Takuji Nishimura, considering the suggestions by    
00012  * Topher Cooper and Marc Rieffel in July-Aug. 1997.             
00013  *
00014  * This library is free software; you can redistribute it and/or 
00015  * modify it under the terms of the GNU Library General Public   
00016  * License as published by the Free Software Foundation; either    
00017  * version 2 of the License, or (at your option) any later         
00018  * version.                                                        
00019  * This library is distributed in the hope that it will be useful, 
00020  * but WITHOUT ANY WARRANTY; without even the implied warranty of  
00021  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.            
00022  * See the GNU Library General Public License for more details.    
00023  * You should have received a copy of the GNU Library General      
00024  * Public License along with this library; if not, write to the    
00025  * Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA    
00026  * 02111-1307  USA                                                 
00027  *
00028  * Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura.       
00029  * When you use this, send an email to: matumoto@math.keio.ac.jp   
00030  * with an appropriate reference to your work.                     
00031  *
00032  * REFERENCE                                                       
00033  * M. Matsumoto and T. Nishimura,                                  
00034  * "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform
00035  * Pseudo-Random Number Generator",                                
00036  * ACM Transactions on Modeling and Computer Simulation,           
00037  * Vol. 8, No. 1, January 1998, pp 3--30.                          
00038  *
00039  * See 
00040  *     http://www.math.keio.ac.jp/~matumoto/emt.html
00041  * and
00042  *     http://www.acm.org/pubs/citations/journals/tomacs/1998-8-1/p3-matsumoto/
00043  *
00044  */
00045 
00046 #ifndef BZ_RAND_MT
00047 #define BZ_RAND_MT
00048 
00049 #include <blitz/blitz.h>
00050 
00051 #include <vector>
00052 #include <string>
00053 #include <sstream> 
00054 #include <iostream>
00055 #include <iterator>
00056 
00057 #ifndef UINT_MAX
00058   #include <limits.h>
00059 #endif
00060 
00061 BZ_NAMESPACE(ranlib)
00062 
00063 class MersenneTwister
00064 {
00065 public:
00066 
00067 #if UINT_MAX < 4294967295U
00068   typedef unsigned long twist_int;  // must be at least 32 bits
00069 #else
00070   typedef unsigned int twist_int;
00071 #endif
00072 
00073 private:
00074 
00075 #if defined(BZ_HAVE_NAMESPACES) && defined(BZ_HAVE_STD)
00076   typedef std::vector<twist_int> State;
00077 #else
00078   typedef vector<twist_int> State;
00079 #endif
00080 
00081   typedef State::iterator Iter;
00082 
00083   struct BitMixer {
00084     enum { K = 0x9908b0df };
00085     BitMixer() : s0(0) {}
00086     inline twist_int low_mask (twist_int s1) const {
00087       return (s1&1u) ? K : 0u;
00088     }
00089     inline twist_int high_mask (twist_int s1) const {
00090       return ((s0&0x80000000)|(s1&0x7fffffff))>>1;
00091     }
00092     inline twist_int operator() (twist_int s1) {
00093       twist_int r = high_mask(s1) ^ low_mask(s1);
00094       s0 = s1;
00095       return r;
00096     }
00097     twist_int s0;
00098   };
00099 
00100 enum { N = 624, PF = 397, reference_seed = 4357 }; 
00101  
00102   void initialize()
00103   {
00104     S.resize(N);
00105     I = S.end();
00106   }
00107  
00108 public: 
00109   MersenneTwister()
00110   {
00111     initialize();
00112     seed();
00113 
00114     // There is a problem: static initialization + templates do not
00115     // mix very well in C++.  If you have a static member
00116     // of a class template, there is no guarantee on its order iin
00117     // static initialization.  This MersenneTwister class is used
00118     // elsewhere as a static member of a template class, and it is
00119     // possible (in fact, I've done so) to create a static initializer
00120     // that will invoke the seed() method of this object before its
00121     // ctor has been called (result: crash). 
00122     // ANSI C++ is stranger than fiction.
00123     // Currently the documentation forbids using RNGs from
00124     // static initializers.  There doesn't seem to be a good
00125     // fix.
00126   }
00127 
00128   MersenneTwister(twist_int initial_seed)
00129   {
00130     initialize();
00131     seed(initial_seed);
00132   }
00133 
00134   void seed (twist_int seed = reference_seed)
00135   {
00136     // seed cannot equal 0
00137     if (seed == 0)
00138       seed = reference_seed;
00139 
00140     enum { Knuth_A = 69069 }; 
00141     twist_int x = seed & 0xFFFFFFFF;
00142     Iter s = S.begin();
00143     twist_int mask = (seed == reference_seed) ? 0 : 0xFFFFFFFF;
00144     for (int j = 0; j < N; ++j) {
00145       // adding j here avoids the risk of all zeros 
00146       // we suppress this term in "compatibility" mode  
00147       *s++ = (x + (mask & j)) & 0xFFFFFFFF; 
00148       x *= Knuth_A;
00149     }
00150 
00151     reload();
00152   }
00153 
00154   void reload (void)
00155   {
00156     // This check is required because it is possible to call random()
00157     // before the constructor.  See the note above about static
00158     // initialization.
00159 
00160     Iter p0 = S.begin();
00161     Iter pM = p0 + PF;
00162     BitMixer twist;
00163     twist (S[0]); // prime the pump
00164     for (Iter pf_end = S.begin()+(N-PF); p0 != pf_end; ++p0, ++pM)
00165       *p0 = *pM ^ twist (p0[1]);
00166     pM = S.begin();
00167     for (Iter s_end = S.begin()+(N-1); p0 != s_end; ++p0, ++pM)
00168       *p0 = *pM ^ twist (p0[1]);
00169     *p0 = *pM ^ twist (S[0]);
00170 
00171     I = S.begin();
00172   }
00173 
00174   inline twist_int random (void)
00175   {
00176     if (I >= S.end()) reload();
00177     twist_int y = *I++;
00178     y ^= (y >> 11);
00179     y ^= (y <<  7) & 0x9D2C5680;
00180     y ^= (y << 15) & 0xEFC60000;
00181     y ^= (y >> 18);
00182     return y;
00183   }
00184 
00185   // functions for getting/setting state
00186   class mt_state {
00187     friend class MersenneTwister;
00188   private:
00189     State S;
00190     int I;
00191   public: 
00192     mt_state() { }
00193     mt_state(State s, int i) : S(s), I(i) { }
00194     mt_state(const std::string& s) {
00195       std::istringstream is(s);
00196       is >> I;
00197       S = State(std::istream_iterator<twist_int>(is),
00198     std::istream_iterator<twist_int>());
00199       assert(!S.empty());
00200     }
00201     operator bool() const { return !S.empty(); }
00202     std::string str() const {
00203       if (S.empty())
00204   return std::string();
00205       std::ostringstream os;
00206       os << I << " ";
00207       std::copy(S.begin(), S.end(),
00208     std::ostream_iterator<twist_int>(os," "));
00209       return os.str();
00210     }
00211   };
00212   
00213   typedef mt_state T_state;
00214   T_state getState() const { return T_state(S, I-S.begin()); }
00215   std::string getStateString() const {
00216     T_state tmp(S, I-S.begin());
00217     return tmp.str();
00218   }
00219   void setState(const T_state& s) {
00220     if (!s) {
00221       std::cerr << "Error: state is empty" << std::endl;
00222       return;
00223     } 
00224     S = s.S;
00225     I = S.begin() + s.I;
00226   }
00227   void setState(const std::string& s) {
00228     T_state tmp(s);
00229     setState(tmp);
00230   }
00231   
00232 private:
00233   State   S;
00234   Iter    I;
00235 };
00236 
00237 
00238 BZ_NAMESPACE_END
00239 
00240 #endif // BZ_RAND_MT

Generated on Tue Dec 4 07:47:59 2007 for blitz by  doxygen 1.4.7