Chris@16: // Copyright 2008 John Maddock
Chris@16: //
Chris@16: // Use, modification and distribution are subject to the
Chris@16: // Boost Software License, Version 1.0.
Chris@16: // (See accompanying file LICENSE_1_0.txt
Chris@16: // or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_CDF_HPP
Chris@16: #define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_CDF_HPP
Chris@16: 
Chris@16: #include <boost/math/policies/error_handling.hpp>
Chris@16: #include <boost/math/distributions/detail/hypergeometric_pdf.hpp>
Chris@16: 
Chris@16: namespace boost{ namespace math{ namespace detail{
Chris@16: 
Chris@16:    template <class T, class Policy>
Chris@16:    T hypergeometric_cdf_imp(unsigned x, unsigned r, unsigned n, unsigned N, bool invert, const Policy& pol)
Chris@16:    {
Chris@16: #ifdef BOOST_MSVC
Chris@16: #  pragma warning(push)
Chris@16: #  pragma warning(disable:4267)
Chris@16: #endif
Chris@16:       BOOST_MATH_STD_USING
Chris@16:       T result = 0;
Chris@16:       T mode = floor(T(r + 1) * T(n + 1) / (N + 2));
Chris@16:       if(x < mode)
Chris@16:       {
Chris@16:          result = hypergeometric_pdf<T>(x, r, n, N, pol);
Chris@16:          T diff = result;
Chris@16:          unsigned lower_limit = static_cast<unsigned>((std::max)(0, (int)(n + r) - (int)(N)));
Chris@16:          while(diff > (invert ? T(1) : result) * tools::epsilon<T>())
Chris@16:          {
Chris@16:             diff = T(x) * T((N + x) - n - r) * diff / (T(1 + n - x) * T(1 + r - x));
Chris@16:             result += diff;
Chris@16:             BOOST_MATH_INSTRUMENT_VARIABLE(x);
Chris@16:             BOOST_MATH_INSTRUMENT_VARIABLE(diff);
Chris@16:             BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:             if(x == lower_limit)
Chris@16:                break;
Chris@16:             --x;
Chris@16:          }
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          invert = !invert;
Chris@16:          unsigned upper_limit = (std::min)(r, n);
Chris@16:          if(x != upper_limit)
Chris@16:          {
Chris@16:             ++x;
Chris@16:             result = hypergeometric_pdf<T>(x, r, n, N, pol);
Chris@16:             T diff = result;
Chris@16:             while((x <= upper_limit) && (diff > (invert ? T(1) : result) * tools::epsilon<T>()))
Chris@16:             {
Chris@16:                diff = T(n - x) * T(r - x) * diff / (T(x + 1) * T((N + x + 1) - n - r));
Chris@16:                result += diff;
Chris@16:                ++x;
Chris@16:                BOOST_MATH_INSTRUMENT_VARIABLE(x);
Chris@16:                BOOST_MATH_INSTRUMENT_VARIABLE(diff);
Chris@16:                BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:             }
Chris@16:          }
Chris@16:       }
Chris@16:       if(invert)
Chris@16:          result = 1 - result;
Chris@16:       return result;
Chris@16: #ifdef BOOST_MSVC
Chris@16: #  pragma warning(pop)
Chris@16: #endif
Chris@16:    }
Chris@16: 
Chris@16:    template <class T, class Policy>
Chris@16:    inline T hypergeometric_cdf(unsigned x, unsigned r, unsigned n, unsigned N, bool invert, const Policy&)
Chris@16:    {
Chris@16:       BOOST_FPU_EXCEPTION_GUARD
Chris@16:       typedef typename tools::promote_args<T>::type result_type;
Chris@16:       typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:       typedef typename policies::normalise<
Chris@16:          Policy, 
Chris@16:          policies::promote_float<false>, 
Chris@16:          policies::promote_double<false>, 
Chris@16:          policies::discrete_quantile<>,
Chris@16:          policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:       value_type result;
Chris@16:       result = detail::hypergeometric_cdf_imp<value_type>(x, r, n, N, invert, forwarding_policy());
Chris@16:       if(result > 1)
Chris@16:       {
Chris@16:          result  = 1;
Chris@16:       }
Chris@16:       if(result < 0)
Chris@16:       {
Chris@16:          result = 0;
Chris@16:       }
Chris@16:       return policies::checked_narrowing_cast<result_type, forwarding_policy>(result, "boost::math::hypergeometric_cdf<%1%>(%1%,%1%,%1%,%1%)");
Chris@16:    }
Chris@16: 
Chris@16: }}} // namespaces
Chris@16: 
Chris@16: #endif
Chris@16: