binomial.hpp 2.9 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/tools/config.hpp>
  11. #include <boost/math/tools/type_traits.hpp>
  12. #include <boost/math/special_functions/math_fwd.hpp>
  13. #include <boost/math/special_functions/factorials.hpp>
  14. #include <boost/math/special_functions/beta.hpp>
  15. #include <boost/math/policies/error_handling.hpp>
  16. namespace boost{ namespace math{
  17. template <class T, class Policy>
  18. BOOST_MATH_GPU_ENABLED T binomial_coefficient(unsigned n, unsigned k, const Policy& pol)
  19. {
  20. static_assert(!boost::math::is_integral<T>::value, "Type T must not be an integral type");
  21. BOOST_MATH_STD_USING
  22. constexpr auto function = "boost::math::binomial_coefficient<%1%>(unsigned, unsigned)";
  23. if(k > n)
  24. return policies::raise_domain_error<T>(function, "The binomial coefficient is undefined for k > n, but got k = %1%.", static_cast<T>(k), pol);
  25. T result; // LCOV_EXCL_LINE
  26. if((k == 0) || (k == n))
  27. return static_cast<T>(1);
  28. if((k == 1) || (k == n-1))
  29. return static_cast<T>(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 = static_cast<T>(k * boost::math::beta(static_cast<T>(k), static_cast<T>(n-k+1), pol));
  42. else
  43. result = static_cast<T>((n - k) * boost::math::beta(static_cast<T>(k+1), static_cast<T>(n-k), pol));
  44. if(result == 0)
  45. return policies::raise_overflow_error<T>(function, nullptr, 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. BOOST_MATH_GPU_ENABLED inline float binomial_coefficient<float, policies::policy<> >(unsigned n, unsigned k, const policies::policy<>&)
  58. {
  59. typedef policies::normalise<
  60. policies::policy<>,
  61. policies::promote_float<true>,
  62. policies::promote_double<false>,
  63. policies::discrete_quantile<>,
  64. policies::assert_undefined<> >::type forwarding_policy;
  65. return policies::checked_narrowing_cast<float, forwarding_policy>(binomial_coefficient<double>(n, k, forwarding_policy()), "boost::math::binomial_coefficient<%1%>(unsigned,unsigned)");
  66. }
  67. template <class T>
  68. BOOST_MATH_GPU_ENABLED inline T binomial_coefficient(unsigned n, unsigned k)
  69. {
  70. return binomial_coefficient<T>(n, k, policies::policy<>());
  71. }
  72. } // namespace math
  73. } // namespace boost
  74. #endif // BOOST_MATH_SF_BINOMIAL_HPP