Updated boost to v1.55.0

[rsem.git] / boost / math / special_functions / binomial.hpp

diff --git a/boost/math/special_functions/binomial.hpp b/boost/math/special_functions/binomial.hpp
new file mode 100644 (file)
index 0000000..16b4f33
--- /dev/null
+++ b/boost/math/special_functions/binomial.hpp
@@ -0,0 +1,81 @@
+//  Copyright John Maddock 2006.
+//  Use, modification and distribution are subject to the
+//  Boost Software License, Version 1.0. (See accompanying file
+//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_SF_BINOMIAL_HPP
+#define BOOST_MATH_SF_BINOMIAL_HPP
+
+#ifdef _MSC_VER
+#pragma once
+#endif
+
+#include <boost/math/special_functions/factorials.hpp>
+#include <boost/math/special_functions/beta.hpp>
+#include <boost/math/policies/error_handling.hpp>
+
+namespace boost{ namespace math{
+
+template <class T, class Policy>
+T binomial_coefficient(unsigned n, unsigned k, const Policy& pol)
+{
+   BOOST_STATIC_ASSERT(!boost::is_integral<T>::value);
+   BOOST_MATH_STD_USING
+   static const char* function = "boost::math::binomial_coefficient<%1%>(unsigned, unsigned)";
+   if(k > n)
+      return policies::raise_domain_error<T>(
+         function, 
+         "The binomial coefficient is undefined for k > n, but got k = %1%.",
+         k, pol);
+   T result;
+   if((k == 0) || (k == n))
+      return 1;
+   if((k == 1) || (k == n-1))
+      return n;
+
+   if(n <= max_factorial<T>::value)
+   {
+      // Use fast table lookup:
+      result = unchecked_factorial<T>(n);
+      result /= unchecked_factorial<T>(n-k);
+      result /= unchecked_factorial<T>(k);
+   }
+   else
+   {
+      // Use the beta function:
+      if(k < n - k)
+         result = k * beta(static_cast<T>(k), static_cast<T>(n-k+1), pol);
+      else
+         result = (n - k) * beta(static_cast<T>(k+1), static_cast<T>(n-k), pol);
+      if(result == 0)
+         return policies::raise_overflow_error<T>(function, 0, pol);
+      result = 1 / result;
+   }
+   // convert to nearest integer:
+   return ceil(result - 0.5f);
+}
+//
+// Type float can only store the first 35 factorials, in order to
+// increase the chance that we can use a table driven implementation
+// we'll promote to double:
+//
+template <>
+inline float binomial_coefficient<float, policies::policy<> >(unsigned n, unsigned k, const policies::policy<>& pol)
+{
+   return policies::checked_narrowing_cast<float, policies::policy<> >(binomial_coefficient<double>(n, k, pol), "boost::math::binomial_coefficient<%1%>(unsigned,unsigned)");
+}
+
+template <class T>
+inline T binomial_coefficient(unsigned n, unsigned k)
+{
+   return binomial_coefficient<T>(n, k, policies::policy<>());
+}
+
+} // namespace math
+} // namespace boost
+
+
+#endif // BOOST_MATH_SF_BINOMIAL_HPP
+
+
+