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