tango.math.random.engines.KissCmwc source code

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