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 // The LLVM C library provides this with conflicting `noexcept` attributes. 101 #if !defined(_LIBCPP_MSVCRT_LIKE) && !defined(__LLVM_LIBC__) 102 extern "C" double lgamma_r(double, int*); 103 #endif 104 105 inline _LIBCPP_HIDE_FROM_ABI double __libcpp_lgamma(double __d) { 106 #if defined(_LIBCPP_MSVCRT_LIKE) || defined(__LLVM_LIBC__) 107 return lgamma(__d); 108 #else 109 int __sign; 110 return lgamma_r(__d, &__sign); 111 #endif 112 } 113 114 template <class _IntType> 115 binomial_distribution<_IntType>::param_type::param_type(result_type __t, double __p) : __t_(__t), __p_(__p) { 116 if (0 < __p_ && __p_ < 1) { 117 __r0_ = static_cast<result_type>((__t_ + 1) * __p_); 118 __pr_ = std::exp( 119 std::__libcpp_lgamma(__t_ + 1.) - std::__libcpp_lgamma(__r0_ + 1.) - std::__libcpp_lgamma(__t_ - __r0_ + 1.) + 120 __r0_ * std::log(__p_) + (__t_ - __r0_) * std::log(1 - __p_)); 121 __odds_ratio_ = __p_ / (1 - __p_); 122 } 123 } 124 125 // Reference: Kemp, C.D. (1986). `A modal method for generating binomial 126 // variables', Commun. Statist. - Theor. Meth. 15(3), 805-813. 127 template <class _IntType> 128 template <class _URNG> 129 _IntType binomial_distribution<_IntType>::operator()(_URNG& __g, const param_type& __pr) { 130 static_assert(__libcpp_random_is_valid_urng<_URNG>::value, ""); 131 if (__pr.__t_ == 0 || __pr.__p_ == 0) 132 return 0; 133 if (__pr.__p_ == 1) 134 return __pr.__t_; 135 uniform_real_distribution<double> __gen; 136 double __u = __gen(__g) - __pr.__pr_; 137 if (__u < 0) 138 return __pr.__r0_; 139 double __pu = __pr.__pr_; 140 double __pd = __pu; 141 result_type __ru = __pr.__r0_; 142 result_type __rd = __ru; 143 while (true) { 144 bool __break = true; 145 if (__rd >= 1) { 146 __pd *= __rd / (__pr.__odds_ratio_ * (__pr.__t_ - __rd + 1)); 147 __u -= __pd; 148 __break = false; 149 if (__u < 0) 150 return __rd - 1; 151 } 152 if (__rd != 0) 153 --__rd; 154 ++__ru; 155 if (__ru <= __pr.__t_) { 156 __pu *= (__pr.__t_ - __ru + 1) * __pr.__odds_ratio_ / __ru; 157 __u -= __pu; 158 __break = false; 159 if (__u < 0) 160 return __ru; 161 } 162 if (__break) 163 return 0; 164 } 165 } 166 167 template <class _CharT, class _Traits, class _IntType> 168 _LIBCPP_HIDE_FROM_ABI basic_ostream<_CharT, _Traits>& 169 operator<<(basic_ostream<_CharT, _Traits>& __os, const binomial_distribution<_IntType>& __x) { 170 __save_flags<_CharT, _Traits> __lx(__os); 171 typedef basic_ostream<_CharT, _Traits> _OStream; 172 __os.flags(_OStream::dec | _OStream::left | _OStream::fixed | _OStream::scientific); 173 _CharT __sp = __os.widen(' '); 174 __os.fill(__sp); 175 return __os << __x.t() << __sp << __x.p(); 176 } 177 178 template <class _CharT, class _Traits, class _IntType> 179 _LIBCPP_HIDE_FROM_ABI basic_istream<_CharT, _Traits>& 180 operator>>(basic_istream<_CharT, _Traits>& __is, binomial_distribution<_IntType>& __x) { 181 typedef binomial_distribution<_IntType> _Eng; 182 typedef typename _Eng::result_type result_type; 183 typedef typename _Eng::param_type param_type; 184 __save_flags<_CharT, _Traits> __lx(__is); 185 typedef basic_istream<_CharT, _Traits> _Istream; 186 __is.flags(_Istream::dec | _Istream::skipws); 187 result_type __t; 188 double __p; 189 __is >> __t >> __p; 190 if (!__is.fail()) 191 __x.param(param_type(__t, __p)); 192 return __is; 193 } 194 195 _LIBCPP_END_NAMESPACE_STD 196 197 _LIBCPP_POP_MACROS 198 199 #endif // _LIBCPP___RANDOM_BINOMIAL_DISTRIBUTION_H 200