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