1
2 /*
3 A C-program for MT19937, with improved initialization 2002/1/26.
4
5 This is an optimized version that amortizes the shift/reload cost,
6 by Eric Landry 2004-03-15.
7
8 Before using, initialize the state by using RandomGen_setSeed(seed) or
9 init_by_array(init_key, key_length).
10
11 Copyright (C) 1997--2004, Makoto Matsumoto, Takuji Nishimura, and
12 Eric Landry; All rights reserved.
13
14 Redistribution and use in source and binary forms, with or without
15 modification, are permitted provided that the following conditions
16 are met:
17
18 1. Redistributions of source code must retain the above copyright
19 notice, this list of conditions and the following disclaimer.
20
21 2. Redistributions in binary form must reproduce the above copyright
22 notice, this list of conditions and the following disclaimer
23 in the documentation and/or other materials provided with the
24 distribution.
25
26 3. The names of its contributors may not be used to endorse or
27 promote products derived from this software without specific
28 prior written permission.
29
30 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
31 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
32 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
33 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
34 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
35 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
36 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
37 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
38 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
39 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
40 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
41
42 Any feedback is very welcome.
43 http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html
44 email: m-mat @ math.sci.hiroshima-u.ac.jp (remove space)
45
46 Reference: M. Matsumoto and T. Nishimura, "Mersenne Twister:
47 A 623-Dimensionally Equidistributed Uniform Pseudo-RandomGen Number
48 Generator", ACM Transactions on Modeling and Computer Simulation,
49 Vol. 8, No. 1, January 1998, pp 3--30.
50 */
51
52 #include "Base.h"
53 #include "RandomGen.h"
54
55 /* Period parameters */
56 #define N RANDOMGEN_N
57 #define M 397
58 #define MATRIX_A 0x9908b0dfUL /* constant vector a */
59 #define UPPER_MASK 0x80000000UL /* most significant w-r bits */
60 #define LOWER_MASK 0x7fffffffUL /* least significant r bits */
61
62 /*static unsigned long self->mt[N];*/ /* the array for the state vector */
63 /*static int self->mti=N+1;*/ /* self->mti==N+1 means self->mt[N] is not initialized */
64
65 /* initializes self->mt[N] with a seed */
66
init_genrand(RandomGen * self,unsigned long s)67 static void init_genrand(RandomGen *self, unsigned long s)
68 {
69 self->mt[0]= s & 0xffffffffUL;
70 for (self->mti=1; self->mti<N; self->mti ++)
71 {
72 self->mt[self->mti] =
73 (1812433253UL * (self->mt[self->mti-1] ^ (self->mt[self->mti-1] >> 30)) + self->mti);
74 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
75 /* In the previous versions, MSBs of the seed affect */
76 /* only MSBs of the array self->mt[]. */
77 /* 2002/01/09 modified by Makoto Matsumoto */
78 self->mt[self->mti] &= 0xffffffffUL;
79 /* for >32 bit machines */
80 }
81 }
82
RandomGen_setSeed(RandomGen * self,unsigned long seed)83 void RandomGen_setSeed(RandomGen *self, unsigned long seed)
84 {
85 init_genrand(self, seed);
86 }
87
88 #include <time.h>
89
RandomGen_chooseRandomSeed(RandomGen * self)90 void RandomGen_chooseRandomSeed(RandomGen *self)
91 {
92 unsigned long seed = 0;
93
94 seed ^= (unsigned long)clock(); // processor time since program start
95 seed ^= (unsigned long)time(NULL); // seconds since 1970
96
97 RandomGen_setSeed(self, seed);
98 }
99
100 /* initialize by an array with array-length */
101 /* init_key is the array for initializing keys */
102 /* key_length is its length */
103 /* slight change for C++, 2004/2/26 */
104 #ifdef EXTRAS
init_by_array(RandomGen * self,unsigned long init_key[],int key_length)105 static void init_by_array(RandomGen *self, unsigned long init_key[],int key_length)
106 {
107 int i, j, k;
108
109 init_genrand(self, 19650218UL);
110 i=1; j=0;
111 k = (N>key_length ? N : key_length);
112
113 for (; k; k--) {
114 self->mt[i] = (self->mt[i] ^ ((self->mt[i-1] ^ (self->mt[i-1] >> 30)) * 1664525UL))
115 + init_key[j] + j; /* non linear */
116 self->mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
117 i ++; j++;
118 if (i>=N) { self->mt[0] = self->mt[N-1]; i=1; }
119 if (j>=key_length) j=0;
120 }
121
122 for (k=N-1; k; k--) {
123 self->mt[i] = (self->mt[i] ^ ((self->mt[i-1] ^ (self->mt[i-1] >> 30)) * 1566083941UL))
124 - i; /* non linear */
125 self->mt[i] &= 0xffffffffUL; /* for WORDSIZE > 32 machines */
126 i ++;
127 if (i>=N) { self->mt[0] = self->mt[N-1]; i=1; }
128 }
129
130 self->mt[0] = 0x80000000UL; /* MSB is 1; assuring non-zero initial array */
131 }
132 #endif
133
134 /* generates a random number on [0,0xffffffff]-interval */
genrand_int32(RandomGen * self)135 static unsigned long genrand_int32(RandomGen *self)
136 {
137 unsigned long y;
138 static unsigned long mag01[2]={0x0UL, MATRIX_A};
139 /* mag01[x] = x * MATRIX_A for x=0,1 */
140
141 if (self->mti >= N) { /* generate N words at one time */
142 int kk;
143
144 if (self->mti == N+1) /* if init_genrand() has not been called, */
145 init_genrand(self, 5489UL); /* a default initial seed is used */
146
147 for (kk=0;kk<N-M;kk++) {
148 y = (self->mt[kk]&UPPER_MASK)|(self->mt[kk+1]&LOWER_MASK);
149 self->mt[kk] = self->mt[kk+M] ^ (y >> 1) ^ mag01[y & 0x1UL];
150 }
151 for (;kk<N-1;kk++) {
152 y = (self->mt[kk]&UPPER_MASK)|(self->mt[kk+1]&LOWER_MASK);
153 self->mt[kk] = self->mt[kk+(M-N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
154 }
155 y = (self->mt[N-1]&UPPER_MASK)|(self->mt[0]&LOWER_MASK);
156 self->mt[N-1] = self->mt[M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
157
158 self->mti = 0;
159 }
160
161 y = self->mt[self->mti ++];
162
163 // Tempering
164 y ^= (y >> 11);
165 y ^= (y << 7) & 0x9d2c5680UL;
166 y ^= (y << 15) & 0xefc60000UL;
167 y ^= (y >> 18);
168
169 return y;
170 }
171
172 /* generates a random number on [0,0x7fffffff]-interval */
173 #ifdef EXTRAS
genrand_int31(RandomGen * self)174 static long genrand_int31(RandomGen *self)
175 {
176 return (long)(genrand_int32(self)>>1);
177 }
178
179 /* generates a random number on [0,1]-real-interval */
genrand_real1(RandomGen * self)180 static double genrand_real1(RandomGen *self)
181 {
182 return genrand_int32(self)*(1.0/4294967295.0);
183 /* divided by 2^32-1 */
184 }
185 #endif
186
187 /* generates a random number on [0,1)-real-interval */
genrand_real2(RandomGen * self)188 static double genrand_real2(RandomGen *self)
189 {
190 return genrand_int32(self)*(1.0/4294967296.0);
191 /* divided by 2^32 */
192 }
193
194 /* generates a random number on (0,1)-real-interval */
195 #ifdef EXTRAS
genrand_real3(RandomGen * self)196 static double genrand_real3(RandomGen *self)
197 {
198 return (((double)genrand_int32(self)) + 0.5)*(1.0/4294967296.0);
199 /* divided by 2^32 */
200 }
201
202 /* generates a random number on [0,1) with 53-bit resolution*/
genrand_res53(RandomGen * self)203 static double genrand_res53(RandomGen *self)
204 {
205 unsigned long a=genrand_int32(self)>>5, b=genrand_int32(self)>>6;
206 return(a*67108864.0+b)*(1.0/9007199254740992.0);
207 }
208 #endif
209
210 /* These real versions are due to Isaku Wada, 2002/01/09 added */
211
212 /*
213 int main(void)
214 {
215 int i;
216 unsigned long init[4]={0x123, 0x234, 0x345, 0x456}, length=4;
217 init_by_array(init, length);
218 printf("1000 outputs of genrand_int32(self)\n");
219 for (i=0; i<1000; i ++) {
220 printf("%10lu ", genrand_int32(self));
221 if (i%5==4) printf("\n");
222 }
223 printf("\n1000 outputs of genrand_real2()\n");
224 for (i=0; i<1000; i ++) {
225 printf("%10.8f ", genrand_real2());
226 if (i%5==4) printf("\n");
227 }
228 return 0;
229 }
230 */
231
232 /* --------------------------------------------------------
233 This stuff added by Steve Dekorte, 2004 07 04
234 */
235
236 #include <time.h>
237 #include <stdlib.h>
238
RandomGen_new(void)239 RandomGen *RandomGen_new(void)
240 {
241 RandomGen *self = (RandomGen *)io_calloc(1, sizeof(RandomGen));
242 unsigned long t1 = (unsigned long)time(NULL);
243 unsigned long t2 = clock();
244 self->mti = RANDOMGEN_N + 1;
245 init_genrand(self, t1 + t2);
246 self->y2 = 0;
247 return self;
248 }
249
RandomGen_free(RandomGen * self)250 void RandomGen_free(RandomGen *self)
251 {
252 io_free(self);
253 }
254
RandomGen_randomDouble(RandomGen * self)255 double RandomGen_randomDouble(RandomGen *self)
256 {
257 return genrand_real2(self);
258 }
259
RandomGen_randomInt(RandomGen * self)260 int RandomGen_randomInt(RandomGen *self)
261 {
262 return (int)genrand_int32(self);
263 }
264
265 #include <math.h>
266 #ifndef M_PI_2 // some windows'es don't define this
267 # define M_PI_2 1.57079632679489661923
268 #endif
269
RandomGen_gaussian(RandomGen * self,double m,double s)270 double RandomGen_gaussian(RandomGen *self, double m, double s)
271 {
272 // http://www.taygeta.com/random/gaussian.html
273 double x1, x2, w, y1; //, y2;
274
275 do {
276 x1 = 2.0 * RandomGen_randomDouble(self) - 1.0;
277 x2 = 2.0 * RandomGen_randomDouble(self) - 1.0;
278 w = x1 * x1 + x2 * x2;
279 } while ( w >= 1.0 );
280
281 w = sqrt( (-2.0 * log( w ) ) / w );
282 y1 = x1 * w;
283 //y2 = x2 * w;
284
285 // The following code resulted in a lot of nans being returned. The
286 // following code *should* also be slower.
287 /*
288 double x1 = 2.0 * RandomGen_randomDouble(self) - 1.0;
289 double x2 = 2.0 * RandomGen_randomDouble(self) - 1.0;
290 double y1 = sqrt( - 2.0 * log(x1) ) * cos( M_PI_2 * x2 );
291 */
292
293 return ( m + y1 * s );
294 }
295
296 /*
297 double RandomGen_gaussian(RandomGen *self, double m, double s)
298 {
299 // mean m, standard deviation s
300 double x1, x2, w, y1;
301
302 if (self->use_last) // use value from previous call
303 {
304 y1 = self->y2;
305 self->use_last = 0;
306 }
307 else
308 {
309 do {
310 x1 = 2.0 * RandomGen_randomDouble(self) - 1.0;
311 x2 = 2.0 * RandomGen_randomDouble(self) - 1.0;
312 w = x1 * x1 + x2 * x2;
313 } while ( w >= 1.0 );
314
315 w = sqrt( (-2.0 * log( w ) ) / w );
316 y1 = x1 * w;
317 self->y2 = x2 * w;
318 self->use_last = 1;
319 }
320
321 return ( m + y1 * s );
322 }
323 */
324
325