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/special_functions/math_fwd.hpp>
0014 #include <boost/math/special_functions/factorials.hpp>
0015 #include <boost/math/special_functions/beta.hpp>
0016 #include <boost/math/policies/error_handling.hpp>
0017 #include <type_traits>
0018
0019 namespace boost{ namespace math{
0020
0021 template <class T, class Policy>
0022 T binomial_coefficient(unsigned n, unsigned k, const Policy& pol)
0023 {
0024    static_assert(!std::is_integral<T>::value, "Type T must not be an integral type");
0025    BOOST_MATH_STD_USING
0026    static const char* function = "boost::math::binomial_coefficient<%1%>(unsigned, unsigned)";
0027    if(k > n)
0028       return policies::raise_domain_error<T>(function, "The binomial coefficient is undefined for k > n, but got k = %1%.", static_cast<T>(k), pol);
0029    T result;  // LCOV_EXCL_LINE
0030    if((k == 0) || (k == n))
0031       return static_cast<T>(1);
0032    if((k == 1) || (k == n-1))
0033       return static_cast<T>(n);
0034
0035    if(n <= max_factorial<T>::value)
0036    {
0037       // Use fast table lookup:
0038       result = unchecked_factorial<T>(n);
0039       result /= unchecked_factorial<T>(n-k);
0040       result /= unchecked_factorial<T>(k);
0041    }
0042    else
0043    {
0044       // Use the beta function:
0045       if(k < n - k)
0046          result = static_cast<T>(k * beta(static_cast<T>(k), static_cast<T>(n-k+1), pol));
0047       else
0048          result = static_cast<T>((n - k) * beta(static_cast<T>(k+1), static_cast<T>(n-k), pol));
0049       if(result == 0)
0050          return policies::raise_overflow_error<T>(function, nullptr, pol);
0051       result = 1 / result;
0052    }
0053    // convert to nearest integer:
0054    return ceil(result - 0.5f);
0055 }
0056 //
0057 // Type float can only store the first 35 factorials, in order to
0058 // increase the chance that we can use a table driven implementation
0059 // we'll promote to double:
0060 //
0061 template <>
0062 inline float binomial_coefficient<float, policies::policy<> >(unsigned n, unsigned k, const policies::policy<>&)
0063 {
0064    typedef policies::normalise<
0065        policies::policy<>,
0066        policies::promote_float<true>,
0067        policies::promote_double<false>,
0068        policies::discrete_quantile<>,
0069        policies::assert_undefined<> >::type forwarding_policy;
0070    return policies::checked_narrowing_cast<float, forwarding_policy>(binomial_coefficient<double>(n, k, forwarding_policy()), "boost::math::binomial_coefficient<%1%>(unsigned,unsigned)");
0071 }
0072
0073 template <class T>
0074 inline T binomial_coefficient(unsigned n, unsigned k)
0075 {
0076    return binomial_coefficient<T>(n, k, policies::policy<>());
0077 }
0078
0079 } // namespace math
0080 } // namespace boost
0081
0082
0083 #endif // BOOST_MATH_SF_BINOMIAL_HPP
0084
0085
0086