ocean.math.random.Twister source code

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 }