1 /*******************************************************************************
2
3 copyright: Copyright (C) 1997--2004, Makoto Matsumoto,
4 Takuji Nishimura, and Eric Landry; All rights reserved
5
6
7 license: BSD style: $(LICENSE)
8
9 version: Jan 2008: Initial release
10
11 author: KeYeR (D interface) keyer@team0xf.com
12
13 *******************************************************************************/
14
15 module tango.math.random.Twister;
16
17
18 version (Win32)
19 private extern(Windows) int QueryPerformanceCounter (ulong *);
20
21 version (Posix)
22 {
23 private import tango.stdc.posix.sys.time;
24 }
25
26 /*******************************************************************************
27
28 Wrapper for the Mersenne twister.
29
30 The Mersenne twister is a pseudorandom number generator linked to
31 CR developed in 1997 by Makoto Matsumoto and Takuji Nishimura that
32 is based on a matrix linear recurrence over a finite binary field
33 F2. It provides for fast generation of very high quality pseudorandom
34 numbers, having been designed specifically to rectify many of the
35 flaws found in older algorithms.
36
37 Mersenne Twister has the following desirable properties:
38 ---
39 1. It was designed to have a period of 2^19937 - 1 (the creators
40 of the algorithm proved this property).
41
42 2. It has a very high order of dimensional equidistribution.
43 This implies that there is negligible serial correlation between
44 successive values in the output sequence.
45
46 3. It passes numerous tests for statistical randomness, including
47 the stringent Diehard tests.
48
49 4. It is fast.
50 ---
51
52 *******************************************************************************/
53
54 struct Twister
55 {
56 public alias natural toInt;
57 public alias fraction toReal;
58
59 private enum : uint // Period parameters
60 {
61 N = 624,
62 M = 397,
63 MATRIX_A = 0x9908b0df, // constant vector a
64 UPPER_MASK = 0x80000000, // most significant w-r bits
65 LOWER_MASK = 0x7fffffff, // least significant r bits
66 }
67
68 private uint[N] mt; // the array for the state vector
69 private uint mti=mt.length+1; // mti==mt.length+1 means mt[] is not initialized
70 private uint vLastRand; // The most recent random uint returned.
71
72 /**********************************************************************
73
74 A global, shared instance, seeded via startup time
75
76 **********************************************************************/
77
78 public static __gshared Twister instance;
79
80 shared static this ()
81 {
82 instance.seed();
83 }
84
85 /**********************************************************************
86
87 Creates and seeds a new generator with the current time
88
89 **********************************************************************/
90
91 static Twister opCall ()
92 {
93 Twister rand;
94 rand.seed();
95 return rand;
96 }
97
98 /**********************************************************************
99
100 Returns X such that 0 <= X < max
101
102 **********************************************************************/
103
104 uint natural (uint max)
105 {
106 return natural() % max;
107 }
108
109 /**********************************************************************
110
111 Returns X such that min <= X < max
112
113 **********************************************************************/
114
115 uint natural (uint min, uint max)
116 {
117 return (natural() % (max-min)) + min;
118 }
119
120 /**********************************************************************
121
122 Returns X such that 0 <= X <= uint.max
123
124 **********************************************************************/
125
126 uint natural (bool pAddEntropy = false)
127 {
128 uint y;
129 static uint mag01[2] =[0, MATRIX_A];
130
131 if (mti >= mt.length) {
132 int kk;
133
134 if (pAddEntropy || mti > mt.length)
135 {
136 seed (5489U, pAddEntropy);
137 }
138
139 for (kk=0;kk<mt.length-M;kk++)
140 {
141 y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
142 mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 1U];
143 }
144 for (;kk<mt.length-1;kk++) {
145 y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
146 mt[kk] = mt[kk+(M-mt.length)] ^ (y >> 1) ^ mag01[y & 1U];
147 }
148 y = (mt[mt.length-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
149 mt[mt.length-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 1U];
150
151 mti = 0;
152 }
153
154 y = mt[mti++];
155
156 y ^= (y >> 11);
157 y ^= (y << 7) & 0x9d2c5680U;
158 y ^= (y << 15) & 0xefc60000U;
159 y ^= (y >> 18);
160
161 vLastRand = y;
162 return y;
163 }
164
165 /**********************************************************************
166
167 generates a random number on [0,1] interval
168
169 **********************************************************************/
170
171 double inclusive ()
172 {
173 return natural()*(1.0/cast(double)uint.max);
174 }
175
176 /**********************************************************************
177
178 generates a random number on (0,1) interval
179
180 **********************************************************************/
181
182 double exclusive ()
183 {
184 return ((cast(double)natural()) + 0.5)*(1.0/(cast(double)uint.max+1.0));
185 }
186
187 /**********************************************************************
188
189 generates a random number on [0,1) interval
190
191 **********************************************************************/
192
193 double fraction ()
194 {
195 return natural()*(1.0/(cast(double)uint.max+1.0));
196 }
197
198 /**********************************************************************
199
200 generates a random number on [0,1) with 53-bit resolution
201
202 **********************************************************************/
203
204 double fractionEx ()
205 {
206 uint a = natural() >> 5, b = natural() >> 6;
207 return(a * 67108864.0 + b) * (1.0 / 9007199254740992.0);
208 }
209
210 /**********************************************************************
211
212 Seed the generator with current time
213
214 **********************************************************************/
215
216 void seed ()
217 {
218 ulong s;
219
220 version (Posix)
221 {
222 timeval tv;
223
224 gettimeofday (&tv, null);
225 s = tv.tv_usec;
226 }
227 version (Win32)
228 QueryPerformanceCounter (&s);
229
230 return seed (cast(uint) s);
231 }
232
233 /**********************************************************************
234
235 initializes the generator with a seed
236
237 **********************************************************************/
238
239 void seed (uint s, bool pAddEntropy = false)
240 {
241 mt[0]= (s + (pAddEntropy ? vLastRand + cast(uint) &this : 0)) & 0xffffffffU;
242 for (mti=1; mti<mt.length; mti++){
243 mt[mti] = (1812433253U * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
244 mt[mti] &= 0xffffffffU;
245 }
246 }
247
248 /**********************************************************************
249
250 initialize by an array with array-length init_key is
251 the array for initializing keys
252
253 **********************************************************************/
254
255 void init (const(uint[]) init_key, bool pAddEntropy = false)
256 {
257 size_t k;
258 int i, j;
259 i=1;
260 j=0;
261
262 seed (19650218U, pAddEntropy);
263 for (k = (mt.length > init_key.length ? mt.length : init_key.length); k; k--) {
264 mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525U))+ init_key[j] + j;
265 mt[i] &= 0xffffffffU;
266 i++;
267 j++;
268
269 if (i >= mt.length){
270 mt[0] = mt[mt.length-1];
271 i=1;
272 }
273
274 if (j >= init_key.length){
275 j=0;
276 }
277 }
278
279 for (k=mt.length-1; k; k--) {
280 mt[i] = (mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941U))- i;
281 mt[i] &= 0xffffffffU;
282 i++;
283
284 if (i>=mt.length){
285 mt[0] = mt[mt.length-1];
286 i=1;
287 }
288 }
289 mt[0] |= 0x80000000U;
290 mti=0;
291 }
292 }
293
294
295 /******************************************************************************
296
297
298 ******************************************************************************/
299
300 debug (Twister)
301 {
302 import tango.io.Stdout;
303 import tango.time.StopWatch;
304
305 void main()
306 {
307 auto dbl = Twister();
308 auto count = 100_000_000;
309 StopWatch w;
310
311 w.start;
312 double v1;
313 for (int i=count; --i;)
314 v1 = dbl.fraction;
315 Stdout.formatln ("{} fraction, {}/s, {:f10}", count, count/w.stop, v1);
316
317 w.start;
318 for (int i=count; --i;)
319 v1 = dbl.fractionEx;
320 Stdout.formatln ("{} fractionEx, {}/s, {:f10}", count, count/w.stop, v1);
321
322 for (int i=count; --i;)
323 {
324 auto v = dbl.fraction;
325 if (v <= 0.0 || v >= 1.0)
326 {
327 Stdout.formatln ("fraction {:f10}", v);
328 break;
329 }
330 v = dbl.fractionEx;
331 if (v <= 0.0 || v >= 1.0)
332 {
333 Stdout.formatln ("fractionEx {:f10}", v);
334 break;
335 }
336 }
337 }
338 }