math/special_functions/binomial.hpp

0001 //  Copyright John Maddock 2006.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_SF_BINOMIAL_HPP
0007 #define BOOST_MATH_SF_BINOMIAL_HPP
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #endif
0012
0013 #include <boost/math/tools/config.hpp>
0014 #include <boost/math/tools/type_traits.hpp>
0015 #include <boost/math/special_functions/math_fwd.hpp>
0016 #include <boost/math/special_functions/factorials.hpp>
0017 #include <boost/math/special_functions/beta.hpp>
0018 #include <boost/math/policies/error_handling.hpp>
0019
0020 namespace boost{ namespace math{
0021
0022 template <class T, class Policy>
0023 BOOST_MATH_GPU_ENABLED T binomial_coefficient(unsigned n, unsigned k, const Policy& pol)
0024 {
0025    static_assert(!boost::math::is_integral<T>::value, "Type T must not be an integral type");
0026    BOOST_MATH_STD_USING
0027    constexpr auto function = "boost::math::binomial_coefficient<%1%>(unsigned, unsigned)";
0028    if(k > n)
0029       return policies::raise_domain_error<T>(function, "The binomial coefficient is undefined for k > n, but got k = %1%.", static_cast<T>(k), pol);
0030    T result;  // LCOV_EXCL_LINE
0031    if((k == 0) || (k == n))
0032       return static_cast<T>(1);
0033    if((k == 1) || (k == n-1))
0034       return static_cast<T>(n);
0035
0036    if(n <= max_factorial<T>::value)
0037    {
0038       // Use fast table lookup:
0039       result = unchecked_factorial<T>(n);
0040       result /= unchecked_factorial<T>(n-k);
0041       result /= unchecked_factorial<T>(k);
0042    }
0043    else
0044    {
0045       // Use the beta function:
0046       if(k < n - k)
0047          result = static_cast<T>(k * boost::math::beta(static_cast<T>(k), static_cast<T>(n-k+1), pol));
0048       else
0049          result = static_cast<T>((n - k) * boost::math::beta(static_cast<T>(k+1), static_cast<T>(n-k), pol));
0050       if(result == 0)
0051          return policies::raise_overflow_error<T>(function, nullptr, pol);
0052       result = 1 / result;
0053    }
0054    // convert to nearest integer:
0055    return ceil(result - 0.5f);
0056 }
0057 //
0058 // Type float can only store the first 35 factorials, in order to
0059 // increase the chance that we can use a table driven implementation
0060 // we'll promote to double:
0061 //
0062 template <>
0063 BOOST_MATH_GPU_ENABLED inline float binomial_coefficient<float, policies::policy<> >(unsigned n, unsigned k, const policies::policy<>&)
0064 {
0065    typedef policies::normalise<
0066        policies::policy<>,
0067        policies::promote_float<true>,
0068        policies::promote_double<false>,
0069        policies::discrete_quantile<>,
0070        policies::assert_undefined<> >::type forwarding_policy;
0071    return policies::checked_narrowing_cast<float, forwarding_policy>(binomial_coefficient<double>(n, k, forwarding_policy()), "boost::math::binomial_coefficient<%1%>(unsigned,unsigned)");
0072 }
0073
0074 template <class T>
0075 BOOST_MATH_GPU_ENABLED inline T binomial_coefficient(unsigned n, unsigned k)
0076 {
0077    return binomial_coefficient<T>(n, k, policies::policy<>());
0078 }
0079
0080 } // namespace math
0081 } // namespace boost
0082
0083
0084 #endif // BOOST_MATH_SF_BINOMIAL_HPP
0085
0086
0087