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)
6 #ifndef BOOST_MATH_SP_FACTORIALS_HPP
7 #define BOOST_MATH_SP_FACTORIALS_HPP
13 #include <boost/math/special_functions/gamma.hpp>
14 #include <boost/math/special_functions/math_fwd.hpp>
15 #include <boost/math/special_functions/detail/unchecked_factorial.hpp>
16 #include <boost/array.hpp>
18 #pragma warning(push) // Temporary until lexical cast fixed.
19 #pragma warning(disable: 4127 4701)
24 #include <boost/config/no_tr1/cmath.hpp>
26 namespace boost { namespace math
29 template <class T, class Policy>
30 inline T factorial(unsigned i, const Policy& pol)
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.
40 BOOST_MATH_STD_USING // Aid ADL for floor.
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);
51 inline T factorial(unsigned i)
53 return factorial<T>(i, policies::policy<>());
56 // Can't have these in a policy enabled world?
58 inline float factorial<float>(unsigned i)
60 if(i <= max_factorial<float>::value)
61 return unchecked_factorial<float>(i);
62 return tools::overflow_error<float>(BOOST_CURRENT_FUNCTION);
66 inline double factorial<double>(unsigned i)
68 if(i <= max_factorial<double>::value)
69 return unchecked_factorial<double>(i);
70 return tools::overflow_error<double>(BOOST_CURRENT_FUNCTION);
73 template <class T, class Policy>
74 T double_factorial(unsigned i, const Policy& pol)
76 BOOST_STATIC_ASSERT(!boost::is_integral<T>::value);
77 BOOST_MATH_STD_USING // ADL lookup of std names
81 if(i < max_factorial<T>::value)
83 unsigned n = (i - 1) / 2;
84 return ceil(unchecked_factorial<T>(i) / (ldexp(T(1), (int)n) * unchecked_factorial<T>(n)) - 0.5f);
87 // Fallthrough: i is too large to use table lookup, try the
88 // gamma function instead.
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);
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);
103 // If we fall through to here then the result is infinite:
105 return policies::raise_overflow_error<T>("boost::math::double_factorial<%1%>(unsigned)", 0, pol);
109 inline T double_factorial(unsigned i)
111 return double_factorial<T>(i, policies::policy<>());
116 template <class T, class Policy>
117 T rising_factorial_imp(T x, int n, const Policy& pol)
119 BOOST_STATIC_ASSERT(!boost::is_integral<T>::value);
123 // For x less than zero, we really have a falling
124 // factorial, modulo a possible change of sign.
126 // Note that the falling factorial isn't defined
127 // for negative n, so we'll get rid of that case
137 T result = ((n&1) ? -1 : 1) * falling_factorial(-x, n, pol);
145 // We don't optimise this for small n, because
146 // tgamma_delta_ratio is alreay optimised for that
149 return 1 / boost::math::tgamma_delta_ratio(x, static_cast<T>(n), pol);
152 template <class T, class Policy>
153 inline T falling_factorial_imp(T x, unsigned n, const Policy& pol)
155 BOOST_STATIC_ASSERT(!boost::is_integral<T>::value);
156 BOOST_MATH_STD_USING // ADL of std names
162 // For x < 0 we really have a rising factorial
163 // modulo a possible change of sign:
165 return (n&1 ? -1 : 1) * rising_factorial(-x, n, pol);
172 // x+1-n will be negative and tgamma_delta_ratio won't
173 // handle it, split the product up into three parts:
176 unsigned n2 = itrunc((T)floor(xp1), pol);
179 T result = boost::math::tgamma_delta_ratio(xp1, -static_cast<T>(n2), pol);
184 result *= falling_factorial(x - 1, n - n2, pol);
188 // Simple case: just the ratio of two
189 // (positive argument) gamma functions.
190 // Note that we don't optimise this for small n,
191 // because tgamma_delta_ratio is alreay optimised
192 // for that use case:
194 return boost::math::tgamma_delta_ratio(x + 1, -static_cast<T>(n), pol);
197 } // namespace detail
200 inline typename tools::promote_args<RT>::type
201 falling_factorial(RT x, unsigned n)
203 typedef typename tools::promote_args<RT>::type result_type;
204 return detail::falling_factorial_imp(
205 static_cast<result_type>(x), n, policies::policy<>());
208 template <class RT, class Policy>
209 inline typename tools::promote_args<RT>::type
210 falling_factorial(RT x, unsigned n, const Policy& pol)
212 typedef typename tools::promote_args<RT>::type result_type;
213 return detail::falling_factorial_imp(
214 static_cast<result_type>(x), n, pol);
218 inline typename tools::promote_args<RT>::type
219 rising_factorial(RT x, int n)
221 typedef typename tools::promote_args<RT>::type result_type;
222 return detail::rising_factorial_imp(
223 static_cast<result_type>(x), n, policies::policy<>());
226 template <class RT, class Policy>
227 inline typename tools::promote_args<RT>::type
228 rising_factorial(RT x, int n, const Policy& pol)
230 typedef typename tools::promote_args<RT>::type result_type;
231 return detail::rising_factorial_imp(
232 static_cast<result_type>(x), n, pol);
238 #endif // BOOST_MATH_SP_FACTORIALS_HPP