1/// <summary>
2/// SimpleRNG is a simple random number generator based on
3/// George Marsaglia's MWC (multiply with carry) generator.
4/// Although it is very simple, it passes Marsaglia's DIEHARD
5/// series of random number generator tests.
6///
7/// Written by John D. Cook
8/// http://www.johndcook.com
9/// Converted to C++ by Pietro Incardona
10/// </summary>
11
12#ifndef SIMPLE_RNG_HPP
13#define SIMPLE_RNG_HPP
14
15#include <string>
16#include <cmath>
17
18class SimpleRNG
19{
20private:
21 unsigned int m_w;
22 unsigned int m_z;
23
24 // This is the heart of the generator.
25 // It uses George Marsaglia's MWC algorithm to produce an unsigned integer.
26 // See http://www.bobwheeler.com/statistics/Password/MarsagliaPost.txt
27 unsigned int GetUint()
28 {
29 m_z = 36969 * (m_z & 65535) + (m_z >> 16);
30 m_w = 18000 * (m_w & 65535) + (m_w >> 16);
31 return (m_z << 16) + m_w;
32 }
33
34public:
35
36 SimpleRNG()
37 {
38 // These values are not magical, just the default values Marsaglia used.
39 // Any pair of unsigned integers should be fine.
40 m_w = 521288629;
41 m_z = 362436069;
42 }
43
44 // The random generator seed can be set three ways:
45 // 1) specifying two non-zero unsigned integers
46 // 2) specifying one non-zero unsigned integer and taking a default value for the second
47 // 3) setting the seed from the system time
48
49 void SetSeed(unsigned int u, unsigned int v)
50 {
51 if (u != 0) m_w = u;
52 if (v != 0) m_z = v;
53 }
54
55 void SetSeed(unsigned int u)
56 {
57 m_w = u;
58 }
59
60 void SetSeedFromSystemTime()
61 {
62 long x = clock();
63 SetSeed((unsigned int)(x >> 16), (unsigned int)(x % 4294967296));
64 }
65
66 // Produce a uniform random sample from the open interval (0, 1).
67 // The method will not return either end point.
68 double GetUniform()
69 {
70 // 0 <= u < 2^32
71 unsigned int u = GetUint();
72 // The magic number below is 1/(2^32 + 2).
73 // The result is strictly between 0 and 1.
74 return (u + 1.0) * 2.328306435454494e-10;
75 }
76
77 // Get normal (Gaussian) random sample with mean 0 and standard deviation 1
78 double GetNormal()
79 {
80 // Use Box-Muller algorithm
81 double u1 = GetUniform();
82 double u2 = GetUniform();
83 double r = std::sqrt( -2.0*std::log(u1) );
84 double theta = 2.0*3.14159265359*u2;
85 return r*std::sin(theta);
86 }
87
88 // Get normal (Gaussian) random sample with specified mean and standard deviation
89 double GetNormal(double mean, double standardDeviation)
90 {
91 if (standardDeviation <= 0.0)
92 {
93 std::cerr << __FILE__ << ":" << __LINE__ << " Shape must be positive. Received." << std::to_string(standardDeviation);
94 throw 100001;
95 }
96 return mean + standardDeviation*GetNormal();
97 }
98
99 // Get exponential random sample with mean 1
100 double GetExponential()
101 {
102 return -std::log( GetUniform() );
103 }
104
105 // Get exponential random sample with specified mean
106 double GetExponential(double mean)
107 {
108 if (mean <= 0.0)
109 {
110 std::cerr << __FILE__ << ":" << __LINE__ << "Mean must be positive. Received " << mean;
111 throw 100002;
112 }
113 return mean*GetExponential();
114 }
115
116 double GetGamma(double shape, double scale)
117 {
118 // Implementation based on "A Simple Method for Generating Gamma Variables"
119 // by George Marsaglia and Wai Wan Tsang. ACM Transactions on Mathematical Software
120 // Vol 26, No 3, September 2000, pages 363-372.
121
122 double d, c, x, xsquared, v, u;
123
124 if (shape >= 1.0)
125 {
126 d = shape - 1.0/3.0;
127 c = 1.0/std::sqrt(9.0*d);
128 for (;;)
129 {
130 do
131 {
132 x = GetNormal();
133 v = 1.0 + c*x;
134 }
135 while (v <= 0.0);
136 v = v*v*v;
137 u = GetUniform();
138 xsquared = x*x;
139 if (u < 1.0 -.0331*xsquared*xsquared || std::log(u) < 0.5*xsquared + d*(1.0 - v + std::log(v)))
140 return scale*d*v;
141 }
142 }
143 else if (shape <= 0.0)
144 {
145 std::cerr << __FILE__ << ":" << __LINE__ << " Shape must be positive. Received" << shape << "\n";
146 throw 100003;
147 }
148 else
149 {
150 double g = GetGamma(shape+1.0, 1.0);
151 double w = GetUniform();
152 return scale*g*std::pow(w, 1.0/shape);
153 }
154 }
155
156 double GetChiSquare(double degreesOfFreedom)
157 {
158 // A chi squared distribution with n degrees of freedom
159 // is a gamma distribution with shape n/2 and scale 2.
160 return GetGamma(0.5 * degreesOfFreedom, 2.0);
161 }
162
163 double GetInverseGamma(double shape, double scale)
164 {
165 // If X is gamma(shape, scale) then
166 // 1/Y is inverse gamma(shape, 1/scale)
167 return 1.0 / GetGamma(shape, 1.0 / scale);
168 }
169
170 double GetWeibull(double shape, double scale)
171 {
172 if (shape <= 0.0 || scale <= 0.0)
173 {
174 std::cerr << __FILE__ << ":" << __LINE__ << " Shape and scale parameters must be positive. Recieved " << shape << " and " << scale;
175 throw 100004;
176 }
177 return scale * std::pow(-std::log(GetUniform()), 1.0 / shape);
178 }
179
180 double GetCauchy(double median, double scale)
181 {
182 if (scale <= 0)
183 {
184 std::cerr << __FILE__ << ":" << __LINE__ << "Scale must be positive. Received " << scale << "\n";
185 throw 100005;
186 }
187
188 double p = GetUniform();
189
190 // Apply inverse of the Cauchy distribution function to a uniform
191 return median + scale*std::tan(3.14159265359*(p - 0.5));
192 }
193
194 double GetStudentT(double degreesOfFreedom)
195 {
196 if (degreesOfFreedom <= 0)
197 {
198 std::cerr << __FILE__ << ":" << __LINE__ << " Degrees of freedom must be positive. Received " << degreesOfFreedom;
199 throw 100006;
200 }
201
202 // See Seminumerical Algorithms by Knuth
203 double y1 = GetNormal();
204 double y2 = GetChiSquare(degreesOfFreedom);
205 return y1 / std::sqrt(y2 / degreesOfFreedom);
206 }
207
208 // The Laplace distribution is also known as the double exponential distribution.
209 double GetLaplace(double mean, double scale)
210 {
211 double u = GetUniform();
212 return (u < 0.5) ?
213 mean + scale*std::log(2.0*u) :
214 mean - scale*std::log(2*(1-u));
215 }
216
217 double GetLogNormal(double mu, double sigma)
218 {
219 return std::exp(GetNormal(mu, sigma));
220 }
221
222 double GetBeta(double a, double b)
223 {
224 if (a <= 0.0 || b <= 0.0)
225 {
226 std::cerr << __FILE__ << ":" << __LINE__ << "Beta parameters must be positive. Received " << a << " and " << b << "\n";
227 throw 100007;
228 }
229
230 // There are more efficient methods for generating beta samples.
231 // However such methods are a little more efficient and much more complicated.
232 // For an explanation of why the following method works, see
233 // http://www.johndcook.com/distribution_chart.html#gamma_beta
234
235 double u = GetGamma(a, 1.0);
236 double v = GetGamma(b, 1.0);
237 return u / (u + v);
238 }
239};
240
241#endif
242