math/special_functions/factorials.hpp

0001 //  Copyright John Maddock 2006, 2010.
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_SP_FACTORIALS_HPP
0007 #define BOOST_MATH_SP_FACTORIALS_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/gamma.hpp>
0015 #include <boost/math/special_functions/detail/unchecked_factorial.hpp>
0016 #include <array>
0017 #ifdef _MSC_VER
0018 #pragma warning(push) // Temporary until lexical cast fixed.
0019 #pragma warning(disable: 4127 4701)
0020 #endif
0021 #ifdef _MSC_VER
0022 #pragma warning(pop)
0023 #endif
0024 #include <type_traits>
0025 #include <cmath>
0026
0027 namespace boost { namespace math
0028 {
0029
0030 template <class T, class Policy>
0031 inline T factorial(unsigned i, const Policy& pol)
0032 {
0033    static_assert(!std::is_integral<T>::value, "Type T must not be an integral type");
0034    // factorial<unsigned int>(n) is not implemented
0035    // because it would overflow integral type T for too small n
0036    // to be useful. Use instead a floating-point type,
0037    // and convert to an unsigned type if essential, for example:
0038    // unsigned int nfac = static_cast<unsigned int>(factorial<double>(n));
0039    // See factorial documentation for more detail.
0040
0041    BOOST_MATH_STD_USING // Aid ADL for floor.
0042
0043    if(i <= max_factorial<T>::value)
0044       return unchecked_factorial<T>(i);
0045    T result = boost::math::tgamma(static_cast<T>(i+1), pol);
0046    if(result > tools::max_value<T>())
0047       return result; // Overflowed value! (But tgamma will have signalled the error already).
0048    return floor(result + 0.5f);
0049 }
0050
0051 template <class T>
0052 inline T factorial(unsigned i)
0053 {
0054    return factorial<T>(i, policies::policy<>());
0055 }
0056 /*
0057 // Can't have these in a policy enabled world?
0058 template<>
0059 inline float factorial<float>(unsigned i)
0060 {
0061    if(i <= max_factorial<float>::value)
0062       return unchecked_factorial<float>(i);
0063    return tools::overflow_error<float>(BOOST_CURRENT_FUNCTION);
0064 }
0065
0066 template<>
0067 inline double factorial<double>(unsigned i)
0068 {
0069    if(i <= max_factorial<double>::value)
0070       return unchecked_factorial<double>(i);
0071    return tools::overflow_error<double>(BOOST_CURRENT_FUNCTION);
0072 }
0073 */
0074 template <class T, class Policy>
0075 T double_factorial(unsigned i, const Policy& pol)
0076 {
0077    static_assert(!std::is_integral<T>::value, "Type T must not be an integral type");
0078    BOOST_MATH_STD_USING  // ADL lookup of std names
0079    if(i & 1)
0080    {
0081       // odd i:
0082       if(i < max_factorial<T>::value)
0083       {
0084          unsigned n = (i - 1) / 2;
0085          return ceil(unchecked_factorial<T>(i) / (ldexp(T(1), (int)n) * unchecked_factorial<T>(n)) - 0.5f);
0086       }
0087       //
0088       // Fallthrough: i is too large to use table lookup, try the
0089       // gamma function instead.
0090       //
0091       T result = boost::math::tgamma(static_cast<T>(i) / 2 + 1, pol) / sqrt(constants::pi<T>());
0092       if(ldexp(tools::max_value<T>(), -static_cast<int>(i+1) / 2) > result)
0093          return ceil(result * ldexp(T(1), static_cast<int>(i+1) / 2) - 0.5f);
0094    }
0095    else
0096    {
0097       // even i:
0098       unsigned n = i / 2;
0099       T result = factorial<T>(n, pol);
0100       if(ldexp(tools::max_value<T>(), -(int)n) > result)
0101          return result * ldexp(T(1), (int)n);
0102    }
0103    //
0104    // If we fall through to here then the result is infinite:
0105    //
0106    return policies::raise_overflow_error<T>("boost::math::double_factorial<%1%>(unsigned)", 0, pol);
0107 }
0108
0109 template <class T>
0110 inline T double_factorial(unsigned i)
0111 {
0112    return double_factorial<T>(i, policies::policy<>());
0113 }
0114
0115 namespace detail{
0116
0117 template <class T, class Policy>
0118 T rising_factorial_imp(T x, int n, const Policy& pol)
0119 {
0120    static_assert(!std::is_integral<T>::value, "Type T must not be an integral type");
0121    if(x < 0)
0122    {
0123       //
0124       // For x less than zero, we really have a falling
0125       // factorial, modulo a possible change of sign.
0126       //
0127       // Note that the falling factorial isn't defined
0128       // for negative n, so we'll get rid of that case
0129       // first:
0130       //
0131       bool inv = false;
0132       if(n < 0)
0133       {
0134          x += n;
0135          n = -n;
0136          inv = true;
0137       }
0138       T result = ((n&1) ? -1 : 1) * falling_factorial(-x, n, pol);
0139       if(inv)
0140          result = 1 / result;
0141       return result;
0142    }
0143    if(n == 0)
0144       return 1;
0145    if(x == 0)
0146    {
0147       if(n < 0)
0148          return static_cast<T>(-boost::math::tgamma_delta_ratio(x + 1, static_cast<T>(-n), pol));
0149       else
0150          return 0;
0151    }
0152    if((x < 1) && (x + n < 0))
0153    {
0154       const auto val = static_cast<T>(boost::math::tgamma_delta_ratio(1 - x, static_cast<T>(-n), pol));
0155       return (n & 1) ? T(-val) : val;
0156    }
0157    //
0158    // We don't optimise this for small n, because
0159    // tgamma_delta_ratio is already optimised for that
0160    // use case:
0161    //
0162    return 1 / static_cast<T>(boost::math::tgamma_delta_ratio(x, static_cast<T>(n), pol));
0163 }
0164
0165 template <class T, class Policy>
0166 inline T falling_factorial_imp(T x, unsigned n, const Policy& pol)
0167 {
0168    static_assert(!std::is_integral<T>::value, "Type T must not be an integral type");
0169    BOOST_MATH_STD_USING // ADL of std names
0170    if(x == 0)
0171       return 0;
0172    if(x < 0)
0173    {
0174       //
0175       // For x < 0 we really have a rising factorial
0176       // modulo a possible change of sign:
0177       //
0178       return (n&1 ? -1 : 1) * rising_factorial(-x, n, pol);
0179    }
0180    if(n == 0)
0181       return 1;
0182    if(x < 0.5f)
0183    {
0184       //
0185       // 1 + x below will throw away digits, so split up calculation:
0186       //
0187       if(n > max_factorial<T>::value - 2)
0188       {
0189          // If the two end of the range are far apart we have a ratio of two very large
0190          // numbers, split the calculation up into two blocks:
0191          T t1 = x * boost::math::falling_factorial(x - 1, max_factorial<T>::value - 2, pol);
0192          T t2 = boost::math::falling_factorial(x - max_factorial<T>::value + 1, n - max_factorial<T>::value + 1, pol);
0193          if(tools::max_value<T>() / fabs(t1) < fabs(t2))
0194             return boost::math::sign(t1) * boost::math::sign(t2) * policies::raise_overflow_error<T>("boost::math::falling_factorial<%1%>", 0, pol);
0195          return t1 * t2;
0196       }
0197       return x * boost::math::falling_factorial(x - 1, n - 1, pol);
0198    }
0199    if(x <= n - 1)
0200    {
0201       //
0202       // x+1-n will be negative and tgamma_delta_ratio won't
0203       // handle it, split the product up into three parts:
0204       //
0205       T xp1 = x + 1;
0206       unsigned n2 = itrunc((T)floor(xp1), pol);
0207       if(n2 == xp1)
0208          return 0;
0209       auto result = static_cast<T>(boost::math::tgamma_delta_ratio(xp1, -static_cast<T>(n2), pol));
0210       x -= n2;
0211       result *= x;
0212       ++n2;
0213       if(n2 < n)
0214          result *= falling_factorial(x - 1, n - n2, pol);
0215       return result;
0216    }
0217    //
0218    // Simple case: just the ratio of two
0219    // (positive argument) gamma functions.
0220    // Note that we don't optimise this for small n,
0221    // because tgamma_delta_ratio is already optimised
0222    // for that use case:
0223    //
0224    return static_cast<T>(boost::math::tgamma_delta_ratio(x + 1, -static_cast<T>(n), pol));
0225 }
0226
0227 } // namespace detail
0228
0229 template <class RT>
0230 inline typename tools::promote_args<RT>::type
0231    falling_factorial(RT x, unsigned n)
0232 {
0233    typedef typename tools::promote_args<RT>::type result_type;
0234    return detail::falling_factorial_imp(
0235       static_cast<result_type>(x), n, policies::policy<>());
0236 }
0237
0238 template <class RT, class Policy>
0239 inline typename tools::promote_args<RT>::type
0240    falling_factorial(RT x, unsigned n, const Policy& pol)
0241 {
0242    typedef typename tools::promote_args<RT>::type result_type;
0243    return detail::falling_factorial_imp(
0244       static_cast<result_type>(x), n, pol);
0245 }
0246
0247 template <class RT>
0248 inline typename tools::promote_args<RT>::type
0249    rising_factorial(RT x, int n)
0250 {
0251    typedef typename tools::promote_args<RT>::type result_type;
0252    return detail::rising_factorial_imp(
0253       static_cast<result_type>(x), n, policies::policy<>());
0254 }
0255
0256 template <class RT, class Policy>
0257 inline typename tools::promote_args<RT>::type
0258    rising_factorial(RT x, int n, const Policy& pol)
0259 {
0260    typedef typename tools::promote_args<RT>::type result_type;
0261    return detail::rising_factorial_imp(
0262       static_cast<result_type>(x), n, pol);
0263 }
0264
0265 } // namespace math
0266 } // namespace boost
0267
0268 #endif // BOOST_MATH_SP_FACTORIALS_HPP
0269