File indexing completed on 2026-05-03 08:13:39
0001
0002
0003
0004
0005
0006
0007
0008
0009 #ifndef _LIBCPP___CXX03___RANDOM_POISSON_DISTRIBUTION_H
0010 #define _LIBCPP___CXX03___RANDOM_POISSON_DISTRIBUTION_H
0011
0012 #include <__cxx03/__config>
0013 #include <__cxx03/__random/clamp_to_integral.h>
0014 #include <__cxx03/__random/exponential_distribution.h>
0015 #include <__cxx03/__random/is_valid.h>
0016 #include <__cxx03/__random/normal_distribution.h>
0017 #include <__cxx03/__random/uniform_real_distribution.h>
0018 #include <__cxx03/cmath>
0019 #include <__cxx03/iosfwd>
0020 #include <__cxx03/limits>
0021
0022 #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
0023 # pragma GCC system_header
0024 #endif
0025
0026 _LIBCPP_PUSH_MACROS
0027 #include <__cxx03/__undef_macros>
0028
0029 _LIBCPP_BEGIN_NAMESPACE_STD
0030
0031 template <class _IntType = int>
0032 class _LIBCPP_TEMPLATE_VIS poisson_distribution {
0033 static_assert(__libcpp_random_is_valid_inttype<_IntType>::value, "IntType must be a supported integer type");
0034
0035 public:
0036
0037 typedef _IntType result_type;
0038
0039 class _LIBCPP_TEMPLATE_VIS param_type {
0040 double __mean_;
0041 double __s_;
0042 double __d_;
0043 double __l_;
0044 double __omega_;
0045 double __c0_;
0046 double __c1_;
0047 double __c2_;
0048 double __c3_;
0049 double __c_;
0050
0051 public:
0052 typedef poisson_distribution distribution_type;
0053
0054 _LIBCPP_HIDE_FROM_ABI explicit param_type(double __mean = 1.0);
0055
0056 _LIBCPP_HIDE_FROM_ABI double mean() const { return __mean_; }
0057
0058 friend _LIBCPP_HIDE_FROM_ABI bool operator==(const param_type& __x, const param_type& __y) {
0059 return __x.__mean_ == __y.__mean_;
0060 }
0061 friend _LIBCPP_HIDE_FROM_ABI bool operator!=(const param_type& __x, const param_type& __y) { return !(__x == __y); }
0062
0063 friend class poisson_distribution;
0064 };
0065
0066 private:
0067 param_type __p_;
0068
0069 public:
0070
0071 #ifndef _LIBCPP_CXX03_LANG
0072 _LIBCPP_HIDE_FROM_ABI poisson_distribution() : poisson_distribution(1.0) {}
0073 _LIBCPP_HIDE_FROM_ABI explicit poisson_distribution(double __mean) : __p_(__mean) {}
0074 #else
0075 _LIBCPP_HIDE_FROM_ABI explicit poisson_distribution(double __mean = 1.0) : __p_(__mean) {}
0076 #endif
0077 _LIBCPP_HIDE_FROM_ABI explicit poisson_distribution(const param_type& __p) : __p_(__p) {}
0078 _LIBCPP_HIDE_FROM_ABI void reset() {}
0079
0080
0081 template <class _URNG>
0082 _LIBCPP_HIDE_FROM_ABI result_type operator()(_URNG& __g) {
0083 return (*this)(__g, __p_);
0084 }
0085 template <class _URNG>
0086 _LIBCPP_HIDE_FROM_ABI result_type operator()(_URNG& __g, const param_type& __p);
0087
0088
0089 _LIBCPP_HIDE_FROM_ABI double mean() const { return __p_.mean(); }
0090
0091 _LIBCPP_HIDE_FROM_ABI param_type param() const { return __p_; }
0092 _LIBCPP_HIDE_FROM_ABI void param(const param_type& __p) { __p_ = __p; }
0093
0094 _LIBCPP_HIDE_FROM_ABI result_type min() const { return 0; }
0095 _LIBCPP_HIDE_FROM_ABI result_type max() const { return numeric_limits<result_type>::max(); }
0096
0097 friend _LIBCPP_HIDE_FROM_ABI bool operator==(const poisson_distribution& __x, const poisson_distribution& __y) {
0098 return __x.__p_ == __y.__p_;
0099 }
0100 friend _LIBCPP_HIDE_FROM_ABI bool operator!=(const poisson_distribution& __x, const poisson_distribution& __y) {
0101 return !(__x == __y);
0102 }
0103 };
0104
0105 template <class _IntType>
0106 poisson_distribution<_IntType>::param_type::param_type(double __mean)
0107
0108
0109
0110 : __mean_(isinf(__mean) ? numeric_limits<double>::max() : __mean) {
0111 if (__mean_ < 10) {
0112 __s_ = 0;
0113 __d_ = 0;
0114 __l_ = std::exp(-__mean_);
0115 __omega_ = 0;
0116 __c3_ = 0;
0117 __c2_ = 0;
0118 __c1_ = 0;
0119 __c0_ = 0;
0120 __c_ = 0;
0121 } else {
0122 __s_ = std::sqrt(__mean_);
0123 __d_ = 6 * __mean_ * __mean_;
0124 __l_ = std::trunc(__mean_ - 1.1484);
0125 __omega_ = .3989423 / __s_;
0126 double __b1 = .4166667E-1 / __mean_;
0127 double __b2 = .3 * __b1 * __b1;
0128 __c3_ = .1428571 * __b1 * __b2;
0129 __c2_ = __b2 - 15. * __c3_;
0130 __c1_ = __b1 - 6. * __b2 + 45. * __c3_;
0131 __c0_ = 1. - __b1 + 3. * __b2 - 15. * __c3_;
0132 __c_ = .1069 / __mean_;
0133 }
0134 }
0135
0136 template <class _IntType>
0137 template <class _URNG>
0138 _IntType poisson_distribution<_IntType>::operator()(_URNG& __urng, const param_type& __pr) {
0139 static_assert(__libcpp_random_is_valid_urng<_URNG>::value, "");
0140 double __tx;
0141 uniform_real_distribution<double> __urd;
0142 if (__pr.__mean_ < 10) {
0143 __tx = 0;
0144 for (double __p = __urd(__urng); __p > __pr.__l_; ++__tx)
0145 __p *= __urd(__urng);
0146 } else {
0147 double __difmuk;
0148 double __g = __pr.__mean_ + __pr.__s_ * normal_distribution<double>()(__urng);
0149 double __u;
0150 if (__g > 0) {
0151 __tx = std::trunc(__g);
0152 if (__tx >= __pr.__l_)
0153 return std::__clamp_to_integral<result_type>(__tx);
0154 __difmuk = __pr.__mean_ - __tx;
0155 __u = __urd(__urng);
0156 if (__pr.__d_ * __u >= __difmuk * __difmuk * __difmuk)
0157 return std::__clamp_to_integral<result_type>(__tx);
0158 }
0159 exponential_distribution<double> __edist;
0160 for (bool __using_exp_dist = false; true; __using_exp_dist = true) {
0161 double __e;
0162 if (__using_exp_dist || __g <= 0) {
0163 double __t;
0164 do {
0165 __e = __edist(__urng);
0166 __u = __urd(__urng);
0167 __u += __u - 1;
0168 __t = 1.8 + (__u < 0 ? -__e : __e);
0169 } while (__t <= -.6744);
0170 __tx = std::trunc(__pr.__mean_ + __pr.__s_ * __t);
0171 __difmuk = __pr.__mean_ - __tx;
0172 __using_exp_dist = true;
0173 }
0174 double __px;
0175 double __py;
0176 if (__tx < 10 && __tx >= 0) {
0177 const double __fac[] = {1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880};
0178 __px = -__pr.__mean_;
0179 __py = std::pow(__pr.__mean_, (double)__tx) / __fac[static_cast<int>(__tx)];
0180 } else {
0181 double __del = .8333333E-1 / __tx;
0182 __del -= 4.8 * __del * __del * __del;
0183 double __v = __difmuk / __tx;
0184 if (std::abs(__v) > 0.25)
0185 __px = __tx * std::log(1 + __v) - __difmuk - __del;
0186 else
0187 __px = __tx * __v * __v *
0188 (((((((.1250060 * __v + -.1384794) * __v + .1421878) * __v + -.1661269) * __v + .2000118) * __v +
0189 -.2500068) *
0190 __v +
0191 .3333333) *
0192 __v +
0193 -.5) -
0194 __del;
0195 __py = .3989423 / std::sqrt(__tx);
0196 }
0197 double __r = (0.5 - __difmuk) / __pr.__s_;
0198 double __r2 = __r * __r;
0199 double __fx = -0.5 * __r2;
0200 double __fy = __pr.__omega_ * (((__pr.__c3_ * __r2 + __pr.__c2_) * __r2 + __pr.__c1_) * __r2 + __pr.__c0_);
0201 if (__using_exp_dist) {
0202 if (__pr.__c_ * std::abs(__u) <= __py * std::exp(__px + __e) - __fy * std::exp(__fx + __e))
0203 break;
0204 } else {
0205 if (__fy - __u * __fy <= __py * std::exp(__px - __fx))
0206 break;
0207 }
0208 }
0209 }
0210 return std::__clamp_to_integral<result_type>(__tx);
0211 }
0212
0213 template <class _CharT, class _Traits, class _IntType>
0214 _LIBCPP_HIDE_FROM_ABI basic_ostream<_CharT, _Traits>&
0215 operator<<(basic_ostream<_CharT, _Traits>& __os, const poisson_distribution<_IntType>& __x) {
0216 __save_flags<_CharT, _Traits> __lx(__os);
0217 typedef basic_ostream<_CharT, _Traits> _OStream;
0218 __os.flags(_OStream::dec | _OStream::left | _OStream::fixed | _OStream::scientific);
0219 return __os << __x.mean();
0220 }
0221
0222 template <class _CharT, class _Traits, class _IntType>
0223 _LIBCPP_HIDE_FROM_ABI basic_istream<_CharT, _Traits>&
0224 operator>>(basic_istream<_CharT, _Traits>& __is, poisson_distribution<_IntType>& __x) {
0225 typedef poisson_distribution<_IntType> _Eng;
0226 typedef typename _Eng::param_type param_type;
0227 __save_flags<_CharT, _Traits> __lx(__is);
0228 typedef basic_istream<_CharT, _Traits> _Istream;
0229 __is.flags(_Istream::dec | _Istream::skipws);
0230 double __mean;
0231 __is >> __mean;
0232 if (!__is.fail())
0233 __x.param(param_type(__mean));
0234 return __is;
0235 }
0236
0237 _LIBCPP_END_NAMESPACE_STD
0238
0239 _LIBCPP_POP_MACROS
0240
0241 #endif
