| // 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> |
| |
| 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 |
| |
| 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 + 0.5f, |
| r + 0.5f, |
| N - n + 0.5f, |
| N - r + 0.5f, |
| N + 0.5f, |
| x + 0.5f, |
| n - x + 0.5f, |
| r - x + 0.5f, |
| N - n - r + x + 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 |
| }; |
| std::sort(sorted_indexes, sorted_indexes + 9, sort_functor<T>(exponents)); |
| |
| 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)); |
| }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_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 |
| |