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

   1 // -*- C++ -*-
   2
   3 // Copyright (C) 2011 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 <testsuite_hooks.h>
  29
  30 namespace __gnu_test
  31 {
  32   // Adapted for libstdc++ from GNU gsl-1.14/randist/test.c
  33   // Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007, 2010
  34   // James Theiler, Brian Gough
  35   template<unsigned long BINS = 100,
  36            unsigned long N = 100000,
  37            typename Distribution, typename Pdf>
  38     void
  39     testDiscreteDist(Distribution& f, Pdf pdf)
  40     {
  41       bool test __attribute__((unused)) = true;
  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 && 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;
  83   }
  84
  85 #ifdef _GLIBCXX_USE_C99_MATH_TR1
  86   inline double
  87   binomial_pdf(int k, double p, int n)
  88   {
  89     if (k < 0 || k > n)
  90       return 0;
  91     else
  92       {
  93         double q;
  94
  95         if (p == 0)
  96           q = (k == 0) ? 1 : 0;
  97         else if (p == 1)
  98           q = (k == n) ? 1 : 0;
  99         else
 100           {
 101             double ln_Cnk = (std::lgamma(n + 1) - std::lgamma(k + 1)
 102                              - std::lgamma(n - k + 1));
 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   geometric_pdf(int k, double p)
 114   {
 115     if (k < 0)
 116       return 0;
 117     else if (k == 0)
 118       return p;
 119     else
 120       return p * std::pow(1 - p, k);
 121   }
 122 } // namespace __gnu_test
 123
 124 #endif // #ifndef _GLIBCXX_TESTSUITE_RANDOM_H