// Copyright 2008 Gautam Sewani | |
// 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_PDF_HPP | |
#define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_PDF_HPP | |
#include <boost/math/constants/constants.hpp> | |
#include <boost/math/special_functions/lanczos.hpp> | |
#include <boost/math/special_functions/gamma.hpp> | |
#include <boost/math/special_functions/pow.hpp> | |
#include <boost/math/special_functions/prime.hpp> | |
#include <boost/math/policies/error_handling.hpp> | |
#ifdef BOOST_MATH_INSTRUMENT | |
#include <typeinfo> | |
#endif | |
namespace boost{ namespace math{ namespace detail{ | |
template <class T, class Func> | |
void bubble_down_one(T* first, T* last, Func f) | |
{ | |
using std::swap; | |
T* next = first; | |
++next; | |
while((next != last) && (!f(*first, *next))) | |
{ | |
swap(*first, *next); | |
++first; | |
++next; | |
} | |
} | |
template <class T> | |
struct sort_functor | |
{ | |
sort_functor(const T* exponents) : m_exponents(exponents){} | |
bool operator()(int i, int j) | |
{ | |
return m_exponents[i] > m_exponents[j]; | |
} | |
private: | |
const T* m_exponents; | |
}; | |
template <class T, class Lanczos, class Policy> | |
T hypergeometric_pdf_lanczos_imp(T /*dummy*/, unsigned x, unsigned r, unsigned n, unsigned N, const Lanczos&, const Policy&) | |
{ | |
BOOST_MATH_STD_USING | |
BOOST_MATH_INSTRUMENT_FPU | |
BOOST_MATH_INSTRUMENT_VARIABLE(x); | |
BOOST_MATH_INSTRUMENT_VARIABLE(r); | |
BOOST_MATH_INSTRUMENT_VARIABLE(n); | |
BOOST_MATH_INSTRUMENT_VARIABLE(N); | |
BOOST_MATH_INSTRUMENT_VARIABLE(typeid(Lanczos).name()); | |
T bases[9] = { | |
T(n) + Lanczos::g() + 0.5f, | |
T(r) + Lanczos::g() + 0.5f, | |
T(N - n) + Lanczos::g() + 0.5f, | |
T(N - r) + Lanczos::g() + 0.5f, | |
1 / (T(N) + Lanczos::g() + 0.5f), | |
1 / (T(x) + Lanczos::g() + 0.5f), | |
1 / (T(n - x) + Lanczos::g() + 0.5f), | |
1 / (T(r - x) + Lanczos::g() + 0.5f), | |
1 / (T(N - n - r + x) + Lanczos::g() + 0.5f) | |
}; | |
T exponents[9] = { | |
n + T(0.5f), | |
r + T(0.5f), | |
N - n + T(0.5f), | |
N - r + T(0.5f), | |
N + T(0.5f), | |
x + T(0.5f), | |
n - x + T(0.5f), | |
r - x + T(0.5f), | |
N - n - r + x + T(0.5f) | |
}; | |
int base_e_factors[9] = { | |
-1, -1, -1, -1, 1, 1, 1, 1, 1 | |
}; | |
int sorted_indexes[9] = { | |
0, 1, 2, 3, 4, 5, 6, 7, 8 | |
}; | |
#ifdef BOOST_MATH_INSTRUMENT | |
BOOST_MATH_INSTRUMENT_FPU | |
for(unsigned i = 0; i < 9; ++i) | |
{ | |
BOOST_MATH_INSTRUMENT_VARIABLE(i); | |
BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]); | |
BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]); | |
BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]); | |
BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]); | |
} | |
#endif | |
std::sort(sorted_indexes, sorted_indexes + 9, sort_functor<T>(exponents)); | |
#ifdef BOOST_MATH_INSTRUMENT | |
BOOST_MATH_INSTRUMENT_FPU | |
for(unsigned i = 0; i < 9; ++i) | |
{ | |
BOOST_MATH_INSTRUMENT_VARIABLE(i); | |
BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]); | |
BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]); | |
BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]); | |
BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]); | |
} | |
#endif | |
do{ | |
exponents[sorted_indexes[0]] -= exponents[sorted_indexes[1]]; | |
bases[sorted_indexes[1]] *= bases[sorted_indexes[0]]; | |
if((bases[sorted_indexes[1]] < tools::min_value<T>()) && (exponents[sorted_indexes[1]] != 0)) | |
{ | |
return 0; | |
} | |
base_e_factors[sorted_indexes[1]] += base_e_factors[sorted_indexes[0]]; | |
bubble_down_one(sorted_indexes, sorted_indexes + 9, sort_functor<T>(exponents)); | |
#ifdef BOOST_MATH_INSTRUMENT | |
for(unsigned i = 0; i < 9; ++i) | |
{ | |
BOOST_MATH_INSTRUMENT_VARIABLE(i); | |
BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]); | |
BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]); | |
BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]); | |
BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]); | |
} | |
#endif | |
}while(exponents[sorted_indexes[1]] > 1); | |
// | |
// Combine equal powers: | |
// | |
int j = 8; | |
while(exponents[sorted_indexes[j]] == 0) --j; | |
while(j) | |
{ | |
while(j && (exponents[sorted_indexes[j-1]] == exponents[sorted_indexes[j]])) | |
{ | |
bases[sorted_indexes[j-1]] *= bases[sorted_indexes[j]]; | |
exponents[sorted_indexes[j]] = 0; | |
base_e_factors[sorted_indexes[j-1]] += base_e_factors[sorted_indexes[j]]; | |
bubble_down_one(sorted_indexes + j, sorted_indexes + 9, sort_functor<T>(exponents)); | |
--j; | |
} | |
--j; | |
#ifdef BOOST_MATH_INSTRUMENT | |
BOOST_MATH_INSTRUMENT_VARIABLE(j); | |
for(unsigned i = 0; i < 9; ++i) | |
{ | |
BOOST_MATH_INSTRUMENT_VARIABLE(i); | |
BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]); | |
BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]); | |
BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]); | |
BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]); | |
} | |
#endif | |
} | |
#ifdef BOOST_MATH_INSTRUMENT | |
BOOST_MATH_INSTRUMENT_FPU | |
for(unsigned i = 0; i < 9; ++i) | |
{ | |
BOOST_MATH_INSTRUMENT_VARIABLE(i); | |
BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]); | |
BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]); | |
BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]); | |
BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]); | |
} | |
#endif | |
T result; | |
BOOST_MATH_INSTRUMENT_VARIABLE(bases[sorted_indexes[0]] * exp(static_cast<T>(base_e_factors[sorted_indexes[0]]))); | |
BOOST_MATH_INSTRUMENT_VARIABLE(exponents[sorted_indexes[0]]); | |
{ | |
BOOST_FPU_EXCEPTION_GUARD | |
result = pow(bases[sorted_indexes[0]] * exp(static_cast<T>(base_e_factors[sorted_indexes[0]])), exponents[sorted_indexes[0]]); | |
} | |
BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
for(unsigned i = 1; (i < 9) && (exponents[sorted_indexes[i]] > 0); ++i) | |
{ | |
BOOST_FPU_EXCEPTION_GUARD | |
if(result < tools::min_value<T>()) | |
return 0; // short circuit further evaluation | |
if(exponents[sorted_indexes[i]] == 1) | |
result *= bases[sorted_indexes[i]] * exp(static_cast<T>(base_e_factors[sorted_indexes[i]])); | |
else if(exponents[sorted_indexes[i]] == 0.5f) | |
result *= sqrt(bases[sorted_indexes[i]] * exp(static_cast<T>(base_e_factors[sorted_indexes[i]]))); | |
else | |
result *= pow(bases[sorted_indexes[i]] * exp(static_cast<T>(base_e_factors[sorted_indexes[i]])), exponents[sorted_indexes[i]]); | |
BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
} | |
result *= Lanczos::lanczos_sum_expG_scaled(static_cast<T>(n + 1)) | |
* Lanczos::lanczos_sum_expG_scaled(static_cast<T>(r + 1)) | |
* Lanczos::lanczos_sum_expG_scaled(static_cast<T>(N - n + 1)) | |
* Lanczos::lanczos_sum_expG_scaled(static_cast<T>(N - r + 1)) | |
/ | |
( Lanczos::lanczos_sum_expG_scaled(static_cast<T>(N + 1)) | |
* Lanczos::lanczos_sum_expG_scaled(static_cast<T>(x + 1)) | |
* Lanczos::lanczos_sum_expG_scaled(static_cast<T>(n - x + 1)) | |
* Lanczos::lanczos_sum_expG_scaled(static_cast<T>(r - x + 1)) | |
* Lanczos::lanczos_sum_expG_scaled(static_cast<T>(N - n - r + x + 1))); | |
BOOST_MATH_INSTRUMENT_VARIABLE(result); | |
return result; | |
} | |
template <class T, class Policy> | |
T hypergeometric_pdf_lanczos_imp(T /*dummy*/, unsigned x, unsigned r, unsigned n, unsigned N, const boost::math::lanczos::undefined_lanczos&, const Policy& pol) | |
{ | |
BOOST_MATH_STD_USING | |
return exp( | |
boost::math::lgamma(T(n + 1), pol) | |
+ boost::math::lgamma(T(r + 1), pol) | |
+ boost::math::lgamma(T(N - n + 1), pol) | |
+ boost::math::lgamma(T(N - r + 1), pol) | |
- boost::math::lgamma(T(N + 1), pol) | |
- boost::math::lgamma(T(x + 1), pol) | |
- boost::math::lgamma(T(n - x + 1), pol) | |
- boost::math::lgamma(T(r - x + 1), pol) | |
- boost::math::lgamma(T(N - n - r + x + 1), pol)); | |
} | |
template <class T> | |
inline T integer_power(const T& x, int ex) | |
{ | |
if(ex < 0) | |
return 1 / integer_power(x, -ex); | |
switch(ex) | |
{ | |
case 0: | |
return 1; | |
case 1: | |
return x; | |
case 2: | |
return x * x; | |
case 3: | |
return x * x * x; | |
case 4: | |
return boost::math::pow<4>(x); | |
case 5: | |
return boost::math::pow<5>(x); | |
case 6: | |
return boost::math::pow<6>(x); | |
case 7: | |
return boost::math::pow<7>(x); | |
case 8: | |
return boost::math::pow<8>(x); | |
} | |
BOOST_MATH_STD_USING | |
#ifdef __SUNPRO_CC | |
return pow(x, T(ex)); | |
#else | |
return pow(x, ex); | |
#endif | |
} | |
template <class T> | |
struct hypergeometric_pdf_prime_loop_result_entry | |
{ | |
T value; | |
const hypergeometric_pdf_prime_loop_result_entry* next; | |
}; | |
#ifdef BOOST_MSVC | |
#pragma warning(push) | |
#pragma warning(disable:4510 4512 4610) | |
#endif | |
struct hypergeometric_pdf_prime_loop_data | |
{ | |
const unsigned x; | |
const unsigned r; | |
const unsigned n; | |
const unsigned N; | |
unsigned prime_index; | |
unsigned current_prime; | |
}; | |
#ifdef BOOST_MSVC | |
#pragma warning(pop) | |
#endif | |
template <class T> | |
T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hypergeometric_pdf_prime_loop_result_entry<T>& result) | |
{ | |
while(data.current_prime <= data.N) | |
{ | |
unsigned base = data.current_prime; | |
int prime_powers = 0; | |
while(base <= data.N) | |
{ | |
prime_powers += data.n / base; | |
prime_powers += data.r / base; | |
prime_powers += (data.N - data.n) / base; | |
prime_powers += (data.N - data.r) / base; | |
prime_powers -= data.N / base; | |
prime_powers -= data.x / base; | |
prime_powers -= (data.n - data.x) / base; | |
prime_powers -= (data.r - data.x) / base; | |
prime_powers -= (data.N - data.n - data.r + data.x) / base; | |
base *= data.current_prime; | |
} | |
if(prime_powers) | |
{ | |
T p = integer_power<T>(data.current_prime, prime_powers); | |
if((p > 1) && (tools::max_value<T>() / p < result.value)) | |
{ | |
// | |
// The next calculation would overflow, use recursion | |
// to sidestep the issue: | |
// | |
hypergeometric_pdf_prime_loop_result_entry<T> t = { p, &result }; | |
data.current_prime = prime(++data.prime_index); | |
return hypergeometric_pdf_prime_loop_imp<T>(data, t); | |
} | |
if((p < 1) && (tools::min_value<T>() / p > result.value)) | |
{ | |
// | |
// The next calculation would underflow, use recursion | |
// to sidestep the issue: | |
// | |
hypergeometric_pdf_prime_loop_result_entry<T> t = { p, &result }; | |
data.current_prime = prime(++data.prime_index); | |
return hypergeometric_pdf_prime_loop_imp<T>(data, t); | |
} | |
result.value *= p; | |
} | |
data.current_prime = prime(++data.prime_index); | |
} | |
// | |
// When we get to here we have run out of prime factors, | |
// the overall result is the product of all the partial | |
// results we have accumulated on the stack so far, these | |
// are in a linked list starting with "data.head" and ending | |
// with "result". | |
// | |
// All that remains is to multiply them together, taking | |
// care not to overflow or underflow. | |
// | |
// Enumerate partial results >= 1 in variable i | |
// and partial results < 1 in variable j: | |
// | |
hypergeometric_pdf_prime_loop_result_entry<T> const *i, *j; | |
i = &result; | |
while(i && i->value < 1) | |
i = i->next; | |
j = &result; | |
while(j && j->value >= 1) | |
j = j->next; | |
T prod = 1; | |
while(i || j) | |
{ | |
while(i && ((prod <= 1) || (j == 0))) | |
{ | |
prod *= i->value; | |
i = i->next; | |
while(i && i->value < 1) | |
i = i->next; | |
} | |
while(j && ((prod >= 1) || (i == 0))) | |
{ | |
prod *= j->value; | |
j = j->next; | |
while(j && j->value >= 1) | |
j = j->next; | |
} | |
} | |
return prod; | |
} | |
template <class T, class Policy> | |
inline T hypergeometric_pdf_prime_imp(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&) | |
{ | |
hypergeometric_pdf_prime_loop_result_entry<T> result = { 1, 0 }; | |
hypergeometric_pdf_prime_loop_data data = { x, r, n, N, 0, prime(0) }; | |
return hypergeometric_pdf_prime_loop_imp<T>(data, result); | |
} | |
template <class T, class Policy> | |
T hypergeometric_pdf_factorial_imp(unsigned x, unsigned r, unsigned n, unsigned N, const Policy&) | |
{ | |
BOOST_MATH_STD_USING | |
BOOST_ASSERT(N < boost::math::max_factorial<T>::value); | |
T result = boost::math::unchecked_factorial<T>(n); | |
T num[3] = { | |
boost::math::unchecked_factorial<T>(r), | |
boost::math::unchecked_factorial<T>(N - n), | |
boost::math::unchecked_factorial<T>(N - r) | |
}; | |
T denom[5] = { | |
boost::math::unchecked_factorial<T>(N), | |
boost::math::unchecked_factorial<T>(x), | |
boost::math::unchecked_factorial<T>(n - x), | |
boost::math::unchecked_factorial<T>(r - x), | |
boost::math::unchecked_factorial<T>(N - n - r + x) | |
}; | |
int i = 0; | |
int j = 0; | |
while((i < 3) || (j < 5)) | |
{ | |
while((j < 5) && ((result >= 1) || (i >= 3))) | |
{ | |
result /= denom[j]; | |
++j; | |
} | |
while((i < 3) && ((result <= 1) || (j >= 5))) | |
{ | |
result *= num[i]; | |
++i; | |
} | |
} | |
return result; | |
} | |
template <class T, class Policy> | |
inline typename tools::promote_args<T>::type | |
hypergeometric_pdf(unsigned x, unsigned r, unsigned n, unsigned 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 lanczos::lanczos<value_type, Policy>::type evaluation_type; | |
typedef typename policies::normalise< | |
Policy, | |
policies::promote_float<false>, | |
policies::promote_double<false>, | |
policies::discrete_quantile<>, | |
policies::assert_undefined<> >::type forwarding_policy; | |
value_type result; | |
if(N <= boost::math::max_factorial<value_type>::value) | |
{ | |
// | |
// If N is small enough then we can evaluate the PDF via the factorials | |
// directly: table lookup of the factorials gives the best performance | |
// of the methods available: | |
// | |
result = detail::hypergeometric_pdf_factorial_imp<value_type>(x, r, n, N, forwarding_policy()); | |
} | |
else if(N <= boost::math::prime(boost::math::max_prime - 1)) | |
{ | |
// | |
// If N is no larger than the largest prime number in our lookup table | |
// (104729) then we can use prime factorisation to evaluate the PDF, | |
// this is slow but accurate: | |
// | |
result = detail::hypergeometric_pdf_prime_imp<value_type>(x, r, n, N, forwarding_policy()); | |
} | |
else | |
{ | |
// | |
// Catch all case - use the lanczos approximation - where available - | |
// to evaluate the ratio of factorials. This is reasonably fast | |
// (almost as quick as using logarithmic evaluation in terms of lgamma) | |
// but only a few digits better in accuracy than using lgamma: | |
// | |
result = detail::hypergeometric_pdf_lanczos_imp(value_type(), x, r, n, N, evaluation_type(), forwarding_policy()); | |
} | |
if(result > 1) | |
{ | |
result = 1; | |
} | |
if(result < 0) | |
{ | |
result = 0; | |
} | |
return policies::checked_narrowing_cast<result_type, forwarding_policy>(result, "boost::math::hypergeometric_pdf<%1%>(%1%,%1%,%1%,%1%)"); | |
} | |
}}} // namespaces | |
#endif | |