distributions/detail/hypergeometric_pdf.hpp

0001 // Copyright 2008 Gautam Sewani
0002 // Copyright 2008 John Maddock
0003 //
0004 // Use, modification and distribution are subject to the
0005 // Boost Software License, Version 1.0.
0006 // (See accompanying file LICENSE_1_0.txt
0007 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0008
0009 #ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_PDF_HPP
0010 #define BOOST_MATH_DISTRIBUTIONS_DETAIL_HG_PDF_HPP
0011
0012 #include <boost/math/constants/constants.hpp>
0013 #include <boost/math/special_functions/lanczos.hpp>
0014 #include <boost/math/special_functions/gamma.hpp>
0015 #include <boost/math/special_functions/pow.hpp>
0016 #include <boost/math/special_functions/prime.hpp>
0017 #include <boost/math/policies/error_handling.hpp>
0018 #include <algorithm>
0019 #include <cstdint>
0020
0021 #ifdef BOOST_MATH_INSTRUMENT
0022 #include <typeinfo>
0023 #endif
0024
0025 namespace boost{ namespace math{ namespace detail{
0026
0027 template <class T, class Func>
0028 void bubble_down_one(T* first, T* last, Func f)
0029 {
0030    using std::swap;
0031    T* next = first;
0032    ++next;
0033    while((next != last) && (!f(*first, *next)))
0034    {
0035       swap(*first, *next);
0036       ++first;
0037       ++next;
0038    }
0039 }
0040
0041 template <class T>
0042 struct sort_functor
0043 {
0044    sort_functor(const T* exponents) : m_exponents(exponents){}
0045    bool operator()(std::size_t i, std::size_t j)
0046    {
0047       return m_exponents[i] > m_exponents[j];
0048    }
0049 private:
0050    const T* m_exponents;
0051 };
0052
0053 template <class T, class Lanczos, class Policy>
0054 T hypergeometric_pdf_lanczos_imp(T /*dummy*/, std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, const Lanczos&, const Policy&)
0055 {
0056    BOOST_MATH_STD_USING
0057
0058    BOOST_MATH_INSTRUMENT_FPU
0059    BOOST_MATH_INSTRUMENT_VARIABLE(x);
0060    BOOST_MATH_INSTRUMENT_VARIABLE(r);
0061    BOOST_MATH_INSTRUMENT_VARIABLE(n);
0062    BOOST_MATH_INSTRUMENT_VARIABLE(N);
0063    BOOST_MATH_INSTRUMENT_VARIABLE(typeid(Lanczos).name());
0064
0065    T bases[9] = {
0066       T(n) + static_cast<T>(Lanczos::g()) + 0.5f,
0067       T(r) + static_cast<T>(Lanczos::g()) + 0.5f,
0068       T(N - n) + static_cast<T>(Lanczos::g()) + 0.5f,
0069       T(N - r) + static_cast<T>(Lanczos::g()) + 0.5f,
0070       1 / (T(N) + static_cast<T>(Lanczos::g()) + 0.5f),
0071       1 / (T(x) + static_cast<T>(Lanczos::g()) + 0.5f),
0072       1 / (T(n - x) + static_cast<T>(Lanczos::g()) + 0.5f),
0073       1 / (T(r - x) + static_cast<T>(Lanczos::g()) + 0.5f),
0074       1 / (T(N - n - r + x) + static_cast<T>(Lanczos::g()) + 0.5f)
0075    };
0076    T exponents[9] = {
0077       n + T(0.5f),
0078       r + T(0.5f),
0079       N - n + T(0.5f),
0080       N - r + T(0.5f),
0081       N + T(0.5f),
0082       x + T(0.5f),
0083       n - x + T(0.5f),
0084       r - x + T(0.5f),
0085       N - n - r + x + T(0.5f)
0086    };
0087    int base_e_factors[9] = {
0088       -1, -1, -1, -1, 1, 1, 1, 1, 1
0089    };
0090    int sorted_indexes[9] = {
0091       0, 1, 2, 3, 4, 5, 6, 7, 8
0092    };
0093 #ifdef BOOST_MATH_INSTRUMENT
0094    BOOST_MATH_INSTRUMENT_FPU
0095    for(unsigned i = 0; i < 9; ++i)
0096    {
0097       BOOST_MATH_INSTRUMENT_VARIABLE(i);
0098       BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]);
0099       BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]);
0100       BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]);
0101       BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]);
0102    }
0103 #endif
0104    std::sort(sorted_indexes, sorted_indexes + 9, sort_functor<T>(exponents));
0105 #ifdef BOOST_MATH_INSTRUMENT
0106    BOOST_MATH_INSTRUMENT_FPU
0107    for(unsigned i = 0; i < 9; ++i)
0108    {
0109       BOOST_MATH_INSTRUMENT_VARIABLE(i);
0110       BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]);
0111       BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]);
0112       BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]);
0113       BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]);
0114    }
0115 #endif
0116
0117    do{
0118       exponents[sorted_indexes[0]] -= exponents[sorted_indexes[1]];
0119       bases[sorted_indexes[1]] *= bases[sorted_indexes[0]];
0120       if((bases[sorted_indexes[1]] < tools::min_value<T>()) && (exponents[sorted_indexes[1]] != 0))
0121       {
0122          return 0;
0123       }
0124       base_e_factors[sorted_indexes[1]] += base_e_factors[sorted_indexes[0]];
0125       bubble_down_one(sorted_indexes, sorted_indexes + 9, sort_functor<T>(exponents));
0126
0127 #ifdef BOOST_MATH_INSTRUMENT
0128       for(unsigned i = 0; i < 9; ++i)
0129       {
0130          BOOST_MATH_INSTRUMENT_VARIABLE(i);
0131          BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]);
0132          BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]);
0133          BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]);
0134          BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]);
0135       }
0136 #endif
0137    }while(exponents[sorted_indexes[1]] > 1);
0138
0139    //
0140    // Combine equal powers:
0141    //
0142    std::size_t j = 8;
0143    while(exponents[sorted_indexes[j]] == 0) --j;
0144    while(j)
0145    {
0146       while(j && (exponents[sorted_indexes[j-1]] == exponents[sorted_indexes[j]]))
0147       {
0148          bases[sorted_indexes[j-1]] *= bases[sorted_indexes[j]];
0149          exponents[sorted_indexes[j]] = 0;
0150          base_e_factors[sorted_indexes[j-1]] += base_e_factors[sorted_indexes[j]];
0151          bubble_down_one(sorted_indexes + j, sorted_indexes + 9, sort_functor<T>(exponents));
0152          --j;
0153       }
0154       --j;
0155
0156 #ifdef BOOST_MATH_INSTRUMENT
0157       BOOST_MATH_INSTRUMENT_VARIABLE(j);
0158       for(unsigned i = 0; i < 9; ++i)
0159       {
0160          BOOST_MATH_INSTRUMENT_VARIABLE(i);
0161          BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]);
0162          BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]);
0163          BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]);
0164          BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]);
0165       }
0166 #endif
0167    }
0168
0169 #ifdef BOOST_MATH_INSTRUMENT
0170    BOOST_MATH_INSTRUMENT_FPU
0171    for(unsigned i = 0; i < 9; ++i)
0172    {
0173       BOOST_MATH_INSTRUMENT_VARIABLE(i);
0174       BOOST_MATH_INSTRUMENT_VARIABLE(bases[i]);
0175       BOOST_MATH_INSTRUMENT_VARIABLE(exponents[i]);
0176       BOOST_MATH_INSTRUMENT_VARIABLE(base_e_factors[i]);
0177       BOOST_MATH_INSTRUMENT_VARIABLE(sorted_indexes[i]);
0178    }
0179 #endif
0180
0181    T result;
0182    BOOST_MATH_INSTRUMENT_VARIABLE(bases[sorted_indexes[0]] * exp(static_cast<T>(base_e_factors[sorted_indexes[0]])));
0183    BOOST_MATH_INSTRUMENT_VARIABLE(exponents[sorted_indexes[0]]);
0184    {
0185       BOOST_FPU_EXCEPTION_GUARD
0186       result = pow(bases[sorted_indexes[0]] * exp(static_cast<T>(base_e_factors[sorted_indexes[0]])), exponents[sorted_indexes[0]]);
0187    }
0188    BOOST_MATH_INSTRUMENT_VARIABLE(result);
0189    for(std::size_t i = 1; (i < 9) && (exponents[sorted_indexes[i]] > 0); ++i)
0190    {
0191       BOOST_FPU_EXCEPTION_GUARD
0192       if(result < tools::min_value<T>())
0193          return 0; // short circuit further evaluation
0194       if(exponents[sorted_indexes[i]] == 1)
0195          result *= bases[sorted_indexes[i]] * exp(static_cast<T>(base_e_factors[sorted_indexes[i]]));
0196       else if(exponents[sorted_indexes[i]] == 0.5f)
0197          result *= sqrt(bases[sorted_indexes[i]] * exp(static_cast<T>(base_e_factors[sorted_indexes[i]])));
0198       else
0199          result *= pow(bases[sorted_indexes[i]] * exp(static_cast<T>(base_e_factors[sorted_indexes[i]])), exponents[sorted_indexes[i]]);
0200
0201       BOOST_MATH_INSTRUMENT_VARIABLE(result);
0202    }
0203
0204    result *= Lanczos::lanczos_sum_expG_scaled(static_cast<T>(n + 1))
0205       * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(r + 1))
0206       * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(N - n + 1))
0207       * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(N - r + 1))
0208       /
0209       ( Lanczos::lanczos_sum_expG_scaled(static_cast<T>(N + 1))
0210          * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(x + 1))
0211          * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(n - x + 1))
0212          * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(r - x + 1))
0213          * Lanczos::lanczos_sum_expG_scaled(static_cast<T>(N - n - r + x + 1)));
0214
0215    BOOST_MATH_INSTRUMENT_VARIABLE(result);
0216    return result;
0217 }
0218
0219 template <class T, class Policy>
0220 T hypergeometric_pdf_lanczos_imp(T /*dummy*/, std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, const boost::math::lanczos::undefined_lanczos&, const Policy& pol)
0221 {
0222    BOOST_MATH_STD_USING
0223    return exp(
0224       boost::math::lgamma(T(n + 1), pol)
0225       + boost::math::lgamma(T(r + 1), pol)
0226       + boost::math::lgamma(T(N - n + 1), pol)
0227       + boost::math::lgamma(T(N - r + 1), pol)
0228       - boost::math::lgamma(T(N + 1), pol)
0229       - boost::math::lgamma(T(x + 1), pol)
0230       - boost::math::lgamma(T(n - x + 1), pol)
0231       - boost::math::lgamma(T(r - x + 1), pol)
0232       - boost::math::lgamma(T(N - n - r + x + 1), pol));
0233 }
0234
0235 template <class T>
0236 inline T integer_power(const T& x, int ex)
0237 {
0238    if(ex < 0)
0239       return 1 / integer_power(x, -ex);
0240    switch(ex)
0241    {
0242    case 0:
0243       return 1;
0244    case 1:
0245       return x;
0246    case 2:
0247       return x * x;
0248    case 3:
0249       return x * x * x;
0250    case 4:
0251       return boost::math::pow<4>(x);
0252    case 5:
0253       return boost::math::pow<5>(x);
0254    case 6:
0255       return boost::math::pow<6>(x);
0256    case 7:
0257       return boost::math::pow<7>(x);
0258    case 8:
0259       return boost::math::pow<8>(x);
0260    }
0261    BOOST_MATH_STD_USING
0262 #ifdef __SUNPRO_CC
0263    return pow(x, T(ex));
0264 #else
0265    return static_cast<T>(pow(x, ex));
0266 #endif
0267 }
0268 template <class T>
0269 struct hypergeometric_pdf_prime_loop_result_entry
0270 {
0271    T value;
0272    const hypergeometric_pdf_prime_loop_result_entry* next;
0273 };
0274
0275 #ifdef _MSC_VER
0276 #pragma warning(push)
0277 #pragma warning(disable:4510 4512 4610)
0278 #endif
0279
0280 struct hypergeometric_pdf_prime_loop_data
0281 {
0282    const std::uint64_t x;
0283    const std::uint64_t r;
0284    const std::uint64_t n;
0285    const std::uint64_t N;
0286    std::size_t prime_index;
0287    std::uint64_t current_prime;
0288 };
0289
0290 #ifdef _MSC_VER
0291 #pragma warning(pop)
0292 #endif
0293
0294 template <class T>
0295 T hypergeometric_pdf_prime_loop_imp(hypergeometric_pdf_prime_loop_data& data, hypergeometric_pdf_prime_loop_result_entry<T>& result)
0296 {
0297    while(data.current_prime <= data.N)
0298    {
0299       std::uint64_t base = data.current_prime;
0300       std::int64_t prime_powers = 0;
0301       while(base <= data.N)
0302       {
0303          prime_powers += data.n / base;
0304          prime_powers += data.r / base;
0305          prime_powers += (data.N - data.n) / base;
0306          prime_powers += (data.N - data.r) / base;
0307          prime_powers -= data.N / base;
0308          prime_powers -= data.x / base;
0309          prime_powers -= (data.n - data.x) / base;
0310          prime_powers -= (data.r - data.x) / base;
0311          prime_powers -= (data.N - data.n - data.r + data.x) / base;
0312          base *= data.current_prime;
0313       }
0314       if(prime_powers)
0315       {
0316          T p = integer_power<T>(static_cast<T>(data.current_prime), prime_powers);
0317          if((p > 1) && (tools::max_value<T>() / p < result.value))
0318          {
0319             //
0320             // The next calculation would overflow, use recursion
0321             // to sidestep the issue:
0322             //
0323             hypergeometric_pdf_prime_loop_result_entry<T> t = { p, &result };
0324             data.current_prime = prime(++data.prime_index);
0325             return hypergeometric_pdf_prime_loop_imp<T>(data, t);
0326          }
0327          if((p < 1) && (tools::min_value<T>() / p > result.value))
0328          {
0329             //
0330             // The next calculation would underflow, use recursion
0331             // to sidestep the issue:
0332             //
0333             hypergeometric_pdf_prime_loop_result_entry<T> t = { p, &result };
0334             data.current_prime = prime(++data.prime_index);
0335             return hypergeometric_pdf_prime_loop_imp<T>(data, t);
0336          }
0337          result.value *= p;
0338       }
0339       data.current_prime = prime(++data.prime_index);
0340    }
0341    //
0342    // When we get to here we have run out of prime factors,
0343    // the overall result is the product of all the partial
0344    // results we have accumulated on the stack so far, these
0345    // are in a linked list starting with "data.head" and ending
0346    // with "result".
0347    //
0348    // All that remains is to multiply them together, taking
0349    // care not to overflow or underflow.
0350    //
0351    // Enumerate partial results >= 1 in variable i
0352    // and partial results < 1 in variable j:
0353    //
0354    hypergeometric_pdf_prime_loop_result_entry<T> const *i, *j;
0355    i = &result;
0356    while(i && i->value < 1)
0357       i = i->next;
0358    j = &result;
0359    while(j && j->value >= 1)
0360       j = j->next;
0361
0362    T prod = 1;
0363
0364    while(i || j)
0365    {
0366       while(i && ((prod <= 1) || (j == 0)))
0367       {
0368          prod *= i->value;
0369          i = i->next;
0370          while(i && i->value < 1)
0371             i = i->next;
0372       }
0373       while(j && ((prod >= 1) || (i == 0)))
0374       {
0375          prod *= j->value;
0376          j = j->next;
0377          while(j && j->value >= 1)
0378             j = j->next;
0379       }
0380    }
0381
0382    return prod;
0383 }
0384
0385 template <class T, class Policy>
0386 inline T hypergeometric_pdf_prime_imp(std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, const Policy&)
0387 {
0388    hypergeometric_pdf_prime_loop_result_entry<T> result = { 1, 0 };
0389    hypergeometric_pdf_prime_loop_data data = { x, r, n, N, 0, prime(0) };
0390    return hypergeometric_pdf_prime_loop_imp<T>(data, result);
0391 }
0392
0393 template <class T, class Policy>
0394 T hypergeometric_pdf_factorial_imp(std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, const Policy&)
0395 {
0396    BOOST_MATH_STD_USING
0397    BOOST_MATH_ASSERT(N <= boost::math::max_factorial<T>::value);
0398    T result = boost::math::unchecked_factorial<T>(n);
0399    T num[3] = {
0400       boost::math::unchecked_factorial<T>(r),
0401       boost::math::unchecked_factorial<T>(N - n),
0402       boost::math::unchecked_factorial<T>(N - r)
0403    };
0404    T denom[5] = {
0405       boost::math::unchecked_factorial<T>(N),
0406       boost::math::unchecked_factorial<T>(x),
0407       boost::math::unchecked_factorial<T>(n - x),
0408       boost::math::unchecked_factorial<T>(r - x),
0409       boost::math::unchecked_factorial<T>(N - n - r + x)
0410    };
0411    std::size_t i = 0;
0412    std::size_t j = 0;
0413    while((i < 3) || (j < 5))
0414    {
0415       while((j < 5) && ((result >= 1) || (i >= 3)))
0416       {
0417          result /= denom[j];
0418          ++j;
0419       }
0420       while((i < 3) && ((result <= 1) || (j >= 5)))
0421       {
0422          result *= num[i];
0423          ++i;
0424       }
0425    }
0426    return result;
0427 }
0428
0429
0430 template <class T, class Policy>
0431 inline typename tools::promote_args<T>::type
0432    hypergeometric_pdf(std::uint64_t x, std::uint64_t r, std::uint64_t n, std::uint64_t N, const Policy&)
0433 {
0434    BOOST_FPU_EXCEPTION_GUARD
0435    typedef typename tools::promote_args<T>::type result_type;
0436    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0437    typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
0438    typedef typename policies::normalise<
0439       Policy,
0440       policies::promote_float<false>,
0441       policies::promote_double<false>,
0442       policies::discrete_quantile<>,
0443       policies::assert_undefined<> >::type forwarding_policy;
0444
0445    value_type result;
0446    if(N <= boost::math::max_factorial<value_type>::value)
0447    {
0448       //
0449       // If N is small enough then we can evaluate the PDF via the factorials
0450       // directly: table lookup of the factorials gives the best performance
0451       // of the methods available:
0452       //
0453       result = detail::hypergeometric_pdf_factorial_imp<value_type>(x, r, n, N, forwarding_policy());
0454    }
0455    else if(N <= boost::math::prime(boost::math::max_prime - 1))
0456    {
0457       //
0458       // If N is no larger than the largest prime number in our lookup table
0459       // (104729) then we can use prime factorisation to evaluate the PDF,
0460       // this is slow but accurate:
0461       //
0462       result = detail::hypergeometric_pdf_prime_imp<value_type>(x, r, n, N, forwarding_policy());
0463    }
0464    else
0465    {
0466       //
0467       // Catch all case - use the lanczos approximation - where available -
0468       // to evaluate the ratio of factorials.  This is reasonably fast
0469       // (almost as quick as using logarithmic evaluation in terms of lgamma)
0470       // but only a few digits better in accuracy than using lgamma:
0471       //
0472       result = detail::hypergeometric_pdf_lanczos_imp(value_type(), x, r, n, N, evaluation_type(), forwarding_policy());
0473    }
0474
0475    if(result > 1)
0476    {
0477       result = 1;
0478    }
0479    if(result < 0)
0480    {
0481       result = 0;
0482    }
0483
0484    return policies::checked_narrowing_cast<result_type, forwarding_policy>(result, "boost::math::hypergeometric_pdf<%1%>(%1%,%1%,%1%,%1%)");
0485 }
0486
0487 }}} // namespaces
0488
0489 #endif
0490