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.transition; 23 import Integer = ocean.text.convert.Integer_tango; 24 import ocean.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 struct Twister 55 { 56 private enum : uint { 57 // Period parameters 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 enum int canCheckpoint=true; 65 enum int canSeed=true; 66 67 private uint[N] mt; // the array for the state vector 68 private uint mti=mt.length+1; // mti==mt.length+1 means mt[] is not initialized 69 70 /// returns a random uint 71 uint next () 72 { 73 uint y; 74 static uint[2] mag01 =[0, MATRIX_A]; 75 76 if (mti >= mt.length) { 77 int kk; 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 return y; 102 } 103 /// returns a random byte 104 ubyte nextB(){ 105 return cast(ubyte)(next() & 0xFF); 106 } 107 /// returns a random long 108 ulong nextL(){ 109 return ((cast(ulong)next)<<32)+cast(ulong)next; 110 } 111 112 /// initializes the generator with a uint as seed 113 void seed (uint s) 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 generator 122 void addEntropy(scope uint delegate() r){ 123 int i, 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 generator 153 void seed(scope uint delegate() r){ 154 seed (19650218UL); 155 addEntropy(r); 156 } 157 /// writes the current status in a string 158 istring toString(){ 159 char[] res=new char[7+(N+1)*9]; 160 int i=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 return assumeUnique(res); 175 } 176 /// reads the current status from a string (that should have been trimmed) 177 /// returns the number of chars read 178 size_t fromString(cstring s){ 179 size_t i; 180 verify(s[0..7]=="Twister","unexpected kind, expected Twister"); 181 i+=7; 182 verify(s[i]=='_',"no separator _ found"); 183 ++i; 184 uint ate; 185 mti=cast(uint)Integer.convert(s[i..i+8],16,&ate); 186 verify(ate==8,"unexpected read size"); 187 i+=8; 188 foreach (ref val;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 return i; 196 } 197 198 } 199