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