blitz Version 0.10
|
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