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



#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_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





