1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
|
/* -*- Mode: C++ -*- */
/* This is public-domain Mersenne Twister code,
* attributed to Michael Brundage. Thanks!
* http://www.qbrundage.com/michaelb/pubs/essays/random_number_generation.html
*/
#undef MT_LEN
#undef MT_IA
class MTRandom {
public:
enum Constants {
MT_LEN = 624,
MT_IA = 397
};
static const uint32_t TEST_SEED1;
static const uint32_t UPPER_MASK;
static const uint32_t LOWER_MASK;
static const uint32_t MATRIX_A;
MTRandom() {
Init(TEST_SEED1);
}
MTRandom(uint32_t seed) {
Init(seed);
}
uint32_t Rand32 () {
uint32_t y;
static unsigned long mag01[2] = {
0 , MATRIX_A
};
if (mt_index_ >= MT_LEN) {
int kk;
for (kk = 0; kk < MT_LEN - MT_IA; kk++) {
y = (mt_buffer_[kk] & UPPER_MASK) | (mt_buffer_[kk + 1] & LOWER_MASK);
mt_buffer_[kk] = mt_buffer_[kk + MT_IA] ^ (y >> 1) ^ mag01[y & 0x1UL];
}
for (;kk < MT_LEN - 1; kk++) {
y = (mt_buffer_[kk] & UPPER_MASK) | (mt_buffer_[kk + 1] & LOWER_MASK);
mt_buffer_[kk] = mt_buffer_[kk + (MT_IA - MT_LEN)] ^ (y >> 1) ^ mag01[y & 0x1UL];
}
y = (mt_buffer_[MT_LEN - 1] & UPPER_MASK) | (mt_buffer_[0] & LOWER_MASK);
mt_buffer_[MT_LEN - 1] = mt_buffer_[MT_IA - 1] ^ (y >> 1) ^ mag01[y & 0x1UL];
mt_index_ = 0;
}
y = mt_buffer_[mt_index_++];
y ^= (y >> 11);
y ^= (y << 7) & 0x9d2c5680UL;
y ^= (y << 15) & 0xefc60000UL;
y ^= (y >> 18);
return y;
}
uint32_t ExpRand32(uint32_t mean) {
double mean_d = mean;
double erand = log (1.0 / (Rand32() / (double)UINT32_MAX));
uint32_t x = (uint32_t) (mean_d * erand + 0.5);
return x;
}
uint64_t Rand64() {
return ((uint64_t)Rand32() << 32) | Rand32();
}
uint64_t ExpRand64(uint64_t mean) {
double mean_d = mean;
double erand = log (1.0 / (Rand64() / (double)UINT32_MAX));
uint64_t x = (uint64_t) (mean_d * erand + 0.5);
return x;
}
template <typename T>
T Rand() {
switch (sizeof(T)) {
case sizeof(uint32_t):
return Rand32();
case sizeof(uint64_t):
return Rand64();
default:
cerr << "Invalid sizeof T" << endl;
abort();
}
}
template <typename T>
T ExpRand(T mean) {
switch (sizeof(T)) {
case sizeof(uint32_t):
return ExpRand32(mean);
case sizeof(uint64_t):
return ExpRand64(mean);
default:
cerr << "Invalid sizeof T" << endl;
abort();
}
}
private:
void Init(uint32_t seed) {
mt_buffer_[0] = seed;
mt_index_ = MT_LEN;
for (int i = 1; i < MT_LEN; i++) {
/* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
/* In the previous versions, MSBs of the seed affect */
/* only MSBs of the array mt[]. */
/* 2002/01/09 modified by Makoto Matsumoto */
mt_buffer_[i] =
(1812433253UL * (mt_buffer_[i-1] ^ (mt_buffer_[i-1] >> 30)) + i);
}
}
int mt_index_;
uint32_t mt_buffer_[MT_LEN];
};
const uint32_t MTRandom::TEST_SEED1 = 5489UL;
const uint32_t MTRandom::UPPER_MASK = 0x80000000;
const uint32_t MTRandom::LOWER_MASK = 0x7FFFFFFF;
const uint32_t MTRandom::MATRIX_A = 0x9908B0DF;
class MTRandom8 {
public:
MTRandom8(MTRandom *rand)
: rand_(rand) {
}
uint8_t Rand8() {
uint32_t r = rand_->Rand32();
// TODO: make this use a single byte at a time?
return (r & 0xff) ^ (r >> 7) ^ (r >> 15) ^ (r >> 21);
}
private:
MTRandom *rand_;
};
|
