xref: /llvm-project/libcxx/test/std/numerics/rand/rand.dist/rand.dist.bern/rand.dist.bern.geo/eval.pass.cpp (revision 09e3a360581dc36d0820d3fb6da9bd7cfed87b5d) 
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 geometric_distribution
15 
16 // template<class _URNG> result_type operator()(_URNG& g);
17 
18 #include <random>
19 #include <cassert>
20 #include <cmath>
21 #include <cstdint>
22 #include <numeric>
23 #include <vector>
24 
25 #include "test_macros.h"
26 
27 template <class T>
28 T sqr(T x) {
29     return x * x;
30 }
31 
32 void test_small_inputs() {
33   std::mt19937 engine;
34   std::geometric_distribution<std::int16_t> distribution(5.45361e-311);
35   typedef std::geometric_distribution<std::int16_t>::result_type result_type;
36   for (int i = 0; i < 1000; ++i) {
37     volatile result_type res = distribution(engine);
38     ((void)res);
39   }
40 }
41 
42 template <class T>
43 void test1() {
44     typedef std::geometric_distribution<T> D;
45     typedef std::mt19937 G;
46     G g;
47     D d(.03125);
48     const int N = 1000000;
49     std::vector<typename D::result_type> u;
50     for (int i = 0; i < N; ++i)
51     {
52         typename D::result_type v = d(g);
53         assert(d.min() <= v && v <= d.max());
54         u.push_back(v);
55     }
56     double mean = std::accumulate(u.begin(), u.end(),
57                                           double(0)) / u.size();
58     double var = 0;
59     double skew = 0;
60     double kurtosis = 0;
61     for (unsigned i = 0; i < u.size(); ++i)
62     {
63         double dbl = (u[i] - mean);
64         double d2 = sqr(dbl);
65         var += d2;
66         skew += dbl * d2;
67         kurtosis += d2 * d2;
68     }
69     var /= u.size();
70     double dev = std::sqrt(var);
71     skew /= u.size() * dev * var;
72     kurtosis /= u.size() * var * var;
73     kurtosis -= 3;
74     double x_mean = (1 - d.p()) / d.p();
75     double x_var = x_mean / d.p();
76     double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
77     double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
78     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
79     assert(std::abs((var - x_var) / x_var) < 0.01);
80     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
81     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
82 }
83 
84 template <class T>
85 void test2() {
86     typedef std::geometric_distribution<T> D;
87     typedef std::mt19937 G;
88     G g;
89     D d(0.05);
90     const int N = 1000000;
91     std::vector<typename D::result_type> u;
92     for (int i = 0; i < N; ++i)
93     {
94         typename D::result_type v = d(g);
95         assert(d.min() <= v && v <= d.max());
96         u.push_back(v);
97     }
98     double mean = std::accumulate(u.begin(), u.end(),
99                                           double(0)) / u.size();
100     double var = 0;
101     double skew = 0;
102     double kurtosis = 0;
103     for (unsigned i = 0; i < u.size(); ++i)
104     {
105         double dbl = (u[i] - mean);
106         double d2 = sqr(dbl);
107         var += d2;
108         skew += dbl * d2;
109         kurtosis += d2 * d2;
110     }
111     var /= u.size();
112     double dev = std::sqrt(var);
113     skew /= u.size() * dev * var;
114     kurtosis /= u.size() * var * var;
115     kurtosis -= 3;
116     double x_mean = (1 - d.p()) / d.p();
117     double x_var = x_mean / d.p();
118     double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
119     double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
120     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
121     assert(std::abs((var - x_var) / x_var) < 0.01);
122     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
123     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
124 }
125 
126 template <class T>
127 void test3() {
128     typedef std::geometric_distribution<T> D;
129     typedef std::minstd_rand G;
130     G g;
131     D d(.25);
132     const int N = 1000000;
133     std::vector<typename D::result_type> u;
134     for (int i = 0; i < N; ++i)
135     {
136         typename D::result_type v = d(g);
137         assert(d.min() <= v && v <= d.max());
138         u.push_back(v);
139     }
140     double mean = std::accumulate(u.begin(), u.end(),
141                                           double(0)) / u.size();
142     double var = 0;
143     double skew = 0;
144     double kurtosis = 0;
145     for (unsigned i = 0; i < u.size(); ++i)
146     {
147         double dbl = (u[i] - mean);
148         double d2 = sqr(dbl);
149         var += d2;
150         skew += dbl * d2;
151         kurtosis += d2 * d2;
152     }
153     var /= u.size();
154     double dev = std::sqrt(var);
155     skew /= u.size() * dev * var;
156     kurtosis /= u.size() * var * var;
157     kurtosis -= 3;
158     double x_mean = (1 - d.p()) / d.p();
159     double x_var = x_mean / d.p();
160     double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
161     double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
162     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
163     assert(std::abs((var - x_var) / x_var) < 0.01);
164     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
165     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
166 }
167 
168 template <class T>
169 void test4() {
170     typedef std::geometric_distribution<T> D;
171     typedef std::mt19937 G;
172     G g;
173     D d(0.5);
174     const int N = 1000000;
175     std::vector<typename D::result_type> u;
176     for (int i = 0; i < N; ++i)
177     {
178         typename D::result_type v = d(g);
179         assert(d.min() <= v && v <= d.max());
180         u.push_back(v);
181     }
182     double mean = std::accumulate(u.begin(), u.end(),
183                                           double(0)) / u.size();
184     double var = 0;
185     double skew = 0;
186     double kurtosis = 0;
187     for (unsigned i = 0; i < u.size(); ++i)
188     {
189         double dbl = (u[i] - mean);
190         double d2 = sqr(dbl);
191         var += d2;
192         skew += dbl * d2;
193         kurtosis += d2 * d2;
194     }
195     var /= u.size();
196     double dev = std::sqrt(var);
197     skew /= u.size() * dev * var;
198     kurtosis /= u.size() * var * var;
199     kurtosis -= 3;
200     double x_mean = (1 - d.p()) / d.p();
201     double x_var = x_mean / d.p();
202     double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
203     double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
204     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
205     assert(std::abs((var - x_var) / x_var) < 0.01);
206     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
207     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
208 }
209 
210 template <class T>
211 void test5() {
212     typedef std::geometric_distribution<T> D;
213     typedef std::mt19937 G;
214     G g;
215     D d(0.75);
216     const int N = 1000000;
217     std::vector<typename D::result_type> u;
218     for (int i = 0; i < N; ++i)
219     {
220         typename D::result_type v = d(g);
221         assert(d.min() <= v && v <= d.max());
222         u.push_back(v);
223     }
224     double mean = std::accumulate(u.begin(), u.end(),
225                                           double(0)) / u.size();
226     double var = 0;
227     double skew = 0;
228     double kurtosis = 0;
229     for (unsigned i = 0; i < u.size(); ++i)
230     {
231         double dbl = (u[i] - mean);
232         double d2 = sqr(dbl);
233         var += d2;
234         skew += dbl * d2;
235         kurtosis += d2 * d2;
236     }
237     var /= u.size();
238     double dev = std::sqrt(var);
239     skew /= u.size() * dev * var;
240     kurtosis /= u.size() * var * var;
241     kurtosis -= 3;
242     double x_mean = (1 - d.p()) / d.p();
243     double x_var = x_mean / d.p();
244     double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
245     double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
246     assert(std::abs((mean - x_mean) / x_mean) < 0.01);
247     assert(std::abs((var - x_var) / x_var) < 0.01);
248     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
249     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.03);
250 }
251 
252 template <class T>
253 void test6() {
254     typedef std::geometric_distribution<T> D;
255     typedef std::mt19937 G;
256     G g;
257     D d(0.96875);
258     const int N = 1000000;
259     std::vector<typename D::result_type> u;
260     for (int i = 0; i < N; ++i)
261     {
262         typename D::result_type v = d(g);
263         assert(d.min() <= v && v <= d.max());
264         u.push_back(v);
265     }
266     double mean = std::accumulate(u.begin(), u.end(),
267                                           double(0)) / u.size();
268     double var = 0;
269     double skew = 0;
270     double kurtosis = 0;
271     for (unsigned i = 0; i < u.size(); ++i)
272     {
273         double dbl = (u[i] - mean);
274         double d2 = sqr(dbl);
275         var += d2;
276         skew += dbl * d2;
277         kurtosis += d2 * d2;
278     }
279     var /= u.size();
280     double dev = std::sqrt(var);
281     skew /= u.size() * dev * var;
282     kurtosis /= u.size() * var * var;
283     kurtosis -= 3;
284     double x_mean = (1 - d.p()) / d.p();
285     double x_var = x_mean / d.p();
286     double x_skew = (2 - d.p()) / std::sqrt((1 - d.p()));
287     double x_kurtosis = 6 + sqr(d.p()) / (1 - d.p());
288     assert(std::abs((mean - x_mean) / x_mean) < 0.02);
289     assert(std::abs((var - x_var) / x_var) < 0.02);
290     assert(std::abs((skew - x_skew) / x_skew) < 0.01);
291     assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.02);
292 }
293 
294 template <class T>
295 void tests() {
296     test1<T>();
297     test2<T>();
298     test3<T>();
299     test4<T>();
300     test5<T>();
301     test6<T>();
302 }
303 
304 int main(int, char**) {
305     test_small_inputs();
306 
307     tests<short>();
308     tests<int>();
309     tests<long>();
310     tests<long long>();
311 
312     tests<unsigned short>();
313     tests<unsigned int>();
314     tests<unsigned long>();
315     tests<unsigned long long>();
316 
317 #if defined(_LIBCPP_VERSION) // extension
318     // TODO: std::geometric_distribution currently doesn't work reliably with small types.
319     // tests<int8_t>();
320     // tests<uint8_t>();
321 #if !defined(TEST_HAS_NO_INT128)
322     tests<__int128_t>();
323     tests<__uint128_t>();
324 #endif
325 #endif
326 
327     return 0;
328 }
329