binomial.hpp 2.9 KB

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