]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 |