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
mean()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
poisson_distribution()77     poisson_distribution() : poisson_distribution(1.0) {}
78     _LIBCPP_INLINE_VISIBILITY
poisson_distribution(double __mean)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
poisson_distribution(const param_type & __p)87     explicit poisson_distribution(const param_type& __p) : __p_(__p) {}
88     _LIBCPP_INLINE_VISIBILITY
reset()89     void reset() {}
90 
91     // generating functions
92     template<class _URNG>
93         _LIBCPP_INLINE_VISIBILITY
operator()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
mean()100     double mean() const {return __p_.mean();}
101 
102     _LIBCPP_INLINE_VISIBILITY
param()103     param_type param() const {return __p_;}
104     _LIBCPP_INLINE_VISIBILITY
param(const param_type & __p)105     void param(const param_type& __p) {__p_ = __p;}
106 
107     _LIBCPP_INLINE_VISIBILITY
min()108     result_type min() const {return 0;}
109     _LIBCPP_INLINE_VISIBILITY
max()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>
param_type(double __mean)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
operator()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 _LIBCPP_HIDE_FROM_ABI 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 _LIBCPP_HIDE_FROM_ABI 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