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