libstdc++-v3/testsuite/util/testsuite_random.h

   1 // -*- C++ -*-
   2
   3 // Copyright (C) 2011-2018 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
  31 namespace __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       double count[BINS], p[BINS];
  43
  44       for (unsigned long i = 0; i < BINS; i++)
  45         count[i] = 0;
  46
  47       for (unsigned long i = 0; i < N; i++)
  48         {
  49           auto r = f();
  50           if (r >= 0 && (unsigned long)r < BINS)
  51             count[r]++;
  52         }
  53
  54       for (unsigned long i = 0; i < BINS; i++)
  55         p[i] = pdf(i);
  56
  57       for (unsigned long i = 0; i < BINS; i++)
  58         {
  59           bool status_i;
  60           double d = std::abs(count[i] - N * p[i]);
  61
  62           if (p[i] != 0)
  63             {
  64               double s = d / std::sqrt(N * p[i]);
  65               status_i = (s > 5) && (d > 1);
  66             }
  67           else
  68             status_i = (count[i] != 0);
  69
  70           VERIFY( !status_i );
  71         }
  72     }
  73
  74   inline double
  75   bernoulli_pdf(int k, double p)
  76   {
  77     if (k == 0)
  78       return 1 - p;
  79     else if (k == 1)
  80       return p;
  81     else
  82       return 0.0;
  83   }
  84
  85 #ifdef _GLIBCXX_USE_C99_MATH_TR1
  86   inline double
  87   binomial_pdf(int k, int n, double p)
  88   {
  89     if (k < 0 || k > n)
  90       return 0.0;
  91     else
  92       {
  93         double q;
  94
  95         if (p == 0.0)
  96           q = (k == 0) ? 1.0 : 0.0;
  97         else if (p == 1.0)
  98           q = (k == n) ? 1.0 : 0.0;
  99         else
 100           {
 101             double ln_Cnk = (std::lgamma(n + 1.0) - std::lgamma(k + 1.0)
 102                              - std::lgamma(n - k + 1.0));
 103             q = ln_Cnk + k * std::log(p) + (n - k) * std::log1p(-p);
 104             q = std::exp(q);
 105           }
 106
 107         return q;
 108       }
 109   }
 110 #endif
 111
 112   inline double
 113   discrete_pdf(int k, std::initializer_list<double> wl)
 114   {
 115     if (!wl.size())
 116       {
 117         static std::initializer_list<double> one = { 1.0 };
 118         wl = one;
 119       }
 120
 121     if (k < 0 || (std::size_t)k >= wl.size())
 122       return 0.0;
 123     else
 124       {
 125         double sum = 0.0;
 126         for (auto it = wl.begin(); it != wl.end(); ++it)
 127           sum += *it;
 128         return wl.begin()[k] / sum;
 129       }
 130   }
 131
 132   inline double
 133   geometric_pdf(int k, double p)
 134   {
 135     if (k < 0)
 136       return 0.0;
 137     else if (k == 0)
 138       return p;
 139     else
 140       return p * std::pow(1 - p, k);
 141   }
 142
 143 #ifdef _GLIBCXX_USE_C99_MATH_TR1
 144   inline double
 145   negative_binomial_pdf(int k, int n, double p)
 146   {
 147     if (k < 0)
 148       return 0.0;
 149     else
 150       {
 151         double f = std::lgamma(k + (double)n);
 152         double a = std::lgamma(n);
 153         double b = std::lgamma(k + 1.0);
 154
 155         return std::exp(f - a - b) * std::pow(p, n) * std::pow(1 - p, k);
 156       }
 157   }
 158
 159   inline double
 160   poisson_pdf(int k, double mu)
 161   {
 162     if (k < 0)
 163       return 0.0;
 164     else
 165       {
 166         double lf = std::lgamma(k + 1.0);
 167         return std::exp(std::log(mu) * k - lf - mu);
 168       }
 169   }
 170 #endif
 171
 172   inline double
 173   uniform_int_pdf(int k, int a, int b)
 174   {
 175     if (k < 0 || k < a || k > b)
 176       return 0.0;
 177     else
 178       return 1.0 / (b - a + 1.0);
 179   }
 180
 181 #ifdef _GLIBCXX_USE_C99_MATH_TR1
 182   inline double
 183   lbincoef(int n, int k)
 184   {
 185     return std::lgamma(double(1 + n))
 186          - std::lgamma(double(1 + k))
 187          - std::lgamma(double(1 + n - k));
 188   }
 189
 190   inline double
 191   hypergeometric_pdf(int k, int N, int K, int n)
 192   {
 193     if (k < 0 || k < std::max(0, n - (N - K)) || k > std::min(K, n))
 194       return 0.0;
 195     else
 196       return lbincoef(K, k) + lbincoef(N - K, n - k) - lbincoef(N, n);
 197   }
 198 #endif
 199
 200 } // namespace __gnu_test
 201
 202 #endif // #ifndef _GLIBCXX_TESTSUITE_RANDOM_H