1 | |
2 | // Copyright John Maddock 2006-7, 2013-14. |
3 | // Copyright Paul A. Bristow 2007, 2013-14. |
4 | // Copyright Nikhar Agrawal 2013-14 |
5 | // Copyright Christopher Kormanyos 2013-14 |
6 | |
7 | // Use, modification and distribution are subject to the |
8 | // Boost Software License, Version 1.0. (See accompanying file |
9 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) |
10 | |
11 | #ifndef BOOST_MATH_SF_GAMMA_HPP |
12 | #define BOOST_MATH_SF_GAMMA_HPP |
13 | |
14 | #ifdef _MSC_VER |
15 | #pragma once |
16 | #endif |
17 | |
18 | #include <boost/config.hpp> |
19 | #include <boost/math/tools/series.hpp> |
20 | #include <boost/math/tools/fraction.hpp> |
21 | #include <boost/math/tools/precision.hpp> |
22 | #include <boost/math/tools/promotion.hpp> |
23 | #include <boost/math/policies/error_handling.hpp> |
24 | #include <boost/math/constants/constants.hpp> |
25 | #include <boost/math/special_functions/math_fwd.hpp> |
26 | #include <boost/math/special_functions/log1p.hpp> |
27 | #include <boost/math/special_functions/trunc.hpp> |
28 | #include <boost/math/special_functions/powm1.hpp> |
29 | #include <boost/math/special_functions/sqrt1pm1.hpp> |
30 | #include <boost/math/special_functions/lanczos.hpp> |
31 | #include <boost/math/special_functions/fpclassify.hpp> |
32 | #include <boost/math/special_functions/detail/igamma_large.hpp> |
33 | #include <boost/math/special_functions/detail/unchecked_factorial.hpp> |
34 | #include <boost/math/special_functions/detail/lgamma_small.hpp> |
35 | #include <boost/math/special_functions/bernoulli.hpp> |
36 | #include <boost/math/special_functions/polygamma.hpp> |
37 | #include <boost/type_traits/is_convertible.hpp> |
38 | #include <boost/assert.hpp> |
39 | #include <boost/mpl/greater.hpp> |
40 | #include <boost/mpl/equal_to.hpp> |
41 | #include <boost/mpl/greater.hpp> |
42 | |
43 | #include <boost/config/no_tr1/cmath.hpp> |
44 | #include <algorithm> |
45 | |
46 | #ifdef BOOST_MSVC |
47 | # pragma warning(push) |
48 | # pragma warning(disable: 4702) // unreachable code (return after domain_error throw). |
49 | # pragma warning(disable: 4127) // conditional expression is constant. |
50 | # pragma warning(disable: 4100) // unreferenced formal parameter. |
51 | // Several variables made comments, |
52 | // but some difficulty as whether referenced on not may depend on macro values. |
53 | // So to be safe, 4100 warnings suppressed. |
54 | // TODO - revisit this? |
55 | #endif |
56 | |
57 | namespace boost{ namespace math{ |
58 | |
59 | namespace detail{ |
60 | |
61 | template <class T> |
62 | inline bool is_odd(T v, const boost::true_type&) |
63 | { |
64 | int i = static_cast<int>(v); |
65 | return i&1; |
66 | } |
67 | template <class T> |
68 | inline bool is_odd(T v, const boost::false_type&) |
69 | { |
70 | // Oh dear can't cast T to int! |
71 | BOOST_MATH_STD_USING |
72 | T modulus = v - 2 * floor(v/2); |
73 | return static_cast<bool>(modulus != 0); |
74 | } |
75 | template <class T> |
76 | inline bool is_odd(T v) |
77 | { |
78 | return is_odd(v, ::boost::is_convertible<T, int>()); |
79 | } |
80 | |
81 | template <class T> |
82 | T sinpx(T z) |
83 | { |
84 | // Ad hoc function calculates x * sin(pi * x), |
85 | // taking extra care near when x is near a whole number. |
86 | BOOST_MATH_STD_USING |
87 | int sign = 1; |
88 | if(z < 0) |
89 | { |
90 | z = -z; |
91 | } |
92 | T fl = floor(z); |
93 | T dist; |
94 | if(is_odd(fl)) |
95 | { |
96 | fl += 1; |
97 | dist = fl - z; |
98 | sign = -sign; |
99 | } |
100 | else |
101 | { |
102 | dist = z - fl; |
103 | } |
104 | BOOST_ASSERT(fl >= 0); |
105 | if(dist > 0.5) |
106 | dist = 1 - dist; |
107 | T result = sin(dist*boost::math::constants::pi<T>()); |
108 | return sign*z*result; |
109 | } // template <class T> T sinpx(T z) |
110 | // |
111 | // tgamma(z), with Lanczos support: |
112 | // |
113 | template <class T, class Policy, class Lanczos> |
114 | T gamma_imp(T z, const Policy& pol, const Lanczos& l) |
115 | { |
116 | BOOST_MATH_STD_USING |
117 | |
118 | T result = 1; |
119 | |
120 | #ifdef BOOST_MATH_INSTRUMENT |
121 | static bool b = false; |
122 | if(!b) |
123 | { |
124 | std::cout << "tgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl; |
125 | b = true; |
126 | } |
127 | #endif |
128 | static const char* function = "boost::math::tgamma<%1%>(%1%)" ; |
129 | |
130 | if(z <= 0) |
131 | { |
132 | if(floor(z) == z) |
133 | return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%." , z, pol); |
134 | if(z <= -20) |
135 | { |
136 | result = gamma_imp(T(-z), pol, l) * sinpx(z); |
137 | BOOST_MATH_INSTRUMENT_VARIABLE(result); |
138 | if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>())) |
139 | return -boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent." , pol); |
140 | result = -boost::math::constants::pi<T>() / result; |
141 | if(result == 0) |
142 | return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent." , pol); |
143 | if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL) |
144 | return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized." , result, pol); |
145 | BOOST_MATH_INSTRUMENT_VARIABLE(result); |
146 | return result; |
147 | } |
148 | |
149 | // shift z to > 1: |
150 | while(z < 0) |
151 | { |
152 | result /= z; |
153 | z += 1; |
154 | } |
155 | } |
156 | BOOST_MATH_INSTRUMENT_VARIABLE(result); |
157 | if((floor(z) == z) && (z < max_factorial<T>::value)) |
158 | { |
159 | result *= unchecked_factorial<T>(itrunc(z, pol) - 1); |
160 | BOOST_MATH_INSTRUMENT_VARIABLE(result); |
161 | } |
162 | else if (z < tools::root_epsilon<T>()) |
163 | { |
164 | if (z < 1 / tools::max_value<T>()) |
165 | result = policies::raise_overflow_error<T>(function, 0, pol); |
166 | result *= 1 / z - constants::euler<T>(); |
167 | } |
168 | else |
169 | { |
170 | result *= Lanczos::lanczos_sum(z); |
171 | T zgh = (z + static_cast<T>(Lanczos::g()) - boost::math::constants::half<T>()); |
172 | T lzgh = log(zgh); |
173 | BOOST_MATH_INSTRUMENT_VARIABLE(result); |
174 | BOOST_MATH_INSTRUMENT_VARIABLE(tools::log_max_value<T>()); |
175 | if(z * lzgh > tools::log_max_value<T>()) |
176 | { |
177 | // we're going to overflow unless this is done with care: |
178 | BOOST_MATH_INSTRUMENT_VARIABLE(zgh); |
179 | if(lzgh * z / 2 > tools::log_max_value<T>()) |
180 | return boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent." , pol); |
181 | T hp = pow(zgh, (z / 2) - T(0.25)); |
182 | BOOST_MATH_INSTRUMENT_VARIABLE(hp); |
183 | result *= hp / exp(zgh); |
184 | BOOST_MATH_INSTRUMENT_VARIABLE(result); |
185 | if(tools::max_value<T>() / hp < result) |
186 | return boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent." , pol); |
187 | result *= hp; |
188 | BOOST_MATH_INSTRUMENT_VARIABLE(result); |
189 | } |
190 | else |
191 | { |
192 | BOOST_MATH_INSTRUMENT_VARIABLE(zgh); |
193 | BOOST_MATH_INSTRUMENT_VARIABLE(pow(zgh, z - boost::math::constants::half<T>())); |
194 | BOOST_MATH_INSTRUMENT_VARIABLE(exp(zgh)); |
195 | result *= pow(zgh, z - boost::math::constants::half<T>()) / exp(zgh); |
196 | BOOST_MATH_INSTRUMENT_VARIABLE(result); |
197 | } |
198 | } |
199 | return result; |
200 | } |
201 | // |
202 | // lgamma(z) with Lanczos support: |
203 | // |
204 | template <class T, class Policy, class Lanczos> |
205 | T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = 0) |
206 | { |
207 | #ifdef BOOST_MATH_INSTRUMENT |
208 | static bool b = false; |
209 | if(!b) |
210 | { |
211 | std::cout << "lgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl; |
212 | b = true; |
213 | } |
214 | #endif |
215 | |
216 | BOOST_MATH_STD_USING |
217 | |
218 | static const char* function = "boost::math::lgamma<%1%>(%1%)" ; |
219 | |
220 | T result = 0; |
221 | int sresult = 1; |
222 | if(z <= -tools::root_epsilon<T>()) |
223 | { |
224 | // reflection formula: |
225 | if(floor(z) == z) |
226 | return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%." , z, pol); |
227 | |
228 | T t = sinpx(z); |
229 | z = -z; |
230 | if(t < 0) |
231 | { |
232 | t = -t; |
233 | } |
234 | else |
235 | { |
236 | sresult = -sresult; |
237 | } |
238 | result = log(boost::math::constants::pi<T>()) - lgamma_imp(z, pol, l) - log(t); |
239 | } |
240 | else if (z < tools::root_epsilon<T>()) |
241 | { |
242 | if (0 == z) |
243 | return policies::raise_pole_error<T>(function, "Evaluation of lgamma at %1%." , z, pol); |
244 | if (fabs(z) < 1 / tools::max_value<T>()) |
245 | result = -log(fabs(z)); |
246 | else |
247 | result = log(fabs(1 / z - constants::euler<T>())); |
248 | if (z < 0) |
249 | sresult = -1; |
250 | } |
251 | else if(z < 15) |
252 | { |
253 | typedef typename policies::precision<T, Policy>::type precision_type; |
254 | typedef boost::integral_constant<int, |
255 | precision_type::value <= 0 ? 0 : |
256 | precision_type::value <= 64 ? 64 : |
257 | precision_type::value <= 113 ? 113 : 0 |
258 | > tag_type; |
259 | |
260 | result = lgamma_small_imp<T>(z, T(z - 1), T(z - 2), tag_type(), pol, l); |
261 | } |
262 | else if((z >= 3) && (z < 100) && (std::numeric_limits<T>::max_exponent >= 1024)) |
263 | { |
264 | // taking the log of tgamma reduces the error, no danger of overflow here: |
265 | result = log(gamma_imp(z, pol, l)); |
266 | } |
267 | else |
268 | { |
269 | // regular evaluation: |
270 | T zgh = static_cast<T>(z + Lanczos::g() - boost::math::constants::half<T>()); |
271 | result = log(zgh) - 1; |
272 | result *= z - 0.5f; |
273 | // |
274 | // Only add on the lanczos sum part if we're going to need it: |
275 | // |
276 | if(result * tools::epsilon<T>() < 20) |
277 | result += log(Lanczos::lanczos_sum_expG_scaled(z)); |
278 | } |
279 | |
280 | if(sign) |
281 | *sign = sresult; |
282 | return result; |
283 | } |
284 | |
285 | // |
286 | // Incomplete gamma functions follow: |
287 | // |
288 | template <class T> |
289 | struct upper_incomplete_gamma_fract |
290 | { |
291 | private: |
292 | T z, a; |
293 | int k; |
294 | public: |
295 | typedef std::pair<T,T> result_type; |
296 | |
297 | upper_incomplete_gamma_fract(T a1, T z1) |
298 | : z(z1-a1+1), a(a1), k(0) |
299 | { |
300 | } |
301 | |
302 | result_type operator()() |
303 | { |
304 | ++k; |
305 | z += 2; |
306 | return result_type(k * (a - k), z); |
307 | } |
308 | }; |
309 | |
310 | template <class T> |
311 | inline T upper_gamma_fraction(T a, T z, T eps) |
312 | { |
313 | // Multiply result by z^a * e^-z to get the full |
314 | // upper incomplete integral. Divide by tgamma(z) |
315 | // to normalise. |
316 | upper_incomplete_gamma_fract<T> f(a, z); |
317 | return 1 / (z - a + 1 + boost::math::tools::continued_fraction_a(f, eps)); |
318 | } |
319 | |
320 | template <class T> |
321 | struct lower_incomplete_gamma_series |
322 | { |
323 | private: |
324 | T a, z, result; |
325 | public: |
326 | typedef T result_type; |
327 | lower_incomplete_gamma_series(T a1, T z1) : a(a1), z(z1), result(1){} |
328 | |
329 | T operator()() |
330 | { |
331 | T r = result; |
332 | a += 1; |
333 | result *= z/a; |
334 | return r; |
335 | } |
336 | }; |
337 | |
338 | template <class T, class Policy> |
339 | inline T lower_gamma_series(T a, T z, const Policy& pol, T init_value = 0) |
340 | { |
341 | // Multiply result by ((z^a) * (e^-z) / a) to get the full |
342 | // lower incomplete integral. Then divide by tgamma(a) |
343 | // to get the normalised value. |
344 | lower_incomplete_gamma_series<T> s(a, z); |
345 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); |
346 | T factor = policies::get_epsilon<T, Policy>(); |
347 | T result = boost::math::tools::sum_series(s, factor, max_iter, init_value); |
348 | policies::check_series_iterations<T>("boost::math::detail::lower_gamma_series<%1%>(%1%)" , max_iter, pol); |
349 | return result; |
350 | } |
351 | |
352 | // |
353 | // Fully generic tgamma and lgamma use Stirling's approximation |
354 | // with Bernoulli numbers. |
355 | // |
356 | template<class T> |
357 | std::size_t highest_bernoulli_index() |
358 | { |
359 | const float digits10_of_type = (std::numeric_limits<T>::is_specialized |
360 | ? static_cast<float>(std::numeric_limits<T>::digits10) |
361 | : static_cast<float>(boost::math::tools::digits<T>() * 0.301F)); |
362 | |
363 | // Find the high index n for Bn to produce the desired precision in Stirling's calculation. |
364 | return static_cast<std::size_t>(18.0F + (0.6F * digits10_of_type)); |
365 | } |
366 | |
367 | template<class T> |
368 | int minimum_argument_for_bernoulli_recursion() |
369 | { |
370 | const float digits10_of_type = (std::numeric_limits<T>::is_specialized |
371 | ? static_cast<float>(std::numeric_limits<T>::digits10) |
372 | : static_cast<float>(boost::math::tools::digits<T>() * 0.301F)); |
373 | |
374 | const float limit = std::ceil(std::pow(1.0f / std::ldexp(1.0f, 1-boost::math::tools::digits<T>()), 1.0f / 20.0f)); |
375 | |
376 | return (int)((std::min)(a: digits10_of_type * 1.7F, b: limit)); |
377 | } |
378 | |
379 | template <class T, class Policy> |
380 | T scaled_tgamma_no_lanczos(const T& z, const Policy& pol, bool islog = false) |
381 | { |
382 | BOOST_MATH_STD_USING |
383 | // |
384 | // Calculates tgamma(z) / (z/e)^z |
385 | // Requires that our argument is large enough for Sterling's approximation to hold. |
386 | // Used internally when combining gamma's of similar magnitude without logarithms. |
387 | // |
388 | BOOST_ASSERT(minimum_argument_for_bernoulli_recursion<T>() <= z); |
389 | |
390 | // Perform the Bernoulli series expansion of Stirling's approximation. |
391 | |
392 | const std::size_t number_of_bernoullis_b2n = policies::get_max_series_iterations<Policy>(); |
393 | |
394 | T one_over_x_pow_two_n_minus_one = 1 / z; |
395 | const T one_over_x2 = one_over_x_pow_two_n_minus_one * one_over_x_pow_two_n_minus_one; |
396 | T sum = (boost::math::bernoulli_b2n<T>(1) / 2) * one_over_x_pow_two_n_minus_one; |
397 | const T target_epsilon_to_break_loop = sum * boost::math::tools::epsilon<T>(); |
398 | const T half_ln_two_pi_over_z = sqrt(boost::math::constants::two_pi<T>() / z); |
399 | T last_term = 2 * sum; |
400 | |
401 | for (std::size_t n = 2U;; ++n) |
402 | { |
403 | one_over_x_pow_two_n_minus_one *= one_over_x2; |
404 | |
405 | const std::size_t n2 = static_cast<std::size_t>(n * 2U); |
406 | |
407 | const T term = (boost::math::bernoulli_b2n<T>(static_cast<int>(n)) * one_over_x_pow_two_n_minus_one) / (n2 * (n2 - 1U)); |
408 | |
409 | if ((n >= 3U) && (abs(term) < target_epsilon_to_break_loop)) |
410 | { |
411 | // We have reached the desired precision in Stirling's expansion. |
412 | // Adding additional terms to the sum of this divergent asymptotic |
413 | // expansion will not improve the result. |
414 | |
415 | // Break from the loop. |
416 | break; |
417 | } |
418 | if (n > number_of_bernoullis_b2n) |
419 | return policies::raise_evaluation_error("scaled_tgamma_no_lanczos<%1%>()" , "Exceeded maximum series iterations without reaching convergence, best approximation was %1%" , T(exp(sum) * half_ln_two_pi_over_z), pol); |
420 | |
421 | sum += term; |
422 | |
423 | // Sanity check for divergence: |
424 | T fterm = fabs(term); |
425 | if(fterm > last_term) |
426 | return policies::raise_evaluation_error("scaled_tgamma_no_lanczos<%1%>()" , "Series became divergent without reaching convergence, best approximation was %1%" , T(exp(sum) * half_ln_two_pi_over_z), pol); |
427 | last_term = fterm; |
428 | } |
429 | |
430 | // Complete Stirling's approximation. |
431 | T scaled_gamma_value = islog ? T(sum + log(half_ln_two_pi_over_z)) : T(exp(sum) * half_ln_two_pi_over_z); |
432 | return scaled_gamma_value; |
433 | } |
434 | |
435 | // Forward declaration of the lgamma_imp template specialization. |
436 | template <class T, class Policy> |
437 | T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign = 0); |
438 | |
439 | template <class T, class Policy> |
440 | T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&) |
441 | { |
442 | BOOST_MATH_STD_USING |
443 | |
444 | static const char* function = "boost::math::tgamma<%1%>(%1%)" ; |
445 | |
446 | // Check if the argument of tgamma is identically zero. |
447 | const bool is_at_zero = (z == 0); |
448 | |
449 | if((boost::math::isnan)(z) || (is_at_zero) || ((boost::math::isinf)(z) && (z < 0))) |
450 | return policies::raise_domain_error<T>(function, "Evaluation of tgamma at %1%." , z, pol); |
451 | |
452 | const bool b_neg = (z < 0); |
453 | |
454 | const bool floor_of_z_is_equal_to_z = (floor(z) == z); |
455 | |
456 | // Special case handling of small factorials: |
457 | if((!b_neg) && floor_of_z_is_equal_to_z && (z < boost::math::max_factorial<T>::value)) |
458 | { |
459 | return boost::math::unchecked_factorial<T>(itrunc(z) - 1); |
460 | } |
461 | |
462 | // Make a local, unsigned copy of the input argument. |
463 | T zz((!b_neg) ? z : -z); |
464 | |
465 | // Special case for ultra-small z: |
466 | if(zz < tools::cbrt_epsilon<T>()) |
467 | { |
468 | const T a0(1); |
469 | const T a1(boost::math::constants::euler<T>()); |
470 | const T six_euler_squared((boost::math::constants::euler<T>() * boost::math::constants::euler<T>()) * 6); |
471 | const T a2((six_euler_squared - boost::math::constants::pi_sqr<T>()) / 12); |
472 | |
473 | const T inverse_tgamma_series = z * ((a2 * z + a1) * z + a0); |
474 | |
475 | return 1 / inverse_tgamma_series; |
476 | } |
477 | |
478 | // Scale the argument up for the calculation of lgamma, |
479 | // and use downward recursion later for the final result. |
480 | const int min_arg_for_recursion = minimum_argument_for_bernoulli_recursion<T>(); |
481 | |
482 | int n_recur; |
483 | |
484 | if(zz < min_arg_for_recursion) |
485 | { |
486 | n_recur = boost::math::itrunc(min_arg_for_recursion - zz) + 1; |
487 | |
488 | zz += n_recur; |
489 | } |
490 | else |
491 | { |
492 | n_recur = 0; |
493 | } |
494 | if (!n_recur) |
495 | { |
496 | if (zz > tools::log_max_value<T>()) |
497 | return policies::raise_overflow_error<T>(function, 0, pol); |
498 | if (log(zz) * zz / 2 > tools::log_max_value<T>()) |
499 | return policies::raise_overflow_error<T>(function, 0, pol); |
500 | } |
501 | T gamma_value = scaled_tgamma_no_lanczos(zz, pol); |
502 | T power_term = pow(zz, zz / 2); |
503 | T exp_term = exp(-zz); |
504 | gamma_value *= (power_term * exp_term); |
505 | if(!n_recur && (tools::max_value<T>() / power_term < gamma_value)) |
506 | return policies::raise_overflow_error<T>(function, 0, pol); |
507 | gamma_value *= power_term; |
508 | |
509 | // Rescale the result using downward recursion if necessary. |
510 | if(n_recur) |
511 | { |
512 | // The order of divides is important, if we keep subtracting 1 from zz |
513 | // we DO NOT get back to z (cancellation error). Further if z < epsilon |
514 | // we would end up dividing by zero. Also in order to prevent spurious |
515 | // overflow with the first division, we must save dividing by |z| till last, |
516 | // so the optimal order of divides is z+1, z+2, z+3...z+n_recur-1,z. |
517 | zz = fabs(z) + 1; |
518 | for(int k = 1; k < n_recur; ++k) |
519 | { |
520 | gamma_value /= zz; |
521 | zz += 1; |
522 | } |
523 | gamma_value /= fabs(z); |
524 | } |
525 | |
526 | // Return the result, accounting for possible negative arguments. |
527 | if(b_neg) |
528 | { |
529 | // Provide special error analysis for: |
530 | // * arguments in the neighborhood of a negative integer |
531 | // * arguments exactly equal to a negative integer. |
532 | |
533 | // Check if the argument of tgamma is exactly equal to a negative integer. |
534 | if(floor_of_z_is_equal_to_z) |
535 | return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%." , z, pol); |
536 | |
537 | gamma_value *= sinpx(z); |
538 | |
539 | BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value); |
540 | |
541 | const bool result_is_too_large_to_represent = ( (abs(gamma_value) < 1) |
542 | && ((tools::max_value<T>() * abs(gamma_value)) < boost::math::constants::pi<T>())); |
543 | |
544 | if(result_is_too_large_to_represent) |
545 | return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent." , pol); |
546 | |
547 | gamma_value = -boost::math::constants::pi<T>() / gamma_value; |
548 | BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value); |
549 | |
550 | if(gamma_value == 0) |
551 | return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent." , pol); |
552 | |
553 | if((boost::math::fpclassify)(gamma_value) == static_cast<int>(FP_SUBNORMAL)) |
554 | return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized." , gamma_value, pol); |
555 | } |
556 | |
557 | return gamma_value; |
558 | } |
559 | |
560 | template <class T, class Policy> |
561 | inline T log_gamma_near_1(const T& z, Policy const& pol) |
562 | { |
563 | // |
564 | // This is for the multiprecision case where there is |
565 | // no lanczos support, use a taylor series at z = 1, |
566 | // see https://www.wolframalpha.com/input/?i=taylor+series+lgamma(x)+at+x+%3D+1 |
567 | // |
568 | BOOST_MATH_STD_USING // ADL of std names |
569 | |
570 | BOOST_ASSERT(fabs(z) < 1); |
571 | |
572 | T result = -constants::euler<T>() * z; |
573 | |
574 | T power_term = z * z / 2; |
575 | int n = 2; |
576 | T term = 0; |
577 | |
578 | do |
579 | { |
580 | term = power_term * boost::math::polygamma(n - 1, T(1)); |
581 | result += term; |
582 | ++n; |
583 | power_term *= z / n; |
584 | } while (fabs(result) * tools::epsilon<T>() < fabs(term)); |
585 | |
586 | return result; |
587 | } |
588 | |
589 | template <class T, class Policy> |
590 | T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign) |
591 | { |
592 | BOOST_MATH_STD_USING |
593 | |
594 | static const char* function = "boost::math::lgamma<%1%>(%1%)" ; |
595 | |
596 | // Check if the argument of lgamma is identically zero. |
597 | const bool is_at_zero = (z == 0); |
598 | |
599 | if(is_at_zero) |
600 | return policies::raise_domain_error<T>(function, "Evaluation of lgamma at zero %1%." , z, pol); |
601 | if((boost::math::isnan)(z)) |
602 | return policies::raise_domain_error<T>(function, "Evaluation of lgamma at %1%." , z, pol); |
603 | if((boost::math::isinf)(z)) |
604 | return policies::raise_overflow_error<T>(function, 0, pol); |
605 | |
606 | const bool b_neg = (z < 0); |
607 | |
608 | const bool floor_of_z_is_equal_to_z = (floor(z) == z); |
609 | |
610 | // Special case handling of small factorials: |
611 | if((!b_neg) && floor_of_z_is_equal_to_z && (z < boost::math::max_factorial<T>::value)) |
612 | { |
613 | if (sign) |
614 | *sign = 1; |
615 | return log(boost::math::unchecked_factorial<T>(itrunc(z) - 1)); |
616 | } |
617 | |
618 | // Make a local, unsigned copy of the input argument. |
619 | T zz((!b_neg) ? z : -z); |
620 | |
621 | const int min_arg_for_recursion = minimum_argument_for_bernoulli_recursion<T>(); |
622 | |
623 | T log_gamma_value; |
624 | |
625 | if (zz < min_arg_for_recursion) |
626 | { |
627 | // Here we simply take the logarithm of tgamma(). This is somewhat |
628 | // inefficient, but simple. The rationale is that the argument here |
629 | // is relatively small and overflow is not expected to be likely. |
630 | if (sign) |
631 | * sign = 1; |
632 | if(fabs(z - 1) < 0.25) |
633 | { |
634 | log_gamma_value = log_gamma_near_1(T(zz - 1), pol); |
635 | } |
636 | else if(fabs(z - 2) < 0.25) |
637 | { |
638 | log_gamma_value = log_gamma_near_1(T(zz - 2), pol) + log(zz - 1); |
639 | } |
640 | else if (z > -tools::root_epsilon<T>()) |
641 | { |
642 | // Reflection formula may fail if z is very close to zero, let the series |
643 | // expansion for tgamma close to zero do the work: |
644 | if (sign) |
645 | *sign = z < 0 ? -1 : 1; |
646 | return log(abs(gamma_imp(z, pol, lanczos::undefined_lanczos()))); |
647 | } |
648 | else |
649 | { |
650 | // No issue with spurious overflow in reflection formula, |
651 | // just fall through to regular code: |
652 | T g = gamma_imp(zz, pol, lanczos::undefined_lanczos()); |
653 | if (sign) |
654 | { |
655 | *sign = g < 0 ? -1 : 1; |
656 | } |
657 | log_gamma_value = log(abs(g)); |
658 | } |
659 | } |
660 | else |
661 | { |
662 | // Perform the Bernoulli series expansion of Stirling's approximation. |
663 | T sum = scaled_tgamma_no_lanczos(zz, pol, true); |
664 | log_gamma_value = zz * (log(zz) - 1) + sum; |
665 | } |
666 | |
667 | int sign_of_result = 1; |
668 | |
669 | if(b_neg) |
670 | { |
671 | // Provide special error analysis if the argument is exactly |
672 | // equal to a negative integer. |
673 | |
674 | // Check if the argument of lgamma is exactly equal to a negative integer. |
675 | if(floor_of_z_is_equal_to_z) |
676 | return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%." , z, pol); |
677 | |
678 | T t = sinpx(z); |
679 | |
680 | if(t < 0) |
681 | { |
682 | t = -t; |
683 | } |
684 | else |
685 | { |
686 | sign_of_result = -sign_of_result; |
687 | } |
688 | |
689 | log_gamma_value = - log_gamma_value |
690 | + log(boost::math::constants::pi<T>()) |
691 | - log(t); |
692 | } |
693 | |
694 | if(sign != static_cast<int*>(0U)) { *sign = sign_of_result; } |
695 | |
696 | return log_gamma_value; |
697 | } |
698 | |
699 | // |
700 | // This helper calculates tgamma(dz+1)-1 without cancellation errors, |
701 | // used by the upper incomplete gamma with z < 1: |
702 | // |
703 | template <class T, class Policy, class Lanczos> |
704 | T tgammap1m1_imp(T dz, Policy const& pol, const Lanczos& l) |
705 | { |
706 | BOOST_MATH_STD_USING |
707 | |
708 | typedef typename policies::precision<T,Policy>::type precision_type; |
709 | |
710 | typedef boost::integral_constant<int, |
711 | precision_type::value <= 0 ? 0 : |
712 | precision_type::value <= 64 ? 64 : |
713 | precision_type::value <= 113 ? 113 : 0 |
714 | > tag_type; |
715 | |
716 | T result; |
717 | if(dz < 0) |
718 | { |
719 | if(dz < -0.5) |
720 | { |
721 | // Best method is simply to subtract 1 from tgamma: |
722 | result = boost::math::tgamma(1+dz, pol) - 1; |
723 | BOOST_MATH_INSTRUMENT_CODE(result); |
724 | } |
725 | else |
726 | { |
727 | // Use expm1 on lgamma: |
728 | result = boost::math::expm1(-boost::math::log1p(dz, pol) |
729 | + lgamma_small_imp<T>(dz+2, dz + 1, dz, tag_type(), pol, l)); |
730 | BOOST_MATH_INSTRUMENT_CODE(result); |
731 | } |
732 | } |
733 | else |
734 | { |
735 | if(dz < 2) |
736 | { |
737 | // Use expm1 on lgamma: |
738 | result = boost::math::expm1(lgamma_small_imp<T>(dz+1, dz, dz-1, tag_type(), pol, l), pol); |
739 | BOOST_MATH_INSTRUMENT_CODE(result); |
740 | } |
741 | else |
742 | { |
743 | // Best method is simply to subtract 1 from tgamma: |
744 | result = boost::math::tgamma(1+dz, pol) - 1; |
745 | BOOST_MATH_INSTRUMENT_CODE(result); |
746 | } |
747 | } |
748 | |
749 | return result; |
750 | } |
751 | |
752 | template <class T, class Policy> |
753 | inline T tgammap1m1_imp(T z, Policy const& pol, |
754 | const ::boost::math::lanczos::undefined_lanczos&) |
755 | { |
756 | BOOST_MATH_STD_USING // ADL of std names |
757 | |
758 | if(fabs(z) < 0.55) |
759 | { |
760 | return boost::math::expm1(log_gamma_near_1(z, pol)); |
761 | } |
762 | return boost::math::expm1(boost::math::lgamma(1 + z, pol)); |
763 | } |
764 | |
765 | // |
766 | // Series representation for upper fraction when z is small: |
767 | // |
768 | template <class T> |
769 | struct small_gamma2_series |
770 | { |
771 | typedef T result_type; |
772 | |
773 | small_gamma2_series(T a_, T x_) : result(-x_), x(-x_), apn(a_+1), n(1){} |
774 | |
775 | T operator()() |
776 | { |
777 | T r = result / (apn); |
778 | result *= x; |
779 | result /= ++n; |
780 | apn += 1; |
781 | return r; |
782 | } |
783 | |
784 | private: |
785 | T result, x, apn; |
786 | int n; |
787 | }; |
788 | // |
789 | // calculate power term prefix (z^a)(e^-z) used in the non-normalised |
790 | // incomplete gammas: |
791 | // |
792 | template <class T, class Policy> |
793 | T full_igamma_prefix(T a, T z, const Policy& pol) |
794 | { |
795 | BOOST_MATH_STD_USING |
796 | |
797 | T prefix; |
798 | if (z > tools::max_value<T>()) |
799 | return 0; |
800 | T alz = a * log(z); |
801 | |
802 | if(z >= 1) |
803 | { |
804 | if((alz < tools::log_max_value<T>()) && (-z > tools::log_min_value<T>())) |
805 | { |
806 | prefix = pow(z, a) * exp(-z); |
807 | } |
808 | else if(a >= 1) |
809 | { |
810 | prefix = pow(z / exp(z/a), a); |
811 | } |
812 | else |
813 | { |
814 | prefix = exp(alz - z); |
815 | } |
816 | } |
817 | else |
818 | { |
819 | if(alz > tools::log_min_value<T>()) |
820 | { |
821 | prefix = pow(z, a) * exp(-z); |
822 | } |
823 | else if(z/a < tools::log_max_value<T>()) |
824 | { |
825 | prefix = pow(z / exp(z/a), a); |
826 | } |
827 | else |
828 | { |
829 | prefix = exp(alz - z); |
830 | } |
831 | } |
832 | // |
833 | // This error handling isn't very good: it happens after the fact |
834 | // rather than before it... |
835 | // |
836 | if((boost::math::fpclassify)(prefix) == (int)FP_INFINITE) |
837 | return policies::raise_overflow_error<T>("boost::math::detail::full_igamma_prefix<%1%>(%1%, %1%)" , "Result of incomplete gamma function is too large to represent." , pol); |
838 | |
839 | return prefix; |
840 | } |
841 | // |
842 | // Compute (z^a)(e^-z)/tgamma(a) |
843 | // most if the error occurs in this function: |
844 | // |
845 | template <class T, class Policy, class Lanczos> |
846 | T regularised_gamma_prefix(T a, T z, const Policy& pol, const Lanczos& l) |
847 | { |
848 | BOOST_MATH_STD_USING |
849 | if (z >= tools::max_value<T>()) |
850 | return 0; |
851 | T agh = a + static_cast<T>(Lanczos::g()) - T(0.5); |
852 | T prefix; |
853 | T d = ((z - a) - static_cast<T>(Lanczos::g()) + T(0.5)) / agh; |
854 | |
855 | if(a < 1) |
856 | { |
857 | // |
858 | // We have to treat a < 1 as a special case because our Lanczos |
859 | // approximations are optimised against the factorials with a > 1, |
860 | // and for high precision types especially (128-bit reals for example) |
861 | // very small values of a can give rather erroneous results for gamma |
862 | // unless we do this: |
863 | // |
864 | // TODO: is this still required? Lanczos approx should be better now? |
865 | // |
866 | if(z <= tools::log_min_value<T>()) |
867 | { |
868 | // Oh dear, have to use logs, should be free of cancellation errors though: |
869 | return exp(a * log(z) - z - lgamma_imp(a, pol, l)); |
870 | } |
871 | else |
872 | { |
873 | // direct calculation, no danger of overflow as gamma(a) < 1/a |
874 | // for small a. |
875 | return pow(z, a) * exp(-z) / gamma_imp(a, pol, l); |
876 | } |
877 | } |
878 | else if((fabs(d*d*a) <= 100) && (a > 150)) |
879 | { |
880 | // special case for large a and a ~ z. |
881 | prefix = a * boost::math::log1pmx(d, pol) + z * static_cast<T>(0.5 - Lanczos::g()) / agh; |
882 | prefix = exp(prefix); |
883 | } |
884 | else |
885 | { |
886 | // |
887 | // general case. |
888 | // direct computation is most accurate, but use various fallbacks |
889 | // for different parts of the problem domain: |
890 | // |
891 | T alz = a * log(z / agh); |
892 | T amz = a - z; |
893 | if(((std::min)(alz, amz) <= tools::log_min_value<T>()) || ((std::max)(alz, amz) >= tools::log_max_value<T>())) |
894 | { |
895 | T amza = amz / a; |
896 | if(((std::min)(alz, amz)/2 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/2 < tools::log_max_value<T>())) |
897 | { |
898 | // compute square root of the result and then square it: |
899 | T sq = pow(z / agh, a / 2) * exp(amz / 2); |
900 | prefix = sq * sq; |
901 | } |
902 | else if(((std::min)(alz, amz)/4 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/4 < tools::log_max_value<T>()) && (z > a)) |
903 | { |
904 | // compute the 4th root of the result then square it twice: |
905 | T sq = pow(z / agh, a / 4) * exp(amz / 4); |
906 | prefix = sq * sq; |
907 | prefix *= prefix; |
908 | } |
909 | else if((amza > tools::log_min_value<T>()) && (amza < tools::log_max_value<T>())) |
910 | { |
911 | prefix = pow((z * exp(amza)) / agh, a); |
912 | } |
913 | else |
914 | { |
915 | prefix = exp(alz + amz); |
916 | } |
917 | } |
918 | else |
919 | { |
920 | prefix = pow(z / agh, a) * exp(amz); |
921 | } |
922 | } |
923 | prefix *= sqrt(agh / boost::math::constants::e<T>()) / Lanczos::lanczos_sum_expG_scaled(a); |
924 | return prefix; |
925 | } |
926 | // |
927 | // And again, without Lanczos support: |
928 | // |
929 | template <class T, class Policy> |
930 | T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined_lanczos& l) |
931 | { |
932 | BOOST_MATH_STD_USING |
933 | |
934 | if((a < 1) && (z < 1)) |
935 | { |
936 | // No overflow possible since the power terms tend to unity as a,z -> 0 |
937 | return pow(z, a) * exp(-z) / boost::math::tgamma(a, pol); |
938 | } |
939 | else if(a > minimum_argument_for_bernoulli_recursion<T>()) |
940 | { |
941 | T scaled_gamma = scaled_tgamma_no_lanczos(a, pol); |
942 | T power_term = pow(z / a, a / 2); |
943 | T a_minus_z = a - z; |
944 | if ((0 == power_term) || (fabs(a_minus_z) > tools::log_max_value<T>())) |
945 | { |
946 | // The result is probably zero, but we need to be sure: |
947 | return exp(a * log(z / a) + a_minus_z - log(scaled_gamma)); |
948 | } |
949 | return (power_term * exp(a_minus_z)) * (power_term / scaled_gamma); |
950 | } |
951 | else |
952 | { |
953 | // |
954 | // Usual case is to calculate the prefix at a+shift and recurse down |
955 | // to the value we want: |
956 | // |
957 | const int min_z = minimum_argument_for_bernoulli_recursion<T>(); |
958 | long shift = 1 + ltrunc(min_z - a); |
959 | T result = regularised_gamma_prefix(T(a + shift), z, pol, l); |
960 | if (result != 0) |
961 | { |
962 | for (long i = 0; i < shift; ++i) |
963 | { |
964 | result /= z; |
965 | result *= a + i; |
966 | } |
967 | return result; |
968 | } |
969 | else |
970 | { |
971 | // |
972 | // We failed, most probably we have z << 1, try again, this time |
973 | // we calculate z^a e^-z / tgamma(a+shift), combining power terms |
974 | // as we go. And again recurse down to the result. |
975 | // |
976 | T scaled_gamma = scaled_tgamma_no_lanczos(T(a + shift), pol); |
977 | T power_term_1 = pow(z / (a + shift), a); |
978 | T power_term_2 = pow(a + shift, -shift); |
979 | T power_term_3 = exp(a + shift - z); |
980 | if ((0 == power_term_1) || (0 == power_term_2) || (0 == power_term_3) || (fabs(a + shift - z) > tools::log_max_value<T>())) |
981 | { |
982 | // We have no test case that gets here, most likely the type T |
983 | // has a high precision but low exponent range: |
984 | return exp(a * log(z) - z - boost::math::lgamma(a, pol)); |
985 | } |
986 | result = power_term_1 * power_term_2 * power_term_3 / scaled_gamma; |
987 | for (long i = 0; i < shift; ++i) |
988 | { |
989 | result *= a + i; |
990 | } |
991 | return result; |
992 | } |
993 | } |
994 | } |
995 | // |
996 | // Upper gamma fraction for very small a: |
997 | // |
998 | template <class T, class Policy> |
999 | inline T tgamma_small_upper_part(T a, T x, const Policy& pol, T* pgam = 0, bool invert = false, T* pderivative = 0) |
1000 | { |
1001 | BOOST_MATH_STD_USING // ADL of std functions. |
1002 | // |
1003 | // Compute the full upper fraction (Q) when a is very small: |
1004 | // |
1005 | T result; |
1006 | result = boost::math::tgamma1pm1(a, pol); |
1007 | if(pgam) |
1008 | *pgam = (result + 1) / a; |
1009 | T p = boost::math::powm1(x, a, pol); |
1010 | result -= p; |
1011 | result /= a; |
1012 | detail::small_gamma2_series<T> s(a, x); |
1013 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>() - 10; |
1014 | p += 1; |
1015 | if(pderivative) |
1016 | *pderivative = p / (*pgam * exp(x)); |
1017 | T init_value = invert ? *pgam : 0; |
1018 | result = -p * tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, (init_value - result) / p); |
1019 | policies::check_series_iterations<T>("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)" , max_iter, pol); |
1020 | if(invert) |
1021 | result = -result; |
1022 | return result; |
1023 | } |
1024 | // |
1025 | // Upper gamma fraction for integer a: |
1026 | // |
1027 | template <class T, class Policy> |
1028 | inline T finite_gamma_q(T a, T x, Policy const& pol, T* pderivative = 0) |
1029 | { |
1030 | // |
1031 | // Calculates normalised Q when a is an integer: |
1032 | // |
1033 | BOOST_MATH_STD_USING |
1034 | T e = exp(-x); |
1035 | T sum = e; |
1036 | if(sum != 0) |
1037 | { |
1038 | T term = sum; |
1039 | for(unsigned n = 1; n < a; ++n) |
1040 | { |
1041 | term /= n; |
1042 | term *= x; |
1043 | sum += term; |
1044 | } |
1045 | } |
1046 | if(pderivative) |
1047 | { |
1048 | *pderivative = e * pow(x, a) / boost::math::unchecked_factorial<T>(itrunc(T(a - 1), pol)); |
1049 | } |
1050 | return sum; |
1051 | } |
1052 | // |
1053 | // Upper gamma fraction for half integer a: |
1054 | // |
1055 | template <class T, class Policy> |
1056 | T finite_half_gamma_q(T a, T x, T* p_derivative, const Policy& pol) |
1057 | { |
1058 | // |
1059 | // Calculates normalised Q when a is a half-integer: |
1060 | // |
1061 | BOOST_MATH_STD_USING |
1062 | T e = boost::math::erfc(sqrt(x), pol); |
1063 | if((e != 0) && (a > 1)) |
1064 | { |
1065 | T term = exp(-x) / sqrt(constants::pi<T>() * x); |
1066 | term *= x; |
1067 | static const T half = T(1) / 2; |
1068 | term /= half; |
1069 | T sum = term; |
1070 | for(unsigned n = 2; n < a; ++n) |
1071 | { |
1072 | term /= n - half; |
1073 | term *= x; |
1074 | sum += term; |
1075 | } |
1076 | e += sum; |
1077 | if(p_derivative) |
1078 | { |
1079 | *p_derivative = 0; |
1080 | } |
1081 | } |
1082 | else if(p_derivative) |
1083 | { |
1084 | // We'll be dividing by x later, so calculate derivative * x: |
1085 | *p_derivative = sqrt(x) * exp(-x) / constants::root_pi<T>(); |
1086 | } |
1087 | return e; |
1088 | } |
1089 | // |
1090 | // Asymptotic approximation for large argument, see: https://dlmf.nist.gov/8.11#E2 |
1091 | // |
1092 | template <class T> |
1093 | struct incomplete_tgamma_large_x_series |
1094 | { |
1095 | typedef T result_type; |
1096 | incomplete_tgamma_large_x_series(const T& a, const T& x) |
1097 | : a_poch(a - 1), z(x), term(1) {} |
1098 | T operator()() |
1099 | { |
1100 | T result = term; |
1101 | term *= a_poch / z; |
1102 | a_poch -= 1; |
1103 | return result; |
1104 | } |
1105 | T a_poch, z, term; |
1106 | }; |
1107 | |
1108 | template <class T, class Policy> |
1109 | T incomplete_tgamma_large_x(const T& a, const T& x, const Policy& pol) |
1110 | { |
1111 | BOOST_MATH_STD_USING |
1112 | incomplete_tgamma_large_x_series<T> s(a, x); |
1113 | boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); |
1114 | T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); |
1115 | boost::math::policies::check_series_iterations<T>("boost::math::tgamma<%1%>(%1%,%1%)" , max_iter, pol); |
1116 | return result; |
1117 | } |
1118 | |
1119 | |
1120 | // |
1121 | // Main incomplete gamma entry point, handles all four incomplete gamma's: |
1122 | // |
1123 | template <class T, class Policy> |
1124 | T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, |
1125 | const Policy& pol, T* p_derivative) |
1126 | { |
1127 | static const char* function = "boost::math::gamma_p<%1%>(%1%, %1%)" ; |
1128 | if(a <= 0) |
1129 | return policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%)." , a, pol); |
1130 | if(x < 0) |
1131 | return policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%)." , x, pol); |
1132 | |
1133 | BOOST_MATH_STD_USING |
1134 | |
1135 | typedef typename lanczos::lanczos<T, Policy>::type lanczos_type; |
1136 | |
1137 | T result = 0; // Just to avoid warning C4701: potentially uninitialized local variable 'result' used |
1138 | |
1139 | if(a >= max_factorial<T>::value && !normalised) |
1140 | { |
1141 | // |
1142 | // When we're computing the non-normalized incomplete gamma |
1143 | // and a is large the result is rather hard to compute unless |
1144 | // we use logs. There are really two options - if x is a long |
1145 | // way from a in value then we can reliably use methods 2 and 4 |
1146 | // below in logarithmic form and go straight to the result. |
1147 | // Otherwise we let the regularized gamma take the strain |
1148 | // (the result is unlikely to underflow in the central region anyway) |
1149 | // and combine with lgamma in the hopes that we get a finite result. |
1150 | // |
1151 | if(invert && (a * 4 < x)) |
1152 | { |
1153 | // This is method 4 below, done in logs: |
1154 | result = a * log(x) - x; |
1155 | if(p_derivative) |
1156 | *p_derivative = exp(result); |
1157 | result += log(upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>())); |
1158 | } |
1159 | else if(!invert && (a > 4 * x)) |
1160 | { |
1161 | // This is method 2 below, done in logs: |
1162 | result = a * log(x) - x; |
1163 | if(p_derivative) |
1164 | *p_derivative = exp(result); |
1165 | T init_value = 0; |
1166 | result += log(detail::lower_gamma_series(a, x, pol, init_value) / a); |
1167 | } |
1168 | else |
1169 | { |
1170 | result = gamma_incomplete_imp(a, x, true, invert, pol, p_derivative); |
1171 | if(result == 0) |
1172 | { |
1173 | if(invert) |
1174 | { |
1175 | // Try http://functions.wolfram.com/06.06.06.0039.01 |
1176 | result = 1 + 1 / (12 * a) + 1 / (288 * a * a); |
1177 | result = log(result) - a + (a - 0.5f) * log(a) + log(boost::math::constants::root_two_pi<T>()); |
1178 | if(p_derivative) |
1179 | *p_derivative = exp(a * log(x) - x); |
1180 | } |
1181 | else |
1182 | { |
1183 | // This is method 2 below, done in logs, we're really outside the |
1184 | // range of this method, but since the result is almost certainly |
1185 | // infinite, we should probably be OK: |
1186 | result = a * log(x) - x; |
1187 | if(p_derivative) |
1188 | *p_derivative = exp(result); |
1189 | T init_value = 0; |
1190 | result += log(detail::lower_gamma_series(a, x, pol, init_value) / a); |
1191 | } |
1192 | } |
1193 | else |
1194 | { |
1195 | result = log(result) + boost::math::lgamma(a, pol); |
1196 | } |
1197 | } |
1198 | if(result > tools::log_max_value<T>()) |
1199 | return policies::raise_overflow_error<T>(function, 0, pol); |
1200 | return exp(result); |
1201 | } |
1202 | |
1203 | BOOST_ASSERT((p_derivative == 0) || (normalised == true)); |
1204 | |
1205 | bool is_int, is_half_int; |
1206 | bool is_small_a = (a < 30) && (a <= x + 1) && (x < tools::log_max_value<T>()); |
1207 | if(is_small_a) |
1208 | { |
1209 | T fa = floor(a); |
1210 | is_int = (fa == a); |
1211 | is_half_int = is_int ? false : (fabs(fa - a) == 0.5f); |
1212 | } |
1213 | else |
1214 | { |
1215 | is_int = is_half_int = false; |
1216 | } |
1217 | |
1218 | int eval_method; |
1219 | |
1220 | if(is_int && (x > 0.6)) |
1221 | { |
1222 | // calculate Q via finite sum: |
1223 | invert = !invert; |
1224 | eval_method = 0; |
1225 | } |
1226 | else if(is_half_int && (x > 0.2)) |
1227 | { |
1228 | // calculate Q via finite sum for half integer a: |
1229 | invert = !invert; |
1230 | eval_method = 1; |
1231 | } |
1232 | else if((x < tools::root_epsilon<T>()) && (a > 1)) |
1233 | { |
1234 | eval_method = 6; |
1235 | } |
1236 | else if ((x > 1000) && ((a < x) || (fabs(a - 50) / x < 1))) |
1237 | { |
1238 | // calculate Q via asymptotic approximation: |
1239 | invert = !invert; |
1240 | eval_method = 7; |
1241 | } |
1242 | else if(x < 0.5) |
1243 | { |
1244 | // |
1245 | // Changeover criterion chosen to give a changeover at Q ~ 0.33 |
1246 | // |
1247 | if(-0.4 / log(x) < a) |
1248 | { |
1249 | eval_method = 2; |
1250 | } |
1251 | else |
1252 | { |
1253 | eval_method = 3; |
1254 | } |
1255 | } |
1256 | else if(x < 1.1) |
1257 | { |
1258 | // |
1259 | // Changover here occurs when P ~ 0.75 or Q ~ 0.25: |
1260 | // |
1261 | if(x * 0.75f < a) |
1262 | { |
1263 | eval_method = 2; |
1264 | } |
1265 | else |
1266 | { |
1267 | eval_method = 3; |
1268 | } |
1269 | } |
1270 | else |
1271 | { |
1272 | // |
1273 | // Begin by testing whether we're in the "bad" zone |
1274 | // where the result will be near 0.5 and the usual |
1275 | // series and continued fractions are slow to converge: |
1276 | // |
1277 | bool use_temme = false; |
1278 | if(normalised && std::numeric_limits<T>::is_specialized && (a > 20)) |
1279 | { |
1280 | T sigma = fabs((x-a)/a); |
1281 | if((a > 200) && (policies::digits<T, Policy>() <= 113)) |
1282 | { |
1283 | // |
1284 | // This limit is chosen so that we use Temme's expansion |
1285 | // only if the result would be larger than about 10^-6. |
1286 | // Below that the regular series and continued fractions |
1287 | // converge OK, and if we use Temme's method we get increasing |
1288 | // errors from the dominant erfc term as it's (inexact) argument |
1289 | // increases in magnitude. |
1290 | // |
1291 | if(20 / a > sigma * sigma) |
1292 | use_temme = true; |
1293 | } |
1294 | else if(policies::digits<T, Policy>() <= 64) |
1295 | { |
1296 | // Note in this zone we can't use Temme's expansion for |
1297 | // types longer than an 80-bit real: |
1298 | // it would require too many terms in the polynomials. |
1299 | if(sigma < 0.4) |
1300 | use_temme = true; |
1301 | } |
1302 | } |
1303 | if(use_temme) |
1304 | { |
1305 | eval_method = 5; |
1306 | } |
1307 | else |
1308 | { |
1309 | // |
1310 | // Regular case where the result will not be too close to 0.5. |
1311 | // |
1312 | // Changeover here occurs at P ~ Q ~ 0.5 |
1313 | // Note that series computation of P is about x2 faster than continued fraction |
1314 | // calculation of Q, so try and use the CF only when really necessary, especially |
1315 | // for small x. |
1316 | // |
1317 | if(x - (1 / (3 * x)) < a) |
1318 | { |
1319 | eval_method = 2; |
1320 | } |
1321 | else |
1322 | { |
1323 | eval_method = 4; |
1324 | invert = !invert; |
1325 | } |
1326 | } |
1327 | } |
1328 | |
1329 | switch(eval_method) |
1330 | { |
1331 | case 0: |
1332 | { |
1333 | result = finite_gamma_q(a, x, pol, p_derivative); |
1334 | if(normalised == false) |
1335 | result *= boost::math::tgamma(a, pol); |
1336 | break; |
1337 | } |
1338 | case 1: |
1339 | { |
1340 | result = finite_half_gamma_q(a, x, p_derivative, pol); |
1341 | if(normalised == false) |
1342 | result *= boost::math::tgamma(a, pol); |
1343 | if(p_derivative && (*p_derivative == 0)) |
1344 | *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type()); |
1345 | break; |
1346 | } |
1347 | case 2: |
1348 | { |
1349 | // Compute P: |
1350 | result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol); |
1351 | if(p_derivative) |
1352 | *p_derivative = result; |
1353 | if(result != 0) |
1354 | { |
1355 | // |
1356 | // If we're going to be inverting the result then we can |
1357 | // reduce the number of series evaluations by quite |
1358 | // a few iterations if we set an initial value for the |
1359 | // series sum based on what we'll end up subtracting it from |
1360 | // at the end. |
1361 | // Have to be careful though that this optimization doesn't |
1362 | // lead to spurious numeric overflow. Note that the |
1363 | // scary/expensive overflow checks below are more often |
1364 | // than not bypassed in practice for "sensible" input |
1365 | // values: |
1366 | // |
1367 | T init_value = 0; |
1368 | bool optimised_invert = false; |
1369 | if(invert) |
1370 | { |
1371 | init_value = (normalised ? 1 : boost::math::tgamma(a, pol)); |
1372 | if(normalised || (result >= 1) || (tools::max_value<T>() * result > init_value)) |
1373 | { |
1374 | init_value /= result; |
1375 | if(normalised || (a < 1) || (tools::max_value<T>() / a > init_value)) |
1376 | { |
1377 | init_value *= -a; |
1378 | optimised_invert = true; |
1379 | } |
1380 | else |
1381 | init_value = 0; |
1382 | } |
1383 | else |
1384 | init_value = 0; |
1385 | } |
1386 | result *= detail::lower_gamma_series(a, x, pol, init_value) / a; |
1387 | if(optimised_invert) |
1388 | { |
1389 | invert = false; |
1390 | result = -result; |
1391 | } |
1392 | } |
1393 | break; |
1394 | } |
1395 | case 3: |
1396 | { |
1397 | // Compute Q: |
1398 | invert = !invert; |
1399 | T g; |
1400 | result = tgamma_small_upper_part(a, x, pol, &g, invert, p_derivative); |
1401 | invert = false; |
1402 | if(normalised) |
1403 | result /= g; |
1404 | break; |
1405 | } |
1406 | case 4: |
1407 | { |
1408 | // Compute Q: |
1409 | result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol); |
1410 | if(p_derivative) |
1411 | *p_derivative = result; |
1412 | if(result != 0) |
1413 | result *= upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>()); |
1414 | break; |
1415 | } |
1416 | case 5: |
1417 | { |
1418 | // |
1419 | // Use compile time dispatch to the appropriate |
1420 | // Temme asymptotic expansion. This may be dead code |
1421 | // if T does not have numeric limits support, or has |
1422 | // too many digits for the most precise version of |
1423 | // these expansions, in that case we'll be calling |
1424 | // an empty function. |
1425 | // |
1426 | typedef typename policies::precision<T, Policy>::type precision_type; |
1427 | |
1428 | typedef boost::integral_constant<int, |
1429 | precision_type::value <= 0 ? 0 : |
1430 | precision_type::value <= 53 ? 53 : |
1431 | precision_type::value <= 64 ? 64 : |
1432 | precision_type::value <= 113 ? 113 : 0 |
1433 | > tag_type; |
1434 | |
1435 | result = igamma_temme_large(a, x, pol, static_cast<tag_type const*>(0)); |
1436 | if(x >= a) |
1437 | invert = !invert; |
1438 | if(p_derivative) |
1439 | *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type()); |
1440 | break; |
1441 | } |
1442 | case 6: |
1443 | { |
1444 | // x is so small that P is necessarily very small too, |
1445 | // use http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/ |
1446 | result = !normalised ? pow(x, a) / (a) : pow(x, a) / boost::math::tgamma(a + 1, pol); |
1447 | result *= 1 - a * x / (a + 1); |
1448 | if (p_derivative) |
1449 | *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type()); |
1450 | break; |
1451 | } |
1452 | case 7: |
1453 | { |
1454 | // x is large, |
1455 | // Compute Q: |
1456 | result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol); |
1457 | if (p_derivative) |
1458 | *p_derivative = result; |
1459 | result /= x; |
1460 | if (result != 0) |
1461 | result *= incomplete_tgamma_large_x(a, x, pol); |
1462 | break; |
1463 | } |
1464 | } |
1465 | |
1466 | if(normalised && (result > 1)) |
1467 | result = 1; |
1468 | if(invert) |
1469 | { |
1470 | T gam = normalised ? 1 : boost::math::tgamma(a, pol); |
1471 | result = gam - result; |
1472 | } |
1473 | if(p_derivative) |
1474 | { |
1475 | // |
1476 | // Need to convert prefix term to derivative: |
1477 | // |
1478 | if((x < 1) && (tools::max_value<T>() * x < *p_derivative)) |
1479 | { |
1480 | // overflow, just return an arbitrarily large value: |
1481 | *p_derivative = tools::max_value<T>() / 2; |
1482 | } |
1483 | |
1484 | *p_derivative /= x; |
1485 | } |
1486 | |
1487 | return result; |
1488 | } |
1489 | |
1490 | // |
1491 | // Ratios of two gamma functions: |
1492 | // |
1493 | template <class T, class Policy, class Lanczos> |
1494 | T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const Lanczos& l) |
1495 | { |
1496 | BOOST_MATH_STD_USING |
1497 | if(z < tools::epsilon<T>()) |
1498 | { |
1499 | // |
1500 | // We get spurious numeric overflow unless we're very careful, this |
1501 | // can occur either inside Lanczos::lanczos_sum(z) or in the |
1502 | // final combination of terms, to avoid this, split the product up |
1503 | // into 2 (or 3) parts: |
1504 | // |
1505 | // G(z) / G(L) = 1 / (z * G(L)) ; z < eps, L = z + delta = delta |
1506 | // z * G(L) = z * G(lim) * (G(L)/G(lim)) ; lim = largest factorial |
1507 | // |
1508 | if(boost::math::max_factorial<T>::value < delta) |
1509 | { |
1510 | T ratio = tgamma_delta_ratio_imp_lanczos(delta, T(boost::math::max_factorial<T>::value - delta), pol, l); |
1511 | ratio *= z; |
1512 | ratio *= boost::math::unchecked_factorial<T>(boost::math::max_factorial<T>::value - 1); |
1513 | return 1 / ratio; |
1514 | } |
1515 | else |
1516 | { |
1517 | return 1 / (z * boost::math::tgamma(z + delta, pol)); |
1518 | } |
1519 | } |
1520 | T zgh = static_cast<T>(z + Lanczos::g() - constants::half<T>()); |
1521 | T result; |
1522 | if(z + delta == z) |
1523 | { |
1524 | if(fabs(delta) < 10) |
1525 | result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol)); |
1526 | else |
1527 | result = 1; |
1528 | } |
1529 | else |
1530 | { |
1531 | if(fabs(delta) < 10) |
1532 | { |
1533 | result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol)); |
1534 | } |
1535 | else |
1536 | { |
1537 | result = pow(zgh / (zgh + delta), z - constants::half<T>()); |
1538 | } |
1539 | // Split the calculation up to avoid spurious overflow: |
1540 | result *= Lanczos::lanczos_sum(z) / Lanczos::lanczos_sum(T(z + delta)); |
1541 | } |
1542 | result *= pow(constants::e<T>() / (zgh + delta), delta); |
1543 | return result; |
1544 | } |
1545 | // |
1546 | // And again without Lanczos support this time: |
1547 | // |
1548 | template <class T, class Policy> |
1549 | T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const lanczos::undefined_lanczos& l) |
1550 | { |
1551 | BOOST_MATH_STD_USING |
1552 | |
1553 | // |
1554 | // We adjust z and delta so that both z and z+delta are large enough for |
1555 | // Sterling's approximation to hold. We can then calculate the ratio |
1556 | // for the adjusted values, and rescale back down to z and z+delta. |
1557 | // |
1558 | // Get the required shifts first: |
1559 | // |
1560 | long numerator_shift = 0; |
1561 | long denominator_shift = 0; |
1562 | const int min_z = minimum_argument_for_bernoulli_recursion<T>(); |
1563 | |
1564 | if (min_z > z) |
1565 | numerator_shift = 1 + ltrunc(min_z - z); |
1566 | if (min_z > z + delta) |
1567 | denominator_shift = 1 + ltrunc(min_z - z - delta); |
1568 | // |
1569 | // If the shifts are zero, then we can just combine scaled tgamma's |
1570 | // and combine the remaining terms: |
1571 | // |
1572 | if (numerator_shift == 0 && denominator_shift == 0) |
1573 | { |
1574 | T scaled_tgamma_num = scaled_tgamma_no_lanczos(z, pol); |
1575 | T scaled_tgamma_denom = scaled_tgamma_no_lanczos(T(z + delta), pol); |
1576 | T result = scaled_tgamma_num / scaled_tgamma_denom; |
1577 | result *= exp(z * boost::math::log1p(-delta / (z + delta), pol)) * pow((delta + z) / constants::e<T>(), -delta); |
1578 | return result; |
1579 | } |
1580 | // |
1581 | // We're going to have to rescale first, get the adjusted z and delta values, |
1582 | // plus the ratio for the adjusted values: |
1583 | // |
1584 | T zz = z + numerator_shift; |
1585 | T dd = delta - (numerator_shift - denominator_shift); |
1586 | T ratio = tgamma_delta_ratio_imp_lanczos(zz, dd, pol, l); |
1587 | // |
1588 | // Use gamma recurrence relations to get back to the original |
1589 | // z and z+delta: |
1590 | // |
1591 | for (long long i = 0; i < numerator_shift; ++i) |
1592 | { |
1593 | ratio /= (z + i); |
1594 | if (i < denominator_shift) |
1595 | ratio *= (z + delta + i); |
1596 | } |
1597 | for (long long i = numerator_shift; i < denominator_shift; ++i) |
1598 | { |
1599 | ratio *= (z + delta + i); |
1600 | } |
1601 | return ratio; |
1602 | } |
1603 | |
1604 | template <class T, class Policy> |
1605 | T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol) |
1606 | { |
1607 | BOOST_MATH_STD_USING |
1608 | |
1609 | if((z <= 0) || (z + delta <= 0)) |
1610 | { |
1611 | // This isn't very sophisticated, or accurate, but it does work: |
1612 | return boost::math::tgamma(z, pol) / boost::math::tgamma(z + delta, pol); |
1613 | } |
1614 | |
1615 | if(floor(delta) == delta) |
1616 | { |
1617 | if(floor(z) == z) |
1618 | { |
1619 | // |
1620 | // Both z and delta are integers, see if we can just use table lookup |
1621 | // of the factorials to get the result: |
1622 | // |
1623 | if((z <= max_factorial<T>::value) && (z + delta <= max_factorial<T>::value)) |
1624 | { |
1625 | return unchecked_factorial<T>((unsigned)itrunc(z, pol) - 1) / unchecked_factorial<T>((unsigned)itrunc(T(z + delta), pol) - 1); |
1626 | } |
1627 | } |
1628 | if(fabs(delta) < 20) |
1629 | { |
1630 | // |
1631 | // delta is a small integer, we can use a finite product: |
1632 | // |
1633 | if(delta == 0) |
1634 | return 1; |
1635 | if(delta < 0) |
1636 | { |
1637 | z -= 1; |
1638 | T result = z; |
1639 | while(0 != (delta += 1)) |
1640 | { |
1641 | z -= 1; |
1642 | result *= z; |
1643 | } |
1644 | return result; |
1645 | } |
1646 | else |
1647 | { |
1648 | T result = 1 / z; |
1649 | while(0 != (delta -= 1)) |
1650 | { |
1651 | z += 1; |
1652 | result /= z; |
1653 | } |
1654 | return result; |
1655 | } |
1656 | } |
1657 | } |
1658 | typedef typename lanczos::lanczos<T, Policy>::type lanczos_type; |
1659 | return tgamma_delta_ratio_imp_lanczos(z, delta, pol, lanczos_type()); |
1660 | } |
1661 | |
1662 | template <class T, class Policy> |
1663 | T tgamma_ratio_imp(T x, T y, const Policy& pol) |
1664 | { |
1665 | BOOST_MATH_STD_USING |
1666 | |
1667 | if((x <= 0) || (boost::math::isinf)(x)) |
1668 | return policies::raise_domain_error<T>("boost::math::tgamma_ratio<%1%>(%1%, %1%)" , "Gamma function ratios only implemented for positive arguments (got a=%1%)." , x, pol); |
1669 | if((y <= 0) || (boost::math::isinf)(y)) |
1670 | return policies::raise_domain_error<T>("boost::math::tgamma_ratio<%1%>(%1%, %1%)" , "Gamma function ratios only implemented for positive arguments (got b=%1%)." , y, pol); |
1671 | |
1672 | if(x <= tools::min_value<T>()) |
1673 | { |
1674 | // Special case for denorms...Ugh. |
1675 | T shift = ldexp(T(1), tools::digits<T>()); |
1676 | return shift * tgamma_ratio_imp(T(x * shift), y, pol); |
1677 | } |
1678 | |
1679 | if((x < max_factorial<T>::value) && (y < max_factorial<T>::value)) |
1680 | { |
1681 | // Rather than subtracting values, lets just call the gamma functions directly: |
1682 | return boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol); |
1683 | } |
1684 | T prefix = 1; |
1685 | if(x < 1) |
1686 | { |
1687 | if(y < 2 * max_factorial<T>::value) |
1688 | { |
1689 | // We need to sidestep on x as well, otherwise we'll underflow |
1690 | // before we get to factor in the prefix term: |
1691 | prefix /= x; |
1692 | x += 1; |
1693 | while(y >= max_factorial<T>::value) |
1694 | { |
1695 | y -= 1; |
1696 | prefix /= y; |
1697 | } |
1698 | return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol); |
1699 | } |
1700 | // |
1701 | // result is almost certainly going to underflow to zero, try logs just in case: |
1702 | // |
1703 | return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol)); |
1704 | } |
1705 | if(y < 1) |
1706 | { |
1707 | if(x < 2 * max_factorial<T>::value) |
1708 | { |
1709 | // We need to sidestep on y as well, otherwise we'll overflow |
1710 | // before we get to factor in the prefix term: |
1711 | prefix *= y; |
1712 | y += 1; |
1713 | while(x >= max_factorial<T>::value) |
1714 | { |
1715 | x -= 1; |
1716 | prefix *= x; |
1717 | } |
1718 | return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol); |
1719 | } |
1720 | // |
1721 | // Result will almost certainly overflow, try logs just in case: |
1722 | // |
1723 | return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol)); |
1724 | } |
1725 | // |
1726 | // Regular case, x and y both large and similar in magnitude: |
1727 | // |
1728 | return boost::math::tgamma_delta_ratio(x, y - x, pol); |
1729 | } |
1730 | |
1731 | template <class T, class Policy> |
1732 | T gamma_p_derivative_imp(T a, T x, const Policy& pol) |
1733 | { |
1734 | BOOST_MATH_STD_USING |
1735 | // |
1736 | // Usual error checks first: |
1737 | // |
1738 | if(a <= 0) |
1739 | return policies::raise_domain_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)" , "Argument a to the incomplete gamma function must be greater than zero (got a=%1%)." , a, pol); |
1740 | if(x < 0) |
1741 | return policies::raise_domain_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)" , "Argument x to the incomplete gamma function must be >= 0 (got x=%1%)." , x, pol); |
1742 | // |
1743 | // Now special cases: |
1744 | // |
1745 | if(x == 0) |
1746 | { |
1747 | return (a > 1) ? 0 : |
1748 | (a == 1) ? 1 : policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)" , 0, pol); |
1749 | } |
1750 | // |
1751 | // Normal case: |
1752 | // |
1753 | typedef typename lanczos::lanczos<T, Policy>::type lanczos_type; |
1754 | T f1 = detail::regularised_gamma_prefix(a, x, pol, lanczos_type()); |
1755 | if((x < 1) && (tools::max_value<T>() * x < f1)) |
1756 | { |
1757 | // overflow: |
1758 | return policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)" , 0, pol); |
1759 | } |
1760 | if(f1 == 0) |
1761 | { |
1762 | // Underflow in calculation, use logs instead: |
1763 | f1 = a * log(x) - x - lgamma(a, pol) - log(x); |
1764 | f1 = exp(f1); |
1765 | } |
1766 | else |
1767 | f1 /= x; |
1768 | |
1769 | return f1; |
1770 | } |
1771 | |
1772 | template <class T, class Policy> |
1773 | inline typename tools::promote_args<T>::type |
1774 | tgamma(T z, const Policy& /* pol */, const boost::true_type) |
1775 | { |
1776 | BOOST_FPU_EXCEPTION_GUARD |
1777 | typedef typename tools::promote_args<T>::type result_type; |
1778 | typedef typename policies::evaluation<result_type, Policy>::type value_type; |
1779 | typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type; |
1780 | typedef typename policies::normalise< |
1781 | Policy, |
1782 | policies::promote_float<false>, |
1783 | policies::promote_double<false>, |
1784 | policies::discrete_quantile<>, |
1785 | policies::assert_undefined<> >::type forwarding_policy; |
1786 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::gamma_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type()), "boost::math::tgamma<%1%>(%1%)" ); |
1787 | } |
1788 | |
1789 | template <class T, class Policy> |
1790 | struct igamma_initializer |
1791 | { |
1792 | struct init |
1793 | { |
1794 | init() |
1795 | { |
1796 | typedef typename policies::precision<T, Policy>::type precision_type; |
1797 | |
1798 | typedef boost::integral_constant<int, |
1799 | precision_type::value <= 0 ? 0 : |
1800 | precision_type::value <= 53 ? 53 : |
1801 | precision_type::value <= 64 ? 64 : |
1802 | precision_type::value <= 113 ? 113 : 0 |
1803 | > tag_type; |
1804 | |
1805 | do_init(tag_type()); |
1806 | } |
1807 | template <int N> |
1808 | static void do_init(const boost::integral_constant<int, N>&) |
1809 | { |
1810 | // If std::numeric_limits<T>::digits is zero, we must not call |
1811 | // our initialization code here as the precision presumably |
1812 | // varies at runtime, and will not have been set yet. Plus the |
1813 | // code requiring initialization isn't called when digits == 0. |
1814 | if(std::numeric_limits<T>::digits) |
1815 | { |
1816 | boost::math::gamma_p(static_cast<T>(400), static_cast<T>(400), Policy()); |
1817 | } |
1818 | } |
1819 | static void do_init(const boost::integral_constant<int, 53>&){} |
1820 | void force_instantiate()const{} |
1821 | }; |
1822 | static const init initializer; |
1823 | static void force_instantiate() |
1824 | { |
1825 | initializer.force_instantiate(); |
1826 | } |
1827 | }; |
1828 | |
1829 | template <class T, class Policy> |
1830 | const typename igamma_initializer<T, Policy>::init igamma_initializer<T, Policy>::initializer; |
1831 | |
1832 | template <class T, class Policy> |
1833 | struct lgamma_initializer |
1834 | { |
1835 | struct init |
1836 | { |
1837 | init() |
1838 | { |
1839 | typedef typename policies::precision<T, Policy>::type precision_type; |
1840 | typedef boost::integral_constant<int, |
1841 | precision_type::value <= 0 ? 0 : |
1842 | precision_type::value <= 64 ? 64 : |
1843 | precision_type::value <= 113 ? 113 : 0 |
1844 | > tag_type; |
1845 | |
1846 | do_init(tag_type()); |
1847 | } |
1848 | static void do_init(const boost::integral_constant<int, 64>&) |
1849 | { |
1850 | boost::math::lgamma(static_cast<T>(2.5), Policy()); |
1851 | boost::math::lgamma(static_cast<T>(1.25), Policy()); |
1852 | boost::math::lgamma(static_cast<T>(1.75), Policy()); |
1853 | } |
1854 | static void do_init(const boost::integral_constant<int, 113>&) |
1855 | { |
1856 | boost::math::lgamma(static_cast<T>(2.5), Policy()); |
1857 | boost::math::lgamma(static_cast<T>(1.25), Policy()); |
1858 | boost::math::lgamma(static_cast<T>(1.5), Policy()); |
1859 | boost::math::lgamma(static_cast<T>(1.75), Policy()); |
1860 | } |
1861 | static void do_init(const boost::integral_constant<int, 0>&) |
1862 | { |
1863 | } |
1864 | void force_instantiate()const{} |
1865 | }; |
1866 | static const init initializer; |
1867 | static void force_instantiate() |
1868 | { |
1869 | initializer.force_instantiate(); |
1870 | } |
1871 | }; |
1872 | |
1873 | template <class T, class Policy> |
1874 | const typename lgamma_initializer<T, Policy>::init lgamma_initializer<T, Policy>::initializer; |
1875 | |
1876 | template <class T1, class T2, class Policy> |
1877 | inline typename tools::promote_args<T1, T2>::type |
1878 | tgamma(T1 a, T2 z, const Policy&, const boost::false_type) |
1879 | { |
1880 | BOOST_FPU_EXCEPTION_GUARD |
1881 | typedef typename tools::promote_args<T1, T2>::type result_type; |
1882 | typedef typename policies::evaluation<result_type, Policy>::type value_type; |
1883 | // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type; |
1884 | typedef typename policies::normalise< |
1885 | Policy, |
1886 | policies::promote_float<false>, |
1887 | policies::promote_double<false>, |
1888 | policies::discrete_quantile<>, |
1889 | policies::assert_undefined<> >::type forwarding_policy; |
1890 | |
1891 | igamma_initializer<value_type, forwarding_policy>::force_instantiate(); |
1892 | |
1893 | return policies::checked_narrowing_cast<result_type, forwarding_policy>( |
1894 | detail::gamma_incomplete_imp(static_cast<value_type>(a), |
1895 | static_cast<value_type>(z), false, true, |
1896 | forwarding_policy(), static_cast<value_type*>(0)), "boost::math::tgamma<%1%>(%1%, %1%)" ); |
1897 | } |
1898 | |
1899 | template <class T1, class T2> |
1900 | inline typename tools::promote_args<T1, T2>::type |
1901 | tgamma(T1 a, T2 z, const boost::false_type& tag) |
1902 | { |
1903 | return tgamma(a, z, policies::policy<>(), tag); |
1904 | } |
1905 | |
1906 | |
1907 | } // namespace detail |
1908 | |
1909 | template <class T> |
1910 | inline typename tools::promote_args<T>::type |
1911 | tgamma(T z) |
1912 | { |
1913 | return tgamma(z, policies::policy<>()); |
1914 | } |
1915 | |
1916 | template <class T, class Policy> |
1917 | inline typename tools::promote_args<T>::type |
1918 | lgamma(T z, int* sign, const Policy&) |
1919 | { |
1920 | BOOST_FPU_EXCEPTION_GUARD |
1921 | typedef typename tools::promote_args<T>::type result_type; |
1922 | typedef typename policies::evaluation<result_type, Policy>::type value_type; |
1923 | typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type; |
1924 | typedef typename policies::normalise< |
1925 | Policy, |
1926 | policies::promote_float<false>, |
1927 | policies::promote_double<false>, |
1928 | policies::discrete_quantile<>, |
1929 | policies::assert_undefined<> >::type forwarding_policy; |
1930 | |
1931 | detail::lgamma_initializer<value_type, forwarding_policy>::force_instantiate(); |
1932 | |
1933 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::lgamma_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type(), sign), "boost::math::lgamma<%1%>(%1%)" ); |
1934 | } |
1935 | |
1936 | template <class T> |
1937 | inline typename tools::promote_args<T>::type |
1938 | lgamma(T z, int* sign) |
1939 | { |
1940 | return lgamma(z, sign, policies::policy<>()); |
1941 | } |
1942 | |
1943 | template <class T, class Policy> |
1944 | inline typename tools::promote_args<T>::type |
1945 | lgamma(T x, const Policy& pol) |
1946 | { |
1947 | return ::boost::math::lgamma(x, 0, pol); |
1948 | } |
1949 | |
1950 | template <class T> |
1951 | inline typename tools::promote_args<T>::type |
1952 | lgamma(T x) |
1953 | { |
1954 | return ::boost::math::lgamma(x, 0, policies::policy<>()); |
1955 | } |
1956 | |
1957 | template <class T, class Policy> |
1958 | inline typename tools::promote_args<T>::type |
1959 | tgamma1pm1(T z, const Policy& /* pol */) |
1960 | { |
1961 | BOOST_FPU_EXCEPTION_GUARD |
1962 | typedef typename tools::promote_args<T>::type result_type; |
1963 | typedef typename policies::evaluation<result_type, Policy>::type value_type; |
1964 | typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type; |
1965 | typedef typename policies::normalise< |
1966 | Policy, |
1967 | policies::promote_float<false>, |
1968 | policies::promote_double<false>, |
1969 | policies::discrete_quantile<>, |
1970 | policies::assert_undefined<> >::type forwarding_policy; |
1971 | |
1972 | return policies::checked_narrowing_cast<typename remove_cv<result_type>::type, forwarding_policy>(detail::tgammap1m1_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type()), "boost::math::tgamma1pm1<%!%>(%1%)" ); |
1973 | } |
1974 | |
1975 | template <class T> |
1976 | inline typename tools::promote_args<T>::type |
1977 | tgamma1pm1(T z) |
1978 | { |
1979 | return tgamma1pm1(z, policies::policy<>()); |
1980 | } |
1981 | |
1982 | // |
1983 | // Full upper incomplete gamma: |
1984 | // |
1985 | template <class T1, class T2> |
1986 | inline typename tools::promote_args<T1, T2>::type |
1987 | tgamma(T1 a, T2 z) |
1988 | { |
1989 | // |
1990 | // Type T2 could be a policy object, or a value, select the |
1991 | // right overload based on T2: |
1992 | // |
1993 | typedef typename policies::is_policy<T2>::type maybe_policy; |
1994 | return detail::tgamma(a, z, maybe_policy()); |
1995 | } |
1996 | template <class T1, class T2, class Policy> |
1997 | inline typename tools::promote_args<T1, T2>::type |
1998 | tgamma(T1 a, T2 z, const Policy& pol) |
1999 | { |
2000 | return detail::tgamma(a, z, pol, boost::false_type()); |
2001 | } |
2002 | // |
2003 | // Full lower incomplete gamma: |
2004 | // |
2005 | template <class T1, class T2, class Policy> |
2006 | inline typename tools::promote_args<T1, T2>::type |
2007 | tgamma_lower(T1 a, T2 z, const Policy&) |
2008 | { |
2009 | BOOST_FPU_EXCEPTION_GUARD |
2010 | typedef typename tools::promote_args<T1, T2>::type result_type; |
2011 | typedef typename policies::evaluation<result_type, Policy>::type value_type; |
2012 | // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type; |
2013 | typedef typename policies::normalise< |
2014 | Policy, |
2015 | policies::promote_float<false>, |
2016 | policies::promote_double<false>, |
2017 | policies::discrete_quantile<>, |
2018 | policies::assert_undefined<> >::type forwarding_policy; |
2019 | |
2020 | detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate(); |
2021 | |
2022 | return policies::checked_narrowing_cast<result_type, forwarding_policy>( |
2023 | detail::gamma_incomplete_imp(static_cast<value_type>(a), |
2024 | static_cast<value_type>(z), false, false, |
2025 | forwarding_policy(), static_cast<value_type*>(0)), "tgamma_lower<%1%>(%1%, %1%)" ); |
2026 | } |
2027 | template <class T1, class T2> |
2028 | inline typename tools::promote_args<T1, T2>::type |
2029 | tgamma_lower(T1 a, T2 z) |
2030 | { |
2031 | return tgamma_lower(a, z, policies::policy<>()); |
2032 | } |
2033 | // |
2034 | // Regularised upper incomplete gamma: |
2035 | // |
2036 | template <class T1, class T2, class Policy> |
2037 | inline typename tools::promote_args<T1, T2>::type |
2038 | gamma_q(T1 a, T2 z, const Policy& /* pol */) |
2039 | { |
2040 | BOOST_FPU_EXCEPTION_GUARD |
2041 | typedef typename tools::promote_args<T1, T2>::type result_type; |
2042 | typedef typename policies::evaluation<result_type, Policy>::type value_type; |
2043 | // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type; |
2044 | typedef typename policies::normalise< |
2045 | Policy, |
2046 | policies::promote_float<false>, |
2047 | policies::promote_double<false>, |
2048 | policies::discrete_quantile<>, |
2049 | policies::assert_undefined<> >::type forwarding_policy; |
2050 | |
2051 | detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate(); |
2052 | |
2053 | return policies::checked_narrowing_cast<result_type, forwarding_policy>( |
2054 | detail::gamma_incomplete_imp(static_cast<value_type>(a), |
2055 | static_cast<value_type>(z), true, true, |
2056 | forwarding_policy(), static_cast<value_type*>(0)), "gamma_q<%1%>(%1%, %1%)" ); |
2057 | } |
2058 | template <class T1, class T2> |
2059 | inline typename tools::promote_args<T1, T2>::type |
2060 | gamma_q(T1 a, T2 z) |
2061 | { |
2062 | return gamma_q(a, z, policies::policy<>()); |
2063 | } |
2064 | // |
2065 | // Regularised lower incomplete gamma: |
2066 | // |
2067 | template <class T1, class T2, class Policy> |
2068 | inline typename tools::promote_args<T1, T2>::type |
2069 | gamma_p(T1 a, T2 z, const Policy&) |
2070 | { |
2071 | BOOST_FPU_EXCEPTION_GUARD |
2072 | typedef typename tools::promote_args<T1, T2>::type result_type; |
2073 | typedef typename policies::evaluation<result_type, Policy>::type value_type; |
2074 | // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type; |
2075 | typedef typename policies::normalise< |
2076 | Policy, |
2077 | policies::promote_float<false>, |
2078 | policies::promote_double<false>, |
2079 | policies::discrete_quantile<>, |
2080 | policies::assert_undefined<> >::type forwarding_policy; |
2081 | |
2082 | detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate(); |
2083 | |
2084 | return policies::checked_narrowing_cast<result_type, forwarding_policy>( |
2085 | detail::gamma_incomplete_imp(static_cast<value_type>(a), |
2086 | static_cast<value_type>(z), true, false, |
2087 | forwarding_policy(), static_cast<value_type*>(0)), "gamma_p<%1%>(%1%, %1%)" ); |
2088 | } |
2089 | template <class T1, class T2> |
2090 | inline typename tools::promote_args<T1, T2>::type |
2091 | gamma_p(T1 a, T2 z) |
2092 | { |
2093 | return gamma_p(a, z, policies::policy<>()); |
2094 | } |
2095 | |
2096 | // ratios of gamma functions: |
2097 | template <class T1, class T2, class Policy> |
2098 | inline typename tools::promote_args<T1, T2>::type |
2099 | tgamma_delta_ratio(T1 z, T2 delta, const Policy& /* pol */) |
2100 | { |
2101 | BOOST_FPU_EXCEPTION_GUARD |
2102 | typedef typename tools::promote_args<T1, T2>::type result_type; |
2103 | typedef typename policies::evaluation<result_type, Policy>::type value_type; |
2104 | typedef typename policies::normalise< |
2105 | Policy, |
2106 | policies::promote_float<false>, |
2107 | policies::promote_double<false>, |
2108 | policies::discrete_quantile<>, |
2109 | policies::assert_undefined<> >::type forwarding_policy; |
2110 | |
2111 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_delta_ratio_imp(static_cast<value_type>(z), static_cast<value_type>(delta), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)" ); |
2112 | } |
2113 | template <class T1, class T2> |
2114 | inline typename tools::promote_args<T1, T2>::type |
2115 | tgamma_delta_ratio(T1 z, T2 delta) |
2116 | { |
2117 | return tgamma_delta_ratio(z, delta, policies::policy<>()); |
2118 | } |
2119 | template <class T1, class T2, class Policy> |
2120 | inline typename tools::promote_args<T1, T2>::type |
2121 | tgamma_ratio(T1 a, T2 b, const Policy&) |
2122 | { |
2123 | typedef typename tools::promote_args<T1, T2>::type result_type; |
2124 | typedef typename policies::evaluation<result_type, Policy>::type value_type; |
2125 | typedef typename policies::normalise< |
2126 | Policy, |
2127 | policies::promote_float<false>, |
2128 | policies::promote_double<false>, |
2129 | policies::discrete_quantile<>, |
2130 | policies::assert_undefined<> >::type forwarding_policy; |
2131 | |
2132 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(b), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)" ); |
2133 | } |
2134 | template <class T1, class T2> |
2135 | inline typename tools::promote_args<T1, T2>::type |
2136 | tgamma_ratio(T1 a, T2 b) |
2137 | { |
2138 | return tgamma_ratio(a, b, policies::policy<>()); |
2139 | } |
2140 | |
2141 | template <class T1, class T2, class Policy> |
2142 | inline typename tools::promote_args<T1, T2>::type |
2143 | gamma_p_derivative(T1 a, T2 x, const Policy&) |
2144 | { |
2145 | BOOST_FPU_EXCEPTION_GUARD |
2146 | typedef typename tools::promote_args<T1, T2>::type result_type; |
2147 | typedef typename policies::evaluation<result_type, Policy>::type value_type; |
2148 | typedef typename policies::normalise< |
2149 | Policy, |
2150 | policies::promote_float<false>, |
2151 | policies::promote_double<false>, |
2152 | policies::discrete_quantile<>, |
2153 | policies::assert_undefined<> >::type forwarding_policy; |
2154 | |
2155 | return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::gamma_p_derivative_imp(static_cast<value_type>(a), static_cast<value_type>(x), forwarding_policy()), "boost::math::gamma_p_derivative<%1%>(%1%, %1%)" ); |
2156 | } |
2157 | template <class T1, class T2> |
2158 | inline typename tools::promote_args<T1, T2>::type |
2159 | gamma_p_derivative(T1 a, T2 x) |
2160 | { |
2161 | return gamma_p_derivative(a, x, policies::policy<>()); |
2162 | } |
2163 | |
2164 | } // namespace math |
2165 | } // namespace boost |
2166 | |
2167 | #ifdef BOOST_MSVC |
2168 | # pragma warning(pop) |
2169 | #endif |
2170 | |
2171 | #include <boost/math/special_functions/detail/igamma_inverse.hpp> |
2172 | #include <boost/math/special_functions/detail/gamma_inva.hpp> |
2173 | #include <boost/math/special_functions/erf.hpp> |
2174 | |
2175 | #endif // BOOST_MATH_SF_GAMMA_HPP |
2176 | |
2177 | |
2178 | |
2179 | |
2180 | |