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
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
|
//---------------------------------------------------------------------
// Algorithmic Conjurings @ http://www.coyotegulch.com
//
// mtprng.cc (libcoyote)
//
// Mersenne Twister -- A pseudorandom Number Generator
//
// ORIGINAL ALGORITHM COPYRIGHT
// ============================
// Copyright (C) 1997, 2002 Makoto Matsumoto and Takuji Nishimura.
// Any feedback is very welcome. For any question, comments, see
// http://www.math.keio.ac.jp/matumoto/emt.html or email
// matumoto@math.keio.ac.jp
//---------------------------------------------------------------------
//
// COPYRIGHT NOTICE, DISCLAIMER, and LICENSE:
//
// This notice applies *only* to this specific expression of this
// algorithm, and does not imply ownership or invention of the
// implemented algorithm.
//
// If you modify this file, you may insert additional notices
// immediately following this sentence.
//
// Copyright 2001, 2002, 2003 Scott Robert Ladd.
// All rights reserved, except as noted herein.
//
// This computer program source file is supplied "AS IS". Scott Robert
// Ladd (hereinafter referred to as "Author") disclaims all warranties,
// expressed or implied, including, without limitation, the warranties
// of merchantability and of fitness for any purpose. The Author
// assumes no liability for direct, indirect, incidental, special,
// exemplary, or consequential damages, which may result from the use
// of this software, even if advised of the possibility of such damage.
//
// The Author hereby grants anyone permission to use, copy, modify, and
// distribute this source code, or portions hereof, for any purpose,
// without fee, subject to the following restrictions:
//
// 1. The origin of this source code must not be misrepresented.
//
// 2. Altered versions must be plainly marked as such and must not
// be misrepresented as being the original source.
//
// 3. This Copyright notice may not be removed or altered from any
// source or altered source distribution.
//
// The Author specifically permits (without fee) and encourages the use
// of this source code for entertainment, education, or decoration. If
// you use this source code in a product, acknowledgment is not required
// but would be appreciated.
//
// Acknowledgement:
// This license is based on the wonderful simple license that
// accompanies libpng.
//
//-----------------------------------------------------------------------
//
// For more information on this software package, please visit
// Scott's web site, Coyote Gulch Productions, at:
//
// http://www.coyotegulch.com
//
//-----------------------------------------------------------------------
#include "mtprng.hh"
using namespace libcoyote;
#if _MSC_VER < 1300
const mtprng::int_type mtprng::MATRIX_A = 0x9908b0dfUL;
const mtprng::int_type mtprng::UPPER_MASK = 0x80000000UL;
const mtprng::int_type mtprng::LOWER_MASK = 0x7fffffffUL;
#endif
//--------------------------------------------------------------------------
// Constructor
mtprng::mtprng(int_type seed)
{
init(seed);
}
//--------------------------------------------------------------------------
// Initializes the generator with "seed"
void mtprng::init(int_type seed)
{
// Save seed for historical purpose
m_seed = seed;
m_mt[0] = seed & 0xffffffffUL;
// Set the seed using values suggested by Matsumoto & Nishimura, using
// a generator by Knuth. See original source for details.
for (m_mti = 1; m_mti < PERIOD_N; ++m_mti)
m_mt[m_mti] = 0xffffffffUL & (1812433253UL * (m_mt[m_mti - 1] ^ (m_mt[m_mti - 1] >> 30)) + m_mti);
}
//--------------------------------------------------------------------------
// Initialize the generator from an array of seeds
void mtprng::init_by_array(const std::vector<int_type> & init_key)
{
// Note: variable names match those in original example
size_t i = 1;
size_t j = 0;
size_t k = static_cast<size_t>((PERIOD_N > init_key.size()) ? PERIOD_N : init_key.size());
init(19650218UL);
for (; k; --k)
{
m_mt[i] = (m_mt[i] ^ ((m_mt[i - 1] ^ (m_mt[i - 1] >> 30)) * 1664525UL)) + init_key[j] + j;
m_mt[i] &= 0xffffffffUL;
++i;
++j;
if (i >= PERIOD_N)
{
m_mt[0] = m_mt[PERIOD_N - 1];
i = 1;
}
if (j >= static_cast<size_t>(init_key.size()))
j = 0;
}
for (k = PERIOD_N - 1; k; --k)
{
m_mt[i] = (m_mt[i] ^ ((m_mt[i - 1] ^ (m_mt[i - 1] >> 30)) * 1566083941UL)) - i;
m_mt[i] &= 0xffffffffUL;
++i;
if (i >= PERIOD_N)
{
m_mt[0] = m_mt[PERIOD_N - 1];
i = 1;
}
}
m_mt[0] = 0x80000000UL;
}
//--------------------------------------------------------------------------
// Obtain the next 32-bit integer in the psuedo-random sequence
mtprng::int_type mtprng::get_rand()
{
// Note: variable names match those in original example
const int_type mag01[2] = { 0, MATRIX_A };
int_type y;
size_t kk;
// Generate N words at a time
if (m_mti >= PERIOD_N)
{
// Fill the m_mt array
for (kk=0; kk < PERIOD_N-PERIOD_M; kk++)
{
y = (m_mt[kk] & UPPER_MASK) | (m_mt[kk+1] & LOWER_MASK);
m_mt[kk] = m_mt[kk + PERIOD_M] ^ (y >> 1) ^ mag01[y & 0x1];
}
for ( ; kk < PERIOD_N-1; kk++)
{
y = (m_mt[kk] & UPPER_MASK) | (m_mt[kk+1] & LOWER_MASK);
m_mt[kk] = m_mt[kk+(PERIOD_M-PERIOD_N)] ^ (y >> 1) ^ mag01[y & 0x1];
}
y = (m_mt[PERIOD_N-1] & UPPER_MASK) | (m_mt[0]&LOWER_MASK);
m_mt[PERIOD_N-1] = m_mt[PERIOD_M-1] ^ (y >> 1) ^ mag01[y & 0x1];
m_mti = 0;
}
// Here is where we actually calculate the number with a series of transformations
y = m_mt[m_mti++];
y ^= (y >> 11);
y ^= (y << 7) & 0x9d2c5680UL;
y ^= (y << 15) & 0xefc60000UL;
y ^= (y >> 18);
return y;
}
//--------------------------------------------------------------------------
// Obtain a psuedorandom real number in the range [0,1), i.e., a number
// greater than or equal to 0 and less than 1, with 53-bit precision.
mtprng::real_type mtprng::get_rand_real53()
{
// privides a granularity of approx. 1.1E-16
int_type a = get_rand() >> 5;
int_type b = get_rand() >> 6;
return real_type(a * 67108864.0 + b) * (1.0 / 9007199254740992.0);
}
|
