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       wl = { 1.0 };
 117
 118     if (k < 0 || (std::size_t)k >= wl.size())
 119       return 0.0;
 120     else
 121       {
 122         double sum = 0.0;
 123         for (auto it = wl.begin(); it != wl.end(); ++it)
 124           sum += *it;
 125         return wl.begin()[k] / sum;
 126       }
 127   }
 128
 129   inline double
 130   geometric_pdf(int k, double p)
 131   {
 132     if (k < 0)
 133       return 0.0;
 134     else if (k == 0)
 135       return p;
 136     else
 137       return p * std::pow(1 - p, k);
 138   }
 139
 140 #ifdef _GLIBCXX_USE_C99_MATH_TR1
 141   inline double
 142   negative_binomial_pdf(int k, int n, double p)
 143   {
 144     if (k < 0)
 145       return 0.0;
 146     else
 147       {
 148         double f = std::lgamma(k + (double)n);
 149         double a = std::lgamma(n);
 150         double b = std::lgamma(k + 1.0);
 151
 152         return std::exp(f - a - b) * std::pow(p, n) * std::pow(1 - p, k);
 153       }
 154   }
 155
 156   inline double
 157   poisson_pdf(int k, double mu)
 158   {
 159     if (k < 0)
 160       return 0.0;
 161     else
 162       {
 163         double lf = std::lgamma(k + 1.0);
 164         return std::exp(std::log(mu) * k - lf - mu);
 165       }
 166   }
 167 #endif
 168
 169   inline double
 170   uniform_int_pdf(int k, int a, int b)
 171   {
 172     if (k < 0 || k < a || k > b)
 173       return 0.0;
 174     else
 175       return 1.0 / (b - a + 1.0);
 176   }
 177
 178 #ifdef _GLIBCXX_USE_C99_MATH_TR1
 179   inline double
 180   lbincoef(int n, int k)
 181   {
 182     return std::lgamma(double(1 + n))
 183          - std::lgamma(double(1 + k))
 184          - std::lgamma(double(1 + n - k));
 185   }
 186
 187   inline double
 188   hypergeometric_pdf(int k, int N, int K, int n)
 189   {
 190     if (k < 0 || k < std::max(0, n - (N - K)) || k > std::min(K, n))
 191       return 0.0;
 192     else
 193       return lbincoef(K, k) + lbincoef(N - K, n - k) - lbincoef(N, n);
 194   }
 195 #endif
 196
 197 } // namespace __gnu_test
 198
 199 #endif // #ifndef _GLIBCXX_TESTSUITE_RANDOM_H