1 /*
2 A C-program for MT19937, with improved initialization 2002/1/26.
3
4 This is an optimized version that amortizes the shift/reload cost,
5 by Eric Landry 2004-03-15.
6
7 Before using, initialize the state by using init_genrand(seed) or
8 init_by_array(init_key, key_length).
9
10 Copyright (C) 1997--2004, Makoto Matsumoto, Takuji Nishimura, and
11 Eric Landry; All rights reserved.
12
13 Redistribution and use in source and binary forms, with or without
14 modification, are permitted provided that the following conditions
15 are met:
16
17 1. Redistributions of source code must retain the above copyright
18 notice, this list of conditions and the following disclaimer.
19
20 2. Redistributions in binary form must reproduce the above copyright
21 notice, this list of conditions and the following disclaimer
22 in the documentation and/or other materials provided with the
23 distribution.
24
25 3. The names of its contributors may not be used to endorse or
26 promote products derived from this software without specific
27 prior written permission.
28
29 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
30 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
31 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
32 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
33 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
34 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
35 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
36 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
37 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
38 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
39 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
40
41 Any feedback is very welcome.
42 http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
43 email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
44
45 Reference: M. Matsumoto and T. Nishimura, "Mersenne Twister:
46 A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number
47 Generator", ACM Transactions on Modeling and Computer Simulation,
48 Vol. 8, No. 1, January 1998, pp 3--30.
49 */
50
51 #include "tm_mt.h"
52
53 /* Period parameters */
54 #define N 624
55 #define M 397
56 #define MATRIX_A 0x9908b0dfUL /* constant vector a */
57 #define UPPER_MASK 0x80000000UL /* most significant w-r bits */
58 #define LOWER_MASK 0x7fffffffUL /* least significant r bits */
59
60 static unsigned long x[N]; /* the array for the state vector */
61 static unsigned long *p0, *p1, *pm;
62
63 /*
64 initialize with a seed
65
66 See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier.
67
68 In the previous versions, MSBs of the seed affect only MSBs of
69 the state.
70
71 2002-01-09 modified by Makoto Matsumoto
72 */
73 void
74 init_genrand(unsigned long s)
75 {
76 int i;
77
78 x[0] = s & 0xffffffffUL;
79 for (i = 1; i < N; ++i) {
80 x[i] = (1812433253UL * (x[i - 1] ^ (x[i - 1] >> 30)) + i)
81 & 0xffffffffUL; /* for >32 bit machines */
82 }
83 p0 = x;
84 p1 = x + 1;
85 pm = x + M;
86 }
87
88 /*
89 initialize by an array with array-length
90
91 init_key is the array for initializing keys
92
93 key_length is its length
94
95 2004-02-26 slight change for C++
96 */
97 void
98 init_by_array(unsigned long init_key[], int key_length)
99 {
100 int i, j, k;
101
102 init_genrand(19650218UL);
103 i = 1;
104 j = 0;
105 for (k = (N > key_length ? N : key_length); k; --k) {
106 /* non linear */
107 x[i] = ((x[i] ^ ((x[i - 1] ^ (x[i - 1] >> 30)) * 1664525UL))
108 + init_key[j] + j) & 0xffffffffUL; /* for WORDSIZE > 32 machines */
109 if (++i >= N) {
110 x[0] = x[N - 1];
111 i = 1;
112 }
113 if (++j >= key_length) {
114 j = 0;
115 }
116 }
117 for (k = N - 1; k; --k) {
118 /* non linear */
119 x[i] = ((x[i] ^ ((x[i - 1] ^ (x[i - 1] >> 30)) * 1566083941UL)) - i)
120 & 0xffffffffUL; /* for WORDSIZE > 32 machines */
121 if (++i >= N) {
122 x[0] = x[N - 1];
123 i = 1;
124 }
125 }
126 x[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
127 }
128
129 /* generates a random number on the interval [0,0xffffffff] */
130 unsigned long
131 genrand_int32(void)
132 {
133 unsigned long y;
134
135 if (!p0) {
136 /* Default seed */
137 init_genrand(5489UL);
138 }
139 /* Twisted feedback */
140 y = *p0 = *pm++ ^ (((*p0 & UPPER_MASK) | (*p1 & LOWER_MASK)) >> 1)
141 ^ (-(*p1 & 1) & MATRIX_A);
142 p0 = p1++;
143 if (pm == x + N) {
144 pm = x;
145 }
146 if (p1 == x + N) {
147 p1 = x;
148 }
149 /* Temper */
150 y ^= y >> 11;
151 y ^= y << 7 & 0x9d2c5680UL;
152 y ^= y << 15 & 0xefc60000UL;
153 y ^= y >> 18;
154 return y;
155 }
156
157 /* generates a random number on the interval [0,0x7fffffff] */
158 long
159 genrand_int31(void)
160 {
161 return (long) (genrand_int32() >> 1);
162 }
163
164 /* generates a random number on the real interval [0,1] */
165 double
166 genrand_real1(void)
167 {
168 return genrand_int32() * (1.0 / 4294967295.0);
169 /* divided by 2^32-1 */
170 }
171
172 /* generates a random number on the real interval [0,1) */
173 double
174 genrand_real2(void)
175 {
176 return genrand_int32() * (1.0 / 4294967296.0);
177 /* divided by 2^32 */
178 }
179
180 /* generates a random number on the real interval (0,1) */
181 double
182 genrand_real3(void)
183 {
184 return (((double) genrand_int32()) + 0.5) * (1.0 / 4294967296.0);
185 /* divided by 2^32 */
186 }
187
188 /* generates a 53-bit random number on the real interval [0,1) */
189 double
190 genrand_res53(void)
191 {
192 unsigned long a = genrand_int32() >> 5, b = genrand_int32() >> 6;
193
194 return (a * 67108864.0 + b) * (1.0 / 9007199254740992.0);
195 }
196
197 /* 2002-01-09 These real versions are due to Isaku Wada */
198