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
24template <class T>
25inline
26T
27sqr(T x)
28{
29    return x * x;
30}
31
32void
33test1()
34{
35    typedef std::geometric_distribution<> D;
36    typedef std::mt19937 G;
37    G g;
38    D d(.03125);
39    const int N = 1000000;
40    std::vector<D::result_type> u;
41    for (int i = 0; i < N; ++i)
42    {
43        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(var);
62    skew /= u.size() * dev * var;
63    kurtosis /= u.size() * var * var;
64    kurtosis -= 3;
65    double x_mean = (1 - d.p()) / d.p();
66    double x_var = x_mean / d.p();
67    double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
68    double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
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.01);
73}
74
75void
76test2()
77{
78    typedef std::geometric_distribution<> D;
79    typedef std::mt19937 G;
80    G g;
81    D d(0.05);
82    const int N = 1000000;
83    std::vector<D::result_type> u;
84    for (int i = 0; i < N; ++i)
85    {
86        D::result_type v = d(g);
87        assert(d.min() <= v && v <= d.max());
88        u.push_back(v);
89    }
90    double mean = std::accumulate(u.begin(), u.end(),
91                                          double(0)) / u.size();
92    double var = 0;
93    double skew = 0;
94    double kurtosis = 0;
95    for (unsigned i = 0; i < u.size(); ++i)
96    {
97        double dbl = (u[i] - mean);
98        double d2 = sqr(dbl);
99        var += d2;
100        skew += dbl * d2;
101        kurtosis += d2 * d2;
102    }
103    var /= u.size();
104    double dev = std::sqrt(var);
105    skew /= u.size() * dev * var;
106    kurtosis /= u.size() * var * var;
107    kurtosis -= 3;
108    double x_mean = (1 - d.p()) / d.p();
109    double x_var = x_mean / d.p();
110    double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
111    double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
112    assert(std::abs((mean - x_mean) / x_mean) < 0.01);
113    assert(std::abs((var - x_var) / x_var) < 0.01);
114    assert(std::abs((skew - x_skew) / x_skew) < 0.01);
115    assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
116}
117
118void
119test3()
120{
121    typedef std::geometric_distribution<> D;
122    typedef std::minstd_rand G;
123    G g;
124    D d(.25);
125    const int N = 1000000;
126    std::vector<D::result_type> u;
127    for (int i = 0; i < N; ++i)
128    {
129        D::result_type v = d(g);
130        assert(d.min() <= v && v <= d.max());
131        u.push_back(v);
132    }
133    double mean = std::accumulate(u.begin(), u.end(),
134                                          double(0)) / u.size();
135    double var = 0;
136    double skew = 0;
137    double kurtosis = 0;
138    for (unsigned i = 0; i < u.size(); ++i)
139    {
140        double dbl = (u[i] - mean);
141        double d2 = sqr(dbl);
142        var += d2;
143        skew += dbl * d2;
144        kurtosis += d2 * d2;
145    }
146    var /= u.size();
147    double dev = std::sqrt(var);
148    skew /= u.size() * dev * var;
149    kurtosis /= u.size() * var * var;
150    kurtosis -= 3;
151    double x_mean = (1 - d.p()) / d.p();
152    double x_var = x_mean / d.p();
153    double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
154    double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
155    assert(std::abs((mean - x_mean) / x_mean) < 0.01);
156    assert(std::abs((var - x_var) / x_var) < 0.01);
157    assert(std::abs((skew - x_skew) / x_skew) < 0.01);
158    assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
159}
160
161void
162test4()
163{
164    typedef std::geometric_distribution<> D;
165    typedef std::mt19937 G;
166    G g;
167    D d(0.5);
168    const int N = 1000000;
169    std::vector<D::result_type> u;
170    for (int i = 0; i < N; ++i)
171    {
172        D::result_type v = d(g);
173        assert(d.min() <= v && v <= d.max());
174        u.push_back(v);
175    }
176    double mean = std::accumulate(u.begin(), u.end(),
177                                          double(0)) / u.size();
178    double var = 0;
179    double skew = 0;
180    double kurtosis = 0;
181    for (unsigned i = 0; i < u.size(); ++i)
182    {
183        double dbl = (u[i] - mean);
184        double d2 = sqr(dbl);
185        var += d2;
186        skew += dbl * d2;
187        kurtosis += d2 * d2;
188    }
189    var /= u.size();
190    double dev = std::sqrt(var);
191    skew /= u.size() * dev * var;
192    kurtosis /= u.size() * var * var;
193    kurtosis -= 3;
194    double x_mean = (1 - d.p()) / d.p();
195    double x_var = x_mean / d.p();
196    double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
197    double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
198    assert(std::abs((mean - x_mean) / x_mean) < 0.01);
199    assert(std::abs((var - x_var) / x_var) < 0.01);
200    assert(std::abs((skew - x_skew) / x_skew) < 0.01);
201    assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
202}
203
204void
205test5()
206{
207    typedef std::geometric_distribution<> D;
208    typedef std::mt19937 G;
209    G g;
210    D d(0.75);
211    const int N = 1000000;
212    std::vector<D::result_type> u;
213    for (int i = 0; i < N; ++i)
214    {
215        D::result_type v = d(g);
216        assert(d.min() <= v && v <= d.max());
217        u.push_back(v);
218    }
219    double mean = std::accumulate(u.begin(), u.end(),
220                                          double(0)) / u.size();
221    double var = 0;
222    double skew = 0;
223    double kurtosis = 0;
224    for (unsigned i = 0; i < u.size(); ++i)
225    {
226        double dbl = (u[i] - mean);
227        double d2 = sqr(dbl);
228        var += d2;
229        skew += dbl * d2;
230        kurtosis += d2 * d2;
231    }
232    var /= u.size();
233    double dev = std::sqrt(var);
234    skew /= u.size() * dev * var;
235    kurtosis /= u.size() * var * var;
236    kurtosis -= 3;
237    double x_mean = (1 - d.p()) / d.p();
238    double x_var = x_mean / d.p();
239    double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
240    double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
241    assert(std::abs((mean - x_mean) / x_mean) < 0.01);
242    assert(std::abs((var - x_var) / x_var) < 0.01);
243    assert(std::abs((skew - x_skew) / x_skew) < 0.01);
244    assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
245}
246
247void
248test6()
249{
250    typedef std::geometric_distribution<> D;
251    typedef std::mt19937 G;
252    G g;
253    D d(0.96875);
254    const int N = 1000000;
255    std::vector<D::result_type> u;
256    for (int i = 0; i < N; ++i)
257    {
258        D::result_type v = d(g);
259        assert(d.min() <= v && v <= d.max());
260        u.push_back(v);
261    }
262    double mean = std::accumulate(u.begin(), u.end(),
263                                          double(0)) / u.size();
264    double var = 0;
265    double skew = 0;
266    double kurtosis = 0;
267    for (unsigned i = 0; i < u.size(); ++i)
268    {
269        double dbl = (u[i] - mean);
270        double d2 = sqr(dbl);
271        var += d2;
272        skew += dbl * d2;
273        kurtosis += d2 * d2;
274    }
275    var /= u.size();
276    double dev = std::sqrt(var);
277    skew /= u.size() * dev * var;
278    kurtosis /= u.size() * var * var;
279    kurtosis -= 3;
280    double x_mean = (1 - d.p()) / d.p();
281    double x_var = x_mean / d.p();
282    double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
283    double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
284    assert(std::abs((mean - x_mean) / x_mean) < 0.01);
285    assert(std::abs((var - x_var) / x_var) < 0.01);
286    assert(std::abs((skew - x_skew) / x_skew) < 0.01);
287    assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
288}
289
290int main()
291{
292    test1();
293    test2();
294    test3();
295    test4();
296    test5();
297    test6();
298}
299