1 //===----------------------------------------------------------------------===//
2 //
3 // The LLVM Compiler Infrastructure
4 //
5 // This file is dual licensed under the MIT and the University of Illinois Open
6 // Source Licenses. See LICENSE.TXT for details.
7 //
8 //===----------------------------------------------------------------------===//
9 //
10 // REQUIRES: long_tests
11
12 // <random>
13
14 // template<class IntType = int>
15 // class geometric_distribution
16
17 // template<class _URNG> result_type operator()(_URNG& g);
18
19 #include <random>
20 #include <numeric>
21 #include <vector>
22 #include <cassert>
23
24 template <class T>
25 inline
26 T
sqr(T x)27 sqr(T x)
28 {
29 return x * x;
30 }
31
32 struct Eng : std::mt19937 {
33 using Base = std::mt19937;
34 using Base::Base;
35 };
36
test_small_inputs()37 void test_small_inputs() {
38 Eng engine;
39 std::geometric_distribution<std::int16_t> distribution(5.45361e-311);
40 for (auto i=0; i < 1000; ++i) {
41 volatile auto res = distribution(engine);
42 ((void)res);
43 }
44 }
45
46 void
test1()47 test1()
48 {
49 typedef std::geometric_distribution<> D;
50 typedef std::mt19937 G;
51 G g;
52 D d(.03125);
53 const int N = 1000000;
54 std::vector<D::result_type> u;
55 for (int i = 0; i < N; ++i)
56 {
57 D::result_type v = d(g);
58 assert(d.min() <= v && v <= d.max());
59 u.push_back(v);
60 }
61 double mean = std::accumulate(u.begin(), u.end(),
62 double(0)) / u.size();
63 double var = 0;
64 double skew = 0;
65 double kurtosis = 0;
66 for (unsigned i = 0; i < u.size(); ++i)
67 {
68 double dbl = (u[i] - mean);
69 double d2 = sqr(dbl);
70 var += d2;
71 skew += dbl * d2;
72 kurtosis += d2 * d2;
73 }
74 var /= u.size();
75 double dev = std::sqrt(var);
76 skew /= u.size() * dev * var;
77 kurtosis /= u.size() * var * var;
78 kurtosis -= 3;
79 double x_mean = (1 - d.p()) / d.p();
80 double x_var = x_mean / d.p();
81 double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
82 double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
83 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
84 assert(std::abs((var - x_var) / x_var) < 0.01);
85 assert(std::abs((skew - x_skew) / x_skew) < 0.01);
86 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
87 }
88
89 void
test2()90 test2()
91 {
92 typedef std::geometric_distribution<> D;
93 typedef std::mt19937 G;
94 G g;
95 D d(0.05);
96 const int N = 1000000;
97 std::vector<D::result_type> u;
98 for (int i = 0; i < N; ++i)
99 {
100 D::result_type v = d(g);
101 assert(d.min() <= v && v <= d.max());
102 u.push_back(v);
103 }
104 double mean = std::accumulate(u.begin(), u.end(),
105 double(0)) / u.size();
106 double var = 0;
107 double skew = 0;
108 double kurtosis = 0;
109 for (unsigned i = 0; i < u.size(); ++i)
110 {
111 double dbl = (u[i] - mean);
112 double d2 = sqr(dbl);
113 var += d2;
114 skew += dbl * d2;
115 kurtosis += d2 * d2;
116 }
117 var /= u.size();
118 double dev = std::sqrt(var);
119 skew /= u.size() * dev * var;
120 kurtosis /= u.size() * var * var;
121 kurtosis -= 3;
122 double x_mean = (1 - d.p()) / d.p();
123 double x_var = x_mean / d.p();
124 double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
125 double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
126 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
127 assert(std::abs((var - x_var) / x_var) < 0.01);
128 assert(std::abs((skew - x_skew) / x_skew) < 0.01);
129 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
130 }
131
132 void
test3()133 test3()
134 {
135 typedef std::geometric_distribution<> D;
136 typedef std::minstd_rand G;
137 G g;
138 D d(.25);
139 const int N = 1000000;
140 std::vector<D::result_type> u;
141 for (int i = 0; i < N; ++i)
142 {
143 D::result_type v = d(g);
144 assert(d.min() <= v && v <= d.max());
145 u.push_back(v);
146 }
147 double mean = std::accumulate(u.begin(), u.end(),
148 double(0)) / u.size();
149 double var = 0;
150 double skew = 0;
151 double kurtosis = 0;
152 for (unsigned i = 0; i < u.size(); ++i)
153 {
154 double dbl = (u[i] - mean);
155 double d2 = sqr(dbl);
156 var += d2;
157 skew += dbl * d2;
158 kurtosis += d2 * d2;
159 }
160 var /= u.size();
161 double dev = std::sqrt(var);
162 skew /= u.size() * dev * var;
163 kurtosis /= u.size() * var * var;
164 kurtosis -= 3;
165 double x_mean = (1 - d.p()) / d.p();
166 double x_var = x_mean / d.p();
167 double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
168 double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
169 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
170 assert(std::abs((var - x_var) / x_var) < 0.01);
171 assert(std::abs((skew - x_skew) / x_skew) < 0.01);
172 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
173 }
174
175 void
test4()176 test4()
177 {
178 typedef std::geometric_distribution<> D;
179 typedef std::mt19937 G;
180 G g;
181 D d(0.5);
182 const int N = 1000000;
183 std::vector<D::result_type> u;
184 for (int i = 0; i < N; ++i)
185 {
186 D::result_type v = d(g);
187 assert(d.min() <= v && v <= d.max());
188 u.push_back(v);
189 }
190 double mean = std::accumulate(u.begin(), u.end(),
191 double(0)) / u.size();
192 double var = 0;
193 double skew = 0;
194 double kurtosis = 0;
195 for (unsigned i = 0; i < u.size(); ++i)
196 {
197 double dbl = (u[i] - mean);
198 double d2 = sqr(dbl);
199 var += d2;
200 skew += dbl * d2;
201 kurtosis += d2 * d2;
202 }
203 var /= u.size();
204 double dev = std::sqrt(var);
205 skew /= u.size() * dev * var;
206 kurtosis /= u.size() * var * var;
207 kurtosis -= 3;
208 double x_mean = (1 - d.p()) / d.p();
209 double x_var = x_mean / d.p();
210 double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
211 double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
212 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
213 assert(std::abs((var - x_var) / x_var) < 0.01);
214 assert(std::abs((skew - x_skew) / x_skew) < 0.01);
215 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
216 }
217
218 void
test5()219 test5()
220 {
221 typedef std::geometric_distribution<> D;
222 typedef std::mt19937 G;
223 G g;
224 D d(0.75);
225 const int N = 1000000;
226 std::vector<D::result_type> u;
227 for (int i = 0; i < N; ++i)
228 {
229 D::result_type v = d(g);
230 assert(d.min() <= v && v <= d.max());
231 u.push_back(v);
232 }
233 double mean = std::accumulate(u.begin(), u.end(),
234 double(0)) / u.size();
235 double var = 0;
236 double skew = 0;
237 double kurtosis = 0;
238 for (unsigned i = 0; i < u.size(); ++i)
239 {
240 double dbl = (u[i] - mean);
241 double d2 = sqr(dbl);
242 var += d2;
243 skew += dbl * d2;
244 kurtosis += d2 * d2;
245 }
246 var /= u.size();
247 double dev = std::sqrt(var);
248 skew /= u.size() * dev * var;
249 kurtosis /= u.size() * var * var;
250 kurtosis -= 3;
251 double x_mean = (1 - d.p()) / d.p();
252 double x_var = x_mean / d.p();
253 double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
254 double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
255 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
256 assert(std::abs((var - x_var) / x_var) < 0.01);
257 assert(std::abs((skew - x_skew) / x_skew) < 0.01);
258 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
259 }
260
261 void
test6()262 test6()
263 {
264 typedef std::geometric_distribution<> D;
265 typedef std::mt19937 G;
266 G g;
267 D d(0.96875);
268 const int N = 1000000;
269 std::vector<D::result_type> u;
270 for (int i = 0; i < N; ++i)
271 {
272 D::result_type v = d(g);
273 assert(d.min() <= v && v <= d.max());
274 u.push_back(v);
275 }
276 double mean = std::accumulate(u.begin(), u.end(),
277 double(0)) / u.size();
278 double var = 0;
279 double skew = 0;
280 double kurtosis = 0;
281 for (unsigned i = 0; i < u.size(); ++i)
282 {
283 double dbl = (u[i] - mean);
284 double d2 = sqr(dbl);
285 var += d2;
286 skew += dbl * d2;
287 kurtosis += d2 * d2;
288 }
289 var /= u.size();
290 double dev = std::sqrt(var);
291 skew /= u.size() * dev * var;
292 kurtosis /= u.size() * var * var;
293 kurtosis -= 3;
294 double x_mean = (1 - d.p()) / d.p();
295 double x_var = x_mean / d.p();
296 double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
297 double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
298 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
299 assert(std::abs((var - x_var) / x_var) < 0.01);
300 assert(std::abs((skew - x_skew) / x_skew) < 0.01);
301 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
302 }
303
main()304 int main()
305 {
306 test1();
307 test2();
308 test3();
309 test4();
310 test5();
311 test6();
312 test_small_inputs();
313 }
314