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 RealType = double>
14// class piecewise_constant_distribution
15
16// template<class _URNG> result_type operator()(_URNG& g, const param_type& parm);
17
18#include <random>
19#include <algorithm>
20#include <cassert>
21#include <cmath>
22#include <cstddef>
23#include <iterator>
24#include <numeric>
25#include <vector>
26
27#include "test_macros.h"
28
29template <class T>
30inline
31T
32sqr(T x)
33{
34 return x*x;
35}
36
37int main(int, char**)
38{
39 {
40 typedef std::piecewise_constant_distribution<> D;
41 typedef D::param_type P;
42 typedef std::mt19937_64 G;
43 G g;
44 double b[] = {10, 14, 16, 17};
45 double p[] = {25, 62.5, 12.5};
46 const std::size_t Np = sizeof(p) / sizeof(p[0]);
47 D d;
48 P pa(b, b+Np+1, p);
49 const int N = 1000000;
50 std::vector<D::result_type> u;
51 for (int i = 0; i < N; ++i)
52 {
53 D::result_type v = d(g, pa);
54 assert(10 <= v && v < 17);
55 u.push_back(v);
56 }
57 std::vector<double> prob(std::begin(p), std::end(p));
58 double s = std::accumulate(prob.begin(), prob.end(), 0.0);
59 for (std::size_t i = 0; i < prob.size(); ++i)
60 prob[i] /= s;
61 std::sort(u.begin(), u.end());
62 for (std::size_t i = 0; i < Np; ++i)
63 {
64 typedef std::vector<D::result_type>::iterator I;
65 I lb = std::lower_bound(u.begin(), u.end(), b[i]);
66 I ub = std::lower_bound(u.begin(), u.end(), b[i+1]);
67 const std::size_t Ni = ub - lb;
68 if (prob[i] == 0)
69 assert(Ni == 0);
70 else
71 {
72 assert(std::abs((double)Ni/N - prob[i]) / prob[i] < .01);
73 double mean = std::accumulate(lb, ub, 0.0) / Ni;
74 double var = 0;
75 double skew = 0;
76 double kurtosis = 0;
77 for (I j = lb; j != ub; ++j)
78 {
79 double dbl = (*j - mean);
80 double d2 = sqr(dbl);
81 var += d2;
82 skew += dbl * d2;
83 kurtosis += d2 * d2;
84 }
85 var /= Ni;
86 double dev = std::sqrt(x: var);
87 skew /= Ni * dev * var;
88 kurtosis /= Ni * var * var;
89 kurtosis -= 3;
90 double x_mean = (b[i+1] + b[i]) / 2;
91 double x_var = sqr(b[i+1] - b[i]) / 12;
92 double x_skew = 0;
93 double x_kurtosis = -6./5;
94 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
95 assert(std::abs((var - x_var) / x_var) < 0.01);
96 assert(std::abs(skew - x_skew) < 0.01);
97 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.01);
98 }
99 }
100 }
101
102 return 0;
103}
104

source code of libcxx/test/std/numerics/rand/rand.dist/rand.dist.samp/rand.dist.samp.pconst/eval_param.pass.cpp