| // 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_SP_FACTORIALS_HPP |
| #define BOOST_MATH_SP_FACTORIALS_HPP |
| |
| #ifdef _MSC_VER |
| #pragma once |
| #endif |
| |
| #include <boost/math/special_functions/gamma.hpp> |
| #include <boost/math/special_functions/math_fwd.hpp> |
| #include <boost/math/special_functions/detail/unchecked_factorial.hpp> |
| #include <boost/array.hpp> |
| #ifdef BOOST_MSVC |
| #pragma warning(push) // Temporary until lexical cast fixed. |
| #pragma warning(disable: 4127 4701) |
| #endif |
| #include <boost/lexical_cast.hpp> |
| #ifdef BOOST_MSVC |
| #pragma warning(pop) |
| #endif |
| #include <boost/config/no_tr1/cmath.hpp> |
| |
| namespace boost { namespace math |
| { |
| |
| template <class T, class Policy> |
| inline T factorial(unsigned i, const Policy& pol) |
| { |
| BOOST_STATIC_ASSERT(!boost::is_integral<T>::value); |
| BOOST_MATH_STD_USING // Aid ADL for floor. |
| |
| if(i <= max_factorial<T>::value) |
| return unchecked_factorial<T>(i); |
| T result = boost::math::tgamma(static_cast<T>(i+1), pol); |
| if(result > tools::max_value<T>()) |
| return result; // Overflowed value! (But tgamma will have signalled the error already). |
| return floor(result + 0.5f); |
| } |
| |
| template <class T> |
| inline T factorial(unsigned i) |
| { |
| return factorial<T>(i, policies::policy<>()); |
| } |
| /* |
| // Can't have these in a policy enabled world? |
| template<> |
| inline float factorial<float>(unsigned i) |
| { |
| if(i <= max_factorial<float>::value) |
| return unchecked_factorial<float>(i); |
| return tools::overflow_error<float>(BOOST_CURRENT_FUNCTION); |
| } |
| |
| template<> |
| inline double factorial<double>(unsigned i) |
| { |
| if(i <= max_factorial<double>::value) |
| return unchecked_factorial<double>(i); |
| return tools::overflow_error<double>(BOOST_CURRENT_FUNCTION); |
| } |
| */ |
| template <class T, class Policy> |
| T double_factorial(unsigned i, const Policy& pol) |
| { |
| BOOST_STATIC_ASSERT(!boost::is_integral<T>::value); |
| BOOST_MATH_STD_USING // ADL lookup of std names |
| if(i & 1) |
| { |
| // odd i: |
| if(i < max_factorial<T>::value) |
| { |
| unsigned n = (i - 1) / 2; |
| return ceil(unchecked_factorial<T>(i) / (ldexp(T(1), (int)n) * unchecked_factorial<T>(n)) - 0.5f); |
| } |
| // |
| // Fallthrough: i is too large to use table lookup, try the |
| // gamma function instead. |
| // |
| T result = boost::math::tgamma(static_cast<T>(i) / 2 + 1, pol) / sqrt(constants::pi<T>()); |
| if(ldexp(tools::max_value<T>(), -static_cast<int>(i+1) / 2) > result) |
| return ceil(result * ldexp(T(1), (i+1) / 2) - 0.5f); |
| } |
| else |
| { |
| // even i: |
| unsigned n = i / 2; |
| T result = factorial<T>(n, pol); |
| if(ldexp(tools::max_value<T>(), -(int)n) > result) |
| return result * ldexp(T(1), (int)n); |
| } |
| // |
| // If we fall through to here then the result is infinite: |
| // |
| return policies::raise_overflow_error<T>("boost::math::double_factorial<%1%>(unsigned)", 0, pol); |
| } |
| |
| template <class T> |
| inline T double_factorial(unsigned i) |
| { |
| return double_factorial<T>(i, policies::policy<>()); |
| } |
| |
| namespace detail{ |
| |
| template <class T, class Policy> |
| T rising_factorial_imp(T x, int n, const Policy& pol) |
| { |
| BOOST_STATIC_ASSERT(!boost::is_integral<T>::value); |
| if(x < 0) |
| { |
| // |
| // For x less than zero, we really have a falling |
| // factorial, modulo a possible change of sign. |
| // |
| // Note that the falling factorial isn't defined |
| // for negative n, so we'll get rid of that case |
| // first: |
| // |
| bool inv = false; |
| if(n < 0) |
| { |
| x += n; |
| n = -n; |
| inv = true; |
| } |
| T result = ((n&1) ? -1 : 1) * falling_factorial(-x, n, pol); |
| if(inv) |
| result = 1 / result; |
| return result; |
| } |
| if(n == 0) |
| return 1; |
| // |
| // We don't optimise this for small n, because |
| // tgamma_delta_ratio is alreay optimised for that |
| // use case: |
| // |
| return 1 / boost::math::tgamma_delta_ratio(x, static_cast<T>(n), pol); |
| } |
| |
| template <class T, class Policy> |
| inline T falling_factorial_imp(T x, unsigned n, const Policy& pol) |
| { |
| BOOST_STATIC_ASSERT(!boost::is_integral<T>::value); |
| BOOST_MATH_STD_USING // ADL of std names |
| if(x == 0) |
| return 0; |
| if(x < 0) |
| { |
| // |
| // For x < 0 we really have a rising factorial |
| // modulo a possible change of sign: |
| // |
| return (n&1 ? -1 : 1) * rising_factorial(-x, n, pol); |
| } |
| if(n == 0) |
| return 1; |
| if(x < n-1) |
| { |
| // |
| // x+1-n will be negative and tgamma_delta_ratio won't |
| // handle it, split the product up into three parts: |
| // |
| T xp1 = x + 1; |
| unsigned n2 = itrunc((T)floor(xp1), pol); |
| if(n2 == xp1) |
| return 0; |
| T result = boost::math::tgamma_delta_ratio(xp1, -static_cast<T>(n2), pol); |
| x -= n2; |
| result *= x; |
| ++n2; |
| if(n2 < n) |
| result *= falling_factorial(x - 1, n - n2, pol); |
| return result; |
| } |
| // |
| // Simple case: just the ratio of two |
| // (positive argument) gamma functions. |
| // Note that we don't optimise this for small n, |
| // because tgamma_delta_ratio is alreay optimised |
| // for that use case: |
| // |
| return boost::math::tgamma_delta_ratio(x + 1, -static_cast<T>(n), pol); |
| } |
| |
| } // namespace detail |
| |
| template <class RT> |
| inline typename tools::promote_args<RT>::type |
| falling_factorial(RT x, unsigned n) |
| { |
| typedef typename tools::promote_args<RT>::type result_type; |
| return detail::falling_factorial_imp( |
| static_cast<result_type>(x), n, policies::policy<>()); |
| } |
| |
| template <class RT, class Policy> |
| inline typename tools::promote_args<RT>::type |
| falling_factorial(RT x, unsigned n, const Policy& pol) |
| { |
| typedef typename tools::promote_args<RT>::type result_type; |
| return detail::falling_factorial_imp( |
| static_cast<result_type>(x), n, pol); |
| } |
| |
| template <class RT> |
| inline typename tools::promote_args<RT>::type |
| rising_factorial(RT x, int n) |
| { |
| typedef typename tools::promote_args<RT>::type result_type; |
| return detail::rising_factorial_imp( |
| static_cast<result_type>(x), n, policies::policy<>()); |
| } |
| |
| template <class RT, class Policy> |
| inline typename tools::promote_args<RT>::type |
| rising_factorial(RT x, int n, const Policy& pol) |
| { |
| typedef typename tools::promote_args<RT>::type result_type; |
| return detail::rising_factorial_imp( |
| static_cast<result_type>(x), n, pol); |
| } |
| |
| } // namespace math |
| } // namespace boost |
| |
| #endif // BOOST_MATH_SP_FACTORIALS_HPP |
| |