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