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_POISSON_DISTRIBUTION_H 10 #define _LIBCPP___RANDOM_POISSON_DISTRIBUTION_H 11 12 #include <__config> 13 #include <__random/clamp_to_integral.h> 14 #include <__random/exponential_distribution.h> 15 #include <__random/is_valid.h> 16 #include <__random/normal_distribution.h> 17 #include <__random/uniform_real_distribution.h> 18 #include <cmath> 19 #include <iosfwd> 20 #include <limits> 21 22 #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER) 23 # pragma GCC system_header 24 #endif 25 26 _LIBCPP_PUSH_MACROS 27 #include <__undef_macros> 28 29 _LIBCPP_BEGIN_NAMESPACE_STD 30 31 template<class _IntType = int> 32 class _LIBCPP_TEMPLATE_VIS poisson_distribution 33 { 34 static_assert(__libcpp_random_is_valid_inttype<_IntType>::value, "IntType must be a supported integer type"); 35 public: 36 // types 37 typedef _IntType result_type; 38 39 class _LIBCPP_TEMPLATE_VIS param_type 40 { 41 double __mean_; 42 double __s_; 43 double __d_; 44 double __l_; 45 double __omega_; 46 double __c0_; 47 double __c1_; 48 double __c2_; 49 double __c3_; 50 double __c_; 51 52 public: 53 typedef poisson_distribution distribution_type; 54 55 _LIBCPP_HIDE_FROM_ABI explicit param_type(double __mean = 1.0); 56 57 _LIBCPP_INLINE_VISIBILITY 58 double mean() const {return __mean_;} 59 60 friend _LIBCPP_INLINE_VISIBILITY 61 bool operator==(const param_type& __x, const param_type& __y) 62 {return __x.__mean_ == __y.__mean_;} 63 friend _LIBCPP_INLINE_VISIBILITY 64 bool operator!=(const param_type& __x, const param_type& __y) 65 {return !(__x == __y);} 66 67 friend class poisson_distribution; 68 }; 69 70 private: 71 param_type __p_; 72 73 public: 74 // constructors and reset functions 75 #ifndef _LIBCPP_CXX03_LANG 76 _LIBCPP_INLINE_VISIBILITY 77 poisson_distribution() : poisson_distribution(1.0) {} 78 _LIBCPP_INLINE_VISIBILITY 79 explicit poisson_distribution(double __mean) 80 : __p_(__mean) {} 81 #else 82 _LIBCPP_INLINE_VISIBILITY 83 explicit poisson_distribution(double __mean = 1.0) 84 : __p_(__mean) {} 85 #endif 86 _LIBCPP_INLINE_VISIBILITY 87 explicit poisson_distribution(const param_type& __p) : __p_(__p) {} 88 _LIBCPP_INLINE_VISIBILITY 89 void reset() {} 90 91 // generating functions 92 template<class _URNG> 93 _LIBCPP_INLINE_VISIBILITY 94 result_type operator()(_URNG& __g) 95 {return (*this)(__g, __p_);} 96 template<class _URNG> 97 _LIBCPP_HIDE_FROM_ABI result_type operator()(_URNG& __g, const param_type& __p); 98 99 // property functions 100 _LIBCPP_INLINE_VISIBILITY 101 double mean() const {return __p_.mean();} 102 103 _LIBCPP_INLINE_VISIBILITY 104 param_type param() const {return __p_;} 105 _LIBCPP_INLINE_VISIBILITY 106 void param(const param_type& __p) {__p_ = __p;} 107 108 _LIBCPP_INLINE_VISIBILITY 109 result_type min() const {return 0;} 110 _LIBCPP_INLINE_VISIBILITY 111 result_type max() const {return numeric_limits<result_type>::max();} 112 113 friend _LIBCPP_INLINE_VISIBILITY 114 bool operator==(const poisson_distribution& __x, 115 const poisson_distribution& __y) 116 {return __x.__p_ == __y.__p_;} 117 friend _LIBCPP_INLINE_VISIBILITY 118 bool operator!=(const poisson_distribution& __x, 119 const poisson_distribution& __y) 120 {return !(__x == __y);} 121 }; 122 123 template<class _IntType> 124 poisson_distribution<_IntType>::param_type::param_type(double __mean) 125 // According to the standard `inf` is a valid input, but it causes the 126 // distribution to hang, so we replace it with the maximum representable 127 // mean. 128 : __mean_(isinf(__mean) ? numeric_limits<double>::max() : __mean) 129 { 130 if (__mean_ < 10) 131 { 132 __s_ = 0; 133 __d_ = 0; 134 __l_ = _VSTD::exp(-__mean_); 135 __omega_ = 0; 136 __c3_ = 0; 137 __c2_ = 0; 138 __c1_ = 0; 139 __c0_ = 0; 140 __c_ = 0; 141 } 142 else 143 { 144 __s_ = _VSTD::sqrt(__mean_); 145 __d_ = 6 * __mean_ * __mean_; 146 __l_ = _VSTD::trunc(__mean_ - 1.1484); 147 __omega_ = .3989423 / __s_; 148 double __b1 = .4166667E-1 / __mean_; 149 double __b2 = .3 * __b1 * __b1; 150 __c3_ = .1428571 * __b1 * __b2; 151 __c2_ = __b2 - 15. * __c3_; 152 __c1_ = __b1 - 6. * __b2 + 45. * __c3_; 153 __c0_ = 1. - __b1 + 3. * __b2 - 15. * __c3_; 154 __c_ = .1069 / __mean_; 155 } 156 } 157 158 template <class _IntType> 159 template<class _URNG> 160 _IntType 161 poisson_distribution<_IntType>::operator()(_URNG& __urng, const param_type& __pr) 162 { 163 static_assert(__libcpp_random_is_valid_urng<_URNG>::value, ""); 164 double __tx; 165 uniform_real_distribution<double> __urd; 166 if (__pr.__mean_ < 10) 167 { 168 __tx = 0; 169 for (double __p = __urd(__urng); __p > __pr.__l_; ++__tx) 170 __p *= __urd(__urng); 171 } 172 else 173 { 174 double __difmuk; 175 double __g = __pr.__mean_ + __pr.__s_ * normal_distribution<double>()(__urng); 176 double __u; 177 if (__g > 0) 178 { 179 __tx = _VSTD::trunc(__g); 180 if (__tx >= __pr.__l_) 181 return _VSTD::__clamp_to_integral<result_type>(__tx); 182 __difmuk = __pr.__mean_ - __tx; 183 __u = __urd(__urng); 184 if (__pr.__d_ * __u >= __difmuk * __difmuk * __difmuk) 185 return _VSTD::__clamp_to_integral<result_type>(__tx); 186 } 187 exponential_distribution<double> __edist; 188 for (bool __using_exp_dist = false; true; __using_exp_dist = true) 189 { 190 double __e; 191 if (__using_exp_dist || __g <= 0) 192 { 193 double __t; 194 do 195 { 196 __e = __edist(__urng); 197 __u = __urd(__urng); 198 __u += __u - 1; 199 __t = 1.8 + (__u < 0 ? -__e : __e); 200 } while (__t <= -.6744); 201 __tx = _VSTD::trunc(__pr.__mean_ + __pr.__s_ * __t); 202 __difmuk = __pr.__mean_ - __tx; 203 __using_exp_dist = true; 204 } 205 double __px; 206 double __py; 207 if (__tx < 10 && __tx >= 0) 208 { 209 const double __fac[] = {1, 1, 2, 6, 24, 120, 720, 5040, 210 40320, 362880}; 211 __px = -__pr.__mean_; 212 __py = _VSTD::pow(__pr.__mean_, (double)__tx) / __fac[static_cast<int>(__tx)]; 213 } 214 else 215 { 216 double __del = .8333333E-1 / __tx; 217 __del -= 4.8 * __del * __del * __del; 218 double __v = __difmuk / __tx; 219 if (_VSTD::abs(__v) > 0.25) 220 __px = __tx * _VSTD::log(1 + __v) - __difmuk - __del; 221 else 222 __px = __tx * __v * __v * (((((((.1250060 * __v + -.1384794) * 223 __v + .1421878) * __v + -.1661269) * __v + .2000118) * 224 __v + -.2500068) * __v + .3333333) * __v + -.5) - __del; 225 __py = .3989423 / _VSTD::sqrt(__tx); 226 } 227 double __r = (0.5 - __difmuk) / __pr.__s_; 228 double __r2 = __r * __r; 229 double __fx = -0.5 * __r2; 230 double __fy = __pr.__omega_ * (((__pr.__c3_ * __r2 + __pr.__c2_) * 231 __r2 + __pr.__c1_) * __r2 + __pr.__c0_); 232 if (__using_exp_dist) 233 { 234 if (__pr.__c_ * _VSTD::abs(__u) <= __py * _VSTD::exp(__px + __e) - 235 __fy * _VSTD::exp(__fx + __e)) 236 break; 237 } 238 else 239 { 240 if (__fy - __u * __fy <= __py * _VSTD::exp(__px - __fx)) 241 break; 242 } 243 } 244 } 245 return _VSTD::__clamp_to_integral<result_type>(__tx); 246 } 247 248 template <class _CharT, class _Traits, class _IntType> 249 _LIBCPP_HIDE_FROM_ABI basic_ostream<_CharT, _Traits>& 250 operator<<(basic_ostream<_CharT, _Traits>& __os, 251 const poisson_distribution<_IntType>& __x) 252 { 253 __save_flags<_CharT, _Traits> __lx(__os); 254 typedef basic_ostream<_CharT, _Traits> _OStream; 255 __os.flags(_OStream::dec | _OStream::left | _OStream::fixed | 256 _OStream::scientific); 257 return __os << __x.mean(); 258 } 259 260 template <class _CharT, class _Traits, class _IntType> 261 _LIBCPP_HIDE_FROM_ABI basic_istream<_CharT, _Traits>& 262 operator>>(basic_istream<_CharT, _Traits>& __is, 263 poisson_distribution<_IntType>& __x) 264 { 265 typedef poisson_distribution<_IntType> _Eng; 266 typedef typename _Eng::param_type param_type; 267 __save_flags<_CharT, _Traits> __lx(__is); 268 typedef basic_istream<_CharT, _Traits> _Istream; 269 __is.flags(_Istream::dec | _Istream::skipws); 270 double __mean; 271 __is >> __mean; 272 if (!__is.fail()) 273 __x.param(param_type(__mean)); 274 return __is; 275 } 276 277 _LIBCPP_END_NAMESPACE_STD 278 279 _LIBCPP_POP_MACROS 280 281 #endif // _LIBCPP___RANDOM_POISSON_DISTRIBUTION_H 282