binomial.hpp 2.4 KB

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