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 poisson_distribution
15
16 // template<class _URNG> result_type operator()(_URNG& g);
17
18 #include <cassert>
19 #include <cmath>
20 #include <cstdint>
21 #include <limits>
22 #include <numeric>
23 #include <random>
24 #include <vector>
25
26 #include "test_macros.h"
27
28 template <class T>
sqr(T x)29 T sqr(T x) {
30 return x * x;
31 }
32
test_bad_ranges()33 void test_bad_ranges() {
34 // Test cases where the mean is around the largest representable integer for
35 // `result_type`. These cases don't generate valid poisson distributions, but
36 // at least they don't blow up.
37 std::mt19937 eng;
38
39 {
40 std::poisson_distribution<std::int16_t> distribution(32710.9);
41 for (int i=0; i < 1000; ++i) {
42 volatile std::int16_t res = distribution(eng);
43 ((void)res);
44 }
45 }
46 {
47 std::poisson_distribution<std::int16_t> distribution(std::numeric_limits<std::int16_t>::max());
48 for (int i=0; i < 1000; ++i) {
49 volatile std::int16_t res = distribution(eng);
50 ((void)res);
51 }
52 }
53 {
54 std::poisson_distribution<std::int16_t> distribution(
55 static_cast<double>(std::numeric_limits<std::int16_t>::max()) + 10);
56 for (int i=0; i < 1000; ++i) {
57 volatile std::int16_t res = distribution(eng);
58 ((void)res);
59 }
60 }
61 {
62 std::poisson_distribution<std::int16_t> distribution(
63 static_cast<double>(std::numeric_limits<std::int16_t>::max()) * 2);
64 for (int i=0; i < 1000; ++i) {
65 volatile std::int16_t res = distribution(eng);
66 ((void)res);
67 }
68 }
69 {
70 // We convert `INF` to `DBL_MAX` otherwise the distribution will hang.
71 std::poisson_distribution<std::int16_t> distribution(std::numeric_limits<double>::infinity());
72 for (int i=0; i < 1000; ++i) {
73 volatile std::int16_t res = distribution(eng);
74 ((void)res);
75 }
76 }
77 {
78 std::poisson_distribution<std::int16_t> distribution(0);
79 for (int i=0; i < 1000; ++i) {
80 volatile std::int16_t res = distribution(eng);
81 ((void)res);
82 }
83 }
84 {
85 // We convert `INF` to `DBL_MAX` otherwise the distribution will hang.
86 std::poisson_distribution<std::int16_t> distribution(-100);
87 for (int i=0; i < 1000; ++i) {
88 volatile std::int16_t res = distribution(eng);
89 ((void)res);
90 }
91 }
92 }
93
94 template <class T>
tests()95 void tests() {
96 {
97 typedef std::poisson_distribution<T> D;
98 typedef std::minstd_rand G;
99 G g;
100 D d(2);
101 const int N = 100000;
102 std::vector<double> u;
103 for (int i = 0; i < N; ++i)
104 {
105 typename D::result_type v = d(g);
106 assert(d.min() <= v && v <= d.max());
107 u.push_back(v);
108 }
109 double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
110 double var = 0;
111 double skew = 0;
112 double kurtosis = 0;
113 for (unsigned i = 0; i < u.size(); ++i)
114 {
115 double dbl = (u[i] - mean);
116 double d2 = sqr(dbl);
117 var += d2;
118 skew += dbl * d2;
119 kurtosis += d2 * d2;
120 }
121 var /= u.size();
122 double dev = std::sqrt(var);
123 skew /= u.size() * dev * var;
124 kurtosis /= u.size() * var * var;
125 kurtosis -= 3;
126 double x_mean = d.mean();
127 double x_var = d.mean();
128 double x_skew = 1 / std::sqrt(x_var);
129 double x_kurtosis = 1 / x_var;
130 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
131 assert(std::abs((var - x_var) / x_var) < 0.01);
132 assert(std::abs((skew - x_skew) / x_skew) < 0.03);
133 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.2);
134 }
135 {
136 typedef std::poisson_distribution<T> D;
137 typedef std::minstd_rand G;
138 G g;
139 D d(0.75);
140 const int N = 100000;
141 std::vector<double> u;
142 for (int i = 0; i < N; ++i)
143 {
144 typename D::result_type v = d(g);
145 assert(d.min() <= v && v <= d.max());
146 u.push_back(v);
147 }
148 double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
149 double var = 0;
150 double skew = 0;
151 double kurtosis = 0;
152 for (unsigned i = 0; i < u.size(); ++i)
153 {
154 double dbl = (u[i] - mean);
155 double d2 = sqr(dbl);
156 var += d2;
157 skew += dbl * d2;
158 kurtosis += d2 * d2;
159 }
160 var /= u.size();
161 double dev = std::sqrt(var);
162 skew /= u.size() * dev * var;
163 kurtosis /= u.size() * var * var;
164 kurtosis -= 3;
165 double x_mean = d.mean();
166 double x_var = d.mean();
167 double x_skew = 1 / std::sqrt(x_var);
168 double x_kurtosis = 1 / x_var;
169 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
170 assert(std::abs((var - x_var) / x_var) < 0.02);
171 assert(std::abs((skew - x_skew) / x_skew) < 0.02);
172 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.09);
173 }
174 {
175 typedef std::poisson_distribution<T> D;
176 typedef std::mt19937 G;
177 G g;
178 D d(20);
179 const int N = 1000000;
180 std::vector<double> u;
181 for (int i = 0; i < N; ++i)
182 {
183 typename D::result_type v = d(g);
184 assert(d.min() <= v && v <= d.max());
185 u.push_back(v);
186 }
187 double mean = std::accumulate(u.begin(), u.end(), 0.0) / u.size();
188 double var = 0;
189 double skew = 0;
190 double kurtosis = 0;
191 for (unsigned i = 0; i < u.size(); ++i)
192 {
193 double dbl = (u[i] - mean);
194 double d2 = sqr(dbl);
195 var += d2;
196 skew += dbl * d2;
197 kurtosis += d2 * d2;
198 }
199 var /= u.size();
200 double dev = std::sqrt(var);
201 skew /= u.size() * dev * var;
202 kurtosis /= u.size() * var * var;
203 kurtosis -= 3;
204 double x_mean = d.mean();
205 double x_var = d.mean();
206 double x_skew = 1 / std::sqrt(x_var);
207 double x_kurtosis = 1 / x_var;
208 assert(std::abs((mean - x_mean) / x_mean) < 0.01);
209 assert(std::abs((var - x_var) / x_var) < 0.01);
210 assert(std::abs((skew - x_skew) / x_skew) < 0.02);
211 assert(std::abs((kurtosis - x_kurtosis) / x_kurtosis) < 0.3);
212 }
213 }
214
main(int,char **)215 int main(int, char**) {
216 test_bad_ranges();
217
218 tests<short>();
219 tests<int>();
220 tests<long>();
221 tests<long long>();
222
223 tests<unsigned short>();
224 tests<unsigned int>();
225 tests<unsigned long>();
226 tests<unsigned long long>();
227
228 #if defined(_LIBCPP_VERSION) // extension
229 // TODO: std::poisson_distribution currently doesn't work reliably with small types.
230 // tests<int8_t>();
231 // tests<uint8_t>();
232 #if !defined(TEST_HAS_NO_INT128)
233 tests<__int128_t>();
234 tests<__uint128_t>();
235 #endif
236 #endif
237
238 return 0;
239 }
240