]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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_CDF_HPP | |
9 | #define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_CDF_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, class Policy> | |
17 | T hypergeometric_cdf_imp(unsigned x, unsigned r, unsigned n, unsigned N, bool invert, const Policy& pol) | |
18 | { | |
19 | #ifdef BOOST_MSVC | |
20 | # pragma warning(push) | |
21 | # pragma warning(disable:4267) | |
22 | #endif | |
23 | BOOST_MATH_STD_USING | |
24 | T result = 0; | |
25 | T mode = floor(T(r + 1) * T(n + 1) / (N + 2)); | |
26 | if(x < mode) | |
27 | { | |
28 | result = hypergeometric_pdf<T>(x, r, n, N, pol); | |
29 | T diff = result; | |
30 | unsigned lower_limit = static_cast<unsigned>((std::max)(0, (int)(n + r) - (int)(N))); | |
31 | while(diff > (invert ? T(1) : result) * tools::epsilon<T>()) | |
32 | { | |
33 | diff = T(x) * T((N + x) - n - r) * diff / (T(1 + n - x) * T(1 + r - x)); | |
34 | result += diff; | |
35 | BOOST_MATH_INSTRUMENT_VARIABLE(x); | |
36 | BOOST_MATH_INSTRUMENT_VARIABLE(diff); | |
37 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
38 | if(x == lower_limit) | |
39 | break; | |
40 | --x; | |
41 | } | |
42 | } | |
43 | else | |
44 | { | |
45 | invert = !invert; | |
46 | unsigned upper_limit = (std::min)(r, n); | |
47 | if(x != upper_limit) | |
48 | { | |
49 | ++x; | |
50 | result = hypergeometric_pdf<T>(x, r, n, N, pol); | |
51 | T diff = result; | |
52 | while((x <= upper_limit) && (diff > (invert ? T(1) : result) * tools::epsilon<T>())) | |
53 | { | |
54 | diff = T(n - x) * T(r - x) * diff / (T(x + 1) * T((N + x + 1) - n - r)); | |
55 | result += diff; | |
56 | ++x; | |
57 | BOOST_MATH_INSTRUMENT_VARIABLE(x); | |
58 | BOOST_MATH_INSTRUMENT_VARIABLE(diff); | |
59 | BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
60 | } | |
61 | } | |
62 | } | |
63 | if(invert) | |
64 | result = 1 - result; | |
65 | return result; | |
66 | #ifdef BOOST_MSVC | |
67 | # pragma warning(pop) | |
68 | #endif | |
69 | } | |
70 | ||
71 | template <class T, class Policy> | |
72 | inline T hypergeometric_cdf(unsigned x, unsigned r, unsigned n, unsigned N, bool invert, const Policy&) | |
73 | { | |
74 | BOOST_FPU_EXCEPTION_GUARD | |
75 | typedef typename tools::promote_args<T>::type result_type; | |
76 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
77 | typedef typename policies::normalise< | |
78 | Policy, | |
79 | policies::promote_float<false>, | |
80 | policies::promote_double<false>, | |
81 | policies::discrete_quantile<>, | |
82 | policies::assert_undefined<> >::type forwarding_policy; | |
83 | ||
84 | value_type result; | |
85 | result = detail::hypergeometric_cdf_imp<value_type>(x, r, n, N, invert, forwarding_policy()); | |
86 | if(result > 1) | |
87 | { | |
88 | result = 1; | |
89 | } | |
90 | if(result < 0) | |
91 | { | |
92 | result = 0; | |
93 | } | |
94 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(result, "boost::math::hypergeometric_cdf<%1%>(%1%,%1%,%1%,%1%)"); | |
95 | } | |
96 | ||
97 | }}} // namespaces | |
98 | ||
99 | #endif | |
100 |