|
132
|
1 /*******************************************************************************
|
|
|
2
|
|
|
3 copyright: Copyright (c) 2004. All rights reserved
|
|
|
4
|
|
|
5 license: BSD style: $(LICENSE)
|
|
|
6
|
|
|
7 version: Initial release: April 2004
|
|
|
8
|
|
|
9 author: Various
|
|
|
10
|
|
|
11 *******************************************************************************/
|
|
|
12
|
|
|
13 module tango.math.Random;
|
|
|
14
|
|
|
15
|
|
|
16 version (Win32)
|
|
|
17 private extern(Windows) int QueryPerformanceCounter (ulong *);
|
|
|
18
|
|
|
19 version (Posix)
|
|
|
20 {
|
|
|
21 private import tango.stdc.posix.sys.time;
|
|
|
22 }
|
|
|
23
|
|
|
24
|
|
|
25 /******************************************************************************
|
|
|
26
|
|
|
27 KISS (via George Marsaglia & Paul Hsieh)
|
|
|
28
|
|
|
29 the idea is to use simple, fast, individually promising
|
|
|
30 generators to get a composite that will be fast, easy to code
|
|
|
31 have a very long period and pass all the tests put to it.
|
|
|
32 The three components of KISS are
|
|
|
33
|
|
|
34 x(n)=a*x(n-1)+1 mod 2^32
|
|
|
35 y(n)=y(n-1)(I+L^13)(I+R^17)(I+L^5),
|
|
|
36 z(n)=2*z(n-1)+z(n-2) +carry mod 2^32
|
|
|
37
|
|
|
38 The y's are a shift register sequence on 32bit binary vectors
|
|
|
39 period 2^32-1; The z's are a simple multiply-with-carry sequence
|
|
|
40 with period 2^63+2^32-1.
|
|
|
41
|
|
|
42 The period of KISS is thus 2^32*(2^32-1)*(2^63+2^32-1) > 2^127
|
|
|
43
|
|
|
44 ******************************************************************************/
|
|
|
45
|
|
|
46 class Random
|
|
|
47 {
|
|
|
48 /**********************************************************************
|
|
|
49
|
|
|
50 Shared instance:
|
|
|
51 ---
|
|
|
52 auto random = Random.shared.next;
|
|
|
53 ---
|
|
|
54
|
|
|
55 **********************************************************************/
|
|
|
56 public static Random shared;
|
|
|
57
|
|
|
58 private uint kiss_k;
|
|
|
59 private uint kiss_m;
|
|
|
60 private uint kiss_x = 1;
|
|
|
61 private uint kiss_y = 2;
|
|
|
62 private uint kiss_z = 4;
|
|
|
63 private uint kiss_w = 8;
|
|
|
64 private uint kiss_carry = 0;
|
|
|
65
|
|
|
66 /**********************************************************************
|
|
|
67
|
|
|
68 Create a static and shared instance:
|
|
|
69 ---
|
|
|
70 auto random = Random.shared.next;
|
|
|
71 ---
|
|
|
72
|
|
|
73 **********************************************************************/
|
|
|
74
|
|
|
75 static this ()
|
|
|
76 {
|
|
|
77 shared = new Random;
|
|
|
78 }
|
|
|
79
|
|
|
80 /**********************************************************************
|
|
|
81
|
|
|
82 Creates and seeds a new generator with the current time
|
|
|
83
|
|
|
84 **********************************************************************/
|
|
|
85
|
|
|
86 this ()
|
|
|
87 {
|
|
|
88 this.seed;
|
|
|
89 }
|
|
|
90
|
|
|
91 /**********************************************************************
|
|
|
92
|
|
|
93 Seed the generator with current time
|
|
|
94
|
|
|
95 **********************************************************************/
|
|
|
96
|
|
|
97 final Random seed ()
|
|
|
98 {
|
|
|
99 ulong s;
|
|
|
100
|
|
|
101 version (Posix)
|
|
|
102 {
|
|
|
103 timeval tv;
|
|
|
104
|
|
|
105 gettimeofday (&tv, null);
|
|
|
106 s = tv.tv_usec;
|
|
|
107 }
|
|
|
108 version (Win32)
|
|
|
109 QueryPerformanceCounter (&s);
|
|
|
110
|
|
|
111 return seed (cast(uint) s);
|
|
|
112 }
|
|
|
113
|
|
|
114 /**********************************************************************
|
|
|
115
|
|
|
116 Seed the generator with a provided value
|
|
|
117
|
|
|
118 **********************************************************************/
|
|
|
119
|
|
|
120 final Random seed (uint seed)
|
|
|
121 {
|
|
|
122 kiss_x = seed | 1;
|
|
|
123 kiss_y = seed | 2;
|
|
|
124 kiss_z = seed | 4;
|
|
|
125 kiss_w = seed | 8;
|
|
|
126 kiss_carry = 0;
|
|
|
127 return this;
|
|
|
128 }
|
|
|
129
|
|
|
130 /**********************************************************************
|
|
|
131
|
|
|
132 Returns X such that 0 <= X <= uint.max
|
|
|
133
|
|
|
134 **********************************************************************/
|
|
|
135
|
|
|
136 final uint next ()
|
|
|
137 {
|
|
|
138 kiss_x = kiss_x * 69069 + 1;
|
|
|
139 kiss_y ^= kiss_y << 13;
|
|
|
140 kiss_y ^= kiss_y >> 17;
|
|
|
141 kiss_y ^= kiss_y << 5;
|
|
|
142 kiss_k = (kiss_z >> 2) + (kiss_w >> 3) + (kiss_carry >> 2);
|
|
|
143 kiss_m = kiss_w + kiss_w + kiss_z + kiss_carry;
|
|
|
144 kiss_z = kiss_w;
|
|
|
145 kiss_w = kiss_m;
|
|
|
146 kiss_carry = kiss_k >> 30;
|
|
|
147 return kiss_x + kiss_y + kiss_w;
|
|
|
148 }
|
|
|
149
|
|
|
150 /**********************************************************************
|
|
|
151
|
|
|
152 Returns X such that 0 <= X < max
|
|
|
153
|
|
|
154 Note that max is exclusive, making it compatible with
|
|
|
155 array indexing
|
|
|
156
|
|
|
157 **********************************************************************/
|
|
|
158
|
|
|
159 final uint next (uint max)
|
|
|
160 {
|
|
|
161 return next() % max;
|
|
|
162 }
|
|
|
163
|
|
|
164 /**********************************************************************
|
|
|
165
|
|
|
166 Returns X such that min <= X < max
|
|
|
167
|
|
|
168 Note that max is exclusive, making it compatible with
|
|
|
169 array indexing
|
|
|
170
|
|
|
171 **********************************************************************/
|
|
|
172
|
|
|
173 final uint next (uint min, uint max)
|
|
|
174 {
|
|
|
175 return next(max-min) + min;
|
|
|
176 }
|
|
|
177 }
|
|
|
178
|
