| // 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 |
| |
| |
| |