tango.math.random.Twister source code

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[2] mag01 =[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 }