// Copyright John Maddock 2006, 2010. | |
// 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); | |
// factorial<unsigned int>(n) is not implemented | |
// because it would overflow integral type T for too small n | |
// to be useful. Use instead a floating-point type, | |
// and convert to an unsigned type if essential, for example: | |
// unsigned int nfac = static_cast<unsigned int>(factorial<double>(n)); | |
// See factorial documentation for more detail. | |
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), static_cast<int>(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 | |