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 {
  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