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_QUANTILE_HPP
Chris@16: #define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_QUANTILE_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>
Chris@16: inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned lbound, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_down>&)
Chris@16: {
Chris@16:    if((p < cum * fudge_factor) && (x != lbound))
Chris@16:    {
Chris@16:       BOOST_MATH_INSTRUMENT_VARIABLE(x-1);
Chris@16:       return --x;
Chris@16:    }
Chris@16:    return x;
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned /*lbound*/, unsigned ubound, const policies::discrete_quantile<policies::integer_round_up>&)
Chris@16: {
Chris@16:    if((cum < p * fudge_factor) && (x != ubound))
Chris@16:    {
Chris@16:       BOOST_MATH_INSTRUMENT_VARIABLE(x+1);
Chris@16:       return ++x;
Chris@16:    }
Chris@16:    return x;
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned lbound, unsigned ubound, const policies::discrete_quantile<policies::integer_round_inwards>&)
Chris@16: {
Chris@16:    if(p >= 0.5)
Chris@16:       return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
Chris@16:    return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline unsigned round_x_from_p(unsigned x, T p, T cum, T fudge_factor, unsigned lbound, unsigned ubound, const policies::discrete_quantile<policies::integer_round_outwards>&)
Chris@16: {
Chris@16:    if(p >= 0.5)
Chris@16:       return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
Chris@16:    return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline unsigned round_x_from_p(unsigned x, T /*p*/, T /*cum*/, T /*fudge_factor*/, unsigned /*lbound*/, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_nearest>&)
Chris@16: {
Chris@16:    return x;
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned lbound, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_down>&)
Chris@16: {
Chris@16:    if((q * fudge_factor > cum) && (x != lbound))
Chris@16:    {
Chris@16:       BOOST_MATH_INSTRUMENT_VARIABLE(x-1);
Chris@16:       return --x;
Chris@16:    }
Chris@16:    return x;
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned /*lbound*/, unsigned ubound, const policies::discrete_quantile<policies::integer_round_up>&)
Chris@16: {
Chris@16:    if((q < cum * fudge_factor) && (x != ubound))
Chris@16:    {
Chris@16:       BOOST_MATH_INSTRUMENT_VARIABLE(x+1);
Chris@16:       return ++x;
Chris@16:    }
Chris@16:    return x;
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned lbound, unsigned ubound, const policies::discrete_quantile<policies::integer_round_inwards>&)
Chris@16: {
Chris@16:    if(q < 0.5)
Chris@16:       return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
Chris@16:    return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline unsigned round_x_from_q(unsigned x, T q, T cum, T fudge_factor, unsigned lbound, unsigned ubound, const policies::discrete_quantile<policies::integer_round_outwards>&)
Chris@16: {
Chris@16:    if(q >= 0.5)
Chris@16:       return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
Chris@16:    return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline unsigned round_x_from_q(unsigned x, T /*q*/, T /*cum*/, T /*fudge_factor*/, unsigned /*lbound*/, unsigned /*ubound*/, const policies::discrete_quantile<policies::integer_round_nearest>&)
Chris@16: {
Chris@16:    return x;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: unsigned hypergeometric_quantile_imp(T p, T q, unsigned r, unsigned n, unsigned N, 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:    typedef typename Policy::discrete_quantile_type discrete_quantile_type;
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    T result;
Chris@16:    T fudge_factor = 1 + tools::epsilon<T>() * ((N <= boost::math::prime(boost::math::max_prime - 1)) ? 50 : 2 * N);
Chris@16:    unsigned base = static_cast<unsigned>((std::max)(0, (int)(n + r) - (int)(N)));
Chris@16:    unsigned lim = (std::min)(r, n);
Chris@16: 
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(p);
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(q);
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(r);
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(n);
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(N);
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(fudge_factor);
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(base);
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(lim);
Chris@16: 
Chris@16:    if(p <= 0.5)
Chris@16:    {
Chris@16:       unsigned x = base;
Chris@16:       result = hypergeometric_pdf<T>(x, r, n, N, pol);
Chris@16:       T diff = result;
Chris@16:       while(result < p)
Chris@16:       {
Chris@16:          diff = (diff > tools::min_value<T>() * 8) 
Chris@16:             ? T(n - x) * T(r - x) * diff / (T(x + 1) * T(N + x + 1 - n - r))
Chris@16:             : hypergeometric_pdf<T>(x + 1, r, n, N, pol);
Chris@16:          if(result + diff / 2 > p)
Chris@16:             break;
Chris@16:          ++x;
Chris@16:          result += diff;
Chris@16: #ifdef BOOST_MATH_INSTRUMENT
Chris@16:          if(diff != 0)
Chris@16:          {
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: #endif
Chris@16:       }
Chris@16:       return round_x_from_p(x, p, result, fudge_factor, base, lim, discrete_quantile_type());
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       unsigned x = lim;
Chris@16:       result = 0;
Chris@16:       T diff = hypergeometric_pdf<T>(x, r, n, N, pol);
Chris@16:       while(result + diff / 2 < q)
Chris@16:       {
Chris@16:          result += diff;
Chris@16:          diff = (diff > tools::min_value<T>() * 8)
Chris@16:             ? x * T(N + x - n - r) * diff / (T(1 + n - x) * T(1 + r - x))
Chris@16:             : hypergeometric_pdf<T>(x - 1, r, n, N, pol);
Chris@16:          --x;
Chris@16: #ifdef BOOST_MATH_INSTRUMENT
Chris@16:          if(diff != 0)
Chris@16:          {
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: #endif
Chris@16:       }
Chris@16:       return round_x_from_q(x, q, result, fudge_factor, base, lim, discrete_quantile_type());
Chris@16:    }
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 unsigned hypergeometric_quantile(T p, T q, unsigned r, unsigned n, unsigned N, 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::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    return detail::hypergeometric_quantile_imp<value_type>(p, q, r, n, N, forwarding_policy());
Chris@16: }
Chris@16: 
Chris@16: }}} // namespaces
Chris@16: 
Chris@16: #endif
Chris@16: