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 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> result_type operator()(_URNG& __g, const param_type& __p); 89 90 // property functions 91 _LIBCPP_INLINE_VISIBILITY 92 result_type t() const {return __p_.t();} 93 _LIBCPP_INLINE_VISIBILITY 94 double p() const {return __p_.p();} 95 96 _LIBCPP_INLINE_VISIBILITY 97 param_type param() const {return __p_;} 98 _LIBCPP_INLINE_VISIBILITY 99 void param(const param_type& __p) {__p_ = __p;} 100 101 _LIBCPP_INLINE_VISIBILITY 102 result_type min() const {return 0;} 103 _LIBCPP_INLINE_VISIBILITY 104 result_type max() const {return t();} 105 106 friend _LIBCPP_INLINE_VISIBILITY 107 bool operator==(const binomial_distribution& __x, 108 const binomial_distribution& __y) 109 {return __x.__p_ == __y.__p_;} 110 friend _LIBCPP_INLINE_VISIBILITY 111 bool operator!=(const binomial_distribution& __x, 112 const binomial_distribution& __y) 113 {return !(__x == __y);} 114 }; 115 116 #ifndef _LIBCPP_MSVCRT_LIKE 117 extern "C" double lgamma_r(double, int *); 118 #endif 119 120 inline _LIBCPP_INLINE_VISIBILITY double __libcpp_lgamma(double __d) { 121 #if defined(_LIBCPP_MSVCRT_LIKE) 122 return lgamma(__d); 123 #else 124 int __sign; 125 return lgamma_r(__d, &__sign); 126 #endif 127 } 128 129 template<class _IntType> 130 binomial_distribution<_IntType>::param_type::param_type(result_type __t, double __p) 131 : __t_(__t), __p_(__p) 132 { 133 if (0 < __p_ && __p_ < 1) 134 { 135 __r0_ = static_cast<result_type>((__t_ + 1) * __p_); 136 __pr_ = _VSTD::exp(__libcpp_lgamma(__t_ + 1.) - 137 __libcpp_lgamma(__r0_ + 1.) - 138 __libcpp_lgamma(__t_ - __r0_ + 1.) + __r0_ * _VSTD::log(__p_) + 139 (__t_ - __r0_) * _VSTD::log(1 - __p_)); 140 __odds_ratio_ = __p_ / (1 - __p_); 141 } 142 } 143 144 // Reference: Kemp, C.D. (1986). `A modal method for generating binomial 145 // variables', Commun. Statist. - Theor. Meth. 15(3), 805-813. 146 template<class _IntType> 147 template<class _URNG> 148 _IntType 149 binomial_distribution<_IntType>::operator()(_URNG& __g, const param_type& __pr) 150 { 151 static_assert(__libcpp_random_is_valid_urng<_URNG>::value, ""); 152 if (__pr.__t_ == 0 || __pr.__p_ == 0) 153 return 0; 154 if (__pr.__p_ == 1) 155 return __pr.__t_; 156 uniform_real_distribution<double> __gen; 157 double __u = __gen(__g) - __pr.__pr_; 158 if (__u < 0) 159 return __pr.__r0_; 160 double __pu = __pr.__pr_; 161 double __pd = __pu; 162 result_type __ru = __pr.__r0_; 163 result_type __rd = __ru; 164 while (true) 165 { 166 bool __break = true; 167 if (__rd >= 1) 168 { 169 __pd *= __rd / (__pr.__odds_ratio_ * (__pr.__t_ - __rd + 1)); 170 __u -= __pd; 171 __break = false; 172 if (__u < 0) 173 return __rd - 1; 174 } 175 if ( __rd != 0 ) 176 --__rd; 177 ++__ru; 178 if (__ru <= __pr.__t_) 179 { 180 __pu *= (__pr.__t_ - __ru + 1) * __pr.__odds_ratio_ / __ru; 181 __u -= __pu; 182 __break = false; 183 if (__u < 0) 184 return __ru; 185 } 186 if (__break) 187 return 0; 188 } 189 } 190 191 template <class _CharT, class _Traits, class _IntType> 192 basic_ostream<_CharT, _Traits>& 193 operator<<(basic_ostream<_CharT, _Traits>& __os, 194 const binomial_distribution<_IntType>& __x) 195 { 196 __save_flags<_CharT, _Traits> __lx(__os); 197 typedef basic_ostream<_CharT, _Traits> _OStream; 198 __os.flags(_OStream::dec | _OStream::left | _OStream::fixed | 199 _OStream::scientific); 200 _CharT __sp = __os.widen(' '); 201 __os.fill(__sp); 202 return __os << __x.t() << __sp << __x.p(); 203 } 204 205 template <class _CharT, class _Traits, class _IntType> 206 basic_istream<_CharT, _Traits>& 207 operator>>(basic_istream<_CharT, _Traits>& __is, 208 binomial_distribution<_IntType>& __x) 209 { 210 typedef binomial_distribution<_IntType> _Eng; 211 typedef typename _Eng::result_type result_type; 212 typedef typename _Eng::param_type param_type; 213 __save_flags<_CharT, _Traits> __lx(__is); 214 typedef basic_istream<_CharT, _Traits> _Istream; 215 __is.flags(_Istream::dec | _Istream::skipws); 216 result_type __t; 217 double __p; 218 __is >> __t >> __p; 219 if (!__is.fail()) 220 __x.param(param_type(__t, __p)); 221 return __is; 222 } 223 224 _LIBCPP_END_NAMESPACE_STD 225 226 _LIBCPP_POP_MACROS 227 228 #endif // _LIBCPP___RANDOM_BINOMIAL_DISTRIBUTION_H 229