00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
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 #if UINT_MAX < 4294967295U
00064 typedef unsigned long twist_int;
00065 #else
00066 typedef unsigned int twist_int;
00067 #endif
00068
00069 class MersenneTwister
00070 {
00071 public:
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 typedef State::size_type SizeType;
00081 typedef State::iterator Iter;
00082
00083
00084 struct BitMixer {
00085 BitMixer() : s0(0), K(0x9908b0df) {};
00086 BitMixer(twist_int k) : s0(0), K(k) {};
00087
00088
00089 inline twist_int low_mask (twist_int s1) const {
00090
00091
00092 return (s1&1u) ? K : 0u;
00093 }
00094
00095
00096
00097
00098 inline twist_int high_mask (twist_int s1) const {
00099 return ( (s0&0x80000000) | (s1&0x7fffffff) ) >>1;
00100 }
00101
00102
00103 inline twist_int operator() (twist_int s1) {
00104
00105
00106
00107
00108 twist_int r = high_mask(s1) ^ low_mask(s1);
00109 s0 = s1;
00110 return r;
00111 }
00112 twist_int s0;
00113 const twist_int K;
00114 };
00115
00116 enum { N = 624,
00117 PF = 397,
00118 reference_seed = 4357 };
00119
00120 void initialize()
00121 {
00122 S.resize(N);
00123 I = S.end();
00124 }
00125
00126 public:
00127 MersenneTwister() : b_(0x9D2C5680), c_(0xEFC60000)
00128 {
00129 initialize();
00130 seed();
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144 }
00145
00146 MersenneTwister(twist_int aa, twist_int bb, twist_int cc) :
00147 twist_(aa), b_(bb), c_(cc)
00148 {
00149 initialize();
00150 seed();
00151 }
00152
00153 MersenneTwister(twist_int initial_seed) : b_(0x9D2C5680), c_(0xEFC60000)
00154 {
00155 initialize();
00156 seed(initial_seed);
00157 }
00158
00159 MersenneTwister(twist_int aa, twist_int bb, twist_int cc,
00160 twist_int initial_seed) : twist_(aa), b_(bb), c_(cc)
00161 {
00162 initialize();
00163 seed(initial_seed);
00164 }
00165
00166
00167
00168
00169 void seed (twist_int seed = reference_seed)
00170 {
00171
00172 if (seed == 0)
00173 seed = reference_seed;
00174
00175 S[0] = seed & 0xFFFFFFFF;
00176 for (SizeType mti=1; mti<S.size(); ++mti) {
00177 S[mti] = (1812433253U * (S[mti-1] ^ (S[mti-1] >> 30)) + mti);
00178 S[mti] &= 0xffffffffU;
00179 }
00180
00181 reload();
00182 }
00183
00184
00185
00186 void seed (State seed_vector)
00187 {
00188 SizeType i, j, k;
00189 seed(19650218U);
00190 i=1; j=0;
00191 const SizeType N=S.size();
00192 const SizeType n=seed_vector.size();
00193 k = (N>n ? N : n);
00194 for (; k; k--) {
00195 S[i] = (S[i] ^ ((S[i-1] ^ (S[i-1] >> 30)) * 1664525U))
00196 + seed_vector[j] + j;
00197 S[i] &= 0xffffffffU;
00198 i++; j++;
00199 if (i>=N) { S[0] = S[N-1]; i=1; }
00200 if (j>=n) j=0;
00201 }
00202 for (k=N-1; k; k--) {
00203 S[i] = (S[i] ^ ((S[i-1] ^ (S[i-1] >> 30)) * 1566083941UL))
00204 - i;
00205 S[i] &= 0xffffffffU;
00206 i++;
00207 if (i>=N) { S[0] = S[N-1]; i=1; }
00208 }
00209
00210 S[0] = 0x80000000U;
00211
00212 reload();
00213 }
00214
00215
00216 void reload (void)
00217 {
00218
00219
00220
00221
00222 Iter p0 = S.begin();
00223 Iter pM = p0 + PF;
00224 twist_ (S[0]);
00225 for (Iter pf_end = S.begin()+(N-PF); p0 != pf_end; ++p0, ++pM)
00226
00227 *p0 = *pM ^ twist_ (p0[1]);
00228
00229
00230 pM = S.begin();
00231 for (Iter s_end = S.begin()+(N-1); p0 != s_end; ++p0, ++pM)
00232 *p0 = *pM ^ twist_ (p0[1]);
00233
00234 *p0 = *pM ^ twist_ (S[0]);
00235
00236 I = S.begin();
00237 }
00238
00239 inline twist_int random (void)
00240 {
00241 if (I >= S.end()) reload();
00242
00243 twist_int y = *I++;
00244
00245 y ^= (y >> 11);
00246 y ^= (y << 7) & b_;
00247 y ^= (y << 15) & c_;
00248 y ^= (y >> 18);
00249 return y;
00250 }
00251
00252
00253 class mt_state {
00254 friend class MersenneTwister;
00255 private:
00256 State S;
00257 State::difference_type I;
00258 public:
00259 mt_state() { }
00260 mt_state(State s, State::difference_type i) : S(s), I(i) { }
00261 mt_state(const std::string& s) {
00262 std::istringstream is(s);
00263 is >> I;
00264 S = State(std::istream_iterator<twist_int>(is),
00265 std::istream_iterator<twist_int>());
00266 assert(!S.empty());
00267 }
00268 operator bool() const { return !S.empty(); }
00269 std::string str() const {
00270 if (S.empty())
00271 return std::string();
00272 std::ostringstream os;
00273 os << I << " ";
00274 std::copy(S.begin(), S.end(),
00275 std::ostream_iterator<twist_int>(os," "));
00276 return os.str();
00277 }
00278 };
00279
00280 typedef mt_state T_state;
00281 T_state getState() const { return T_state(S, I-S.begin()); }
00282 std::string getStateString() const {
00283 T_state tmp(S, I-S.begin());
00284 return tmp.str();
00285 }
00286 void setState(const T_state& s) {
00287 if (!s) {
00288 std::cerr << "Error: state is empty" << std::endl;
00289 return;
00290 }
00291 S = s.S;
00292 I = S.begin() + s.I;
00293 }
00294 void setState(const std::string& s) {
00295 T_state tmp(s);
00296 setState(tmp);
00297 }
00298
00299 private:
00300 BitMixer twist_;
00301 const twist_int b_,c_;
00302
00303 State S;
00304 Iter I;
00305 };
00306
00307
00310 class MersenneTwisterCreator
00311 {
00312 public:
00313 static MersenneTwister create(unsigned int i) {
00314
00315 i = i % n;
00316 return MersenneTwister(a_[i], b_[i], c_[i]);
00317 };
00318
00319 private:
00320 static const unsigned int n=48;
00321 static const twist_int a_[n];
00322 static const twist_int b_[n];
00323 static const twist_int c_[n];
00324 };
00325
00326 BZ_NAMESPACE_END
00327
00328 #endif // BZ_RAND_MT