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 static_assert(__libcpp_random_is_valid_inttype<_IntType>::value, "IntType must be a supported integer type"); 30 31public: 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 58private: 59 param_type __p_; 60 61public: 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 101extern "C" double lgamma_r(double, int*); 102#endif 103 104inline _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 113template <class _IntType> 114binomial_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. 126template <class _IntType> 127template <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 166template <class _CharT, class _Traits, class _IntType> 167_LIBCPP_HIDE_FROM_ABI basic_ostream<_CharT, _Traits>& 168operator<<(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 177template <class _CharT, class _Traits, class _IntType> 178_LIBCPP_HIDE_FROM_ABI basic_istream<_CharT, _Traits>& 179operator>>(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