libcxx/include/__random/poisson_distribution.h

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