Chris@16: //  Copyright John Maddock 2006, 2010.
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_MATH_SP_FACTORIALS_HPP
Chris@16: #define BOOST_MATH_SP_FACTORIALS_HPP
Chris@16: 
Chris@16: #ifdef _MSC_VER
Chris@16: #pragma once
Chris@16: #endif
Chris@16: 
Chris@101: #include <boost/math/special_functions/math_fwd.hpp>
Chris@16: #include <boost/math/special_functions/gamma.hpp>
Chris@16: #include <boost/math/special_functions/detail/unchecked_factorial.hpp>
Chris@16: #include <boost/array.hpp>
Chris@16: #ifdef BOOST_MSVC
Chris@16: #pragma warning(push) // Temporary until lexical cast fixed.
Chris@16: #pragma warning(disable: 4127 4701)
Chris@16: #endif
Chris@16: #ifdef BOOST_MSVC
Chris@16: #pragma warning(pop)
Chris@16: #endif
Chris@16: #include <boost/config/no_tr1/cmath.hpp>
Chris@16: 
Chris@16: namespace boost { namespace math
Chris@16: {
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline T factorial(unsigned i, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_STATIC_ASSERT(!boost::is_integral<T>::value);
Chris@16:    // factorial<unsigned int>(n) is not implemented
Chris@16:    // because it would overflow integral type T for too small n
Chris@16:    // to be useful. Use instead a floating-point type,
Chris@16:    // and convert to an unsigned type if essential, for example:
Chris@16:    // unsigned int nfac = static_cast<unsigned int>(factorial<double>(n));
Chris@16:    // See factorial documentation for more detail.
Chris@16: 
Chris@16:    BOOST_MATH_STD_USING // Aid ADL for floor.
Chris@16: 
Chris@16:    if(i <= max_factorial<T>::value)
Chris@16:       return unchecked_factorial<T>(i);
Chris@16:    T result = boost::math::tgamma(static_cast<T>(i+1), pol);
Chris@16:    if(result > tools::max_value<T>())
Chris@16:       return result; // Overflowed value! (But tgamma will have signalled the error already).
Chris@16:    return floor(result + 0.5f);
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T factorial(unsigned i)
Chris@16: {
Chris@16:    return factorial<T>(i, policies::policy<>());
Chris@16: }
Chris@16: /*
Chris@16: // Can't have these in a policy enabled world?
Chris@16: template<>
Chris@16: inline float factorial<float>(unsigned i)
Chris@16: {
Chris@16:    if(i <= max_factorial<float>::value)
Chris@16:       return unchecked_factorial<float>(i);
Chris@16:    return tools::overflow_error<float>(BOOST_CURRENT_FUNCTION);
Chris@16: }
Chris@16: 
Chris@16: template<>
Chris@16: inline double factorial<double>(unsigned i)
Chris@16: {
Chris@16:    if(i <= max_factorial<double>::value)
Chris@16:       return unchecked_factorial<double>(i);
Chris@16:    return tools::overflow_error<double>(BOOST_CURRENT_FUNCTION);
Chris@16: }
Chris@16: */
Chris@16: template <class T, class Policy>
Chris@16: T double_factorial(unsigned i, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_STATIC_ASSERT(!boost::is_integral<T>::value);
Chris@16:    BOOST_MATH_STD_USING  // ADL lookup of std names
Chris@16:    if(i & 1)
Chris@16:    {
Chris@16:       // odd i:
Chris@16:       if(i < max_factorial<T>::value)
Chris@16:       {
Chris@16:          unsigned n = (i - 1) / 2;
Chris@16:          return ceil(unchecked_factorial<T>(i) / (ldexp(T(1), (int)n) * unchecked_factorial<T>(n)) - 0.5f);
Chris@16:       }
Chris@16:       //
Chris@16:       // Fallthrough: i is too large to use table lookup, try the
Chris@16:       // gamma function instead.
Chris@16:       //
Chris@16:       T result = boost::math::tgamma(static_cast<T>(i) / 2 + 1, pol) / sqrt(constants::pi<T>());
Chris@16:       if(ldexp(tools::max_value<T>(), -static_cast<int>(i+1) / 2) > result)
Chris@16:          return ceil(result * ldexp(T(1), static_cast<int>(i+1) / 2) - 0.5f);
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       // even i:
Chris@16:       unsigned n = i / 2;
Chris@16:       T result = factorial<T>(n, pol);
Chris@16:       if(ldexp(tools::max_value<T>(), -(int)n) > result)
Chris@16:          return result * ldexp(T(1), (int)n);
Chris@16:    }
Chris@16:    //
Chris@16:    // If we fall through to here then the result is infinite:
Chris@16:    //
Chris@16:    return policies::raise_overflow_error<T>("boost::math::double_factorial<%1%>(unsigned)", 0, pol);
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T double_factorial(unsigned i)
Chris@16: {
Chris@16:    return double_factorial<T>(i, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: namespace detail{
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: T rising_factorial_imp(T x, int n, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_STATIC_ASSERT(!boost::is_integral<T>::value);
Chris@16:    if(x < 0)
Chris@16:    {
Chris@16:       //
Chris@16:       // For x less than zero, we really have a falling
Chris@16:       // factorial, modulo a possible change of sign.
Chris@16:       //
Chris@16:       // Note that the falling factorial isn't defined
Chris@16:       // for negative n, so we'll get rid of that case
Chris@16:       // first:
Chris@16:       //
Chris@16:       bool inv = false;
Chris@16:       if(n < 0)
Chris@16:       {
Chris@16:          x += n;
Chris@16:          n = -n;
Chris@16:          inv = true;
Chris@16:       }
Chris@16:       T result = ((n&1) ? -1 : 1) * falling_factorial(-x, n, pol);
Chris@16:       if(inv)
Chris@16:          result = 1 / result;
Chris@16:       return result;
Chris@16:    }
Chris@16:    if(n == 0)
Chris@16:       return 1;
Chris@101:    if(x == 0)
Chris@101:    {
Chris@101:       if(n < 0)
Chris@101:          return -boost::math::tgamma_delta_ratio(x + 1, static_cast<T>(-n), pol);
Chris@101:       else
Chris@101:          return 0;
Chris@101:    }
Chris@101:    if((x < 1) && (x + n < 0))
Chris@101:    {
Chris@101:       T val = boost::math::tgamma_delta_ratio(1 - x, static_cast<T>(-n), pol);
Chris@101:       return (n & 1) ? T(-val) : val;
Chris@101:    }
Chris@16:    //
Chris@16:    // We don't optimise this for small n, because
Chris@16:    // tgamma_delta_ratio is alreay optimised for that
Chris@16:    // use case:
Chris@16:    //
Chris@16:    return 1 / boost::math::tgamma_delta_ratio(x, static_cast<T>(n), pol);
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline T falling_factorial_imp(T x, unsigned n, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_STATIC_ASSERT(!boost::is_integral<T>::value);
Chris@16:    BOOST_MATH_STD_USING // ADL of std names
Chris@101:    if((x == 0) && (n >= 0))
Chris@16:       return 0;
Chris@16:    if(x < 0)
Chris@16:    {
Chris@16:       //
Chris@16:       // For x < 0 we really have a rising factorial
Chris@16:       // modulo a possible change of sign:
Chris@16:       //
Chris@16:       return (n&1 ? -1 : 1) * rising_factorial(-x, n, pol);
Chris@16:    }
Chris@16:    if(n == 0)
Chris@16:       return 1;
Chris@101:    if(x < 0.5f)
Chris@101:    {
Chris@101:       //
Chris@101:       // 1 + x below will throw away digits, so split up calculation:
Chris@101:       //
Chris@101:       if(n > max_factorial<T>::value - 2)
Chris@101:       {
Chris@101:          // If the two end of the range are far apart we have a ratio of two very large
Chris@101:          // numbers, split the calculation up into two blocks:
Chris@101:          T t1 = x * boost::math::falling_factorial(x - 1, max_factorial<T>::value - 2);
Chris@101:          T t2 = boost::math::falling_factorial(x - max_factorial<T>::value + 1, n - max_factorial<T>::value + 1);
Chris@101:          if(tools::max_value<T>() / fabs(t1) < fabs(t2))
Chris@101:             return boost::math::sign(t1) * boost::math::sign(t2) * policies::raise_overflow_error<T>("boost::math::falling_factorial<%1%>", 0, pol);
Chris@101:          return t1 * t2;
Chris@101:       }
Chris@101:       return x * boost::math::falling_factorial(x - 1, n - 1);
Chris@101:    }
Chris@101:    if(x <= n - 1)
Chris@16:    {
Chris@16:       //
Chris@16:       // x+1-n will be negative and tgamma_delta_ratio won't
Chris@16:       // handle it, split the product up into three parts:
Chris@16:       //
Chris@16:       T xp1 = x + 1;
Chris@16:       unsigned n2 = itrunc((T)floor(xp1), pol);
Chris@16:       if(n2 == xp1)
Chris@16:          return 0;
Chris@16:       T result = boost::math::tgamma_delta_ratio(xp1, -static_cast<T>(n2), pol);
Chris@16:       x -= n2;
Chris@16:       result *= x;
Chris@16:       ++n2;
Chris@16:       if(n2 < n)
Chris@16:          result *= falling_factorial(x - 1, n - n2, pol);
Chris@16:       return result;
Chris@16:    }
Chris@16:    //
Chris@16:    // Simple case: just the ratio of two
Chris@16:    // (positive argument) gamma functions.
Chris@16:    // Note that we don't optimise this for small n,
Chris@16:    // because tgamma_delta_ratio is alreay optimised
Chris@16:    // for that use case:
Chris@16:    //
Chris@16:    return boost::math::tgamma_delta_ratio(x + 1, -static_cast<T>(n), pol);
Chris@16: }
Chris@16: 
Chris@16: } // namespace detail
Chris@16: 
Chris@16: template <class RT>
Chris@16: inline typename tools::promote_args<RT>::type
Chris@16:    falling_factorial(RT x, unsigned n)
Chris@16: {
Chris@16:    typedef typename tools::promote_args<RT>::type result_type;
Chris@16:    return detail::falling_factorial_imp(
Chris@16:       static_cast<result_type>(x), n, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class RT, class Policy>
Chris@16: inline typename tools::promote_args<RT>::type
Chris@16:    falling_factorial(RT x, unsigned n, const Policy& pol)
Chris@16: {
Chris@16:    typedef typename tools::promote_args<RT>::type result_type;
Chris@16:    return detail::falling_factorial_imp(
Chris@16:       static_cast<result_type>(x), n, pol);
Chris@16: }
Chris@16: 
Chris@16: template <class RT>
Chris@16: inline typename tools::promote_args<RT>::type
Chris@16:    rising_factorial(RT x, int n)
Chris@16: {
Chris@16:    typedef typename tools::promote_args<RT>::type result_type;
Chris@16:    return detail::rising_factorial_imp(
Chris@16:       static_cast<result_type>(x), n, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class RT, class Policy>
Chris@16: inline typename tools::promote_args<RT>::type
Chris@16:    rising_factorial(RT x, int n, const Policy& pol)
Chris@16: {
Chris@16:    typedef typename tools::promote_args<RT>::type result_type;
Chris@16:    return detail::rising_factorial_imp(
Chris@16:       static_cast<result_type>(x), n, pol);
Chris@16: }
Chris@16: 
Chris@16: } // namespace math
Chris@16: } // namespace boost
Chris@16: 
Chris@16: #endif // BOOST_MATH_SP_FACTORIALS_HPP
Chris@16: