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 27template<class _IntType = int> 28class _LIBCPP_TEMPLATE_VIS binomial_distribution 29{ 30 static_assert(__libcpp_random_is_valid_inttype<_IntType>::value, "IntType must be a supported integer type"); 31public: 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 62private: 63 param_type __p_; 64 65public: 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 117extern "C" double lgamma_r(double, int *); 118#endif 119 120inline _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 129template<class _IntType> 130binomial_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(std::__libcpp_lgamma(__t_ + 1.) - 137 std::__libcpp_lgamma(__r0_ + 1.) - 138 std::__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. 146template<class _IntType> 147template<class _URNG> 148_IntType 149binomial_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 191template <class _CharT, class _Traits, class _IntType> 192_LIBCPP_HIDE_FROM_ABI basic_ostream<_CharT, _Traits>& 193operator<<(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 205template <class _CharT, class _Traits, class _IntType> 206_LIBCPP_HIDE_FROM_ABI basic_istream<_CharT, _Traits>& 207operator>>(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