1//===----------------------------------------------------------------------===//
2//
3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4// See https://llvm.org/LICENSE.txt for license information.
5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6//
7//===----------------------------------------------------------------------===//
8//
9// REQUIRES: long_tests
10
11// <random>
12
13// template<class IntType = int>
14// class binomial_distribution
15
16// template<class _URNG> result_type operator()(_URNG& g);
17
18#include <cassert>
19#include <cstdint>
20#include <numeric>
21#include <random>
22#include <type_traits>
23#include <vector>
24
25#include "test_macros.h"
26
27template <class T>
28T sqr(T x) {
29 return x * x;
30}
31
32template <class T>
33void test1() {
34 typedef std::binomial_distribution<T> D;
35 typedef std::mt19937_64 G;
36 G g;
37 D d(5, .75);
38 const int N = 1000000;
39 std::vector<typename D::result_type> u;
40 for (int i = 0; i < N; ++i)
41 {
42 typename D::result_type v = d(g);
43 assert(d.min() <= v && v <= d.max());
44 u.push_back(v);
45 }
46 double mean = std::accumulate(u.begin(), u.end(),
47 double(0)) / u.size();
48 double var = 0;
49 double skew = 0;
50 double kurtosis = 0;
51 for (unsigned i = 0; i < u.size(); ++i)
52 {
53 double dbl = (u[i] - mean);
54 double d2 = sqr(dbl);
55 var += d2;
56 skew += dbl * d2;
57 kurtosis += d2 * d2;
58 }
59 var /= u.size();
60 double dev = std::sqrt(x: var);
61 skew /= u.size() * dev * var;
62 kurtosis /= u.size() * var * var;
63 kurtosis -= 3;
64 double x_mean = d.t() * d.p();
65 double x_var = x_mean*(1-d.p());
66 double x_skew = (1-2*d.p()) / std::sqrt(x: x_var);
67 double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
68 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
69 assert(std::abs((var - x_var) / x_var) < 0.01);
70 assert(std::abs((skew - x_skew) / x_skew) < 0.01);
71 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.04);
72}
73
74template <class T>
75void test2() {
76 typedef std::binomial_distribution<T> D;
77 typedef std::mt19937 G;
78 G g;
79 D d(30, .03125);
80 const int N = 100000;
81 std::vector<typename D::result_type> u;
82 for (int i = 0; i < N; ++i)
83 {
84 typename D::result_type v = d(g);
85 assert(d.min() <= v && v <= d.max());
86 u.push_back(v);
87 }
88 double mean = std::accumulate(u.begin(), u.end(),
89 double(0)) / u.size();
90 double var = 0;
91 double skew = 0;
92 double kurtosis = 0;
93 for (unsigned i = 0; i < u.size(); ++i)
94 {
95 double dbl = (u[i] - mean);
96 double d2 = sqr(x: dbl);
97 var += d2;
98 skew += dbl * d2;
99 kurtosis += d2 * d2;
100 }
101 var /= u.size();
102 double dev = std::sqrt(x: var);
103 skew /= u.size() * dev * var;
104 kurtosis /= u.size() * var * var;
105 kurtosis -= 3;
106 double x_mean = d.t() * d.p();
107 double x_var = x_mean*(1-d.p());
108 double x_skew = (1-2*d.p()) / std::sqrt(x: x_var);
109 double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
110 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
111 assert(std::abs((var - x_var) / x_var) < 0.01);
112 assert(std::abs((skew - x_skew) / x_skew) < 0.01);
113 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
114}
115
116template <class T>
117void test3() {
118 typedef std::binomial_distribution<T> D;
119 typedef std::mt19937 G;
120 G g;
121 D d(40, .25);
122 const int N = 100000;
123 std::vector<typename D::result_type> u;
124 for (int i = 0; i < N; ++i)
125 {
126 typename D::result_type v = d(g);
127 assert(d.min() <= v && v <= d.max());
128 u.push_back(v);
129 }
130 double mean = std::accumulate(u.begin(), u.end(),
131 double(0)) / u.size();
132 double var = 0;
133 double skew = 0;
134 double kurtosis = 0;
135 for (unsigned i = 0; i < u.size(); ++i)
136 {
137 double dbl = (u[i] - mean);
138 double d2 = sqr(x: dbl);
139 var += d2;
140 skew += dbl * d2;
141 kurtosis += d2 * d2;
142 }
143 var /= u.size();
144 double dev = std::sqrt(x: var);
145 skew /= u.size() * dev * var;
146 kurtosis /= u.size() * var * var;
147 kurtosis -= 3;
148 double x_mean = d.t() * d.p();
149 double x_var = x_mean*(1-d.p());
150 double x_skew = (1-2*d.p()) / std::sqrt(x: x_var);
151 double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
152 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
153 assert(std::abs((var - x_var) / x_var) < 0.01);
154 assert(std::abs((skew - x_skew) / x_skew) < 0.03);
155 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.3);
156}
157
158template <class T>
159void test4() {
160 typedef std::binomial_distribution<T> D;
161 typedef std::mt19937 G;
162 G g;
163 D d(40, 0);
164 const int N = 100000;
165 std::vector<typename D::result_type> u;
166 for (int i = 0; i < N; ++i)
167 {
168 typename D::result_type v = d(g);
169 assert(d.min() <= v && v <= d.max());
170 u.push_back(v);
171 }
172 double mean = std::accumulate(u.begin(), u.end(),
173 double(0)) / u.size();
174 double var = 0;
175 double skew = 0;
176 double kurtosis = 0;
177 for (unsigned i = 0; i < u.size(); ++i)
178 {
179 double dbl = (u[i] - mean);
180 double d2 = sqr(x: dbl);
181 var += d2;
182 skew += dbl * d2;
183 kurtosis += d2 * d2;
184 }
185 var /= u.size();
186 double dev = std::sqrt(x: var);
187 // In this case:
188 // skew computes to 0./0. == nan
189 // kurtosis computes to 0./0. == nan
190 // x_skew == inf
191 // x_kurtosis == inf
192 skew /= u.size() * dev * var;
193 kurtosis /= u.size() * var * var;
194 kurtosis -= 3;
195 double x_mean = d.t() * d.p();
196 double x_var = x_mean*(1-d.p());
197 double x_skew = (1-2*d.p()) / std::sqrt(x: x_var);
198 double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
199 assert(mean == x_mean);
200 assert(var == x_var);
201 // assert(skew == x_skew);
202 (void)skew; (void)x_skew;
203 // assert(kurtosis == x_kurtosis);
204 (void)kurtosis; (void)x_kurtosis;
205}
206
207template <class T>
208void test5() {
209 typedef std::binomial_distribution<T> D;
210 typedef std::mt19937 G;
211 G g;
212 D d(40, 1);
213 const int N = 100000;
214 std::vector<typename D::result_type> u;
215 for (int i = 0; i < N; ++i)
216 {
217 typename D::result_type v = d(g);
218 assert(d.min() <= v && v <= d.max());
219 u.push_back(v);
220 }
221 double mean = std::accumulate(u.begin(), u.end(),
222 double(0)) / u.size();
223 double var = 0;
224 double skew = 0;
225 double kurtosis = 0;
226 for (unsigned i = 0; i < u.size(); ++i)
227 {
228 double dbl = (u[i] - mean);
229 double d2 = sqr(x: dbl);
230 var += d2;
231 skew += dbl * d2;
232 kurtosis += d2 * d2;
233 }
234 var /= u.size();
235 double dev = std::sqrt(x: var);
236 // In this case:
237 // skew computes to 0./0. == nan
238 // kurtosis computes to 0./0. == nan
239 // x_skew == -inf
240 // x_kurtosis == inf
241 skew /= u.size() * dev * var;
242 kurtosis /= u.size() * var * var;
243 kurtosis -= 3;
244 double x_mean = d.t() * d.p();
245 double x_var = x_mean*(1-d.p());
246 double x_skew = (1-2*d.p()) / std::sqrt(x: x_var);
247 double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
248 assert(mean == x_mean);
249 assert(var == x_var);
250 // assert(skew == x_skew);
251 (void)skew; (void)x_skew;
252 // assert(kurtosis == x_kurtosis);
253 (void)kurtosis; (void)x_kurtosis;
254}
255
256template <class T>
257void test6() {
258 typedef std::binomial_distribution<T> D;
259 typedef std::mt19937 G;
260 G g;
261 D d(127, 0.5);
262 const int N = 100000;
263 std::vector<typename D::result_type> u;
264 for (int i = 0; i < N; ++i)
265 {
266 typename D::result_type v = d(g);
267 assert(d.min() <= v && v <= d.max());
268 u.push_back(v);
269 }
270 double mean = std::accumulate(u.begin(), u.end(),
271 double(0)) / u.size();
272 double var = 0;
273 double skew = 0;
274 double kurtosis = 0;
275 for (unsigned i = 0; i < u.size(); ++i)
276 {
277 double dbl = (u[i] - mean);
278 double d2 = sqr(x: dbl);
279 var += d2;
280 skew += dbl * d2;
281 kurtosis += d2 * d2;
282 }
283 var /= u.size();
284 double dev = std::sqrt(x: var);
285 skew /= u.size() * dev * var;
286 kurtosis /= u.size() * var * var;
287 kurtosis -= 3;
288 double x_mean = d.t() * d.p();
289 double x_var = x_mean*(1-d.p());
290 double x_skew = (1-2*d.p()) / std::sqrt(x: x_var);
291 double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
292 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
293 assert(std::abs((var - x_var) / x_var) < 0.01);
294 assert(std::abs(skew - x_skew) < 0.02);
295 assert(std::abs(kurtosis - x_kurtosis) < 0.01);
296}
297
298template <class T>
299void test7() {
300 typedef std::binomial_distribution<T> D;
301 typedef std::mt19937 G;
302 G g;
303 D d(1, 0.5);
304 const int N = 100000;
305 std::vector<typename D::result_type> u;
306 for (int i = 0; i < N; ++i)
307 {
308 typename D::result_type v = d(g);
309 assert(d.min() <= v && v <= d.max());
310 u.push_back(v);
311 }
312 double mean = std::accumulate(u.begin(), u.end(),
313 double(0)) / u.size();
314 double var = 0;
315 double skew = 0;
316 double kurtosis = 0;
317 for (unsigned i = 0; i < u.size(); ++i)
318 {
319 double dbl = (u[i] - mean);
320 double d2 = sqr(x: dbl);
321 var += d2;
322 skew += dbl * d2;
323 kurtosis += d2 * d2;
324 }
325 var /= u.size();
326 double dev = std::sqrt(x: var);
327 skew /= u.size() * dev * var;
328 kurtosis /= u.size() * var * var;
329 kurtosis -= 3;
330 double x_mean = d.t() * d.p();
331 double x_var = x_mean*(1-d.p());
332 double x_skew = (1-2*d.p()) / std::sqrt(x: x_var);
333 double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
334 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
335 assert(std::abs((var - x_var) / x_var) < 0.01);
336 assert(std::abs(skew - x_skew) < 0.01);
337 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
338}
339
340template <class T>
341void test8() {
342 const int N = 100000;
343 std::mt19937 gen1;
344 std::mt19937 gen2;
345
346 using UnsignedT = typename std::make_unsigned<T>::type;
347 std::binomial_distribution<T> dist1(5, 0.1);
348 std::binomial_distribution<UnsignedT> dist2(5, 0.1);
349
350 for (int i = 0; i < N; ++i) {
351 T r1 = dist1(gen1);
352 UnsignedT r2 = dist2(gen2);
353 assert(r1 >= 0);
354 assert(static_cast<UnsignedT>(r1) == r2);
355 }
356}
357
358template <class T>
359void test9() {
360 typedef std::binomial_distribution<T> D;
361 typedef std::mt19937 G;
362 G g;
363 D d(0, 0.005);
364 const int N = 100000;
365 std::vector<typename D::result_type> u;
366 for (int i = 0; i < N; ++i)
367 {
368 typename D::result_type v = d(g);
369 assert(d.min() <= v && v <= d.max());
370 u.push_back(v);
371 }
372 double mean = std::accumulate(u.begin(), u.end(),
373 double(0)) / u.size();
374 double var = 0;
375 double skew = 0;
376 double kurtosis = 0;
377 for (unsigned i = 0; i < u.size(); ++i)
378 {
379 double dbl = (u[i] - mean);
380 double d2 = sqr(x: dbl);
381 var += d2;
382 skew += dbl * d2;
383 kurtosis += d2 * d2;
384 }
385 var /= u.size();
386 double dev = std::sqrt(x: var);
387 // In this case:
388 // skew computes to 0./0. == nan
389 // kurtosis computes to 0./0. == nan
390 // x_skew == inf
391 // x_kurtosis == inf
392 skew /= u.size() * dev * var;
393 kurtosis /= u.size() * var * var;
394 kurtosis -= 3;
395 double x_mean = d.t() * d.p();
396 double x_var = x_mean*(1-d.p());
397 double x_skew = (1-2*d.p()) / std::sqrt(x: x_var);
398 double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
399 assert(mean == x_mean);
400 assert(var == x_var);
401 // assert(skew == x_skew);
402 (void)skew; (void)x_skew;
403 // assert(kurtosis == x_kurtosis);
404 (void)kurtosis; (void)x_kurtosis;
405}
406
407template <class T>
408void test10() {
409 typedef std::binomial_distribution<T> D;
410 typedef std::mt19937 G;
411 G g;
412 D d(0, 0);
413 const int N = 100000;
414 std::vector<typename D::result_type> u;
415 for (int i = 0; i < N; ++i)
416 {
417 typename D::result_type v = d(g);
418 assert(d.min() <= v && v <= d.max());
419 u.push_back(v);
420 }
421 double mean = std::accumulate(u.begin(), u.end(),
422 double(0)) / u.size();
423 double var = 0;
424 double skew = 0;
425 double kurtosis = 0;
426 for (unsigned i = 0; i < u.size(); ++i)
427 {
428 double dbl = (u[i] - mean);
429 double d2 = sqr(x: dbl);
430 var += d2;
431 skew += dbl * d2;
432 kurtosis += d2 * d2;
433 }
434 var /= u.size();
435 double dev = std::sqrt(x: var);
436 // In this case:
437 // skew computes to 0./0. == nan
438 // kurtosis computes to 0./0. == nan
439 // x_skew == inf
440 // x_kurtosis == inf
441 skew /= u.size() * dev * var;
442 kurtosis /= u.size() * var * var;
443 kurtosis -= 3;
444 double x_mean = d.t() * d.p();
445 double x_var = x_mean*(1-d.p());
446 double x_skew = (1-2*d.p()) / std::sqrt(x: x_var);
447 double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
448 assert(mean == x_mean);
449 assert(var == x_var);
450 // assert(skew == x_skew);
451 (void)skew; (void)x_skew;
452 // assert(kurtosis == x_kurtosis);
453 (void)kurtosis; (void)x_kurtosis;
454}
455
456template <class T>
457void test11() {
458 typedef std::binomial_distribution<T> D;
459 typedef std::mt19937 G;
460 G g;
461 D d(0, 1);
462 const int N = 100000;
463 std::vector<typename D::result_type> u;
464 for (int i = 0; i < N; ++i)
465 {
466 typename D::result_type v = d(g);
467 assert(d.min() <= v && v <= d.max());
468 u.push_back(v);
469 }
470 double mean = std::accumulate(u.begin(), u.end(),
471 double(0)) / u.size();
472 double var = 0;
473 double skew = 0;
474 double kurtosis = 0;
475 for (unsigned i = 0; i < u.size(); ++i)
476 {
477 double dbl = (u[i] - mean);
478 double d2 = sqr(x: dbl);
479 var += d2;
480 skew += dbl * d2;
481 kurtosis += d2 * d2;
482 }
483 var /= u.size();
484 double dev = std::sqrt(x: var);
485 // In this case:
486 // skew computes to 0./0. == nan
487 // kurtosis computes to 0./0. == nan
488 // x_skew == -inf
489 // x_kurtosis == inf
490 skew /= u.size() * dev * var;
491 kurtosis /= u.size() * var * var;
492 kurtosis -= 3;
493 double x_mean = d.t() * d.p();
494 double x_var = x_mean*(1-d.p());
495 double x_skew = (1-2*d.p()) / std::sqrt(x: x_var);
496 double x_kurtosis = (1-6*d.p()*(1-d.p())) / x_var;
497 assert(mean == x_mean);
498 assert(var == x_var);
499 // assert(skew == x_skew);
500 (void)skew; (void)x_skew;
501 // assert(kurtosis == x_kurtosis);
502 (void)kurtosis; (void)x_kurtosis;
503}
504
505template <class T>
506void tests() {
507 test1<T>();
508 test2<T>();
509 test3<T>();
510 test4<T>();
511 test5<T>();
512 test6<T>();
513 test7<T>();
514 test8<T>();
515 test9<T>();
516 test10<T>();
517 test11<T>();
518}
519
520int main(int, char**) {
521 tests<short>();
522 tests<int>();
523 tests<long>();
524 tests<long long>();
525
526 tests<unsigned short>();
527 tests<unsigned int>();
528 tests<unsigned long>();
529 tests<unsigned long long>();
530
531#if defined(_LIBCPP_VERSION) // extension
532 tests<std::int8_t>();
533 tests<std::uint8_t>();
534#if !defined(TEST_HAS_NO_INT128)
535 tests<__int128_t>();
536 tests<__uint128_t>();
537#endif
538#endif
539
540 return 0;
541}
542

source code of libcxx/test/std/numerics/rand/rand.dist/rand.dist.bern/rand.dist.bern.bin/eval.pass.cpp