ceph/src/boost/libs/math/include/boost/math/special_functions/factorials.hpp

   1 //  Copyright John Maddock 2006, 2010.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_SP_FACTORIALS_HPP
   7 #define BOOST_MATH_SP_FACTORIALS_HPP
   8
   9 #ifdef _MSC_VER
  10 #pragma once
  11 #endif
  12
  13 #include <boost/math/special_functions/math_fwd.hpp>
  14 #include <boost/math/special_functions/gamma.hpp>
  15 #include <boost/math/special_functions/detail/unchecked_factorial.hpp>
  16 #include <boost/array.hpp>
  17 #ifdef BOOST_MSVC
  18 #pragma warning(push) // Temporary until lexical cast fixed.
  19 #pragma warning(disable: 4127 4701)
  20 #endif
  21 #ifdef BOOST_MSVC
  22 #pragma warning(pop)
  23 #endif
  24 #include <boost/config/no_tr1/cmath.hpp>
  25
  26 namespace boost { namespace math
  27 {
  28
  29 template <class T, class Policy>
  30 inline T factorial(unsigned i, const Policy& pol)
  31 {
  32    BOOST_STATIC_ASSERT(!boost::is_integral<T>::value);
  33    // factorial<unsigned int>(n) is not implemented
  34    // because it would overflow integral type T for too small n
  35    // to be useful. Use instead a floating-point type,
  36    // and convert to an unsigned type if essential, for example:
  37    // unsigned int nfac = static_cast<unsigned int>(factorial<double>(n));
  38    // See factorial documentation for more detail.
  39
  40    BOOST_MATH_STD_USING // Aid ADL for floor.
  41
  42    if(i <= max_factorial<T>::value)
  43       return unchecked_factorial<T>(i);
  44    T result = boost::math::tgamma(static_cast<T>(i+1), pol);
  45    if(result > tools::max_value<T>())
  46       return result; // Overflowed value! (But tgamma will have signalled the error already).
  47    return floor(result + 0.5f);
  48 }
  49
  50 template <class T>
  51 inline T factorial(unsigned i)
  52 {
  53    return factorial<T>(i, policies::policy<>());
  54 }
  55 /*
  56 // Can't have these in a policy enabled world?
  57 template<>
  58 inline float factorial<float>(unsigned i)
  59 {
  60    if(i <= max_factorial<float>::value)
  61       return unchecked_factorial<float>(i);
  62    return tools::overflow_error<float>(BOOST_CURRENT_FUNCTION);
  63 }
  64
  65 template<>
  66 inline double factorial<double>(unsigned i)
  67 {
  68    if(i <= max_factorial<double>::value)
  69       return unchecked_factorial<double>(i);
  70    return tools::overflow_error<double>(BOOST_CURRENT_FUNCTION);
  71 }
  72 */
  73 template <class T, class Policy>
  74 T double_factorial(unsigned i, const Policy& pol)
  75 {
  76    BOOST_STATIC_ASSERT(!boost::is_integral<T>::value);
  77    BOOST_MATH_STD_USING  // ADL lookup of std names
  78    if(i & 1)
  79    {
  80       // odd i:
  81       if(i < max_factorial<T>::value)
  82       {
  83          unsigned n = (i - 1) / 2;
  84          return ceil(unchecked_factorial<T>(i) / (ldexp(T(1), (int)n) * unchecked_factorial<T>(n)) - 0.5f);
  85       }
  86       //
  87       // Fallthrough: i is too large to use table lookup, try the
  88       // gamma function instead.
  89       //
  90       T result = boost::math::tgamma(static_cast<T>(i) / 2 + 1, pol) / sqrt(constants::pi<T>());
  91       if(ldexp(tools::max_value<T>(), -static_cast<int>(i+1) / 2) > result)
  92          return ceil(result * ldexp(T(1), static_cast<int>(i+1) / 2) - 0.5f);
  93    }
  94    else
  95    {
  96       // even i:
  97       unsigned n = i / 2;
  98       T result = factorial<T>(n, pol);
  99       if(ldexp(tools::max_value<T>(), -(int)n) > result)
 100          return result * ldexp(T(1), (int)n);
 101    }
 102    //
 103    // If we fall through to here then the result is infinite:
 104    //
 105    return policies::raise_overflow_error<T>("boost::math::double_factorial<%1%>(unsigned)", 0, pol);
 106 }
 107
 108 template <class T>
 109 inline T double_factorial(unsigned i)
 110 {
 111    return double_factorial<T>(i, policies::policy<>());
 112 }
 113
 114 namespace detail{
 115
 116 template <class T, class Policy>
 117 T rising_factorial_imp(T x, int n, const Policy& pol)
 118 {
 119    BOOST_STATIC_ASSERT(!boost::is_integral<T>::value);
 120    if(x < 0)
 121    {
 122       //
 123       // For x less than zero, we really have a falling
 124       // factorial, modulo a possible change of sign.
 125       //
 126       // Note that the falling factorial isn't defined
 127       // for negative n, so we'll get rid of that case
 128       // first:
 129       //
 130       bool inv = false;
 131       if(n < 0)
 132       {
 133          x += n;
 134          n = -n;
 135          inv = true;
 136       }
 137       T result = ((n&1) ? -1 : 1) * falling_factorial(-x, n, pol);
 138       if(inv)
 139          result = 1 / result;
 140       return result;
 141    }
 142    if(n == 0)
 143       return 1;
 144    if(x == 0)
 145    {
 146       if(n < 0)
 147          return -boost::math::tgamma_delta_ratio(x + 1, static_cast<T>(-n), pol);
 148       else
 149          return 0;
 150    }
 151    if((x < 1) && (x + n < 0))
 152    {
 153       T val = boost::math::tgamma_delta_ratio(1 - x, static_cast<T>(-n), pol);
 154       return (n & 1) ? T(-val) : val;
 155    }
 156    //
 157    // We don't optimise this for small n, because
 158    // tgamma_delta_ratio is alreay optimised for that
 159    // use case:
 160    //
 161    return 1 / boost::math::tgamma_delta_ratio(x, static_cast<T>(n), pol);
 162 }
 163
 164 template <class T, class Policy>
 165 inline T falling_factorial_imp(T x, unsigned n, const Policy& pol)
 166 {
 167    BOOST_STATIC_ASSERT(!boost::is_integral<T>::value);
 168    BOOST_MATH_STD_USING // ADL of std names
 169    if((x == 0) && (n >= 0))
 170       return 0;
 171    if(x < 0)
 172    {
 173       //
 174       // For x < 0 we really have a rising factorial
 175       // modulo a possible change of sign:
 176       //
 177       return (n&1 ? -1 : 1) * rising_factorial(-x, n, pol);
 178    }
 179    if(n == 0)
 180       return 1;
 181    if(x < 0.5f)
 182    {
 183       //
 184       // 1 + x below will throw away digits, so split up calculation:
 185       //
 186       if(n > max_factorial<T>::value - 2)
 187       {
 188          // If the two end of the range are far apart we have a ratio of two very large
 189          // numbers, split the calculation up into two blocks:
 190          T t1 = x * boost::math::falling_factorial(x - 1, max_factorial<T>::value - 2);
 191          T t2 = boost::math::falling_factorial(x - max_factorial<T>::value + 1, n - max_factorial<T>::value + 1);
 192          if(tools::max_value<T>() / fabs(t1) < fabs(t2))
 193             return boost::math::sign(t1) * boost::math::sign(t2) * policies::raise_overflow_error<T>("boost::math::falling_factorial<%1%>", 0, pol);
 194          return t1 * t2;
 195       }
 196       return x * boost::math::falling_factorial(x - 1, n - 1);
 197    }
 198    if(x <= n - 1)
 199    {
 200       //
 201       // x+1-n will be negative and tgamma_delta_ratio won't
 202       // handle it, split the product up into three parts:
 203       //
 204       T xp1 = x + 1;
 205       unsigned n2 = itrunc((T)floor(xp1), pol);
 206       if(n2 == xp1)
 207          return 0;
 208       T result = boost::math::tgamma_delta_ratio(xp1, -static_cast<T>(n2), pol);
 209       x -= n2;
 210       result *= x;
 211       ++n2;
 212       if(n2 < n)
 213          result *= falling_factorial(x - 1, n - n2, pol);
 214       return result;
 215    }
 216    //
 217    // Simple case: just the ratio of two
 218    // (positive argument) gamma functions.
 219    // Note that we don't optimise this for small n,
 220    // because tgamma_delta_ratio is alreay optimised
 221    // for that use case:
 222    //
 223    return boost::math::tgamma_delta_ratio(x + 1, -static_cast<T>(n), pol);
 224 }
 225
 226 } // namespace detail
 227
 228 template <class RT>
 229 inline typename tools::promote_args<RT>::type
 230    falling_factorial(RT x, unsigned n)
 231 {
 232    typedef typename tools::promote_args<RT>::type result_type;
 233    return detail::falling_factorial_imp(
 234       static_cast<result_type>(x), n, policies::policy<>());
 235 }
 236
 237 template <class RT, class Policy>
 238 inline typename tools::promote_args<RT>::type
 239    falling_factorial(RT x, unsigned n, const Policy& pol)
 240 {
 241    typedef typename tools::promote_args<RT>::type result_type;
 242    return detail::falling_factorial_imp(
 243       static_cast<result_type>(x), n, pol);
 244 }
 245
 246 template <class RT>
 247 inline typename tools::promote_args<RT>::type
 248    rising_factorial(RT x, int n)
 249 {
 250    typedef typename tools::promote_args<RT>::type result_type;
 251    return detail::rising_factorial_imp(
 252       static_cast<result_type>(x), n, policies::policy<>());
 253 }
 254
 255 template <class RT, class Policy>
 256 inline typename tools::promote_args<RT>::type
 257    rising_factorial(RT x, int n, const Policy& pol)
 258 {
 259    typedef typename tools::promote_args<RT>::type result_type;
 260    return detail::rising_factorial_imp(
 261       static_cast<result_type>(x), n, pol);
 262 }
 263
 264 } // namespace math
 265 } // namespace boost
 266
 267 #endif // BOOST_MATH_SP_FACTORIALS_HPP
 268