distributions/detail/hypergeometric_cdf.hpp

0001 // Copyright 2008 John Maddock
0002 //
0003 // Use, modification and distribution are subject to the
0004 // Boost Software License, Version 1.0.
0005 // (See accompanying file LICENSE_1_0.txt
0006 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0007
0008 #ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_CDF_HPP
0009 #define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_CDF_HPP
0010
0011 #include <boost/math/policies/error_handling.hpp>
0012 #include <boost/math/distributions/detail/hypergeometric_pdf.hpp>
0013 #include <cstdint>
0014
0015 namespace boost{ namespace math{ namespace detail{
0016
0017    template <class T, class Policy>
0018    T hypergeometric_cdf_imp(std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, bool invert, const Policy& pol)
0019    {
0020 #ifdef _MSC_VER
0021 #  pragma warning(push)
0022 #  pragma warning(disable:4267)
0023 #endif
0024       BOOST_MATH_STD_USING
0025       T result = 0;
0026       T mode = floor(T(r + 1) * T(n + 1) / (N + 2));
0027       if(x < mode)
0028       {
0029          result = hypergeometric_pdf<T>(x, r, n, N, pol);
0030          T diff = result;
0031          const auto lower_limit = static_cast<std::uint64_t>((std::max)(INT64_C(0), static_cast<std::int64_t>(n + r) - static_cast<std::int64_t>(N)));
0032          while(diff > (invert ? T(1) : result) * tools::epsilon<T>())
0033          {
0034             diff = T(x) * T((N + x) - n - r) * diff / (T(1 + n - x) * T(1 + r - x));
0035             result += diff;
0036             BOOST_MATH_INSTRUMENT_VARIABLE(x);
0037             BOOST_MATH_INSTRUMENT_VARIABLE(diff);
0038             BOOST_MATH_INSTRUMENT_VARIABLE(result);
0039             if(x == lower_limit)
0040                break;
0041             --x;
0042          }
0043       }
0044       else
0045       {
0046          invert = !invert;
0047          const auto upper_limit = (std::min)(r, n);
0048          if(x != upper_limit)
0049          {
0050             ++x;
0051             result = hypergeometric_pdf<T>(x, r, n, N, pol);
0052             T diff = result;
0053             while((x <= upper_limit) && (diff > (invert ? T(1) : result) * tools::epsilon<T>()))
0054             {
0055                diff = T(n - x) * T(r - x) * diff / (T(x + 1) * T((N + x + 1) - n - r));
0056                result += diff;
0057                ++x;
0058                BOOST_MATH_INSTRUMENT_VARIABLE(x);
0059                BOOST_MATH_INSTRUMENT_VARIABLE(diff);
0060                BOOST_MATH_INSTRUMENT_VARIABLE(result);
0061             }
0062          }
0063       }
0064       if(invert)
0065          result = 1 - result;
0066       return result;
0067 #ifdef _MSC_VER
0068 #  pragma warning(pop)
0069 #endif
0070    }
0071
0072    template <class T, class Policy>
0073    inline T hypergeometric_cdf(std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, bool invert, const Policy&)
0074    {
0075       BOOST_FPU_EXCEPTION_GUARD
0076       typedef typename tools::promote_args<T>::type result_type;
0077       typedef typename policies::evaluation<result_type, Policy>::type value_type;
0078       typedef typename policies::normalise<
0079          Policy,
0080          policies::promote_float<false>,
0081          policies::promote_double<false>,
0082          policies::discrete_quantile<>,
0083          policies::assert_undefined<> >::type forwarding_policy;
0084
0085       value_type result;
0086       result = detail::hypergeometric_cdf_imp<value_type>(x, r, n, N, invert, forwarding_policy());
0087       if(result > 1)
0088       {
0089          result  = 1;
0090       }
0091       if(result < 0)
0092       {
0093          result = 0;
0094       }
0095       return policies::checked_narrowing_cast<result_type, forwarding_policy>(result, "boost::math::hypergeometric_cdf<%1%>(%1%,%1%,%1%,%1%)");
0096    }
0097
0098 }}} // namespaces
0099
0100 #endif
0101