shylie/llvm-project
1
0
Fork 0
Code Issues Pull Requests Actions 6 Packages Projects Releases Wiki Activity
llvm-project/third-party/boost-math/include/boost/math/distributions/detail/hypergeometric_quantile.hpp
Louis Dionne 585da50d7d
[third-party] Add a snapshot of Boost.Math 1.89 standalone (#141508) 
This PR adds the code of Boost.Math as of version 1.89 into the
third-party directory, as discussed in a recent RFC [1].

The goal is for this code to be used as a back-end for the C++17
Math Special Functions.

As explained in third-paty/README.md, this code is cleared for
usage inside libc++ for the Math Special functions, however
the LLVM Foundation should be consulted before using this
code anywhere else in the LLVM project, due to the fact
that it is under the Boost Software License (as opposed
to the usual LLVM license). See the RFC [1] for more details.

[1]: https://discourse.llvm.org/t/rfc-libc-taking-a-dependency-on-boost-math-for-the-c-17-math-special-functions
2025-10-27 14:43:57 -07:00

246 lines
9.8 KiB
C++
Raw Blame History

// Copyright 2008 John Maddock
//
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)
#ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_QUANTILE_HPP
#define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_QUANTILE_HPP
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/distributions/detail/hypergeometric_pdf.hpp>
namespace boost{ namespace math{ namespace detail{
template <class T>
inline std::uint64_t round_x_from_p(std::uint64_t x, T p, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t /*ubound*/, const policies::discrete_quantile<policies::integer_round_down>&)
{
if((p < cum * fudge_factor) && (x != lbound))
{
BOOST_MATH_INSTRUMENT_VARIABLE(x-1);
return --x;
}
return x;
}
template <class T>
inline std::uint64_t round_x_from_p(std::uint64_t x, T p, T cum, T fudge_factor, std::uint64_t /*lbound*/, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_up>&)
{
if((cum < p * fudge_factor) && (x != ubound))
{
BOOST_MATH_INSTRUMENT_VARIABLE(x+1);
return ++x;
}
return x;
}
template <class T>
inline std::uint64_t round_x_from_p(std::uint64_t x, T p, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_inwards>&)
{
if(p >= 0.5)
return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
}
template <class T>
inline std::uint64_t round_x_from_p(std::uint64_t x, T p, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_outwards>&)
{
if(p >= 0.5)
return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
return round_x_from_p(x, p, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
}
template <class T>
inline std::uint64_t round_x_from_p(std::uint64_t x, T /*p*/, T /*cum*/, T /*fudge_factor*/, std::uint64_t /*lbound*/, std::uint64_t /*ubound*/, const policies::discrete_quantile<policies::integer_round_nearest>&)
{
return x;
}
template <class T>
inline std::uint64_t round_x_from_q(std::uint64_t x, T q, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t /*ubound*/, const policies::discrete_quantile<policies::integer_round_down>&)
{
if((q * fudge_factor > cum) && (x != lbound))
{
BOOST_MATH_INSTRUMENT_VARIABLE(x-1);
return --x;
}
return x;
}
template <class T>
inline std::uint64_t round_x_from_q(std::uint64_t x, T q, T cum, T fudge_factor, std::uint64_t /*lbound*/, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_up>&)
{
if((q < cum * fudge_factor) && (x != ubound))
{
BOOST_MATH_INSTRUMENT_VARIABLE(x+1);
return ++x;
}
return x;
}
template <class T>
inline std::uint64_t round_x_from_q(std::uint64_t x, T q, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_inwards>&)
{
if(q < 0.5)
return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
}
template <class T>
inline std::uint64_t round_x_from_q(std::uint64_t x, T q, T cum, T fudge_factor, std::uint64_t lbound, std::uint64_t ubound, const policies::discrete_quantile<policies::integer_round_outwards>&)
{
if(q >= 0.5)
return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_down>());
return round_x_from_q(x, q, cum, fudge_factor, lbound, ubound, policies::discrete_quantile<policies::integer_round_up>());
}
template <class T>
inline std::uint64_t round_x_from_q(std::uint64_t x, T /*q*/, T /*cum*/, T /*fudge_factor*/, std::uint64_t /*lbound*/, std::uint64_t /*ubound*/, const policies::discrete_quantile<policies::integer_round_nearest>&)
{
return x;
}
template <class T, class Policy>
std::uint64_t hypergeometric_quantile_imp(T p, T q, std::uint64_t r, std::uint64_t n, std::uint64_t N, const Policy& pol)
{
#ifdef _MSC_VER
# pragma warning(push)
# pragma warning(disable:4267)
#endif
typedef typename Policy::discrete_quantile_type discrete_quantile_type;
BOOST_MATH_STD_USING
BOOST_FPU_EXCEPTION_GUARD
T result;
T fudge_factor = 1 + tools::epsilon<T>() * ((N <= boost::math::prime(boost::math::max_prime - 1)) ? 50 : 2 * N);
std::uint64_t base = static_cast<std::uint64_t>((std::max)(0, static_cast<int>(n + r) - static_cast<int>(N)));
std::uint64_t lim = (std::min)(r, n);
BOOST_MATH_INSTRUMENT_VARIABLE(p);
BOOST_MATH_INSTRUMENT_VARIABLE(q);
BOOST_MATH_INSTRUMENT_VARIABLE(r);
BOOST_MATH_INSTRUMENT_VARIABLE(n);
BOOST_MATH_INSTRUMENT_VARIABLE(N);
BOOST_MATH_INSTRUMENT_VARIABLE(fudge_factor);
BOOST_MATH_INSTRUMENT_VARIABLE(base);
BOOST_MATH_INSTRUMENT_VARIABLE(lim);
if(p <= 0.5)
{
std::uint64_t x = base;
result = hypergeometric_pdf<T>(x, r, n, N, pol);
T diff = result;
if (diff == 0)
{
++x;
// We want to skip through x values as fast as we can until we start getting non-zero values,
// otherwise we're just making lots of expensive PDF calls:
T log_pdf = boost::math::lgamma(static_cast<T>(n + 1), pol)
+ boost::math::lgamma(static_cast<T>(r + 1), pol)
+ boost::math::lgamma(static_cast<T>(N - n + 1), pol)
+ boost::math::lgamma(static_cast<T>(N - r + 1), pol)
- boost::math::lgamma(static_cast<T>(N + 1), pol)
- boost::math::lgamma(static_cast<T>(x + 1), pol)
- boost::math::lgamma(static_cast<T>(n - x + 1), pol)
- boost::math::lgamma(static_cast<T>(r - x + 1), pol)
- boost::math::lgamma(static_cast<T>(N - n - r + x + 1), pol);
while (log_pdf < tools::log_min_value<T>())
{
log_pdf += -log(static_cast<T>(x + 1)) + log(static_cast<T>(n - x)) + log(static_cast<T>(r - x)) - log(static_cast<T>(N - n - r + x + 1));
++x;
}
// By the time we get here, log_pdf may be fairly inaccurate due to
// roundoff errors, get a fresh PDF calculation before proceeding:
diff = hypergeometric_pdf<T>(x, r, n, N, pol);
}
while(result < p)
{
diff = (diff > tools::min_value<T>() * 8)
? T(n - x) * T(r - x) * diff / (T(x + 1) * T(N + x + 1 - n - r))
: hypergeometric_pdf<T>(x + 1, r, n, N, pol);
if(result + diff / 2 > p)
break;
++x;
result += diff;
#ifdef BOOST_MATH_INSTRUMENT
if(diff != 0)
{
BOOST_MATH_INSTRUMENT_VARIABLE(x);
BOOST_MATH_INSTRUMENT_VARIABLE(diff);
BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
#endif
}
return round_x_from_p(x, p, result, fudge_factor, base, lim, discrete_quantile_type());
}
else
{
std::uint64_t x = lim;
result = 0;
T diff = hypergeometric_pdf<T>(x, r, n, N, pol);
if (diff == 0)
{
// We want to skip through x values as fast as we can until we start getting non-zero values,
// otherwise we're just making lots of expensive PDF calls:
--x;
T log_pdf = boost::math::lgamma(static_cast<T>(n + 1), pol)
+ boost::math::lgamma(static_cast<T>(r + 1), pol)
+ boost::math::lgamma(static_cast<T>(N - n + 1), pol)
+ boost::math::lgamma(static_cast<T>(N - r + 1), pol)
- boost::math::lgamma(static_cast<T>(N + 1), pol)
- boost::math::lgamma(static_cast<T>(x + 1), pol)
- boost::math::lgamma(static_cast<T>(n - x + 1), pol)
- boost::math::lgamma(static_cast<T>(r - x + 1), pol)
- boost::math::lgamma(static_cast<T>(N - n - r + x + 1), pol);
while (log_pdf < tools::log_min_value<T>())
{
log_pdf += log(static_cast<T>(x)) - log(static_cast<T>(n - x + 1)) - log(static_cast<T>(r - x + 1)) + log(static_cast<T>(N - n - r + x));
--x;
}
// By the time we get here, log_pdf may be fairly inaccurate due to
// roundoff errors, get a fresh PDF calculation before proceeding:
diff = hypergeometric_pdf<T>(x, r, n, N, pol);
}
while(result + diff / 2 < q)
{
result += diff;
diff = (diff > tools::min_value<T>() * 8)
? x * T(N + x - n - r) * diff / (T(1 + n - x) * T(1 + r - x))
: hypergeometric_pdf<T>(x - 1, r, n, N, pol);
--x;
#ifdef BOOST_MATH_INSTRUMENT
if(diff != 0)
{
BOOST_MATH_INSTRUMENT_VARIABLE(x);
BOOST_MATH_INSTRUMENT_VARIABLE(diff);
BOOST_MATH_INSTRUMENT_VARIABLE(result);
}
#endif
}
return round_x_from_q(x, q, result, fudge_factor, base, lim, discrete_quantile_type());
}
#ifdef _MSC_VER
# pragma warning(pop)
#endif
}
template <class T, class Policy>
inline std::uint64_t hypergeometric_quantile(T p, T q, std::uint64_t r, std::uint64_t n, std::uint64_t N, const Policy&)
{
BOOST_FPU_EXCEPTION_GUARD
typedef typename tools::promote_args<T>::type result_type;
typedef typename policies::evaluation<result_type, Policy>::type value_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::assert_undefined<> >::type forwarding_policy;
return detail::hypergeometric_quantile_imp<value_type>(p, q, r, n, N, forwarding_policy());
}
}}} // namespaces
#endif
Reference in New Issue View Git Blame Copy Permalink