Chris@16: //  Copyright John Maddock 2006.
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_MATH_SF_BINOMIAL_HPP
Chris@16: #define BOOST_MATH_SF_BINOMIAL_HPP
Chris@16: 
Chris@16: #ifdef _MSC_VER
Chris@16: #pragma once
Chris@16: #endif
Chris@16: 
Chris@101: #include <boost/math/special_functions/math_fwd.hpp>
Chris@16: #include <boost/math/special_functions/factorials.hpp>
Chris@16: #include <boost/math/special_functions/beta.hpp>
Chris@16: #include <boost/math/policies/error_handling.hpp>
Chris@16: 
Chris@16: namespace boost{ namespace math{
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: T binomial_coefficient(unsigned n, unsigned k, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_STATIC_ASSERT(!boost::is_integral<T>::value);
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    static const char* function = "boost::math::binomial_coefficient<%1%>(unsigned, unsigned)";
Chris@16:    if(k > n)
Chris@16:       return policies::raise_domain_error<T>(
Chris@16:          function, 
Chris@16:          "The binomial coefficient is undefined for k > n, but got k = %1%.",
Chris@101:          static_cast<T>(k), pol);
Chris@16:    T result;
Chris@16:    if((k == 0) || (k == n))
Chris@101:       return static_cast<T>(1);
Chris@16:    if((k == 1) || (k == n-1))
Chris@101:       return static_cast<T>(n);
Chris@16: 
Chris@16:    if(n <= max_factorial<T>::value)
Chris@16:    {
Chris@16:       // Use fast table lookup:
Chris@16:       result = unchecked_factorial<T>(n);
Chris@16:       result /= unchecked_factorial<T>(n-k);
Chris@16:       result /= unchecked_factorial<T>(k);
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       // Use the beta function:
Chris@16:       if(k < n - k)
Chris@16:          result = k * beta(static_cast<T>(k), static_cast<T>(n-k+1), pol);
Chris@16:       else
Chris@16:          result = (n - k) * beta(static_cast<T>(k+1), static_cast<T>(n-k), pol);
Chris@16:       if(result == 0)
Chris@16:          return policies::raise_overflow_error<T>(function, 0, pol);
Chris@16:       result = 1 / result;
Chris@16:    }
Chris@16:    // convert to nearest integer:
Chris@16:    return ceil(result - 0.5f);
Chris@16: }
Chris@16: //
Chris@16: // Type float can only store the first 35 factorials, in order to
Chris@16: // increase the chance that we can use a table driven implementation
Chris@16: // we'll promote to double:
Chris@16: //
Chris@16: template <>
Chris@16: inline float binomial_coefficient<float, policies::policy<> >(unsigned n, unsigned k, const policies::policy<>& pol)
Chris@16: {
Chris@16:    return policies::checked_narrowing_cast<float, policies::policy<> >(binomial_coefficient<double>(n, k, pol), "boost::math::binomial_coefficient<%1%>(unsigned,unsigned)");
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T binomial_coefficient(unsigned n, unsigned k)
Chris@16: {
Chris@16:    return binomial_coefficient<T>(n, k, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: } // namespace math
Chris@16: } // namespace boost
Chris@16: 
Chris@16: 
Chris@16: #endif // BOOST_MATH_SF_BINOMIAL_HPP
Chris@16: 
Chris@16: 
Chris@16: