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