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