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