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
32int main()
33{
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 (int i = 0; i < u.size(); ++i)
53        {
54            double d = (u[i] - mean);
55            double d2 = sqr(d);
56            var += d2;
57            skew += d * 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    {
75        typedef std::geometric_distribution<> D;
76        typedef std::mt19937 G;
77        G g;
78        D d(0.05);
79        const int N = 1000000;
80        std::vector<D::result_type> u;
81        for (int i = 0; i < N; ++i)
82        {
83            D::result_type v = d(g);
84            assert(d.min() <= v && v <= d.max());
85            u.push_back(v);
86        }
87        double mean = std::accumulate(u.begin(), u.end(),
88                                              double(0)) / u.size();
89        double var = 0;
90        double skew = 0;
91        double kurtosis = 0;
92        for (int i = 0; i < u.size(); ++i)
93        {
94            double d = (u[i] - mean);
95            double d2 = sqr(d);
96            var += d2;
97            skew += d * d2;
98            kurtosis += d2 * d2;
99        }
100        var /= u.size();
101        double dev = std::sqrt(var);
102        skew /= u.size() * dev * var;
103        kurtosis /= u.size() * var * var;
104        kurtosis -= 3;
105        double x_mean = (1 - d.p()) / d.p();
106        double x_var = x_mean / d.p();
107        double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
108        double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
109        assert(std::abs((mean - x_mean) / x_mean) < 0.01);
110        assert(std::abs((var - x_var) / x_var) < 0.01);
111        assert(std::abs((skew - x_skew) / x_skew) < 0.01);
112        assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
113    }
114    {
115        typedef std::geometric_distribution<> D;
116        typedef std::minstd_rand G;
117        G g;
118        D d(.25);
119        const int N = 1000000;
120        std::vector<D::result_type> u;
121        for (int i = 0; i < N; ++i)
122        {
123            D::result_type v = d(g);
124            assert(d.min() <= v && v <= d.max());
125            u.push_back(v);
126        }
127        double mean = std::accumulate(u.begin(), u.end(),
128                                              double(0)) / u.size();
129        double var = 0;
130        double skew = 0;
131        double kurtosis = 0;
132        for (int i = 0; i < u.size(); ++i)
133        {
134            double d = (u[i] - mean);
135            double d2 = sqr(d);
136            var += d2;
137            skew += d * d2;
138            kurtosis += d2 * d2;
139        }
140        var /= u.size();
141        double dev = std::sqrt(var);
142        skew /= u.size() * dev * var;
143        kurtosis /= u.size() * var * var;
144        kurtosis -= 3;
145        double x_mean = (1 - d.p()) / d.p();
146        double x_var = x_mean / d.p();
147        double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
148        double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
149        assert(std::abs((mean - x_mean) / x_mean) < 0.01);
150        assert(std::abs((var - x_var) / x_var) < 0.01);
151        assert(std::abs((skew - x_skew) / x_skew) < 0.01);
152        assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
153    }
154    {
155        typedef std::geometric_distribution<> D;
156        typedef std::mt19937 G;
157        G g;
158        D d(0.5);
159        const int N = 1000000;
160        std::vector<D::result_type> u;
161        for (int i = 0; i < N; ++i)
162        {
163            D::result_type v = d(g);
164            assert(d.min() <= v && v <= d.max());
165            u.push_back(v);
166        }
167        double mean = std::accumulate(u.begin(), u.end(),
168                                              double(0)) / u.size();
169        double var = 0;
170        double skew = 0;
171        double kurtosis = 0;
172        for (int i = 0; i < u.size(); ++i)
173        {
174            double d = (u[i] - mean);
175            double d2 = sqr(d);
176            var += d2;
177            skew += d * d2;
178            kurtosis += d2 * d2;
179        }
180        var /= u.size();
181        double dev = std::sqrt(var);
182        skew /= u.size() * dev * var;
183        kurtosis /= u.size() * var * var;
184        kurtosis -= 3;
185        double x_mean = (1 - d.p()) / d.p();
186        double x_var = x_mean / d.p();
187        double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
188        double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
189        assert(std::abs((mean - x_mean) / x_mean) < 0.01);
190        assert(std::abs((var - x_var) / x_var) < 0.01);
191        assert(std::abs((skew - x_skew) / x_skew) < 0.01);
192        assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
193    }
194    {
195        typedef std::geometric_distribution<> D;
196        typedef std::mt19937 G;
197        G g;
198        D d(0.75);
199        const int N = 1000000;
200        std::vector<D::result_type> u;
201        for (int i = 0; i < N; ++i)
202        {
203            D::result_type v = d(g);
204            assert(d.min() <= v && v <= d.max());
205            u.push_back(v);
206        }
207        double mean = std::accumulate(u.begin(), u.end(),
208                                              double(0)) / u.size();
209        double var = 0;
210        double skew = 0;
211        double kurtosis = 0;
212        for (int i = 0; i < u.size(); ++i)
213        {
214            double d = (u[i] - mean);
215            double d2 = sqr(d);
216            var += d2;
217            skew += d * d2;
218            kurtosis += d2 * d2;
219        }
220        var /= u.size();
221        double dev = std::sqrt(var);
222        skew /= u.size() * dev * var;
223        kurtosis /= u.size() * var * var;
224        kurtosis -= 3;
225        double x_mean = (1 - d.p()) / d.p();
226        double x_var = x_mean / d.p();
227        double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
228        double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
229        assert(std::abs((mean - x_mean) / x_mean) < 0.01);
230        assert(std::abs((var - x_var) / x_var) < 0.01);
231        assert(std::abs((skew - x_skew) / x_skew) < 0.01);
232        assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
233    }
234    {
235        typedef std::geometric_distribution<> D;
236        typedef std::mt19937 G;
237        G g;
238        D d(0.96875);
239        const int N = 1000000;
240        std::vector<D::result_type> u;
241        for (int i = 0; i < N; ++i)
242        {
243            D::result_type v = d(g);
244            assert(d.min() <= v && v <= d.max());
245            u.push_back(v);
246        }
247        double mean = std::accumulate(u.begin(), u.end(),
248                                              double(0)) / u.size();
249        double var = 0;
250        double skew = 0;
251        double kurtosis = 0;
252        for (int i = 0; i < u.size(); ++i)
253        {
254            double d = (u[i] - mean);
255            double d2 = sqr(d);
256            var += d2;
257            skew += d * d2;
258            kurtosis += d2 * d2;
259        }
260        var /= u.size();
261        double dev = std::sqrt(var);
262        skew /= u.size() * dev * var;
263        kurtosis /= u.size() * var * var;
264        kurtosis -= 3;
265        double x_mean = (1 - d.p()) / d.p();
266        double x_var = x_mean / d.p();
267        double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
268        double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
269        assert(std::abs((mean - x_mean) / x_mean) < 0.01);
270        assert(std::abs((var - x_var) / x_var) < 0.01);
271        assert(std::abs((skew - x_skew) / x_skew) < 0.01);
272        assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
273    }
274}
275