240 lines · cpp
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-exception6//7//===----------------------------------------------------------------------===//8//9// REQUIRES: long_tests10 11// <random>12 13// template<class IntType = int>14// class poisson_distribution15 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 28template <class T>29T sqr(T x) {30 return x * x;31}32 33void test_bad_ranges() {34 // Test cases where the mean is around the largest representable integer for35 // `result_type`. These cases don't generate valid poisson distributions, but36 // 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 94template <class T>95void 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 215int 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) // extension229 // 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#endif236#endif237 238 return 0;239}240