102 lines · plain
1// Copyright 2008 John Maddock2//3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0.5// (See accompanying file LICENSE_1_0.txt6// or copy at http://www.boost.org/LICENSE_1_0.txt)7 8#ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_CDF_HPP9#define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_CDF_HPP10 11#include <boost/math/policies/error_handling.hpp>12#include <boost/math/distributions/detail/hypergeometric_pdf.hpp>13#include <cstdint>14 15namespace boost{ namespace math{ namespace detail{16 17 template <class T, class Policy>18 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)19 {20#ifdef _MSC_VER21# pragma warning(push)22# pragma warning(disable:4267)23#endif24 BOOST_MATH_STD_USING25 T result = 0;26 T mode = floor(T(r + 1) * T(n + 1) / (N + 2));27 if(x < mode)28 {29 result = hypergeometric_pdf<T>(x, r, n, N, pol);30 T diff = result;31 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)));32 while(diff > (invert ? T(1) : result) * tools::epsilon<T>())33 {34 diff = T(x) * T((N + x) - n - r) * diff / (T(1 + n - x) * T(1 + r - x));35 result += diff;36 BOOST_MATH_INSTRUMENT_VARIABLE(x);37 BOOST_MATH_INSTRUMENT_VARIABLE(diff);38 BOOST_MATH_INSTRUMENT_VARIABLE(result);39 if(x == lower_limit)40 break;41 --x;42 }43 }44 else45 {46 invert = !invert;47 const auto upper_limit = (std::min)(r, n);48 if(x != upper_limit)49 {50 ++x;51 result = hypergeometric_pdf<T>(x, r, n, N, pol);52 T diff = result;53 while((x <= upper_limit) && (diff > (invert ? T(1) : result) * tools::epsilon<T>()))54 {55 diff = T(n - x) * T(r - x) * diff / (T(x + 1) * T((N + x + 1) - n - r));56 result += diff;57 ++x;58 BOOST_MATH_INSTRUMENT_VARIABLE(x);59 BOOST_MATH_INSTRUMENT_VARIABLE(diff);60 BOOST_MATH_INSTRUMENT_VARIABLE(result);61 }62 }63 }64 if(invert)65 result = 1 - result;66 return result;67#ifdef _MSC_VER68# pragma warning(pop)69#endif70 }71 72 template <class T, class Policy>73 inline T hypergeometric_cdf(std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, bool invert, const Policy&)74 {75 BOOST_FPU_EXCEPTION_GUARD76 typedef typename tools::promote_args<T>::type result_type;77 typedef typename policies::evaluation<result_type, Policy>::type value_type;78 typedef typename policies::normalise<79 Policy,80 policies::promote_float<false>,81 policies::promote_double<false>,82 policies::discrete_quantile<>,83 policies::assert_undefined<> >::type forwarding_policy;84 85 value_type result;86 result = detail::hypergeometric_cdf_imp<value_type>(x, r, n, N, invert, forwarding_policy());87 if(result > 1)88 {89 result = 1;90 }91 if(result < 0)92 {93 result = 0;94 }95 return policies::checked_narrowing_cast<result_type, forwarding_policy>(result, "boost::math::hypergeometric_cdf<%1%>(%1%,%1%,%1%,%1%)");96 }97 98}}} // namespaces99 100#endif101 102