binomial.hpp 2.7 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788
  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<>&)
  59. {
  60. typedef policies::normalise<
  61. policies::policy<>,
  62. policies::promote_float<true>,
  63. policies::promote_double<false>,
  64. policies::discrete_quantile<>,
  65. policies::assert_undefined<> >::type forwarding_policy;
  66. return policies::checked_narrowing_cast<float, forwarding_policy>(binomial_coefficient<double>(n, k, forwarding_policy()), "boost::math::binomial_coefficient<%1%>(unsigned,unsigned)");
  67. }
  68. template <class T>
  69. inline T binomial_coefficient(unsigned n, unsigned k)
  70. {
  71. return binomial_coefficient<T>(n, k, policies::policy<>());
  72. }
  73. } // namespace math
  74. } // namespace boost
  75. #endif // BOOST_MATH_SF_BINOMIAL_HPP