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)
5
6#ifndef BOOST_MATH_SF_BINOMIAL_HPP
7#define BOOST_MATH_SF_BINOMIAL_HPP
8
9#ifdef _MSC_VER
10#pragma once
11#endif
12
13#include <boost/math/special_functions/math_fwd.hpp>
14#include <boost/math/special_functions/factorials.hpp>
15#include <boost/math/special_functions/beta.hpp>
16#include <boost/math/policies/error_handling.hpp>
17
18namespace boost{ namespace math{
19
20template <class T, class Policy>
21T binomial_coefficient(unsigned n, unsigned k, const Policy& pol)
22{
23 BOOST_STATIC_ASSERT(!boost::is_integral<T>::value);
24 BOOST_MATH_STD_USING
25 static const char* function = "boost::math::binomial_coefficient<%1%>(unsigned, unsigned)";
26 if(k > n)
27 return policies::raise_domain_error<T>(
28 function,
29 "The binomial coefficient is undefined for k > n, but got k = %1%.",
30 static_cast<T>(k), pol);
31 T result;
32 if((k == 0) || (k == n))
33 return static_cast<T>(1);
34 if((k == 1) || (k == n-1))
35 return static_cast<T>(n);
36
37 if(n <= max_factorial<T>::value)
38 {
39 // Use fast table lookup:
40 result = unchecked_factorial<T>(n);
41 result /= unchecked_factorial<T>(n-k);
42 result /= unchecked_factorial<T>(k);
43 }
44 else
45 {
46 // Use the beta function:
47 if(k < n - k)
48 result = k * beta(static_cast<T>(k), static_cast<T>(n-k+1), pol);
49 else
50 result = (n - k) * beta(static_cast<T>(k+1), static_cast<T>(n-k), pol);
51 if(result == 0)
52 return policies::raise_overflow_error<T>(function, 0, pol);
53 result = 1 / result;
54 }
55 // convert to nearest integer:
56 return ceil(result - 0.5f);
57}
58//
59// Type float can only store the first 35 factorials, in order to
60// increase the chance that we can use a table driven implementation
61// we'll promote to double:
62//
63template <>
64inline float binomial_coefficient<float, policies::policy<> >(unsigned n, unsigned k, const policies::policy<>& pol)
65{
66 return policies::checked_narrowing_cast<float, policies::policy<> >(val: binomial_coefficient<double>(n, k, pol), function: "boost::math::binomial_coefficient<%1%>(unsigned,unsigned)");
67}
68
69template <class T>
70inline T binomial_coefficient(unsigned n, unsigned k)
71{
72 return binomial_coefficient<T>(n, k, policies::policy<>());
73}
74
75} // namespace math
76} // namespace boost
77
78
79#endif // BOOST_MATH_SF_BINOMIAL_HPP
80
81
82
83

source code of include/boost/math/special_functions/binomial.hpp