1 /******************************************************************************* 2 3 Copyright: 4 Copyright (c) 2008. Fawzi Mohamed 5 Some parts copyright (c) 2009-2016 dunnhumby Germany GmbH. 6 All rights reserved. 7 8 License: 9 Tango Dual License: 3-Clause BSD License / Academic Free License v3.0. 10 See LICENSE_TANGO.txt for details. 11 12 Version: Initial release: July 2008 13 14 Authors: Fawzi Mohamed 15 16 *******************************************************************************/ 17 18 module ocean.math.random.engines.KissCmwc; 19 import ocean.meta.types.Qualifiers; 20 import ocean.core.TypeConvert: assumeUnique; 21 import Integer = ocean.text.convert.Integer_tango; 22 import ocean.core.Verify; 23 24 /+ CMWC and KISS random number generators combined, for extra security wrt. plain CMWC and 25 + Marisaglia, Journal of Modern Applied Statistical Methods (2003), vol.2,No.1,p 2-13 26 + a simple safe and fast RNG that passes all statistical tests, has a large seed, and is reasonably fast 27 + This is the engine, *never* use it directly, always use it though a RandomG class 28 +/ 29 struct KissCmwc(uint cmwc_r=1024U,ulong cmwc_a=987769338UL){ 30 uint[cmwc_r] cmwc_q=void; 31 uint cmwc_i=cmwc_r-1u; 32 uint cmwc_c=123; 33 private uint kiss_x = 123456789; 34 private uint kiss_y = 362436000; 35 private uint kiss_z = 521288629; 36 private uint kiss_c = 7654321; 37 uint nBytes = 0; 38 uint restB = 0; 39 40 enum int canCheckpoint=true; 41 enum int canSeed=true; 42 43 void skip(uint n){ 44 for (int i=n;i!=n;--i){ 45 next; 46 } 47 } 48 ubyte nextB(){ 49 if (nBytes>0) { 50 ubyte res=cast(ubyte)(restB & 0xFF); 51 restB >>= 8; 52 --nBytes; 53 return res; 54 } else { 55 restB=next; 56 ubyte res=cast(ubyte)(restB & 0xFF); 57 restB >>= 8; 58 nBytes=3; 59 return res; 60 } 61 } 62 uint next(){ 63 enum uint rMask=cmwc_r-1u; 64 static assert((rMask&cmwc_r)==0,"cmwc_r is supposed to be a power of 2"); // put a more stringent test? 65 enum uint m=0xFFFF_FFFE; 66 cmwc_i=(cmwc_i+1)&rMask; 67 ulong t=cmwc_a*cmwc_q[cmwc_i]+cmwc_c; 68 cmwc_c=cast(uint)(t>>32); 69 uint x=cast(uint)t+cmwc_c; 70 if (x<cmwc_c) { 71 ++x; ++cmwc_c; 72 } 73 enum ulong a = 698769069UL; 74 ulong k_t; 75 kiss_x = 69069*kiss_x+12345; 76 kiss_y ^= (kiss_y<<13); kiss_y ^= (kiss_y>>17); kiss_y ^= (kiss_y<<5); 77 k_t = a*kiss_z+kiss_c; kiss_c = cast(uint)(k_t>>32); 78 kiss_z=cast(uint)k_t; 79 return (cmwc_q[cmwc_i]=m-x)+kiss_x+kiss_y+kiss_z; // xor to avoid overflow? 80 } 81 82 ulong nextL(){ 83 return ((cast(ulong)next)<<32)+cast(ulong)next; 84 } 85 86 void seed(scope uint delegate() rSeed){ 87 kiss_x = rSeed(); 88 for (int i=0;i<100;++i){ 89 kiss_y=rSeed(); 90 if (kiss_y!=0) break; 91 } 92 if (kiss_y==0) kiss_y=362436000; 93 kiss_z=rSeed(); 94 /* Don’t really need to seed c as well (is reset after a next), 95 but doing it allows to completely restore a given internal state */ 96 kiss_c = rSeed() % 698769069; /* Should be less than 698769069 */ 97 98 cmwc_i=cmwc_r-1u; // randomize also this? 99 for (int ii=0;ii<10;++ii){ 100 for (uint i=0;i<cmwc_r;++i){ 101 cmwc_q[i]=rSeed(); 102 } 103 uint nV=rSeed(); 104 uint to=(uint.max/cmwc_a)*cmwc_a; 105 for (int i=0;i<10;++i){ 106 if (nV<to) break; 107 nV=rSeed(); 108 } 109 cmwc_c=cast(uint)(nV%cmwc_a); 110 nBytes = 0; 111 restB=0; 112 if (cmwc_c==0){ 113 for (uint i=0;i<cmwc_r;++i) 114 if (cmwc_q[i]!=0) return; 115 } else if (cmwc_c==cmwc_a-1){ 116 for (uint i=0;i<cmwc_r;++i) 117 if (cmwc_q[i]!=0xFFFF_FFFF) return; 118 } else return; 119 } 120 cmwc_c=1; 121 kiss_x = rSeed(); 122 for (int i=0;i<100;++i){ 123 kiss_y=rSeed(); 124 if (kiss_y!=0) break; 125 } 126 if (kiss_y==0) kiss_y=362436000; 127 kiss_z=rSeed(); 128 /* Don’t really need to seed c as well (is reset after a next), 129 but doing it allows to completely restore a given internal state */ 130 kiss_c = rSeed() % 698769069; /* Should be less than 698769069 */ 131 } 132 /// writes the current status in a string 133 istring toString(){ 134 char[] res=new char[11+16+(cmwc_r+9)*9]; 135 int i=0; 136 res[i..i+11]="CMWC+KISS99"; 137 i+=11; 138 Integer.format(res[i..i+16],cmwc_a,"x16"); 139 i+=16; 140 res[i]='_'; 141 ++i; 142 Integer.format(res[i..i+8],cmwc_r,"x8"); 143 i+=8; 144 foreach (val;cmwc_q){ 145 res[i]='_'; 146 ++i; 147 Integer.format(res[i..i+8],val,"x8"); 148 i+=8; 149 } 150 foreach (val;[cmwc_i,cmwc_c,nBytes,restB,kiss_x,kiss_y,kiss_z,kiss_c]){ 151 res[i]='_'; 152 ++i; 153 Integer.format(res[i..i+8],val,"x8"); 154 i+=8; 155 } 156 verify(i==res.length,"unexpected size"); 157 return assumeUnique(res); 158 } 159 /// reads the current status from a string (that should have been trimmed) 160 /// returns the number of chars read 161 size_t fromString(cstring s){ 162 char[16] tmpC; 163 size_t i=0; 164 verify(s[i..i+11]=="CMWC+KISS99","unexpected kind, expected CMWC+KISS99"); 165 i+=11; 166 verify(s[i..i+16]==Integer.format(tmpC,cmwc_a,"x16"),"unexpected cmwc_a"); 167 i+=16; 168 verify(s[i]=='_',"unexpected format, expected _"); 169 i++; 170 verify(s[i..i+8]==Integer.format(tmpC,cmwc_r,"x8"),"unexpected cmwc_r"); 171 i+=8; 172 foreach (ref val;cmwc_q){ 173 verify(s[i]=='_',"no separator _ found"); 174 ++i; 175 uint ate; 176 val=cast(uint)Integer.convert(s[i..i+8],16,&ate); 177 verify(ate==8,"unexpected read size"); 178 i+=8; 179 } 180 foreach (val;[&cmwc_i,&cmwc_c,&nBytes,&restB,&kiss_x,&kiss_y,&kiss_z,&kiss_c]){ 181 verify(s[i]=='_',"no separator _ found"); 182 ++i; 183 uint ate; 184 *val=cast(uint)Integer.convert(s[i..i+8],16,&ate); 185 verify(ate==8,"unexpected read size"); 186 i+=8; 187 } 188 return i; 189 } 190 } 191 192 /// some variations of the CMWC part, the first has a period of ~10^39461 193 /// the first number (r) is basically the size of the seed and storage (and all bit patterns 194 /// of that size are guarenteed to crop up in the period), the period is (2^32-1)^r*a 195 alias KissCmwc!(4096U,18782UL) KissCmwc_4096_1; 196 alias KissCmwc!(2048U,1030770UL) KissCmwc_2048_1; 197 alias KissCmwc!(2048U,1047570UL) KissCmwc_2048_2; 198 alias KissCmwc!(1024U,5555698UL) KissCmwc_1024_1; 199 alias KissCmwc!(1024U,987769338UL) KissCmwc_1024_2; 200 alias KissCmwc!(512U,123462658UL) KissCmwc_512_1; 201 alias KissCmwc!(512U,123484214UL) KissCmwc_512_2; 202 alias KissCmwc!(256U,987662290UL) KissCmwc_256_1; 203 alias KissCmwc!(256U,987665442UL) KissCmwc_256_2; 204 alias KissCmwc!(128U,987688302UL) KissCmwc_128_1; 205 alias KissCmwc!(128U,987689614UL) KissCmwc_128_2; 206 alias KissCmwc!(64U ,987651206UL) KissCmwc_64_1; 207 alias KissCmwc!(64U ,987657110UL) KissCmwc_64_2; 208 alias KissCmwc!(32U ,987655670UL) KissCmwc_32_1; 209 alias KissCmwc!(32U ,987655878UL) KissCmwc_32_2; 210 alias KissCmwc!(16U ,987651178UL) KissCmwc_16_1; 211 alias KissCmwc!(16U ,987651182UL) KissCmwc_16_2; 212 alias KissCmwc!(8U ,987651386UL) KissCmwc_8_1; 213 alias KissCmwc!(8U ,987651670UL) KissCmwc_8_2; 214 alias KissCmwc!(4U ,987654366UL) KissCmwc_4_1; 215 alias KissCmwc!(4U ,987654978UL) KissCmwc_4_2; 216 alias KissCmwc_1024_2 KissCmwc_default;