1 /*******************************************************************************
2
3 With gratitude to Dr Jurgen A Doornik. See his paper entitled
4 "Conversion of high-period random numbers to floating point"
5
6 Copyright:
7 Copyright (c) 2008 Tango contributors.
8 Some parts copyright (c) 2009-2016 dunnhumby Germany GmbH.
9 All rights reserved.
10
11 License:
12 Tango Dual License: 3-Clause BSD License / Academic Free License v3.0.
13 See LICENSE_TANGO.txt for details.
14
15 Version: Initial release: May 2008
16
17 Authors: Various
18
19 *******************************************************************************/
20
21 module ocean.math.random.Kiss;
22
23 import core.sys.posix.sys.time : timeval, gettimeofday;
24
25 /******************************************************************************
26
27 KISS (from George Marsaglia)
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
39 The y's are a shift register sequence on 32bit binary vectors
40 period 2^32-1; The z's are a simple multiply-with-carry sequence
41 with period 2^63+2^32-1. The period of KISS is thus
42 ---
43 2^32*(2^32-1)*(2^63+2^32-1) > 2^127
44 ---
45
46 Note that this should be passed by reference, unless you really
47 intend to provide a local copy to a callee
48
49 ******************************************************************************/
50
51 struct Kiss
52 {
53 ///
54 public alias natural toInt;
55 ///
56 public alias fraction toReal;
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 private enum double M_RAN_INVM32 = 2.32830643653869628906e-010,
67 M_RAN_INVM52 = 2.22044604925031308085e-016;
68
69 /**********************************************************************
70
71 A global, shared instance, seeded via startup time
72
73 **********************************************************************/
74
75 public static Kiss instance;
76
77 static this ()
78 {
79 instance.seed;
80 }
81
82 /**********************************************************************
83
84 Creates and seeds a new generator with the current time
85
86 **********************************************************************/
87
88 static Kiss opCall ()
89 {
90 Kiss rand;
91 rand.seed;
92 return rand;
93 }
94
95 /**********************************************************************
96
97 Seed the generator with current time
98
99 **********************************************************************/
100
101 void seed ()
102 {
103 ulong s;
104
105 timeval tv;
106
107 gettimeofday (&tv, null);
108 s = tv.tv_usec;
109
110 return seed (cast(uint) s);
111 }
112
113 /**********************************************************************
114
115 Seed the generator with a provided value
116
117 **********************************************************************/
118
119 void seed (uint seed)
120 {
121 kiss_x = seed | 1;
122 kiss_y = seed | 2;
123 kiss_z = seed | 4;
124 kiss_w = seed | 8;
125 kiss_carry = 0;
126 }
127
128 /**********************************************************************
129
130 Returns X such that 0 <= X <= uint.max
131
132 **********************************************************************/
133
134 uint natural ()
135 {
136 kiss_x = kiss_x * 69069 + 1;
137 kiss_y ^= kiss_y << 13;
138 kiss_y ^= kiss_y >> 17;
139 kiss_y ^= kiss_y << 5;
140 kiss_k = (kiss_z >> 2) + (kiss_w >> 3) + (kiss_carry >> 2);
141 kiss_m = kiss_w + kiss_w + kiss_z + kiss_carry;
142 kiss_z = kiss_w;
143 kiss_w = kiss_m;
144 kiss_carry = kiss_k >> 30;
145 return kiss_x + kiss_y + kiss_w;
146 }
147
148 /**********************************************************************
149
150 Returns X such that 0 <= X < max
151
152 Note that max is exclusive, making it compatible with
153 array indexing
154
155 **********************************************************************/
156
157 uint natural (uint max)
158 {
159 return natural % max;
160 }
161
162 /**********************************************************************
163
164 Returns X such that min <= X < max
165
166 Note that max is exclusive, making it compatible with
167 array indexing
168
169 **********************************************************************/
170
171 uint natural (uint min, uint max)
172 {
173 return (natural % (max-min)) + min;
174 }
175
176 /**********************************************************************
177
178 Returns a value in the range [0, 1$(RPAREN) using 32 bits
179 of precision (with thanks to Dr Jurgen A Doornik)
180
181 **********************************************************************/
182
183 double fraction ()
184 {
185 return ((cast(int) natural) * M_RAN_INVM32 + (0.5 + M_RAN_INVM32 / 2));
186 }
187
188 /**********************************************************************
189
190 Returns a value in the range [0, 1$(RPAREN) using 52 bits
191 of precision (with thanks to Dr Jurgen A Doornik)
192
193 **********************************************************************/
194
195 double fractionEx ()
196 {
197 return ((cast(int) natural) * M_RAN_INVM32 + (0.5 + M_RAN_INVM52 / 2) +
198 ((cast(int) natural) & 0x000FFFFF) * M_RAN_INVM52);
199 }
200 }
201
202
203
204 /******************************************************************************
205
206
207 ******************************************************************************/
208
209 debug (Kiss)
210 {
211 import ocean.io.Stdout;
212 import ocean.time.StopWatch;
213
214 void main()
215 {
216 auto dbl = Kiss();
217 auto count = 100_000_000;
218 StopWatch w;
219
220 w.start;
221 double v1;
222 for (int i=count; --i;)
223 v1 = dbl.fraction;
224 Stdout.formatln ("{} fraction, {}/s, {:f10}", count, count/w.stop, v1);
225
226 w.start;
227 for (int i=count; --i;)
228 v1 = dbl.fractionEx;
229 Stdout.formatln ("{} fractionEx, {}/s, {:f10}", count, count/w.stop, v1);
230
231 for (int i=count; --i;)
232 {
233 auto v = dbl.fraction;
234 if (v <= 0.0 || v >= 1.0)
235 {
236 Stdout.formatln ("fraction {:f10}", v);
237 break;
238 }
239 v = dbl.fractionEx;
240 if (v <= 0.0 || v >= 1.0)
241 {
242 Stdout.formatln ("fractionEx {:f10}", v);
243 break;
244 }
245 }
246 }
247 }