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