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

source code of libcxx/test/std/numerics/rand/rand.dist/rand.dist.bern/rand.dist.bern.negbin/eval.pass.cpp