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