## root / include / mersennetwister.h @ aeae20a1

History | View | Annotate | Download (14.6 KB)

1 | 01873262 | Georg Kunz | ```
// MersenneTwister.h
``` |
---|---|---|---|

2 | ```
// Mersenne Twister random number generator -- a C++ class MTRand
``` |
||

3 | ```
// Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
``` |
||

4 | ```
// Richard J. Wagner v1.0 15 May 2003 rjwagner@writeme.com
``` |
||

5 | |||

6 | ```
// The Mersenne Twister is an algorithm for generating random numbers. It
``` |
||

7 | ```
// was designed with consideration of the flaws in various other generators.
``` |
||

8 | ```
// The period, 2^19937-1, and the order of equidistribution, 623 dimensions,
``` |
||

9 | ```
// are far greater. The generator is also fast; it avoids multiplication and
``` |
||

10 | ```
// division, and it benefits from caches and pipelines. For more information
``` |
||

11 | ```
// see the inventors' web page at http://www.math.keio.ac.jp/~matumoto/emt.html
``` |
||

12 | |||

13 | ```
// Reference
``` |
||

14 | ```
// M. Matsumoto and T. Nishimura, "Mersenne Twister: A 623-Dimensionally
``` |
||

15 | ```
// Equidistributed Uniform Pseudo-Random Number Generator", ACM Transactions on
``` |
||

16 | ```
// Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, pp 3-30.
``` |
||

17 | |||

18 | ```
// Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
``` |
||

19 | ```
// Copyright (C) 2000 - 2003, Richard J. Wagner
``` |
||

20 | ```
// All rights reserved.
``` |
||

21 | ```
//
``` |
||

22 | ```
// Redistribution and use in source and binary forms, with or without
``` |
||

23 | ```
// modification, are permitted provided that the following conditions
``` |
||

24 | ```
// are met:
``` |
||

25 | ```
//
``` |
||

26 | ```
// 1. Redistributions of source code must retain the above copyright
``` |
||

27 | ```
// notice, this list of conditions and the following disclaimer.
``` |
||

28 | ```
//
``` |
||

29 | ```
// 2. Redistributions in binary form must reproduce the above copyright
``` |
||

30 | ```
// notice, this list of conditions and the following disclaimer in the
``` |
||

31 | ```
// documentation and/or other materials provided with the distribution.
``` |
||

32 | ```
//
``` |
||

33 | ```
// 3. The names of its contributors may not be used to endorse or promote
``` |
||

34 | ```
// products derived from this software without specific prior written
``` |
||

35 | ```
// permission.
``` |
||

36 | ```
//
``` |
||

37 | ```
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
``` |
||

38 | ```
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
``` |
||

39 | ```
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
``` |
||

40 | ```
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
``` |
||

41 | ```
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
``` |
||

42 | ```
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
``` |
||

43 | ```
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
``` |
||

44 | ```
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
``` |
||

45 | ```
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
``` |
||

46 | ```
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
``` |
||

47 | ```
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
``` |
||

48 | |||

49 | ```
// The original code included the following notice:
``` |
||

50 | ```
//
``` |
||

51 | ```
// When you use this, send an email to: matumoto@math.keio.ac.jp
``` |
||

52 | ```
// with an appropriate reference to your work.
``` |
||

53 | ```
//
``` |
||

54 | ```
// It would be nice to CC: rjwagner@writeme.com and Cokus@math.washington.edu
``` |
||

55 | ```
// when you write.
``` |
||

56 | |||

57 | ```
#ifndef MERSENNETWISTER_H
``` |
||

58 | ```
#define MERSENNETWISTER_H
``` |
||

59 | |||

60 | ```
// Not thread safe (unless auto-initialization is avoided and each thread has
``` |
||

61 | ```
// its own MTRand object)
``` |
||

62 | |||

63 | #include <iostream> |
||

64 | #include <limits.h> |
||

65 | #include <stdio.h> |
||

66 | #include <time.h> |
||

67 | #include <math.h> |
||

68 | |||

69 | NAMESPACE_BEGIN |
||

70 | |||

71 | class MTRand { |
||

72 | ```
// Data
``` |
||

73 | ```
public:
``` |
||

74 | typedef unsigned long uint32; // unsigned integer type, at least 32 bits |
||

75 | |||

76 | enum Dummy1 { N = 624 }; // length of state vector |
||

77 | enum Dummy2 { SAVE = N + 1 }; // length of array for save() |
||

78 | ```
// Note: DummyX names needed by buggy gcc 4.0.1 on OS/X (Andras)
``` |
||

79 | |||

80 | ```
protected:
``` |
||

81 | enum Dummy3 { M = 397 }; // period parameter |
||

82 | |||

83 | ```
uint32 state[N]; // internal state
``` |
||

84 | ```
uint32 *pNext; // next value to get from state
``` |
||

85 | int left; // number of values left before reload needed |
||

86 | |||

87 | |||

88 | ```
//Methods
``` |
||

89 | ```
public:
``` |
||

90 | MTRand( const uint32& oneSeed ); // initialize with a simple uint32 |
||

91 | MTRand( uint32 *const bigSeed, uint32 const seedLength = N ); // or an array |
||

92 | ```
MTRand(); // auto-initialize with /dev/urandom or time() and clock()
``` |
||

93 | |||

94 | ```
// Do NOT use for CRYPTOGRAPHY without securely hashing several returned
``` |
||

95 | ```
// values together, otherwise the generator state can be learned after
``` |
||

96 | ```
// reading 624 consecutive values.
``` |
||

97 | |||

98 | ```
// Access to 32-bit random numbers
``` |
||

99 | double rand(); // real number in [0,1] |
||

100 | double rand( const double& n ); // real number in [0,n] |
||

101 | double randExc(); // real number in [0,1) |
||

102 | double randExc( const double& n ); // real number in [0,n) |
||

103 | double randDblExc(); // real number in (0,1) |
||

104 | double randDblExc( const double& n ); // real number in (0,n) |
||

105 | ```
uint32 randInt(); // integer in [0,2^32-1]
``` |
||

106 | uint32 randInt( const uint32& n ); // integer in [0,n] for n < 2^32 |
||

107 | double operator()() { return rand(); } // same as rand() |
||

108 | |||

109 | ```
// Access to 53-bit random numbers (capacity of IEEE double precision)
``` |
||

110 | double rand53(); // real number in [0,1) |
||

111 | |||

112 | ```
// Access to nonuniform random number distributions
``` |
||

113 | double randNorm( const double& mean = 0.0, const double& variance = 0.0 ); |
||

114 | |||

115 | ```
// Re-seeding functions with same behavior as initializers
``` |
||

116 | void seed( const uint32 oneSeed ); |
||

117 | void seed( uint32 *const bigSeed, const uint32 seedLength = N ); |
||

118 | ```
void seed();
``` |
||

119 | |||

120 | ```
// Saving and loading generator state
``` |
||

121 | void save( uint32* saveArray ) const; // to array of size SAVE |
||

122 | void load( uint32 *const loadArray ); // from such array |
||

123 | ```
friend std::ostream& operator<<( std::ostream& os, const MTRand& mtrand );
``` |
||

124 | friend std::istream& operator>>( std::istream& is, MTRand& mtrand ); |
||

125 | |||

126 | ```
protected:
``` |
||

127 | void initialize( const uint32 oneSeed ); |
||

128 | ```
void reload();
``` |
||

129 | uint32 hiBit( const uint32& u ) const { return u & 0x80000000UL; } |
||

130 | uint32 loBit( const uint32& u ) const { return u & 0x00000001UL; } |
||

131 | uint32 loBits( const uint32& u ) const { return u & 0x7fffffffUL; } |
||

132 | uint32 mixBits( const uint32& u, const uint32& v ) const |
||

133 | ```
{ return hiBit(u) | loBits(v); }
``` |
||

134 | uint32 twist( const uint32& m, const uint32& s0, const uint32& s1 ) const |
||

135 | { return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfUL); } |
||

136 | ```
static uint32 hash( time_t t, clock_t c );
``` |
||

137 | }; |
||

138 | |||

139 | |||

140 | inline MTRand::MTRand( const uint32& oneSeed ) |
||

141 | { seed(oneSeed); } |
||

142 | |||

143 | inline MTRand::MTRand( uint32 *const bigSeed, const uint32 seedLength ) |
||

144 | { seed(bigSeed,seedLength); } |
||

145 | |||

146 | ```
inline MTRand::MTRand()
``` |
||

147 | { seed(); } |
||

148 | |||

149 | inline double MTRand::rand() |
||

150 | { return double(randInt()) * (1.0/4294967295.0); } |
||

151 | |||

152 | inline double MTRand::rand( const double& n ) |
||

153 | ```
{ return rand() * n; }
``` |
||

154 | |||

155 | inline double MTRand::randExc() |
||

156 | { return double(randInt()) * (1.0/4294967296.0); } |
||

157 | |||

158 | inline double MTRand::randExc( const double& n ) |
||

159 | ```
{ return randExc() * n; }
``` |
||

160 | |||

161 | inline double MTRand::randDblExc() |
||

162 | { return ( double(randInt()) + 0.5 ) * (1.0/4294967296.0); } |
||

163 | |||

164 | inline double MTRand::randDblExc( const double& n ) |
||

165 | ```
{ return randDblExc() * n; }
``` |
||

166 | |||

167 | inline double MTRand::rand53() |
||

168 | { |
||

169 | uint32 a = randInt() >> 5, b = randInt() >> 6; |
||

170 | return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0); // by Isaku Wada |
||

171 | } |
||

172 | |||

173 | inline double MTRand::randNorm( const double& mean, const double& variance ) |
||

174 | { |
||

175 | ```
// Return a real number from a normal (Gaussian) distribution with given
``` |
||

176 | ```
// mean and variance by Box-Muller method
``` |
||

177 | double r = sqrt( -2.0 * log( 1.0-randDblExc()) ) * variance; |
||

178 | double phi = 2.0 * 3.14159265358979323846264338328 * randExc(); |
||

179 | ```
return mean + r * cos(phi);
``` |
||

180 | } |
||

181 | |||

182 | ```
inline MTRand::uint32 MTRand::randInt()
``` |
||

183 | { |
||

184 | ```
// Pull a 32-bit integer from the generator state
``` |
||

185 | ```
// Every other access function simply transforms the numbers extracted here
``` |
||

186 | |||

187 | if( left == 0 ) reload(); |
||

188 | --left; |
||

189 | |||

190 | ```
register uint32 s1;
``` |
||

191 | s1 = *pNext++; |
||

192 | ```
s1 ^= (s1 >> 11);
``` |
||

193 | s1 ^= (s1 << 7) & 0x9d2c5680UL; |
||

194 | s1 ^= (s1 << 15) & 0xefc60000UL; |
||

195 | return ( s1 ^ (s1 >> 18) ); |
||

196 | } |
||

197 | |||

198 | inline MTRand::uint32 MTRand::randInt( const uint32& n ) |
||

199 | { |
||

200 | ```
// Find which bits are used in n
``` |
||

201 | ```
// Optimized by Magnus Jonsson (magnus@smartelectronix.com)
``` |
||

202 | uint32 used = n; |
||

203 | ```
used |= used >> 1;
``` |
||

204 | ```
used |= used >> 2;
``` |
||

205 | ```
used |= used >> 4;
``` |
||

206 | ```
used |= used >> 8;
``` |
||

207 | ```
used |= used >> 16;
``` |
||

208 | |||

209 | ```
// Draw numbers until one is found in [0,n]
``` |
||

210 | uint32 i; |
||

211 | ```
do
``` |
||

212 | ```
i = randInt() & used; // toss unused bits to shorten search
``` |
||

213 | ```
while( i > n );
``` |
||

214 | ```
return i;
``` |
||

215 | } |
||

216 | |||

217 | |||

218 | inline void MTRand::seed( const uint32 oneSeed ) |
||

219 | { |
||

220 | ```
// Seed the generator with a simple uint32
``` |
||

221 | initialize(oneSeed); |
||

222 | reload(); |
||

223 | } |
||

224 | |||

225 | |||

226 | inline void MTRand::seed( uint32 *const bigSeed, const uint32 seedLength ) |
||

227 | { |
||

228 | ```
// Seed the generator with an array of uint32's
``` |
||

229 | ```
// There are 2^19937-1 possible initial states. This function allows
``` |
||

230 | ```
// all of those to be accessed by providing at least 19937 bits (with a
``` |
||

231 | ```
// default seed length of N = 624 uint32's). Any bits above the lower 32
``` |
||

232 | ```
// in each element are discarded.
``` |
||

233 | ```
// Just call seed() if you want to get array from /dev/urandom
``` |
||

234 | ```
initialize(19650218UL);
``` |
||

235 | register int i = 1; |
||

236 | register uint32 j = 0; |
||

237 | register int k = ( N > seedLength ? N : seedLength ); |
||

238 | ```
for( ; k; --k )
``` |
||

239 | { |
||

240 | state[i] = |
||

241 | state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1664525UL ); |
||

242 | ```
state[i] += ( bigSeed[j] & 0xffffffffUL ) + j;
``` |
||

243 | ```
state[i] &= 0xffffffffUL;
``` |
||

244 | ++i; ++j; |
||

245 | if( i >= N ) { state[0] = state[N-1]; i = 1; } |
||

246 | if( j >= seedLength ) j = 0; |
||

247 | } |
||

248 | for( k = N - 1; k; --k ) |
||

249 | { |
||

250 | state[i] = |
||

251 | state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL ); |
||

252 | state[i] -= i; |
||

253 | ```
state[i] &= 0xffffffffUL;
``` |
||

254 | ++i; |
||

255 | if( i >= N ) { state[0] = state[N-1]; i = 1; } |
||

256 | } |
||

257 | state[0] = 0x80000000UL; // MSB is 1, assuring non-zero initial array |
||

258 | reload(); |
||

259 | } |
||

260 | |||

261 | |||

262 | inline void MTRand::seed() |
||

263 | { |
||

264 | ```
// Seed the generator with an array from /dev/urandom if available
``` |
||

265 | ```
// Otherwise use a hash of time() and clock() values
``` |
||

266 | |||

267 | ```
// First try getting an array from /dev/urandom
``` |
||

268 | FILE* urandom = fopen( "/dev/urandom", "rb" ); |
||

269 | ```
if( urandom )
``` |
||

270 | { |
||

271 | uint32 bigSeed[N]; |
||

272 | ```
register uint32 *s = bigSeed;
``` |
||

273 | register int i = N; |
||

274 | register bool success = true; |
||

275 | ```
while( success && i-- )
``` |
||

276 | success = fread( s++, sizeof(uint32), 1, urandom ); |
||

277 | fclose(urandom); |
||

278 | if( success ) { seed( bigSeed, N ); return; } |
||

279 | } |
||

280 | |||

281 | ```
// Was not successful, so use time() and clock() instead
``` |
||

282 | ```
seed( hash( time(NULL), clock() ) );
``` |
||

283 | } |
||

284 | |||

285 | |||

286 | inline void MTRand::initialize( const uint32 seed ) |
||

287 | { |
||

288 | ```
// Initialize generator state with seed
``` |
||

289 | ```
// See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier.
``` |
||

290 | ```
// In previous versions, most significant bits (MSBs) of the seed affect
``` |
||

291 | ```
// only MSBs of the state array. Modified 9 Jan 2002 by Makoto Matsumoto.
``` |
||

292 | ```
register uint32 *s = state;
``` |
||

293 | ```
register uint32 *r = state;
``` |
||

294 | register int i = 1; |
||

295 | ```
*s++ = seed & 0xffffffffUL;
``` |
||

296 | ```
for( ; i < N; ++i )
``` |
||

297 | { |
||

298 | *s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL; |
||

299 | r++; |
||

300 | } |
||

301 | } |
||

302 | |||

303 | |||

304 | inline void MTRand::reload() |
||

305 | { |
||

306 | ```
// Generate N new values in state
``` |
||

307 | ```
// Made clearer and faster by Matthew Bellew (matthew.bellew@home.com)
``` |
||

308 | ```
register uint32 *p = state;
``` |
||

309 | register int i; |
||

310 | ```
for( i = N - M; i--; ++p )
``` |
||

311 | *p = twist( p[M], p[0], p[1] ); |
||

312 | ```
for( i = M; --i; ++p )
``` |
||

313 | *p = twist( p[M-N], p[0], p[1] ); |
||

314 | *p = twist( p[M-N], p[0], state[0] ); |
||

315 | |||

316 | left = N, pNext = state; |
||

317 | } |
||

318 | |||

319 | |||

320 | ```
inline MTRand::uint32 MTRand::hash( time_t t, clock_t c )
``` |
||

321 | { |
||

322 | ```
// Get a uint32 from t and c
``` |
||

323 | ```
// Better than uint32(x) in case x is floating point in [0,1]
``` |
||

324 | ```
// Based on code by Lawrence Kirby (fred@genesis.demon.co.uk)
``` |
||

325 | |||

326 | static uint32 differ = 0; // guarantee time-based seeds will change |
||

327 | |||

328 | ```
uint32 h1 = 0;
``` |
||

329 | unsigned char *p = (unsigned char *) &t; |
||

330 | for( size_t i = 0; i < sizeof(t); ++i ) |
||

331 | { |
||

332 | ```
h1 *= UCHAR_MAX + 2U;
``` |
||

333 | h1 += p[i]; |
||

334 | } |
||

335 | ```
uint32 h2 = 0;
``` |
||

336 | p = (unsigned char *) &c; |
||

337 | for( size_t j = 0; j < sizeof(c); ++j ) |
||

338 | { |
||

339 | ```
h2 *= UCHAR_MAX + 2U;
``` |
||

340 | h2 += p[j]; |
||

341 | } |
||

342 | ```
return ( h1 + differ++ ) ^ h2;
``` |
||

343 | } |
||

344 | |||

345 | |||

346 | inline void MTRand::save( uint32* saveArray ) const |
||

347 | { |
||

348 | ```
register uint32 *sa = saveArray;
``` |
||

349 | register const uint32 *s = state; |
||

350 | register int i = N; |
||

351 | ```
for( ; i--; *sa++ = *s++ ) {}
``` |
||

352 | *sa = left; |
||

353 | } |
||

354 | |||

355 | |||

356 | inline void MTRand::load( uint32 *const loadArray ) |
||

357 | { |
||

358 | ```
register uint32 *s = state;
``` |
||

359 | ```
register uint32 *la = loadArray;
``` |
||

360 | register int i = N; |
||

361 | ```
for( ; i--; *s++ = *la++ ) {}
``` |
||

362 | left = *la; |
||

363 | pNext = &state[N-left]; |
||

364 | } |
||

365 | |||

366 | |||

367 | inline std::ostream& operator<<( std::ostream& os, const MTRand& mtrand ) |
||

368 | { |
||

369 | register const MTRand::uint32 *s = mtrand.state; |
||

370 | register int i = mtrand.N; |
||

371 | for( ; i--; os << *s++ << "\t" ) {} |
||

372 | ```
return os << mtrand.left;
``` |
||

373 | } |
||

374 | |||

375 | |||

376 | ```
inline std::istream& operator>>( std::istream& is, MTRand& mtrand )
``` |
||

377 | { |
||

378 | ```
register MTRand::uint32 *s = mtrand.state;
``` |
||

379 | register int i = mtrand.N; |
||

380 | ```
for( ; i--; is >> *s++ ) {}
``` |
||

381 | is >> mtrand.left; |
||

382 | mtrand.pNext = &mtrand.state[mtrand.N-mtrand.left]; |
||

383 | ```
return is;
``` |
||

384 | } |
||

385 | |||

386 | NAMESPACE_END |
||

387 | |||

388 | |||

389 | #endif // MERSENNETWISTER_H |
||

390 | |||

391 | ```
// Change log:
``` |
||

392 | ```
//
``` |
||

393 | ```
// v0.1 - First release on 15 May 2000
``` |
||

394 | ```
// - Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
``` |
||

395 | ```
// - Translated from C to C++
``` |
||

396 | ```
// - Made completely ANSI compliant
``` |
||

397 | ```
// - Designed convenient interface for initialization, seeding, and
``` |
||

398 | ```
// obtaining numbers in default or user-defined ranges
``` |
||

399 | ```
// - Added automatic seeding from /dev/urandom or time() and clock()
``` |
||

400 | ```
// - Provided functions for saving and loading generator state
``` |
||

401 | ```
//
``` |
||

402 | ```
// v0.2 - Fixed bug which reloaded generator one step too late
``` |
||

403 | ```
//
``` |
||

404 | ```
// v0.3 - Switched to clearer, faster reload() code from Matthew Bellew
``` |
||

405 | ```
//
``` |
||

406 | ```
// v0.4 - Removed trailing newline in saved generator format to be consistent
``` |
||

407 | ```
// with output format of built-in types
``` |
||

408 | ```
//
``` |
||

409 | ```
// v0.5 - Improved portability by replacing static const int's with enum's and
``` |
||

410 | ```
// clarifying return values in seed(); suggested by Eric Heimburg
``` |
||

411 | ```
// - Removed MAXINT constant; use 0xffffffffUL instead
``` |
||

412 | ```
//
``` |
||

413 | ```
// v0.6 - Eliminated seed overflow when uint32 is larger than 32 bits
``` |
||

414 | ```
// - Changed integer [0,n] generator to give better uniformity
``` |
||

415 | ```
//
``` |
||

416 | ```
// v0.7 - Fixed operator precedence ambiguity in reload()
``` |
||

417 | ```
// - Added access for real numbers in (0,1) and (0,n)
``` |
||

418 | ```
//
``` |
||

419 | ```
// v0.8 - Included time.h header to properly support time_t and clock_t
``` |
||

420 | ```
//
``` |
||

421 | ```
// v1.0 - Revised seeding to match 26 Jan 2002 update of Nishimura and Matsumoto
``` |
||

422 | ```
// - Allowed for seeding with arrays of any length
``` |
||

423 | ```
// - Added access for real numbers in [0,1) with 53-bit resolution
``` |
||

424 | ```
// - Added access for real numbers from normal (Gaussian) distributions
``` |
||

425 | ```
// - Increased overall speed by optimizing twist()
``` |
||

426 | ```
// - Doubled speed of integer [0,n] generation
``` |
||

427 | ```
// - Fixed out-of-range number generation on 64-bit machines
``` |
||

428 | ```
// - Improved portability by substituting literal constants for long enum's
``` |
||

429 | `// - Changed license from GNU LGPL to BSD` |