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
26namespace boost { namespace math
27{
28
29template <class T, class Policy>
30inline 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
50template <class T>
51inline T factorial(unsigned i)
52{
53 return factorial<T>(i, policies::policy<>());
54}
55/*
56// Can't have these in a policy enabled world?
57template<>
58inline 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
65template<>
66inline 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*/
73template <class T, class Policy>
74T 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
108template <class T>
109inline T double_factorial(unsigned i)
110{
111 return double_factorial<T>(i, policies::policy<>());
112}
113
114namespace detail{
115
116template <class T, class Policy>
117T 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 already optimised for that
159 // use case:
160 //
161 return 1 / boost::math::tgamma_delta_ratio(x, static_cast<T>(n), pol);
162}
163
164template <class T, class Policy>
165inline 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)
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 already 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
228template <class RT>
229inline 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
237template <class RT, class Policy>
238inline 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
246template <class RT>
247inline 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
255template <class RT, class Policy>
256inline 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
269

source code of include/boost/math/special_functions/factorials.hpp