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 moduleocean.math.random.engines.Twister;
22 importocean.meta.types.Qualifiers;
23 importocean.core.TypeConvert: assumeUnique;
24 importInteger = ocean.text.convert.Integer_tango;
25 importocean.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 structTwister56 {
57 privateenum : uint {
58 // Period parameters59 N = 624,
60 M = 397,
61 MATRIX_A = 0x9908b0df, // constant vector a62 UPPER_MASK = 0x80000000, // most significant w-r bits63 LOWER_MASK = 0x7fffffff, // least significant r bits64 }
65 enumintcanCheckpoint=true;
66 enumintcanSeed=true;
67 68 privateuint[N] mt; // the array for the state vector69 privateuintmti=mt.length+1; // mti==mt.length+1 means mt[] is not initialized70 71 /// returns a random uint72 uintnext ()
73 {
74 uinty;
75 staticuint[2] mag01 =[0, MATRIX_A];
76 77 if (mti >= mt.length) {
78 intkk;
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 returny;
103 }
104 /// returns a random byte105 ubytenextB(){
106 returncast(ubyte)(next() & 0xFF);
107 }
108 /// returns a random long109 ulongnextL(){
110 return ((cast(ulong)next)<<32)+cast(ulong)next;
111 }
112 113 /// initializes the generator with a uint as seed114 voidseed (uints)
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 generator123 voidaddEntropy(scopeuintdelegate() r){
124 inti, 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 generator154 voidseed(scopeuintdelegate() r){
155 seed (19650218UL);
156 addEntropy(r);
157 }
158 /// writes the current status in a string159 istringtoString(){
160 char[] res=newchar[7+(N+1)*9];
161 inti=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 returnassumeUnique(res);
176 }
177 /// reads the current status from a string (that should have been trimmed)178 /// returns the number of chars read179 size_tfromString(cstrings){
180 size_ti;
181 verify(s[0..7]=="Twister","unexpected kind, expected Twister");
182 i+=7;
183 verify(s[i]=='_',"no separator _ found");
184 ++i;
185 uintate;
186 mti=cast(uint)Integer.convert(s[i..i+8],16,&ate);
187 verify(ate==8,"unexpected read size");
188 i+=8;
189 foreach (refval;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 returni;
197 }
198 199 }
200