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_BINOMIAL_DISTRIBUTION_H
10#define _LIBCPP___RANDOM_BINOMIAL_DISTRIBUTION_H
11
12#include <__config>
13#include <__random/is_valid.h>
14#include <__random/uniform_real_distribution.h>
15#include <cmath>
16#include <iosfwd>
17
18#if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
19#  pragma GCC system_header
20#endif
21
22_LIBCPP_PUSH_MACROS
23#include <__undef_macros>
24
25_LIBCPP_BEGIN_NAMESPACE_STD
26
27template<class _IntType = int>
28class _LIBCPP_TEMPLATE_VIS binomial_distribution
29{
30    static_assert(__libcpp_random_is_valid_inttype<_IntType>::value, "IntType must be a supported integer type");
31public:
32    // types
33    typedef _IntType result_type;
34
35    class _LIBCPP_TEMPLATE_VIS param_type
36    {
37        result_type __t_;
38        double __p_;
39        double __pr_;
40        double __odds_ratio_;
41        result_type __r0_;
42    public:
43        typedef binomial_distribution distribution_type;
44
45        explicit param_type(result_type __t = 1, double __p = 0.5);
46
47        _LIBCPP_INLINE_VISIBILITY
48        result_type t() const {return __t_;}
49        _LIBCPP_INLINE_VISIBILITY
50        double p() const {return __p_;}
51
52        friend _LIBCPP_INLINE_VISIBILITY
53            bool operator==(const param_type& __x, const param_type& __y)
54            {return __x.__t_ == __y.__t_ && __x.__p_ == __y.__p_;}
55        friend _LIBCPP_INLINE_VISIBILITY
56            bool operator!=(const param_type& __x, const param_type& __y)
57            {return !(__x == __y);}
58
59        friend class binomial_distribution;
60    };
61
62private:
63    param_type __p_;
64
65public:
66    // constructors and reset functions
67#ifndef _LIBCPP_CXX03_LANG
68    _LIBCPP_INLINE_VISIBILITY
69    binomial_distribution() : binomial_distribution(1) {}
70    _LIBCPP_INLINE_VISIBILITY
71    explicit binomial_distribution(result_type __t, double __p = 0.5)
72        : __p_(param_type(__t, __p)) {}
73#else
74    _LIBCPP_INLINE_VISIBILITY
75    explicit binomial_distribution(result_type __t = 1, double __p = 0.5)
76        : __p_(param_type(__t, __p)) {}
77#endif
78    _LIBCPP_INLINE_VISIBILITY
79    explicit binomial_distribution(const param_type& __p) : __p_(__p) {}
80    _LIBCPP_INLINE_VISIBILITY
81    void reset() {}
82
83    // generating functions
84    template<class _URNG>
85        _LIBCPP_INLINE_VISIBILITY
86        result_type operator()(_URNG& __g)
87        {return (*this)(__g, __p_);}
88    template<class _URNG> result_type operator()(_URNG& __g, const param_type& __p);
89
90    // property functions
91    _LIBCPP_INLINE_VISIBILITY
92    result_type t() const {return __p_.t();}
93    _LIBCPP_INLINE_VISIBILITY
94    double p() const {return __p_.p();}
95
96    _LIBCPP_INLINE_VISIBILITY
97    param_type param() const {return __p_;}
98    _LIBCPP_INLINE_VISIBILITY
99    void param(const param_type& __p) {__p_ = __p;}
100
101    _LIBCPP_INLINE_VISIBILITY
102    result_type min() const {return 0;}
103    _LIBCPP_INLINE_VISIBILITY
104    result_type max() const {return t();}
105
106    friend _LIBCPP_INLINE_VISIBILITY
107        bool operator==(const binomial_distribution& __x,
108                        const binomial_distribution& __y)
109        {return __x.__p_ == __y.__p_;}
110    friend _LIBCPP_INLINE_VISIBILITY
111        bool operator!=(const binomial_distribution& __x,
112                        const binomial_distribution& __y)
113        {return !(__x == __y);}
114};
115
116#ifndef _LIBCPP_MSVCRT_LIKE
117extern "C" double lgamma_r(double, int *);
118#endif
119
120inline _LIBCPP_INLINE_VISIBILITY double __libcpp_lgamma(double __d) {
121#if defined(_LIBCPP_MSVCRT_LIKE)
122  return lgamma(__d);
123#else
124  int __sign;
125  return lgamma_r(__d, &__sign);
126#endif
127}
128
129template<class _IntType>
130binomial_distribution<_IntType>::param_type::param_type(result_type __t, double __p)
131    : __t_(__t), __p_(__p)
132{
133    if (0 < __p_ && __p_ < 1)
134    {
135        __r0_ = static_cast<result_type>((__t_ + 1) * __p_);
136        __pr_ = _VSTD::exp(std::__libcpp_lgamma(__t_ + 1.) -
137                           std::__libcpp_lgamma(__r0_ + 1.) -
138                           std::__libcpp_lgamma(__t_ - __r0_ + 1.) + __r0_ * _VSTD::log(__p_) +
139                           (__t_ - __r0_) * _VSTD::log(1 - __p_));
140        __odds_ratio_ = __p_ / (1 - __p_);
141    }
142}
143
144// Reference: Kemp, C.D. (1986). `A modal method for generating binomial
145//           variables', Commun. Statist. - Theor. Meth. 15(3), 805-813.
146template<class _IntType>
147template<class _URNG>
148_IntType
149binomial_distribution<_IntType>::operator()(_URNG& __g, const param_type& __pr)
150{
151    static_assert(__libcpp_random_is_valid_urng<_URNG>::value, "");
152    if (__pr.__t_ == 0 || __pr.__p_ == 0)
153        return 0;
154    if (__pr.__p_ == 1)
155        return __pr.__t_;
156    uniform_real_distribution<double> __gen;
157    double __u = __gen(__g) - __pr.__pr_;
158    if (__u < 0)
159        return __pr.__r0_;
160    double __pu = __pr.__pr_;
161    double __pd = __pu;
162    result_type __ru = __pr.__r0_;
163    result_type __rd = __ru;
164    while (true)
165    {
166        bool __break = true;
167        if (__rd >= 1)
168        {
169            __pd *= __rd / (__pr.__odds_ratio_ * (__pr.__t_ - __rd + 1));
170            __u -= __pd;
171            __break = false;
172            if (__u < 0)
173                return __rd - 1;
174        }
175        if ( __rd != 0 )
176            --__rd;
177        ++__ru;
178        if (__ru <= __pr.__t_)
179        {
180            __pu *= (__pr.__t_ - __ru + 1) * __pr.__odds_ratio_ / __ru;
181            __u -= __pu;
182            __break = false;
183            if (__u < 0)
184                return __ru;
185        }
186        if (__break)
187            return 0;
188    }
189}
190
191template <class _CharT, class _Traits, class _IntType>
192_LIBCPP_HIDE_FROM_ABI basic_ostream<_CharT, _Traits>&
193operator<<(basic_ostream<_CharT, _Traits>& __os,
194           const binomial_distribution<_IntType>& __x)
195{
196    __save_flags<_CharT, _Traits> __lx(__os);
197    typedef basic_ostream<_CharT, _Traits> _OStream;
198    __os.flags(_OStream::dec | _OStream::left | _OStream::fixed |
199               _OStream::scientific);
200    _CharT __sp = __os.widen(' ');
201    __os.fill(__sp);
202    return __os << __x.t() << __sp << __x.p();
203}
204
205template <class _CharT, class _Traits, class _IntType>
206_LIBCPP_HIDE_FROM_ABI basic_istream<_CharT, _Traits>&
207operator>>(basic_istream<_CharT, _Traits>& __is,
208           binomial_distribution<_IntType>& __x)
209{
210    typedef binomial_distribution<_IntType> _Eng;
211    typedef typename _Eng::result_type result_type;
212    typedef typename _Eng::param_type param_type;
213    __save_flags<_CharT, _Traits> __lx(__is);
214    typedef basic_istream<_CharT, _Traits> _Istream;
215    __is.flags(_Istream::dec | _Istream::skipws);
216    result_type __t;
217    double __p;
218    __is >> __t >> __p;
219    if (!__is.fail())
220        __x.param(param_type(__t, __p));
221    return __is;
222}
223
224_LIBCPP_END_NAMESPACE_STD
225
226_LIBCPP_POP_MACROS
227
228#endif // _LIBCPP___RANDOM_BINOMIAL_DISTRIBUTION_H
229