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.transition;
23 importInteger = ocean.text.convert.Integer_tango;
24 importocean.core.Verify;
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 structTwister55 {
56 privateenum : uint {
57 // Period parameters58 N = 624,
59 M = 397,
60 MATRIX_A = 0x9908b0df, // constant vector a61 UPPER_MASK = 0x80000000, // most significant w-r bits62 LOWER_MASK = 0x7fffffff, // least significant r bits63 }
64 enumintcanCheckpoint=true;
65 enumintcanSeed=true;
66 67 privateuint[N] mt; // the array for the state vector68 privateuintmti=mt.length+1; // mti==mt.length+1 means mt[] is not initialized69 70 /// returns a random uint71 uintnext ()
72 {
73 uinty;
74 staticuint[2] mag01 =[0, MATRIX_A];
75 76 if (mti >= mt.length) {
77 intkk;
78 79 for (kk=0;kk<mt.length-M;kk++)
80 {
81 y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
82 mt[kk] = mt[kk+M] ^ (y >> 1) ^ mag01[y & 1U];
83 }
84 for (;kk<mt.length-1;kk++) {
85 y = (mt[kk]&UPPER_MASK)|(mt[kk+1]&LOWER_MASK);
86 mt[kk] = mt[kk+(M-mt.length)] ^ (y >> 1) ^ mag01[y & 1U];
87 }
88 y = (mt[mt.length-1]&UPPER_MASK)|(mt[0]&LOWER_MASK);
89 mt[mt.length-1] = mt[M-1] ^ (y >> 1) ^ mag01[y & 1U];
90 91 mti = 0;
92 }
93 94 y = mt[mti++];
95 96 y ^= (y >> 11);
97 y ^= (y << 7) & 0x9d2c5680UL;
98 y ^= (y << 15) & 0xefc60000UL;
99 y ^= (y >> 18);
100 101 returny;
102 }
103 /// returns a random byte104 ubytenextB(){
105 returncast(ubyte)(next() & 0xFF);
106 }
107 /// returns a random long108 ulongnextL(){
109 return ((cast(ulong)next)<<32)+cast(ulong)next;
110 }
111 112 /// initializes the generator with a uint as seed113 voidseed (uints)
114 {
115 mt[0]= s & 0xffff_ffffU; // this is very suspicious, was the previous line incorrectly translated from C???116 for (mti=1; mti<mt.length; mti++){
117 mt[mti] = cast(uint)(1812433253UL * (mt[mti-1] ^ (mt[mti-1] >> 30)) + mti);
118 mt[mti] &= 0xffff_ffffUL; // this is very suspicious, was the previous line incorrectly translated from C???119 }
120 }
121 /// adds entropy to the generator122 voidaddEntropy(scopeuintdelegate() r){
123 inti, j, k;
124 i=1;
125 j=0;
126 127 for (k = mt.length; k; k--) {
128 mt[i] = cast(uint)((mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1664525UL))+ r() + j);
129 mt[i] &= 0xffff_ffffUL; // this is very suspicious, was the previous line incorrectly translated from C???130 i++;
131 j++;
132 133 if (i >= mt.length){
134 mt[0] = mt[mt.length-1];
135 i=1;
136 }
137 }
138 139 for (k=mt.length-1; k; k--) {
140 mt[i] = cast(uint)((mt[i] ^ ((mt[i-1] ^ (mt[i-1] >> 30)) * 1566083941UL))- i);
141 mt[i] &= 0xffffffffUL; // this is very suspicious, was the previous line incorrectly translated from C???142 i++;
143 144 if (i>=mt.length){
145 mt[0] = mt[mt.length-1];
146 i=1;
147 }
148 }
149 mt[0] |= 0x80000000UL;
150 mti=0;
151 }
152 /// seeds the generator153 voidseed(scopeuintdelegate() r){
154 seed (19650218UL);
155 addEntropy(r);
156 }
157 /// writes the current status in a string158 istringtoString(){
159 char[] res=newchar[7+(N+1)*9];
160 inti=0;
161 res[i..i+7]="Twister";
162 i+=7;
163 res[i]='_';
164 ++i;
165 Integer.format(res[i..i+8],mti,"x8");
166 i+=8;
167 foreach (val;mt){
168 res[i]='_';
169 ++i;
170 Integer.format(res[i..i+8],val,"x8");
171 i+=8;
172 }
173 verify(i==res.length,"unexpected size");
174 returnassumeUnique(res);
175 }
176 /// reads the current status from a string (that should have been trimmed)177 /// returns the number of chars read178 size_tfromString(cstrings){
179 size_ti;
180 verify(s[0..7]=="Twister","unexpected kind, expected Twister");
181 i+=7;
182 verify(s[i]=='_',"no separator _ found");
183 ++i;
184 uintate;
185 mti=cast(uint)Integer.convert(s[i..i+8],16,&ate);
186 verify(ate==8,"unexpected read size");
187 i+=8;
188 foreach (refval;mt){
189 verify(s[i]=='_',"no separator _ found");
190 ++i;
191 val=cast(uint)Integer.convert(s[i..i+8],16,&ate);
192 verify(ate==8,"unexpected read size");
193 i+=8;
194 }
195 returni;
196 }
197 198 }
199