1 // Copyright 2008 John Maddock
2 //
3 // Use, modification and distribution are subject to the
4 // Boost Software License, Version 1.0.
5 // (See accompanying file LICENSE_1_0.txt
6 // or copy at http://www.boost.org/LICENSE_1_0.txt)
7
8 #ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_QUANTILE_HPP
9 #define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_QUANTILE_HPP
10
11 #include <boost/math/policies/error_handling.hpp>
12 #include <boost/math/distributions/detail/hypergeometric_pdf.hpp>
13
14 namespace boost{ namespace math{ namespace detail{
15
16 template <class T>
round_x_from_p(unsigned x,T p,T cum,T fudge_factor,unsigned lbound,unsigned,const policies::discrete_quantile<policies::integer_round_down> &)17 inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned lbound, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_down>&)
18 {
19 if((p < cum * fudge_factor) && (x != lbound))
20 {
21 BOOST_MATH_INSTRUMENT_VARIABLE(x-1);
22 return --x;
23 }
24 return x;
25 }
26
27 template <class T>
round_x_from_p(unsigned x,T p,T cum,T fudge_factor,unsigned,unsigned ubound,const policies::discrete_quantile<policies::integer_round_up> &)28 inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned /*lbound*/, unsigned ubound, const policies::discrete_quantile<policies::integer_round_up>&)
29 {
30 if((cum < p * fudge_factor) && (x != ubound))
31 {
32 BOOST_MATH_INSTRUMENT_VARIABLE(x+1);
33 return ++x;
34 }
35 return x;
36 }
37
38 template <class T>
round_x_from_p(unsigned x,T p,T cum,T fudge_factor,unsigned lbound,unsigned ubound,const policies::discrete_quantile<policies::integer_round_inwards> &)39 inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned lbound, unsigned ubound, const policies::discrete_quantile<policies::integer_round_inwards>&)
40 {
41 if(p >= 0.5)
42 return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
43 return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
44 }
45
46 template <class T>
round_x_from_p(unsigned x,T p,T cum,T fudge_factor,unsigned lbound,unsigned ubound,const policies::discrete_quantile<policies::integer_round_outwards> &)47 inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned lbound, unsigned ubound, const policies::discrete_quantile<policies::integer_round_outwards>&)
48 {
49 if(p >= 0.5)
50 return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
51 return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
52 }
53
54 template <class T>
round_x_from_p(unsigned x,T,T,T,unsigned,unsigned,const policies::discrete_quantile<policies::integer_round_nearest> &)55 inline unsigned round_x_from_p(unsigned x, T /*p*/, T /*cum*/, T /*fudge_factor*/, unsigned /*lbound*/, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_nearest>&)
56 {
57 return x;
58 }
59
60 template <class T>
round_x_from_q(unsigned x,T q,T cum,T fudge_factor,unsigned lbound,unsigned,const policies::discrete_quantile<policies::integer_round_down> &)61 inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned lbound, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_down>&)
62 {
63 if((q * fudge_factor > cum) && (x != lbound))
64 {
65 BOOST_MATH_INSTRUMENT_VARIABLE(x-1);
66 return --x;
67 }
68 return x;
69 }
70
71 template <class T>
round_x_from_q(unsigned x,T q,T cum,T fudge_factor,unsigned,unsigned ubound,const policies::discrete_quantile<policies::integer_round_up> &)72 inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned /*lbound*/, unsigned ubound, const policies::discrete_quantile<policies::integer_round_up>&)
73 {
74 if((q < cum * fudge_factor) && (x != ubound))
75 {
76 BOOST_MATH_INSTRUMENT_VARIABLE(x+1);
77 return ++x;
78 }
79 return x;
80 }
81
82 template <class T>
round_x_from_q(unsigned x,T q,T cum,T fudge_factor,unsigned lbound,unsigned ubound,const policies::discrete_quantile<policies::integer_round_inwards> &)83 inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned lbound, unsigned ubound, const policies::discrete_quantile<policies::integer_round_inwards>&)
84 {
85 if(q < 0.5)
86 return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
87 return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
88 }
89
90 template <class T>
round_x_from_q(unsigned x,T q,T cum,T fudge_factor,unsigned lbound,unsigned ubound,const policies::discrete_quantile<policies::integer_round_outwards> &)91 inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned lbound, unsigned ubound, const policies::discrete_quantile<policies::integer_round_outwards>&)
92 {
93 if(q >= 0.5)
94 return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
95 return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
96 }
97
98 template <class T>
round_x_from_q(unsigned x,T,T,T,unsigned,unsigned,const policies::discrete_quantile<policies::integer_round_nearest> &)99 inline unsigned round_x_from_q(unsigned x, T /*q*/, T /*cum*/, T /*fudge_factor*/, unsigned /*lbound*/, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_nearest>&)
100 {
101 return x;
102 }
103
104 template <class T, class Policy>
hypergeometric_quantile_imp(T p,T q,unsigned r,unsigned n,unsigned N,const Policy & pol)105 unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned N, const Policy& pol)
106 {
107 #ifdef BOOST_MSVC
108 # pragma warning(push)
109 # pragma warning(disable:4267)
110 #endif
111 typedef typename Policy::discrete_quantile_type discrete_quantile_type;
112 BOOST_MATH_STD_USING
113 BOOST_FPU_EXCEPTION_GUARD
114 T result;
115 T fudge_factor = 1 + tools::epsilon<T>() * ((N <= boost::math::prime(boost::math::max_prime - 1)) ? 50 : 2 * N);
116 unsigned base = static_cast<unsigned>((std::max)(0, (int)(n + r) - (int)(N)));
117 unsigned lim = (std::min)(r, n);
118
119 BOOST_MATH_INSTRUMENT_VARIABLE(p);
120 BOOST_MATH_INSTRUMENT_VARIABLE(q);
121 BOOST_MATH_INSTRUMENT_VARIABLE(r);
122 BOOST_MATH_INSTRUMENT_VARIABLE(n);
123 BOOST_MATH_INSTRUMENT_VARIABLE(N);
124 BOOST_MATH_INSTRUMENT_VARIABLE(fudge_factor);
125 BOOST_MATH_INSTRUMENT_VARIABLE(base);
126 BOOST_MATH_INSTRUMENT_VARIABLE(lim);
127
128 if(p <= 0.5)
129 {
130 unsigned x = base;
131 result = hypergeometric_pdf<T>(x, r, n, N, pol);
132 T diff = result;
133 while(result < p)
134 {
135 diff = (diff > tools::min_value<T>() * 8)
136 ? T(n - x) * T(r - x) * diff / (T(x + 1) * T(N + x + 1 - n - r))
137 : hypergeometric_pdf<T>(x + 1, r, n, N, pol);
138 if(result + diff / 2 > p)
139 break;
140 ++x;
141 result += diff;
142 #ifdef BOOST_MATH_INSTRUMENT
143 if(diff != 0)
144 {
145 BOOST_MATH_INSTRUMENT_VARIABLE(x);
146 BOOST_MATH_INSTRUMENT_VARIABLE(diff);
147 BOOST_MATH_INSTRUMENT_VARIABLE(result);
148 }
149 #endif
150 }
151 return round_x_from_p(x, p, result, fudge_factor, base, lim, discrete_quantile_type());
152 }
153 else
154 {
155 unsigned x = lim;
156 result = 0;
157 T diff = hypergeometric_pdf<T>(x, r, n, N, pol);
158 while(result + diff / 2 < q)
159 {
160 result += diff;
161 diff = (diff > tools::min_value<T>() * 8)
162 ? x * T(N + x - n - r) * diff / (T(1 + n - x) * T(1 + r - x))
163 : hypergeometric_pdf<T>(x - 1, r, n, N, pol);
164 --x;
165 #ifdef BOOST_MATH_INSTRUMENT
166 if(diff != 0)
167 {
168 BOOST_MATH_INSTRUMENT_VARIABLE(x);
169 BOOST_MATH_INSTRUMENT_VARIABLE(diff);
170 BOOST_MATH_INSTRUMENT_VARIABLE(result);
171 }
172 #endif
173 }
174 return round_x_from_q(x, q, result, fudge_factor, base, lim, discrete_quantile_type());
175 }
176 #ifdef BOOST_MSVC
177 # pragma warning(pop)
178 #endif
179 }
180
181 template <class T, class Policy>
hypergeometric_quantile(T p,T q,unsigned r,unsigned n,unsigned N,const Policy &)182 inline unsigned hypergeometric_quantile(T p, T q, unsigned r, unsigned n, unsigned N, const Policy&)
183 {
184 BOOST_FPU_EXCEPTION_GUARD
185 typedef typename tools::promote_args<T>::type result_type;
186 typedef typename policies::evaluation<result_type, Policy>::type value_type;
187 typedef typename policies::normalise<
188 Policy,
189 policies::promote_float<false>,
190 policies::promote_double<false>,
191 policies::assert_undefined<> >::type forwarding_policy;
192
193 return detail::hypergeometric_quantile_imp<value_type>(p, q, r, n, N, forwarding_policy());
194 }
195
196 }}} // namespaces
197
198 #endif
199
200