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

   1 // -*- C++ -*-
   2
   3 // Copyright (C) 2011-2025 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 <system_error>
  30 #include <testsuite_hooks.h>
  31
  32 namespace __gnu_test
  33 {
  34   // Adapted for libstdc++ from GNU gsl-1.14/randist/test.c
  35   // Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007, 2010
  36   // James Theiler, Brian Gough
  37   template<unsigned long BINS = 100,
  38            unsigned long N = 100000,
  39            typename Distribution, typename Pdf>
  40     void
  41     testDiscreteDist(Distribution& f, Pdf pdf)
  42     {
  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_FUNCS
  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       {
 118         static std::initializer_list<double> one = { 1.0 };
 119         wl = one;
 120       }
 121
 122     if (k < 0 || (std::size_t)k >= wl.size())
 123       return 0.0;
 124     else
 125       {
 126         double sum = 0.0;
 127         for (auto it = wl.begin(); it != wl.end(); ++it)
 128           sum += *it;
 129         return wl.begin()[k] / sum;
 130       }
 131   }
 132
 133   inline double
 134   geometric_pdf(int k, double p)
 135   {
 136     if (k < 0)
 137       return 0.0;
 138     else if (k == 0)
 139       return p;
 140     else
 141       return p * std::pow(1 - p, k);
 142   }
 143
 144 #ifdef _GLIBCXX_USE_C99_MATH_FUNCS
 145   inline double
 146   negative_binomial_pdf(int k, int n, double p)
 147   {
 148     if (k < 0)
 149       return 0.0;
 150     else
 151       {
 152         double f = std::lgamma(k + (double)n);
 153         double a = std::lgamma(n);
 154         double b = std::lgamma(k + 1.0);
 155
 156         return std::exp(f - a - b) * std::pow(p, n) * std::pow(1 - p, k);
 157       }
 158   }
 159
 160   inline double
 161   poisson_pdf(int k, double mu)
 162   {
 163     if (k < 0)
 164       return 0.0;
 165     else
 166       {
 167         double lf = std::lgamma(k + 1.0);
 168         return std::exp(std::log(mu) * k - lf - mu);
 169       }
 170   }
 171 #endif
 172
 173   inline double
 174   uniform_int_pdf(int k, int a, int b)
 175   {
 176     if (k < 0 || k < a || k > b)
 177       return 0.0;
 178     else
 179       return 1.0 / (b - a + 1.0);
 180   }
 181
 182 #ifdef _GLIBCXX_USE_C99_MATH_FUNCS
 183   inline double
 184   lbincoef(int n, int k)
 185   {
 186     return std::lgamma(double(1 + n))
 187          - std::lgamma(double(1 + k))
 188          - std::lgamma(double(1 + n - k));
 189   }
 190
 191   inline double
 192   hypergeometric_pdf(int k, int N, int K, int n)
 193   {
 194     if (k < 0 || k < std::max(0, n - (N - K)) || k > std::min(K, n))
 195       return 0.0;
 196     else
 197       return lbincoef(K, k) + lbincoef(N - K, n - k) - lbincoef(N, n);
 198   }
 199 #endif
 200
 201   // Check whether TOKEN can construct a std::random_device successfully.
 202   inline bool
 203   random_device_available(const std::string& token) noexcept
 204   {
 205     try {
 206       std::random_device dev(token);
 207       return true;
 208     } catch (const std::system_error& /* See PR libstdc++/105081 */) {
 209       return false;
 210     }
 211   }
 212
 213 } // namespace __gnu_test
 214
 215 #endif // #ifndef _GLIBCXX_TESTSUITE_RANDOM_H