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