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

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