libcxx/include/__random/binomial_distribution.h

   1 //===----------------------------------------------------------------------===//
   2 //
   3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
   4 // See https://llvm.org/LICENSE.txt for license information.
   5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
   6 //
   7 //===----------------------------------------------------------------------===//
   8
   9 #ifndef _LIBCPP___RANDOM_BINOMIAL_DISTRIBUTION_H
  10 #define _LIBCPP___RANDOM_BINOMIAL_DISTRIBUTION_H
  11
  12 #include <__config>
  13 #include <__random/uniform_real_distribution.h>
  14 #include <cmath>
  15 #include <iosfwd>
  16
  17 #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
  18 #pragma GCC system_header
  19 #endif
  20
  21 _LIBCPP_PUSH_MACROS
  22 #include <__undef_macros>
  23
  24 _LIBCPP_BEGIN_NAMESPACE_STD
  25
  26 template<class _IntType = int>
  27 class _LIBCPP_TEMPLATE_VIS binomial_distribution
  28 {
  29 public:
  30     // types
  31     typedef _IntType result_type;
  32
  33     class _LIBCPP_TEMPLATE_VIS param_type
  34     {
  35         result_type __t_;
  36         double __p_;
  37         double __pr_;
  38         double __odds_ratio_;
  39         result_type __r0_;
  40     public:
  41         typedef binomial_distribution distribution_type;
  42
  43         explicit param_type(result_type __t = 1, double __p = 0.5);
  44
  45         _LIBCPP_INLINE_VISIBILITY
  46         result_type t() const {return __t_;}
  47         _LIBCPP_INLINE_VISIBILITY
  48         double p() const {return __p_;}
  49
  50         friend _LIBCPP_INLINE_VISIBILITY
  51             bool operator==(const param_type& __x, const param_type& __y)
  52             {return __x.__t_ == __y.__t_ && __x.__p_ == __y.__p_;}
  53         friend _LIBCPP_INLINE_VISIBILITY
  54             bool operator!=(const param_type& __x, const param_type& __y)
  55             {return !(__x == __y);}
  56
  57         friend class binomial_distribution;
  58     };
  59
  60 private:
  61     param_type __p_;
  62
  63 public:
  64     // constructors and reset functions
  65 #ifndef _LIBCPP_CXX03_LANG
  66     _LIBCPP_INLINE_VISIBILITY
  67     binomial_distribution() : binomial_distribution(1) {}
  68     _LIBCPP_INLINE_VISIBILITY
  69     explicit binomial_distribution(result_type __t, double __p = 0.5)
  70         : __p_(param_type(__t, __p)) {}
  71 #else
  72     _LIBCPP_INLINE_VISIBILITY
  73     explicit binomial_distribution(result_type __t = 1, double __p = 0.5)
  74         : __p_(param_type(__t, __p)) {}
  75 #endif
  76     _LIBCPP_INLINE_VISIBILITY
  77     explicit binomial_distribution(const param_type& __p) : __p_(__p) {}
  78     _LIBCPP_INLINE_VISIBILITY
  79     void reset() {}
  80
  81     // generating functions
  82     template<class _URNG>
  83         _LIBCPP_INLINE_VISIBILITY
  84         result_type operator()(_URNG& __g)
  85         {return (*this)(__g, __p_);}
  86     template<class _URNG> result_type operator()(_URNG& __g, const param_type& __p);
  87
  88     // property functions
  89     _LIBCPP_INLINE_VISIBILITY
  90     result_type t() const {return __p_.t();}
  91     _LIBCPP_INLINE_VISIBILITY
  92     double p() const {return __p_.p();}
  93
  94     _LIBCPP_INLINE_VISIBILITY
  95     param_type param() const {return __p_;}
  96     _LIBCPP_INLINE_VISIBILITY
  97     void param(const param_type& __p) {__p_ = __p;}
  98
  99     _LIBCPP_INLINE_VISIBILITY
 100     result_type min() const {return 0;}
 101     _LIBCPP_INLINE_VISIBILITY
 102     result_type max() const {return t();}
 103
 104     friend _LIBCPP_INLINE_VISIBILITY
 105         bool operator==(const binomial_distribution& __x,
 106                         const binomial_distribution& __y)
 107         {return __x.__p_ == __y.__p_;}
 108     friend _LIBCPP_INLINE_VISIBILITY
 109         bool operator!=(const binomial_distribution& __x,
 110                         const binomial_distribution& __y)
 111         {return !(__x == __y);}
 112 };
 113
 114 #ifndef _LIBCPP_MSVCRT_LIKE
 115 extern "C" double lgamma_r(double, int *);
 116 #endif
 117
 118 inline _LIBCPP_INLINE_VISIBILITY double __libcpp_lgamma(double __d) {
 119 #if defined(_LIBCPP_MSVCRT_LIKE)
 120   return lgamma(__d);
 121 #else
 122   int __sign;
 123   return lgamma_r(__d, &__sign);
 124 #endif
 125 }
 126
 127 template<class _IntType>
 128 binomial_distribution<_IntType>::param_type::param_type(result_type __t, double __p)
 129     : __t_(__t), __p_(__p)
 130 {
 131     if (0 < __p_ && __p_ < 1)
 132     {
 133         __r0_ = static_cast<result_type>((__t_ + 1) * __p_);
 134         __pr_ = _VSTD::exp(__libcpp_lgamma(__t_ + 1.) -
 135                            __libcpp_lgamma(__r0_ + 1.) -
 136                            __libcpp_lgamma(__t_ - __r0_ + 1.) + __r0_ * _VSTD::log(__p_) +
 137                            (__t_ - __r0_) * _VSTD::log(1 - __p_));
 138         __odds_ratio_ = __p_ / (1 - __p_);
 139     }
 140 }
 141
 142 // Reference: Kemp, C.D. (1986). `A modal method for generating binomial
 143 //           variables', Commun. Statist. - Theor. Meth. 15(3), 805-813.
 144 template<class _IntType>
 145 template<class _URNG>
 146 _IntType
 147 binomial_distribution<_IntType>::operator()(_URNG& __g, const param_type& __pr)
 148 {
 149     if (__pr.__t_ == 0 || __pr.__p_ == 0)
 150         return 0;
 151     if (__pr.__p_ == 1)
 152         return __pr.__t_;
 153     uniform_real_distribution<double> __gen;
 154     double __u = __gen(__g) - __pr.__pr_;
 155     if (__u < 0)
 156         return __pr.__r0_;
 157     double __pu = __pr.__pr_;
 158     double __pd = __pu;
 159     result_type __ru = __pr.__r0_;
 160     result_type __rd = __ru;
 161     while (true)
 162     {
 163         bool __break = true;
 164         if (__rd >= 1)
 165         {
 166             __pd *= __rd / (__pr.__odds_ratio_ * (__pr.__t_ - __rd + 1));
 167             __u -= __pd;
 168             __break = false;
 169             if (__u < 0)
 170                 return __rd - 1;
 171         }
 172         if ( __rd != 0 )
 173             --__rd;
 174         ++__ru;
 175         if (__ru <= __pr.__t_)
 176         {
 177             __pu *= (__pr.__t_ - __ru + 1) * __pr.__odds_ratio_ / __ru;
 178             __u -= __pu;
 179             __break = false;
 180             if (__u < 0)
 181                 return __ru;
 182         }
 183         if (__break)
 184             return 0;
 185     }
 186 }
 187
 188 template <class _CharT, class _Traits, class _IntType>
 189 basic_ostream<_CharT, _Traits>&
 190 operator<<(basic_ostream<_CharT, _Traits>& __os,
 191            const binomial_distribution<_IntType>& __x)
 192 {
 193     __save_flags<_CharT, _Traits> __lx(__os);
 194     typedef basic_ostream<_CharT, _Traits> _OStream;
 195     __os.flags(_OStream::dec | _OStream::left | _OStream::fixed |
 196                _OStream::scientific);
 197     _CharT __sp = __os.widen(' ');
 198     __os.fill(__sp);
 199     return __os << __x.t() << __sp << __x.p();
 200 }
 201
 202 template <class _CharT, class _Traits, class _IntType>
 203 basic_istream<_CharT, _Traits>&
 204 operator>>(basic_istream<_CharT, _Traits>& __is,
 205            binomial_distribution<_IntType>& __x)
 206 {
 207     typedef binomial_distribution<_IntType> _Eng;
 208     typedef typename _Eng::result_type result_type;
 209     typedef typename _Eng::param_type param_type;
 210     __save_flags<_CharT, _Traits> __lx(__is);
 211     typedef basic_istream<_CharT, _Traits> _Istream;
 212     __is.flags(_Istream::dec | _Istream::skipws);
 213     result_type __t;
 214     double __p;
 215     __is >> __t >> __p;
 216     if (!__is.fail())
 217         __x.param(param_type(__t, __p));
 218     return __is;
 219 }
 220
 221 _LIBCPP_END_NAMESPACE_STD
 222
 223 _LIBCPP_POP_MACROS
 224
 225 #endif // _LIBCPP___RANDOM_BINOMIAL_DISTRIBUTION_H