contrib/llvm-project/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/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