]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | ||
18 | namespace boost{ namespace math{ | |
19 | ||
20 | template <class T, class Policy> | |
21 | T 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 | // | |
63 | template <> | |
20effc67 | 64 | inline float binomial_coefficient<float, policies::policy<> >(unsigned n, unsigned k, const policies::policy<>&) |
7c673cae | 65 | { |
20effc67 TL |
66 | typedef policies::normalise< |
67 | policies::policy<>, | |
68 | policies::promote_float<true>, | |
69 | policies::promote_double<false>, | |
70 | policies::discrete_quantile<>, | |
71 | policies::assert_undefined<> >::type forwarding_policy; | |
72 | return policies::checked_narrowing_cast<float, forwarding_policy>(binomial_coefficient<double>(n, k, forwarding_policy()), "boost::math::binomial_coefficient<%1%>(unsigned,unsigned)"); | |
7c673cae FG |
73 | } |
74 | ||
75 | template <class T> | |
76 | inline T binomial_coefficient(unsigned n, unsigned k) | |
77 | { | |
78 | return binomial_coefficient<T>(n, k, policies::policy<>()); | |
79 | } | |
80 | ||
81 | } // namespace math | |
82 | } // namespace boost | |
83 | ||
84 | ||
85 | #endif // BOOST_MATH_SF_BINOMIAL_HPP | |
86 | ||
87 | ||
88 |