ceph/src/boost/libs/math/include/boost/math/special_functions/binomial.hpp

   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
   6 #ifndef BOOST_MATH_SF_BINOMIAL_HPP
   7 #define BOOST_MATH_SF_BINOMIAL_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 #include <boost/math/special_functions/math_fwd.hpp>
  14 #include <boost/math/special_functions/factorials.hpp>
  15 #include <boost/math/special_functions/beta.hpp>
  16 #include <boost/math/policies/error_handling.hpp>
  17
  18 namespace boost{ namespace math{
  19
  20 template <class T, class Policy>
  21 T binomial_coefficient(unsigned n, unsigned k, const Policy& pol)
  22 {
  23    BOOST_STATIC_ASSERT(!boost::is_integral<T>::value);
  24    BOOST_MATH_STD_USING
  25    static const char* function = "boost::math::binomial_coefficient<%1%>(unsigned, unsigned)";
  26    if(k > n)
  27       return policies::raise_domain_error<T>(
  28          function,
  29          "The binomial coefficient is undefined for k > n, but got k = %1%.",
  30          static_cast<T>(k), pol);
  31    T result;
  32    if((k == 0) || (k == n))
  33       return static_cast<T>(1);
  34    if((k == 1) || (k == n-1))
  35       return static_cast<T>(n);
  36
  37    if(n <= max_factorial<T>::value)
  38    {
  39       // Use fast table lookup:
  40       result = unchecked_factorial<T>(n);
  41       result /= unchecked_factorial<T>(n-k);
  42       result /= unchecked_factorial<T>(k);
  43    }
  44    else
  45    {
  46       // Use the beta function:
  47       if(k < n - k)
  48          result = k * beta(static_cast<T>(k), static_cast<T>(n-k+1), pol);
  49       else
  50          result = (n - k) * beta(static_cast<T>(k+1), static_cast<T>(n-k), pol);
  51       if(result == 0)
  52          return policies::raise_overflow_error<T>(function, 0, pol);
  53       result = 1 / result;
  54    }
  55    // convert to nearest integer:
  56    return ceil(result - 0.5f);
  57 }
  58 //
  59 // Type float can only store the first 35 factorials, in order to
  60 // increase the chance that we can use a table driven implementation
  61 // we'll promote to double:
  62 //
  63 template <>
  64 inline float binomial_coefficient<float, policies::policy<> >(unsigned n, unsigned k, const policies::policy<>& pol)
  65 {
  66    return policies::checked_narrowing_cast<float, policies::policy<> >(binomial_coefficient<double>(n, k, pol), "boost::math::binomial_coefficient<%1%>(unsigned,unsigned)");
  67 }
  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
  75 } // namespace math
  76 } // namespace boost
  77
  78
  79 #endif // BOOST_MATH_SF_BINOMIAL_HPP
  80
  81
  82