tango.math.random.engines.KISS 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.KISS;
8 private import Integer = tango.text.convert.Integer;
9 
10 /+ Kiss99 random number generator, by Marisaglia
11 + a simple RNG that passes all statistical tests
12 + This is the engine, *never* use it directly, always use it though a RandomG class
13 +/
14 struct Kiss99{
15     private uint kiss_x = 123456789;
16     private uint kiss_y = 362436000;
17     private uint kiss_z = 521288629;
18     private uint kiss_c = 7654321;
19     private uint nBytes = 0;
20     private uint restB  = 0;
21     
22     enum int canCheckpoint=true;
23     enum int canSeed=true;
24     
25     void skip(uint n){
26         for (int i=n;i!=n;--i){
27             next();
28         }
29     }
30     ubyte nextB(){
31         if (nBytes>0) {
32             ubyte res=cast(ubyte)(restB & 0xFF);
33             restB >>= 8;
34             --nBytes;
35             return res;
36         } else {
37             restB=next();
38             ubyte res=cast(ubyte)(restB & 0xFF);
39             restB >>= 8;
40             nBytes=3;
41             return res;
42         }
43     }
44     uint next(){
45         enum ulong a = 698769069UL;
46         ulong t;
47         kiss_x = 69069*kiss_x+12345;
48         kiss_y ^= (kiss_y<<13); kiss_y ^= (kiss_y>>17); kiss_y ^= (kiss_y<<5);
49         t = a*kiss_z+kiss_c; kiss_c = cast(uint)(t>>32);
50         kiss_z=cast(uint)t;
51         return kiss_x+kiss_y+kiss_z;
52     }
53     ulong nextL(){
54         return ((cast(ulong)next())<<32)+cast(ulong)next();
55     }
56     
57     void seed(scope uint delegate() r){
58         kiss_x = r();
59         for (int i=0;i<100;++i){
60             kiss_y=r();
61             if (kiss_y!=0) break;
62         }
63         if (kiss_y==0) kiss_y=362436000;
64         kiss_z=r();
65         /* Don’t really need to seed c as well (is reset after a next),
66            but doing it allows to completely restore a given internal state */
67         kiss_c = r() % 698769069; /* Should be less than 698769069 */
68         nBytes = 0;
69         restB=0;
70     }
71     /// writes the current status in a string
72     immutable(char)[] toString(){
73         char[] res=new char[6+6*9];
74         int i=0;
75         res[i..i+6]="KISS99";
76         i+=6;
77         foreach (val;[kiss_x,kiss_y,kiss_z,kiss_c,nBytes,restB]){
78             res[i]='_';
79             ++i;
80             Integer.format(res[i..i+8],val,cast(char[])"x8");
81             i+=8;
82         }
83         assert(i==res.length,"unexpected size");
84         return cast(immutable(char)[])res;
85     }
86     /// reads the current status from a string (that should have been trimmed)
87     /// returns the number of chars read
88     size_t fromString(const(char[]) s){
89         size_t i=0;
90         assert(s[i..i+4]=="KISS","unexpected kind, expected KISS");
91         assert(s[i+4..i+7]=="99_","unexpected version, expected 99");
92         i+=6;
93         foreach (val;[&kiss_x,&kiss_y,&kiss_z,&kiss_c,&nBytes,&restB]){
94             assert(s[i]=='_',"no separator _ found");
95             ++i;
96             size_t ate;
97             *val=cast(uint)Integer.convert(s[i..i+8],16,&ate);
98             assert(ate==8,"unexpected read size");
99             i+=8;
100         }
101         return i;
102     }
103 }