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 class MersenneTwister
00064 {
00065 public:
00066
00067 #if UINT_MAX < 4294967295U
00068 typedef unsigned long twist_int;
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
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
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
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
00146
00147 *s++ = (x + (mask & j)) & 0xFFFFFFFF;
00148 x *= Knuth_A;
00149 }
00150
00151 reload();
00152 }
00153
00154 void reload (void)
00155 {
00156
00157
00158
00159
00160 Iter p0 = S.begin();
00161 Iter pM = p0 + PF;
00162 BitMixer twist;
00163 twist (S[0]);
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
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