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