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/is_valid.h> 14 #include <__random/uniform_real_distribution.h> 15 #include <cmath> 16 #include <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 <__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 199