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 module ocean.math.random.engines.CMWC; 18 import ocean.meta.types.Qualifiers; 19 import ocean.core.TypeConvert: assumeUnique; 20 import Integer=ocean.text.convert.Integer_tango; 21 import ocean.core.Verify; 22 23 /+ CMWC a random number generator, 24 + Marisaglia, Journal of Modern Applied Statistical Methods (2003), vol.2,No.1,p 2-13 25 + a simple and fast RNG that passes all statistical tests, has a large seed, and is very fast 26 + By default ComplimentaryMultiplyWithCarry with r=1024, a=987769338, b=2^32-1, period a*b^r>10^9873 27 + This is the engine, *never* use it directly, always use it though a RandomG class 28 +/ 29 struct CMWC(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 uint nBytes = 0; 34 uint restB = 0; 35 36 enum int canCheckpoint=true; 37 enum int canSeed=true; 38 39 void skip(uint n){ 40 for (int i=n;i!=n;--i){ 41 next; 42 } 43 } 44 ubyte nextB(){ 45 if (nBytes>0) { 46 ubyte res=cast(ubyte)(restB & 0xFF); 47 restB >>= 8; 48 --nBytes; 49 return res; 50 } else { 51 restB=next; 52 ubyte res=cast(ubyte)(restB & 0xFF); 53 restB >>= 8; 54 nBytes=3; 55 return res; 56 } 57 } 58 uint next(){ 59 enum uint rMask=cmwc_r-1u; 60 static assert((rMask&cmwc_r)==0,"cmwc_r is supposed to be a power of 2"); // put a more stringent test? 61 enum uint m=0xFFFF_FFFE; 62 cmwc_i=(cmwc_i+1)&rMask; 63 ulong t=cmwc_a*cmwc_q[cmwc_i]+cmwc_c; 64 cmwc_c=cast(uint)(t>>32); 65 uint x=cast(uint)t+cmwc_c; 66 if (x<cmwc_c) { 67 ++x; ++cmwc_c; 68 } 69 return (cmwc_q[cmwc_i]=m-x); 70 } 71 72 ulong nextL(){ 73 return ((cast(ulong)next)<<32)+cast(ulong)next; 74 } 75 76 void seed(scope uint delegate() rSeed){ 77 cmwc_i=cmwc_r-1u; // randomize also this? 78 for (int ii=0;ii<10;++ii){ 79 for (uint i=0;i<cmwc_r;++i){ 80 cmwc_q[i]=rSeed(); 81 } 82 uint nV=rSeed(); 83 uint to=(uint.max/cmwc_a)*cmwc_a; 84 for (int i=0;i<10;++i){ 85 if (nV<to) break; 86 nV=rSeed(); 87 } 88 cmwc_c=cast(uint)(nV%cmwc_a); 89 nBytes = 0; 90 restB=0; 91 if (cmwc_c==0){ 92 for (uint i=0;i<cmwc_r;++i) 93 if (cmwc_q[i]!=0) return; 94 } else if (cmwc_c==cmwc_a-1){ 95 for (uint i=0;i<cmwc_r;++i) 96 if (cmwc_q[i]!=0xFFFF_FFFF) return; 97 } else return; 98 } 99 cmwc_c=1; 100 } 101 /// writes the current status in a string 102 istring toString(){ 103 char[] res=new char[4+16+(cmwc_r+5)*9]; 104 int i=0; 105 res[i..i+4]="CMWC"; 106 i+=4; 107 Integer.format(res[i..i+16],cmwc_a,"x16"); 108 i+=16; 109 res[i]='_'; 110 ++i; 111 Integer.format(res[i..i+8],cmwc_r,"x8"); 112 i+=8; 113 foreach (val;cmwc_q){ 114 res[i]='_'; 115 ++i; 116 Integer.format(res[i..i+8],val,"x8"); 117 i+=8; 118 } 119 foreach (val;[cmwc_i,cmwc_c,nBytes,restB]){ 120 res[i]='_'; 121 ++i; 122 Integer.format(res[i..i+8],val,"x8"); 123 i+=8; 124 } 125 verify(i==res.length,"unexpected size"); 126 return assumeUnique(res); 127 } 128 /// reads the current status from a string (that should have been trimmed) 129 /// returns the number of chars read 130 size_t fromString(cstring s){ 131 char[16] tmpC; 132 size_t i=0; 133 verify(s[i..i+4]=="CMWC","unexpected kind, expected CMWC"); 134 i+=4; 135 verify(s[i..i+16]==Integer.format(tmpC,cmwc_a,"x16"),"unexpected cmwc_a"); 136 i+=16; 137 verify(s[i]=='_',"unexpected format, expected _"); 138 i++; 139 verify(s[i..i+8]==Integer.format(tmpC,cmwc_r,"x8"),"unexpected cmwc_r"); 140 i+=8; 141 foreach (ref val;cmwc_q){ 142 verify(s[i]=='_',"no separator _ found"); 143 ++i; 144 uint ate; 145 val=cast(uint)Integer.convert(s[i..i+8],16,&ate); 146 verify(ate==8,"unexpected read size"); 147 i+=8; 148 } 149 foreach (val;[&cmwc_i,&cmwc_c,&nBytes,&restB]){ 150 verify(s[i]=='_',"no separator _ found"); 151 ++i; 152 uint ate; 153 *val=cast(uint)Integer.convert(s[i..i+8],16,&ate); 154 verify(ate==8,"unexpected read size"); 155 i+=8; 156 } 157 return i; 158 } 159 } 160 161 /// some variations, the first has a period of ~10^39461, the first number (r) is basically the size of the seed 162 /// (and all bit patterns of that size are guarenteed to crop up in the period), the period is 2^(32*r)*a 163 alias CMWC!(4096U,18782UL) CMWC_4096_1; 164 alias CMWC!(2048U,1030770UL) CMWC_2048_1; 165 alias CMWC!(2048U,1047570UL) CMWC_2048_2; 166 alias CMWC!(1024U,5555698UL) CMWC_1024_1; 167 alias CMWC!(1024U,987769338UL) CMWC_1024_2; 168 alias CMWC!(512U,123462658UL) CMWC_512_1; 169 alias CMWC!(512U,123484214UL) CMWC_512_2; 170 alias CMWC!(256U,987662290UL) CMWC_256_1; 171 alias CMWC!(256U,987665442UL) CMWC_256_2; 172 alias CMWC!(128U,987688302UL) CMWC_128_1; 173 alias CMWC!(128U,987689614UL) CMWC_128_2; 174 alias CMWC!(64U ,987651206UL) CMWC_64_1; 175 alias CMWC!(64U ,987657110UL) CMWC_64_2; 176 alias CMWC!(32U ,987655670UL) CMWC_32_1; 177 alias CMWC!(32U ,987655878UL) CMWC_32_2; 178 alias CMWC!(16U ,987651178UL) CMWC_16_1; 179 alias CMWC!(16U ,987651182UL) CMWC_16_2; 180 alias CMWC!(8U ,987651386UL) CMWC_8_1; 181 alias CMWC!(8U ,987651670UL) CMWC_8_2; 182 alias CMWC!(4U ,987654366UL) CMWC_4_1; 183 alias CMWC!(4U ,987654978UL) CMWC_4_2; 184 alias CMWC_1024_2 CMWC_default;