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
|
/* A C-program for MT19937: Integer version (1999/10/28) */
/* genrand() generates one pseudorandom unsigned integer (32bit) */
/* which is uniformly distributed among 0 to 2^32-1 for each */
/* call. sgenrand(seed) sets initial values to the working area */
/* of 624 words. Before genrand(), sgenrand(seed) must be */
/* called once. (seed is any 32-bit integer.) */
/* Coded by Takuji Nishimura, considering the suggestions by */
/* Topher Cooper and Marc Rieffel in July-Aug. 1997. */
/* This library is free software under the Artistic license: */
/* see the file COPYING distributed together with this code. */
/* For the verification of the code, its output sequence file */
/* mt19937int.out is attached (2001/4/2) */
/* Copyright (C) 1997, 1999 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 */
/* REFERENCE */
/* M. Matsumoto and T. Nishimura, */
/* "Mersenne Twister: A 623-Dimensionally Equidistributed Uniform */
/* Pseudo-Random Number Generator", */
/* ACM Transactions on Modeling and Computer Simulation, */
/* Vol. 8, No. 1, January 1998, pp 3--30. */
#include<stdio.h>
#include "mt19937int.h"
/* Period parameters */
#define N 624
#define M 397
#define MATRIX_A 0x9908b0df /* constant vector a */
#define UPPER_MASK 0x80000000 /* most significant w-r bits */
#define LOWER_MASK 0x7fffffff /* least significant r bits */
/* Tempering parameters */
#define TEMPERING_MASK_B 0x9d2c5680
#define TEMPERING_MASK_C 0xefc60000
#define TEMPERING_SHIFT_U(y) (y >> 11)
#define TEMPERING_SHIFT_S(y) (y << 7)
#define TEMPERING_SHIFT_T(y) (y << 15)
#define TEMPERING_SHIFT_L(y) (y >> 18)
static unsigned long mt[N]; /* the array for the state vector */
static int mti=N+1; /* mti==N+1 means mt[N] is not initialized */
/* Initializing the array with a seed */
void sgenrand(signed long seed)
{
int i;
for (i=0;i<N;i++) {
mt[i] = seed & 0xffff0000;
seed = 69069 * seed + 1;
mt[i] |= (seed & 0xffff0000) >> 16;
seed = 69069 * seed + 1;
}
mti = N;
}
/* Initialization by "sgenrand()" is an example. Theoretically, */
/* there are 2^19937-1 possible states as an intial state. */
/* This function allows to choose any of 2^19937-1 ones. */
/* Essential bits in "seed_array[]" is following 19937 bits: */
/* (seed_array[0]&UPPER_MASK), seed_array[1], ..., seed_array[N-1]. */
/* (seed_array[0]&LOWER_MASK) is discarded. */
/* Theoretically, */
/* (seed_array[0]&UPPER_MASK), seed_array[1], ..., seed_array[N-1] */
/* can take any values except all zeros. */
void lsgenrand(unsigned long *seed_array)
/* the length of seed_array[] must be at least N */
{
int i;
for (i=0;i<N;i++)
mt[i] = seed_array[i];
mti=N;
}
unsigned long genrand()
{
unsigned long y;
static unsigned long mag01[2]={0x0, MATRIX_A};
/* mag01[x] = x * MATRIX_A for x=0,1 */
if (mti >= N) { /* generate N words at one time */
int kk;
if (mti == N+1) /* if sgenrand() has not been called, */
sgenrand(4357); /* a default initial seed is used */
for (kk=0;kk<N-M;kk++) {
y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1];
}
for (;kk<N-1;kk++) {
y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
mt[kk] = mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1];
}
y = (mt[N-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
mt[N-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1];
mti = 0;
}
y = mt[mti++];
y ^= TEMPERING_SHIFT_U(y);
y ^= TEMPERING_SHIFT_S(y) & TEMPERING_MASK_B;
y ^= TEMPERING_SHIFT_T(y) & TEMPERING_MASK_C;
y ^= TEMPERING_SHIFT_L(y);
return y;
}
/* This main() outputs first 1000 generated numbers. */
/*
main()
{
int i;
sgenrand(4357);
for (i=0; i<1000; i++) {
printf("%10lu ", genrand());
if (i%5==4) printf("\n");
}
}
*/
|
