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:
16             KeYeR (D interface) keyer@team0xf.com
17             fawzi (converted to engine)
18 
19 *******************************************************************************/
20 
21 module ocean.math.random.engines.Twister;
22 import ocean.meta.types.Qualifiers;
23 import ocean.core.TypeConvert: assumeUnique;
24 import Integer = ocean.text.convert.Integer_tango;
25 import ocean.core.Verify;
26 
27 /*******************************************************************************
28 
29         Wrapper for the Mersenne twister.
30 
31         The Mersenne twister is a pseudorandom number generator linked to
32         CR developed in 1997 by Makoto Matsumoto and Takuji Nishimura that
33         is based on a matrix linear recurrence over a finite binary field
34         F2. It provides for fast generation of very high quality pseudorandom
35         numbers, having been designed specifically to rectify many of the
36         flaws found in older algorithms.
37 
38         Mersenne Twister has the following desirable properties:
39         ---
40         1. It was designed to have a period of 2^19937 - 1 (the creators
41            of the algorithm proved this property).
42 
43         2. It has a very high order of dimensional equidistribution.
44            This implies that there is negligible serial correlation between
45            successive values in the output sequence.
46 
47         3. It passes numerous tests for statistical randomness, including
48            the stringent Diehard tests.
49 
50         4. It is fast.
51         ---
52 
53 *******************************************************************************/
54 
55 struct Twister
56 {
57     private enum : uint {
58         // Period parameters
59         N          = 624,
60         M          = 397,
61         MATRIX_A   = 0x9908b0df,        // constant vector a
62         UPPER_MASK = 0x80000000,        //  most significant w-r bits
63         LOWER_MASK = 0x7fffffff,        // least significant r bits
64     }
65     enum int canCheckpoint=true;
66     enum int canSeed=true;
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 
71     /// returns a random uint
72     uint next ()
73     {
74         uint y;
75         static uint[2] mag01 =[0, MATRIX_A];
76 
77         if (mti >= mt.length) {
78             int kk;
79 
80             for (kk=0;kk<mt.length-M;kk++)
81             {
82                 y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
83                 mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 1U];
84             }
85             for (;kk<mt.length-1;kk++) {
86                 y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
87                 mt[kk] = mt[kk+(M-mt.length)] ^ (y >> 1) ^ mag01[y & 1U];
88             }
89             y = (mt[mt.length-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
90             mt[mt.length-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 1U];
91 
92             mti = 0;
93         }
94 
95         y = mt[mti++];
96 
97         y ^= (y >> 11);
98         y ^= (y << 7)  &  0x9d2c5680UL;
99         y ^= (y << 15) &  0xefc60000UL;
100         y ^= (y >> 18);
101 
102         return y;
103     }
104     /// returns a random byte
105     ubyte nextB(){
106         return cast(ubyte)(next() & 0xFF);
107     }
108     /// returns a random long
109     ulong nextL(){
110         return ((cast(ulong)next)<<32)+cast(ulong)next;
111     }
112 
113     /// initializes the generator with a uint as seed
114     void seed (uint s)
115     {
116         mt[0]= s & 0xffff_ffffU;  // this is very suspicious, was the previous line incorrectly translated from C???
117         for (mti=1; mti<mt.length; mti++){
118             mt[mti] = cast(uint)(1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
119             mt[mti] &= 0xffff_ffffUL; // this is very suspicious, was the previous line incorrectly translated from C???
120         }
121     }
122     /// adds entropy to the generator
123     void addEntropy(scope uint delegate() r){
124         int i, j, k;
125         i=1;
126         j=0;
127 
128         for (k = mt.length; k; k--)   {
129             mt[i] = cast(uint)((mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))+ r() + j);
130             mt[i] &=  0xffff_ffffUL;  // this is very suspicious, was the previous line incorrectly translated from C???
131             i++;
132             j++;
133 
134             if (i >= mt.length){
135                     mt[0] = mt[mt.length-1];
136                     i=1;
137             }
138         }
139 
140         for (k=mt.length-1; k; k--)     {
141             mt[i] = cast(uint)((mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))- i);
142             mt[i] &=  0xffffffffUL;  // this is very suspicious, was the previous line incorrectly translated from C???
143             i++;
144 
145             if (i>=mt.length){
146                     mt[0] = mt[mt.length-1];
147                     i=1;
148             }
149         }
150         mt[0] |=  0x80000000UL;
151         mti=0;
152     }
153     /// seeds the generator
154     void seed(scope uint delegate() r){
155         seed (19650218UL);
156         addEntropy(r);
157     }
158     /// writes the current status in a string
159     istring toString(){
160         char[] res=new char[7+(N+1)*9];
161         int i=0;
162         res[i..i+7]="Twister";
163         i+=7;
164         res[i]='_';
165         ++i;
166         Integer.format(res[i..i+8],mti,"x8");
167         i+=8;
168         foreach (val;mt){
169             res[i]='_';
170             ++i;
171             Integer.format(res[i..i+8],val,"x8");
172             i+=8;
173         }
174         verify(i==res.length,"unexpected size");
175         return assumeUnique(res);
176     }
177     /// reads the current status from a string (that should have been trimmed)
178     /// returns the number of chars read
179     size_t fromString(cstring s){
180         size_t i;
181         verify(s[0..7]=="Twister","unexpected kind, expected Twister");
182         i+=7;
183         verify(s[i]=='_',"no separator _ found");
184         ++i;
185         uint ate;
186         mti=cast(uint)Integer.convert(s[i..i+8],16,&ate);
187         verify(ate==8,"unexpected read size");
188         i+=8;
189         foreach (ref val;mt){
190             verify(s[i]=='_',"no separator _ found");
191             ++i;
192             val=cast(uint)Integer.convert(s[i..i+8],16,&ate);
193             verify(ate==8,"unexpected read size");
194             i+=8;
195         }
196         return i;
197     }
198 
199 }
200