00001 /* 00002 * lib/numerics/rng/mt19937.h 00003 * 00004 * Implementation of fast Mersenne Twister PRNG with period 2^19937-1. 00005 * 00006 * Converted to class 2004 by Wolfgang Wieser ] wwieser (a) gmx <*> de [ based on 00007 * mt19937ar-cok.c by Takuji Nishimura and Makoto Matsumoto; see mt19937.cc 00008 * for more legal information. 00009 * 00010 * This file may be distributed and/or modified under the terms of the 00011 * GNU General Public License version 2 as published by the Free Software 00012 * Foundation. (See COPYING.GPL for details.) 00013 * 00014 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE 00015 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. 00016 * 00017 */ 00018 00019 #ifndef _LIB_NUMERICS_RNG_MT19937_H_ 00020 #define _LIB_NUMERICS_RNG_MT19937_H_ 1 00021 00035 #include <lib/numerics/rng/rng.h> 00036 00037 00038 namespace NUM // numerics 00039 { 00040 00073 class RandomNumberGenerator_MersenneTwister19937 : 00074 public RandomNumberGenerator_Base 00075 { 00076 public: 00077 static const int N=624; 00078 static const int M=397; 00079 static const uint32 MATRIX_A=0x9908b0dfU; /* constant vector a */ 00080 static const uint32 UMASK=0x80000000U; /* most significant w-r bits */ 00081 static const uint32 LMASK=0x7fffffffU; /* least significant r bits */ 00082 00083 private: 00084 uint32 state[N]; /* the array for the state vector */ 00085 uint32 *next; 00086 int left; 00087 00088 inline uint32 MIXBITS(uint32 u,uint32 v) 00089 { return( (u & UMASK) | (v & LMASK) ); } 00090 inline uint32 TWIST(uint32 u,uint32 v) 00091 { return((MIXBITS(u,v) >> 1) ^ (v&1U ? MATRIX_A : 0U)); } 00092 00093 void next_state(); 00094 00096 inline uint32 _rand32() 00097 { 00098 if(!--left) next_state(); 00099 uint32 y=*next++; 00100 /* Tempering */ 00101 y ^= (y >> 11); 00102 y ^= (y << 7) & 0x9d2c5680U; 00103 y ^= (y << 15) & 0xefc60000U; 00104 y ^= (y >> 18); 00105 return(y); 00106 } 00107 00109 RandomNumberGenerator_MersenneTwister19937( 00110 const RandomNumberGenerator_MersenneTwister19937 &); 00112 void operator=(const RandomNumberGenerator_MersenneTwister19937 &); 00113 public: 00115 RandomNumberGenerator_MersenneTwister19937(uint32 s) : 00116 RandomNumberGenerator_Base() 00117 { seed(s); } 00119 RandomNumberGenerator_MersenneTwister19937(const uint32 *init_key, 00120 int nelem) : RandomNumberGenerator_Base() 00121 { seed(init_key,nelem); } 00123 ~RandomNumberGenerator_MersenneTwister19937() {} 00124 00131 void seed(uint32 seed); 00132 00143 void seed(const uint32 *init_key,int nelem); 00144 00146 uint32 rand_uint32() // not inline (overriding virtual) 00147 { return(_rand32()); } 00148 00150 int32 rand_int32() // not inline (overriding virtual) 00151 { return((int32)(_rand32()>>1)); } 00152 00154 dbl rand_dbl32o() // not inline (overriding virtual) 00155 { return(_rand32()*FPConst<dbl>::div_1_pow_2_32); } 00156 00158 dbl rand_dbl32() // not inline (overriding virtual) 00159 { return(_rand32()*FPConst<dbl>::div_1_pow_2_32m1); } 00160 00162 dbl rand_dbl53o() // not inline (overriding virtual) 00163 { 00164 uint32 a=_rand32()>>5, b=_rand32()>>6; 00165 return( (a*67108864.0+b)*FPConst<dbl>::div_1_pow_2_53 ); 00166 } 00167 00169 dbl rand_dbl53() // not inline (overriding virtual) 00170 { 00171 uint32 a=_rand32()>>5, b=_rand32()>>6; 00172 return( (a*67108864.0+b)*FPConst<dbl>::div_1_pow_2_53m1 ); 00173 } 00174 }; 00175 00176 } // end of namespace NUM 00177 00178 #endif /* _LIB_NUMERICS_RNG_MT19937_H_ */