blitz Version 0.9
|
00001 /* A C-program for MT19937: Integer version (1998/4/6) */ 00002 /* genrand() generates one pseudorandom unsigned integer (32bit) */ 00003 /* which is uniformly distributed among 0 to 2^32-1 for each */ 00004 /* call. sgenrand(seed) set initial values to the working area */ 00005 /* of 624 words. Before genrand(), sgenrand(seed) must be */ 00006 /* called once. (seed is any 32-bit integer except for 0). */ 00007 /* Coded by Takuji Nishimura, considering the suggestions by */ 00008 /* Topher Cooper and Marc Rieffel in July-Aug. 1997. */ 00009 00010 /* This library is free software; you can redistribute it and/or */ 00011 /* modify it under the terms of the GNU Library General Public */ 00012 /* License as published by the Free Software Foundation; either */ 00013 /* version 2 of the License, or (at your option) any later */ 00014 /* version. */ 00015 /* This library is distributed in the hope that it will be useful, */ 00016 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */ 00017 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. */ 00018 /* See the GNU Library General Public License for more details. */ 00019 /* You should have received a copy of the GNU Library General */ 00020 /* Public License along with this library; if not, write to the */ 00021 /* Free Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA */ 00022 /* 02111-1307 USA */ 00023 00024 /* Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura. */ 00025 /* When you use this, send an email to: matumoto@math.keio.ac.jp */ 00026 /* with an appropriate reference to your work. */ 00027 00028 /* REFERENCE */ 00029 /* M. Matsumoto and T. Nishimura, */ 00030 /* "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform */ 00031 /* Pseudo-Random Number Generator", */ 00032 /* ACM Transactions on Modeling and Computer Simulation, */ 00033 /* Vol. 8, No. 1, January 1998, pp 3--30. */ 00034 00035 // See http://www.math.keio.ac.jp/~matumoto/emt.html 00036 00037 // 1999-01-25 adapted to STL-like idiom 00038 // allan@stokes.ca (Allan Stokes) www.stokes.ca 00039 00040 #ifndef BZ_RAND_MT 00041 #define BZ_RAND_MT 00042 00043 #ifndef BZ_BLITZ_H 00044 #include <blitz/blitz.h> 00045 #endif 00046 00047 #include <vector> 00048 00049 BZ_NAMESPACE(blitz) 00050 00051 // decomposition issues: 00052 // machine representation of integer types 00053 // output buffer option verses inline post-conditioning 00054 00055 class MersenneTwister 00056 { 00057 private: 00058 typedef unsigned int twist_int; // must be at least 32 bits 00059 // larger might be faster 00060 typedef vector<twist_int> State; 00061 typedef State::iterator Iter; 00062 00063 struct BitMixer { 00064 enum { K = 0x9908b0df }; 00065 BitMixer() : s0(0) {} 00066 inline friend twist_int low_mask (twist_int s1) { 00067 return (s1&1u) ? K : 0u; 00068 } 00069 inline twist_int high_mask (twist_int s1) const { 00070 return ((s0&0x80000000)|(s1&0x7fffffff))>>1; 00071 } 00072 inline twist_int operator() (twist_int s1) { 00073 twist_int r = high_mask(s1) ^ low_mask(s1); 00074 s0 = s1; 00075 return r; 00076 } 00077 twist_int s0; 00078 }; 00079 00080 enum { N = 624, PF = 397, reference_seed = 4357 }; 00081 00082 public: 00083 MersenneTwister () {} // S empty will trigger auto-seed 00084 00085 void seed (twist_int seed = reference_seed) 00086 { 00087 if (!S.size()) S.resize(N); 00088 enum { Knuth_A = 69069 }; 00089 twist_int x = seed & 0xFFFFFFFF; 00090 Iter s = S.begin(); 00091 twist_int mask = (seed == reference_seed) ? 0 : 0xFFFFFFFF; 00092 for (int j = 0; j < N; ++j) { 00093 // adding j here avoids the risk of all zeros 00094 // we suppress this term in "compatibility" mode 00095 *s++ = (x + (mask & j)) & 0xFFFFFFFF; 00096 x *= Knuth_A; 00097 } 00098 reload(); 00099 } 00100 00101 void reload (void) 00102 { 00103 if (!S.size()) seed (); // auto-seed detection 00104 00105 Iter p0 = S.begin(); 00106 Iter pM = p0 + PF; 00107 BitMixer twist; 00108 twist (S[0]); // prime the pump 00109 for (Iter pf_end = S.begin()+(N-PF); p0 != pf_end; ++p0, ++pM) 00110 *p0 = *pM ^ twist (p0[1]); 00111 pM = S.begin(); 00112 for (Iter s_end = S.begin()+(N-1); p0 != s_end; ++p0, ++pM) 00113 *p0 = *pM ^ twist (p0[1]); 00114 *p0 = *pM ^ twist (S[0]); 00115 00116 I = S.begin(); 00117 } 00118 00119 inline twist_int random (void) 00120 { 00121 if (I >= S.end()) reload(); 00122 twist_int y = *I++; 00123 y ^= (y >> 11); 00124 y ^= (y << 7) & 0x9D2C5680; 00125 y ^= (y << 15) & 0xEFC60000; 00126 y ^= (y >> 18); 00127 return y; 00128 } 00129 00130 private: 00131 State S; 00132 Iter I; 00133 }; 00134 00135 00136 // This version returns a double in the range [0,1). 00137 00138 class MersenneTwisterDouble { 00139 00140 public: 00141 MersenneTwisterDouble() 00142 { 00143 // f = 1/(2^32); 00144 f = (1.0 / 65536) / 65536; 00145 } 00146 00147 void randomize(unsigned int seed) 00148 { 00149 gen_.seed(seed); 00150 } 00151 00152 double random() 00153 { 00154 unsigned long y1 = gen_.random(); 00155 unsigned long y2 = gen_.random(); 00156 00157 return ((y1 * f) * y2 * f); 00158 } 00159 00160 private: 00161 MersenneTwister gen_; 00162 double f; 00163 }; 00164 00165 BZ_NAMESPACE_END 00166 00167 #endif // BZ_RAND_MT