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 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> result_type operator()(_URNG& __g, const param_type& __p); 97 98 // property functions 99 _LIBCPP_INLINE_VISIBILITY 100 double mean() const {return __p_.mean();} 101 102 _LIBCPP_INLINE_VISIBILITY 103 param_type param() const {return __p_;} 104 _LIBCPP_INLINE_VISIBILITY 105 void param(const param_type& __p) {__p_ = __p;} 106 107 _LIBCPP_INLINE_VISIBILITY 108 result_type min() const {return 0;} 109 _LIBCPP_INLINE_VISIBILITY 110 result_type max() const {return numeric_limits<result_type>::max();} 111 112 friend _LIBCPP_INLINE_VISIBILITY 113 bool operator==(const poisson_distribution& __x, 114 const poisson_distribution& __y) 115 {return __x.__p_ == __y.__p_;} 116 friend _LIBCPP_INLINE_VISIBILITY 117 bool operator!=(const poisson_distribution& __x, 118 const poisson_distribution& __y) 119 {return !(__x == __y);} 120 }; 121 122 template<class _IntType> 123 poisson_distribution<_IntType>::param_type::param_type(double __mean) 124 // According to the standard `inf` is a valid input, but it causes the 125 // distribution to hang, so we replace it with the maximum representable 126 // mean. 127 : __mean_(isinf(__mean) ? numeric_limits<double>::max() : __mean) 128 { 129 if (__mean_ < 10) 130 { 131 __s_ = 0; 132 __d_ = 0; 133 __l_ = _VSTD::exp(-__mean_); 134 __omega_ = 0; 135 __c3_ = 0; 136 __c2_ = 0; 137 __c1_ = 0; 138 __c0_ = 0; 139 __c_ = 0; 140 } 141 else 142 { 143 __s_ = _VSTD::sqrt(__mean_); 144 __d_ = 6 * __mean_ * __mean_; 145 __l_ = _VSTD::trunc(__mean_ - 1.1484); 146 __omega_ = .3989423 / __s_; 147 double __b1_ = .4166667E-1 / __mean_; 148 double __b2_ = .3 * __b1_ * __b1_; 149 __c3_ = .1428571 * __b1_ * __b2_; 150 __c2_ = __b2_ - 15. * __c3_; 151 __c1_ = __b1_ - 6. * __b2_ + 45. * __c3_; 152 __c0_ = 1. - __b1_ + 3. * __b2_ - 15. * __c3_; 153 __c_ = .1069 / __mean_; 154 } 155 } 156 157 template <class _IntType> 158 template<class _URNG> 159 _IntType 160 poisson_distribution<_IntType>::operator()(_URNG& __urng, const param_type& __pr) 161 { 162 static_assert(__libcpp_random_is_valid_urng<_URNG>::value, ""); 163 double __tx; 164 uniform_real_distribution<double> __urd; 165 if (__pr.__mean_ < 10) 166 { 167 __tx = 0; 168 for (double __p = __urd(__urng); __p > __pr.__l_; ++__tx) 169 __p *= __urd(__urng); 170 } 171 else 172 { 173 double __difmuk; 174 double __g = __pr.__mean_ + __pr.__s_ * normal_distribution<double>()(__urng); 175 double __u; 176 if (__g > 0) 177 { 178 __tx = _VSTD::trunc(__g); 179 if (__tx >= __pr.__l_) 180 return _VSTD::__clamp_to_integral<result_type>(__tx); 181 __difmuk = __pr.__mean_ - __tx; 182 __u = __urd(__urng); 183 if (__pr.__d_ * __u >= __difmuk * __difmuk * __difmuk) 184 return _VSTD::__clamp_to_integral<result_type>(__tx); 185 } 186 exponential_distribution<double> __edist; 187 for (bool __using_exp_dist = false; true; __using_exp_dist = true) 188 { 189 double __e; 190 if (__using_exp_dist || __g <= 0) 191 { 192 double __t; 193 do 194 { 195 __e = __edist(__urng); 196 __u = __urd(__urng); 197 __u += __u - 1; 198 __t = 1.8 + (__u < 0 ? -__e : __e); 199 } while (__t <= -.6744); 200 __tx = _VSTD::trunc(__pr.__mean_ + __pr.__s_ * __t); 201 __difmuk = __pr.__mean_ - __tx; 202 __using_exp_dist = true; 203 } 204 double __px; 205 double __py; 206 if (__tx < 10 && __tx >= 0) 207 { 208 const double __fac[] = {1, 1, 2, 6, 24, 120, 720, 5040, 209 40320, 362880}; 210 __px = -__pr.__mean_; 211 __py = _VSTD::pow(__pr.__mean_, (double)__tx) / __fac[static_cast<int>(__tx)]; 212 } 213 else 214 { 215 double __del = .8333333E-1 / __tx; 216 __del -= 4.8 * __del * __del * __del; 217 double __v = __difmuk / __tx; 218 if (_VSTD::abs(__v) > 0.25) 219 __px = __tx * _VSTD::log(1 + __v) - __difmuk - __del; 220 else 221 __px = __tx * __v * __v * (((((((.1250060 * __v + -.1384794) * 222 __v + .1421878) * __v + -.1661269) * __v + .2000118) * 223 __v + -.2500068) * __v + .3333333) * __v + -.5) - __del; 224 __py = .3989423 / _VSTD::sqrt(__tx); 225 } 226 double __r = (0.5 - __difmuk) / __pr.__s_; 227 double __r2 = __r * __r; 228 double __fx = -0.5 * __r2; 229 double __fy = __pr.__omega_ * (((__pr.__c3_ * __r2 + __pr.__c2_) * 230 __r2 + __pr.__c1_) * __r2 + __pr.__c0_); 231 if (__using_exp_dist) 232 { 233 if (__pr.__c_ * _VSTD::abs(__u) <= __py * _VSTD::exp(__px + __e) - 234 __fy * _VSTD::exp(__fx + __e)) 235 break; 236 } 237 else 238 { 239 if (__fy - __u * __fy <= __py * _VSTD::exp(__px - __fx)) 240 break; 241 } 242 } 243 } 244 return _VSTD::__clamp_to_integral<result_type>(__tx); 245 } 246 247 template <class _CharT, class _Traits, class _IntType> 248 basic_ostream<_CharT, _Traits>& 249 operator<<(basic_ostream<_CharT, _Traits>& __os, 250 const poisson_distribution<_IntType>& __x) 251 { 252 __save_flags<_CharT, _Traits> __lx(__os); 253 typedef basic_ostream<_CharT, _Traits> _OStream; 254 __os.flags(_OStream::dec | _OStream::left | _OStream::fixed | 255 _OStream::scientific); 256 return __os << __x.mean(); 257 } 258 259 template <class _CharT, class _Traits, class _IntType> 260 basic_istream<_CharT, _Traits>& 261 operator>>(basic_istream<_CharT, _Traits>& __is, 262 poisson_distribution<_IntType>& __x) 263 { 264 typedef poisson_distribution<_IntType> _Eng; 265 typedef typename _Eng::param_type param_type; 266 __save_flags<_CharT, _Traits> __lx(__is); 267 typedef basic_istream<_CharT, _Traits> _Istream; 268 __is.flags(_Istream::dec | _Istream::skipws); 269 double __mean; 270 __is >> __mean; 271 if (!__is.fail()) 272 __x.param(param_type(__mean)); 273 return __is; 274 } 275 276 _LIBCPP_END_NAMESPACE_STD 277 278 _LIBCPP_POP_MACROS 279 280 #endif // _LIBCPP___RANDOM_POISSON_DISTRIBUTION_H 281