1 // Copyright John Maddock 2006.
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
6 #ifndef BOOST_MATH_SF_BINOMIAL_HPP
7 #define BOOST_MATH_SF_BINOMIAL_HPP
13 #include <boost/math/special_functions/factorials.hpp>
14 #include <boost/math/special_functions/beta.hpp>
15 #include <boost/math/policies/error_handling.hpp>
17 namespace boost{ namespace math{
19 template <class T, class Policy>
20 T binomial_coefficient(unsigned n, unsigned k, const Policy& pol)
22 BOOST_STATIC_ASSERT(!boost::is_integral<T>::value);
24 static const char* function = "boost::math::binomial_coefficient<%1%>(unsigned, unsigned)";
26 return policies::raise_domain_error<T>(
28 "The binomial coefficient is undefined for k > n, but got k = %1%.",
31 if((k == 0) || (k == n))
33 if((k == 1) || (k == n-1))
36 if(n <= max_factorial<T>::value)
38 // Use fast table lookup:
39 result = unchecked_factorial<T>(n);
40 result /= unchecked_factorial<T>(n-k);
41 result /= unchecked_factorial<T>(k);
45 // Use the beta function:
47 result = k * beta(static_cast<T>(k), static_cast<T>(n-k+1), pol);
49 result = (n - k) * beta(static_cast<T>(k+1), static_cast<T>(n-k), pol);
51 return policies::raise_overflow_error<T>(function, 0, pol);
54 // convert to nearest integer:
55 return ceil(result - 0.5f);
58 // Type float can only store the first 35 factorials, in order to
59 // increase the chance that we can use a table driven implementation
60 // we'll promote to double:
63 inline float binomial_coefficient<float, policies::policy<> >(unsigned n, unsigned k, const policies::policy<>& pol)
65 return policies::checked_narrowing_cast<float, policies::policy<> >(binomial_coefficient<double>(n, k, pol), "boost::math::binomial_coefficient<%1%>(unsigned,unsigned)");
69 inline T binomial_coefficient(unsigned n, unsigned k)
71 return binomial_coefficient<T>(n, k, policies::policy<>());
78 #endif // BOOST_MATH_SF_BINOMIAL_HPP