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#include <type_traits>
18
19namespace boost{ namespace math{
20
21template <class T, class Policy>
22T binomial_coefficient(unsigned n, unsigned k, const Policy& pol)
23{
24 static_assert(!std::is_integral<T>::value, "Type T must not be an integral type");
25 BOOST_MATH_STD_USING
26 static const char* function = "boost::math::binomial_coefficient<%1%>(unsigned, unsigned)";
27 if(k > n)
28 return policies::raise_domain_error<T>(function, "The binomial coefficient is undefined for k > n, but got k = %1%.", static_cast<T>(k), pol);
29 T result; // LCOV_EXCL_LINE
30 if((k == 0) || (k == n))
31 return static_cast<T>(1);
32 if((k == 1) || (k == n-1))
33 return static_cast<T>(n);
34
35 if(n <= max_factorial<T>::value)
36 {
37 // Use fast table lookup:
38 result = unchecked_factorial<T>(n);
39 result /= unchecked_factorial<T>(n-k);
40 result /= unchecked_factorial<T>(k);
41 }
42 else
43 {
44 // Use the beta function:
45 if(k < n - k)
46 result = static_cast<T>(k * beta(static_cast<T>(k), static_cast<T>(n-k+1), pol));
47 else
48 result = static_cast<T>((n - k) * beta(static_cast<T>(k+1), static_cast<T>(n-k), pol));
49 if(result == 0)
50 return policies::raise_overflow_error<T>(function, nullptr, pol);
51 result = 1 / result;
52 }
53 // convert to nearest integer:
54 return ceil(result - 0.5f);
55}
56//
57// Type float can only store the first 35 factorials, in order to
58// increase the chance that we can use a table driven implementation
59// we'll promote to double:
60//
61template <>
62inline float binomial_coefficient<float, policies::policy<> >(unsigned n, unsigned k, const policies::policy<>&)
63{
64 typedef policies::normalise<
65 policies::policy<>,
66 policies::promote_float<true>,
67 policies::promote_double<false>,
68 policies::discrete_quantile<>,
69 policies::assert_undefined<> >::type forwarding_policy;
70 return policies::checked_narrowing_cast<float, forwarding_policy>(val: binomial_coefficient<double>(n, k, pol: forwarding_policy()), function: "boost::math::binomial_coefficient<%1%>(unsigned,unsigned)");
71}
72
73template <class T>
74inline T binomial_coefficient(unsigned n, unsigned k)
75{
76 return binomial_coefficient<T>(n, k, policies::policy<>());
77}
78
79} // namespace math
80} // namespace boost
81
82
83#endif // BOOST_MATH_SF_BINOMIAL_HPP
84
85
86
87

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