3 // Copyright (C) 2011 Free Software Foundation, Inc.
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
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.
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/>.
21 * @file testsuite_random.h
24 #ifndef _GLIBCXX_TESTSUITE_RANDOM_H
25 #define _GLIBCXX_TESTSUITE_RANDOM_H
28 #include <initializer_list>
29 #include <testsuite_hooks.h>
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>
40 testDiscreteDist(Distribution& f, Pdf pdf)
42 bool test __attribute__((unused)) = true;
43 double count[BINS], p[BINS];
45 for (unsigned long i = 0; i < BINS; i++)
48 for (unsigned long i = 0; i < N; i++)
51 if (r >= 0 && (unsigned long)r < BINS)
55 for (unsigned long i = 0; i < BINS; i++)
58 for (unsigned long i = 0; i < BINS; i++)
61 double d = std::abs(count[i] - N * p[i]);
65 double s = d / std::sqrt(N * p[i]);
66 status_i = (s > 5) && (d > 1);
69 status_i = (count[i] != 0);
76 bernoulli_pdf(int k, double p)
86 #ifdef _GLIBCXX_USE_C99_MATH_TR1
88 binomial_pdf(int k, int n, double p)
97 q = (k == 0) ? 1.0 : 0.0;
99 q = (k == n) ? 1.0 : 0.0;
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);
114 discrete_pdf(int k, std::initializer_list<double> wl)
119 if (k < 0 || (std::size_t)k >= wl.size())
124 for (auto it = wl.begin(); it != wl.end(); ++it)
126 return wl.begin()[k] / sum;
131 geometric_pdf(int k, double p)
138 return p * std::pow(1 - p, k);
141 #ifdef _GLIBCXX_USE_C99_MATH_TR1
143 negative_binomial_pdf(int k, int n, double p)
149 double f = std::lgamma(k + (double)n);
150 double a = std::lgamma(n);
151 double b = std::lgamma(k + 1.0);
153 return std::exp(f - a - b) * std::pow(p, n) * std::pow(1 - p, k);
158 poisson_pdf(int k, double mu)
164 double lf = std::lgamma(k + 1.0);
165 return std::exp(std::log(mu) * k - lf - mu);
171 uniform_int_pdf(int k, int a, int b)
173 if (k < 0 || k < a || k > b)
176 return 1.0 / (b - a + 1.0);
179 } // namespace __gnu_test
181 #endif // #ifndef _GLIBCXX_TESTSUITE_RANDOM_H