1// -*- C++ -*- 2 3// Copyright (C) 2011-2015 Free Software Foundation, Inc. 4// 5// This file is part of the GNU ISO C++ Library. This library is free 6// software; you can redistribute it and/or modify it under the terms 7// of the GNU General Public License as published by the Free Software 8// Foundation; either version 3, or (at your option) any later 9// version. 10 11// This library is distributed in the hope that it will be useful, but 12// WITHOUT ANY WARRANTY; without even the implied warranty of 13// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 14// General Public License for more details. 15 16// You should have received a copy of the GNU General Public License along 17// with this library; see the file COPYING3. If not see 18// <http://www.gnu.org/licenses/>. 19 20/** 21 * @file testsuite_random.h 22 */ 23 24#ifndef _GLIBCXX_TESTSUITE_RANDOM_H 25#define _GLIBCXX_TESTSUITE_RANDOM_H 26 27#include <cmath> 28#include <initializer_list> 29#include <testsuite_hooks.h> 30 31namespace __gnu_test 32{ 33 // Adapted for libstdc++ from GNU gsl-1.14/randist/test.c 34 // Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007, 2010 35 // James Theiler, Brian Gough 36 template<unsigned long BINS = 100, 37 unsigned long N = 100000, 38 typename Distribution, typename Pdf> 39 void 40 testDiscreteDist(Distribution& f, Pdf pdf) 41 { 42 bool test __attribute__((unused)) = true; 43 double count[BINS], p[BINS]; 44 45 for (unsigned long i = 0; i < BINS; i++) 46 count[i] = 0; 47 48 for (unsigned long i = 0; i < N; i++) 49 { 50 auto r = f(); 51 if (r >= 0 && (unsigned long)r < BINS) 52 count[r]++; 53 } 54 55 for (unsigned long i = 0; i < BINS; i++) 56 p[i] = pdf(i); 57 58 for (unsigned long i = 0; i < BINS; i++) 59 { 60 bool status_i; 61 double d = std::abs(count[i] - N * p[i]); 62 63 if (p[i] != 0) 64 { 65 double s = d / std::sqrt(N * p[i]); 66 status_i = (s > 5) && (d > 1); 67 } 68 else 69 status_i = (count[i] != 0); 70 71 VERIFY( !status_i ); 72 } 73 } 74 75 inline double 76 bernoulli_pdf(int k, double p) 77 { 78 if (k == 0) 79 return 1 - p; 80 else if (k == 1) 81 return p; 82 else 83 return 0.0; 84 } 85 86#ifdef _GLIBCXX_USE_C99_MATH_TR1 87 inline double 88 binomial_pdf(int k, int n, double p) 89 { 90 if (k < 0 || k > n) 91 return 0.0; 92 else 93 { 94 double q; 95 96 if (p == 0.0) 97 q = (k == 0) ? 1.0 : 0.0; 98 else if (p == 1.0) 99 q = (k == n) ? 1.0 : 0.0; 100 else 101 { 102 double ln_Cnk = (std::lgamma(n + 1.0) - std::lgamma(k + 1.0) 103 - std::lgamma(n - k + 1.0)); 104 q = ln_Cnk + k * std::log(p) + (n - k) * std::log1p(-p); 105 q = std::exp(q); 106 } 107 108 return q; 109 } 110 } 111#endif 112 113 inline double 114 discrete_pdf(int k, std::initializer_list<double> wl) 115 { 116 if (!wl.size()) 117 wl = { 1.0 }; 118 119 if (k < 0 || (std::size_t)k >= wl.size()) 120 return 0.0; 121 else 122 { 123 double sum = 0.0; 124 for (auto it = wl.begin(); it != wl.end(); ++it) 125 sum += *it; 126 return wl.begin()[k] / sum; 127 } 128 } 129 130 inline double 131 geometric_pdf(int k, double p) 132 { 133 if (k < 0) 134 return 0.0; 135 else if (k == 0) 136 return p; 137 else 138 return p * std::pow(1 - p, k); 139 } 140 141#ifdef _GLIBCXX_USE_C99_MATH_TR1 142 inline double 143 negative_binomial_pdf(int k, int n, double p) 144 { 145 if (k < 0) 146 return 0.0; 147 else 148 { 149 double f = std::lgamma(k + (double)n); 150 double a = std::lgamma(n); 151 double b = std::lgamma(k + 1.0); 152 153 return std::exp(f - a - b) * std::pow(p, n) * std::pow(1 - p, k); 154 } 155 } 156 157 inline double 158 poisson_pdf(int k, double mu) 159 { 160 if (k < 0) 161 return 0.0; 162 else 163 { 164 double lf = std::lgamma(k + 1.0); 165 return std::exp(std::log(mu) * k - lf - mu); 166 } 167 } 168#endif 169 170 inline double 171 uniform_int_pdf(int k, int a, int b) 172 { 173 if (k < 0 || k < a || k > b) 174 return 0.0; 175 else 176 return 1.0 / (b - a + 1.0); 177 } 178 179#ifdef _GLIBCXX_USE_C99_MATH_TR1 180 inline double 181 lbincoef(int n, int k) 182 { 183 return std::lgamma(double(1 + n)) 184 - std::lgamma(double(1 + k)) 185 - std::lgamma(double(1 + n - k)); 186 } 187 188 inline double 189 hypergeometric_pdf(int k, int N, int K, int n) 190 { 191 if (k < 0 || k < std::max(0, n - (N - K)) || k > std::min(K, n)) 192 return 0.0; 193 else 194 return lbincoef(K, k) + lbincoef(N - K, n - k) - lbincoef(N, n); 195 } 196#endif 197 198} // namespace __gnu_test 199 200#endif // #ifndef _GLIBCXX_TESTSUITE_RANDOM_H 201