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

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