binomial.hpp 2.4 KB

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