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