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