88 lines · plain
1// Copyright John Maddock 2006.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_SF_BINOMIAL_HPP7#define BOOST_MATH_SF_BINOMIAL_HPP8 9#ifdef _MSC_VER10#pragma once11#endif12 13#include <boost/math/tools/config.hpp>14#include <boost/math/tools/type_traits.hpp>15#include <boost/math/special_functions/math_fwd.hpp>16#include <boost/math/special_functions/factorials.hpp>17#include <boost/math/special_functions/beta.hpp>18#include <boost/math/policies/error_handling.hpp>19 20namespace boost{ namespace math{21 22template <class T, class Policy>23BOOST_MATH_GPU_ENABLED T binomial_coefficient(unsigned n, unsigned k, const Policy& pol)24{25 static_assert(!boost::math::is_integral<T>::value, "Type T must not be an integral type");26 BOOST_MATH_STD_USING27 constexpr auto function = "boost::math::binomial_coefficient<%1%>(unsigned, unsigned)";28 if(k > n)29 return policies::raise_domain_error<T>(function, "The binomial coefficient is undefined for k > n, but got k = %1%.", static_cast<T>(k), pol);30 T result; // LCOV_EXCL_LINE31 if((k == 0) || (k == n))32 return static_cast<T>(1);33 if((k == 1) || (k == n-1))34 return static_cast<T>(n);35 36 if(n <= max_factorial<T>::value)37 {38 // Use fast table lookup:39 result = unchecked_factorial<T>(n);40 result /= unchecked_factorial<T>(n-k);41 result /= unchecked_factorial<T>(k);42 }43 else44 {45 // Use the beta function:46 if(k < n - k)47 result = static_cast<T>(k * boost::math::beta(static_cast<T>(k), static_cast<T>(n-k+1), pol));48 else49 result = static_cast<T>((n - k) * boost::math::beta(static_cast<T>(k+1), static_cast<T>(n-k), pol));50 if(result == 0)51 return policies::raise_overflow_error<T>(function, nullptr, pol);52 result = 1 / result;53 }54 // convert to nearest integer:55 return ceil(result - 0.5f);56}57//58// Type float can only store the first 35 factorials, in order to59// increase the chance that we can use a table driven implementation60// we'll promote to double:61//62template <>63BOOST_MATH_GPU_ENABLED inline float binomial_coefficient<float, policies::policy<> >(unsigned n, unsigned k, const policies::policy<>&)64{65 typedef policies::normalise<66 policies::policy<>,67 policies::promote_float<true>,68 policies::promote_double<false>,69 policies::discrete_quantile<>,70 policies::assert_undefined<> >::type forwarding_policy;71 return policies::checked_narrowing_cast<float, forwarding_policy>(binomial_coefficient<double>(n, k, forwarding_policy()), "boost::math::binomial_coefficient<%1%>(unsigned,unsigned)");72}73 74template <class T>75BOOST_MATH_GPU_ENABLED inline T binomial_coefficient(unsigned n, unsigned k)76{77 return binomial_coefficient<T>(n, k, policies::policy<>());78}79 80} // namespace math81} // namespace boost82 83 84#endif // BOOST_MATH_SF_BINOMIAL_HPP85 86 87 88