Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:40:10

0001 //  Copyright John Maddock 2006.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005 
0006 #ifndef BOOST_MATH_SF_BINOMIAL_HPP
0007 #define BOOST_MATH_SF_BINOMIAL_HPP
0008 
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #endif
0012 
0013 #include <boost/math/special_functions/math_fwd.hpp>
0014 #include <boost/math/special_functions/factorials.hpp>
0015 #include <boost/math/special_functions/beta.hpp>
0016 #include <boost/math/policies/error_handling.hpp>
0017 #include <type_traits>
0018 
0019 namespace boost{ namespace math{
0020 
0021 template <class T, class Policy>
0022 T binomial_coefficient(unsigned n, unsigned k, const Policy& pol)
0023 {
0024    static_assert(!std::is_integral<T>::value, "Type T must not be an integral type");
0025    BOOST_MATH_STD_USING
0026    static const char* function = "boost::math::binomial_coefficient<%1%>(unsigned, unsigned)";
0027    if(k > n)
0028       return policies::raise_domain_error<T>(
0029          function,
0030          "The binomial coefficient is undefined for k > n, but got k = %1%.",
0031          static_cast<T>(k), pol);
0032    T result;
0033    if((k == 0) || (k == n))
0034       return static_cast<T>(1);
0035    if((k == 1) || (k == n-1))
0036       return static_cast<T>(n);
0037 
0038    if(n <= max_factorial<T>::value)
0039    {
0040       // Use fast table lookup:
0041       result = unchecked_factorial<T>(n);
0042       result /= unchecked_factorial<T>(n-k);
0043       result /= unchecked_factorial<T>(k);
0044    }
0045    else
0046    {
0047       // Use the beta function:
0048       if(k < n - k)
0049          result = static_cast<T>(k * beta(static_cast<T>(k), static_cast<T>(n-k+1), pol));
0050       else
0051          result = static_cast<T>((n - k) * beta(static_cast<T>(k+1), static_cast<T>(n-k), pol));
0052       if(result == 0)
0053          return policies::raise_overflow_error<T>(function, nullptr, pol);
0054       result = 1 / result;
0055    }
0056    // convert to nearest integer:
0057    return ceil(result - 0.5f);
0058 }
0059 //
0060 // Type float can only store the first 35 factorials, in order to
0061 // increase the chance that we can use a table driven implementation
0062 // we'll promote to double:
0063 //
0064 template <>
0065 inline float binomial_coefficient<float, policies::policy<> >(unsigned n, unsigned k, const policies::policy<>&)
0066 {
0067    typedef policies::normalise<
0068        policies::policy<>,
0069        policies::promote_float<true>,
0070        policies::promote_double<false>,
0071        policies::discrete_quantile<>,
0072        policies::assert_undefined<> >::type forwarding_policy;
0073    return policies::checked_narrowing_cast<float, forwarding_policy>(binomial_coefficient<double>(n, k, forwarding_policy()), "boost::math::binomial_coefficient<%1%>(unsigned,unsigned)");
0074 }
0075 
0076 template <class T>
0077 inline T binomial_coefficient(unsigned n, unsigned k)
0078 {
0079    return binomial_coefficient<T>(n, k, policies::policy<>());
0080 }
0081 
0082 } // namespace math
0083 } // namespace boost
0084 
0085 
0086 #endif // BOOST_MATH_SF_BINOMIAL_HPP
0087 
0088 
0089