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