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>
hypergeometric_cdf(unsigned x,unsigned r,unsigned n,unsigned N,bool invert,const 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 
101