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
57namespace boost{ namespace math{
58
59namespace detail{
60
61template <class T>
62inline bool is_odd(T v, const boost::true_type&)
63{
64 int i = static_cast<int>(v);
65 return i&1;
66}
67template <class T>
68inline 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}
75template <class T>
76inline bool is_odd(T v)
77{
78 return is_odd(v, ::boost::is_convertible<T, int>());
79}
80
81template <class T>
82T 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//
113template <class T, class Policy, class Lanczos>
114T 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//
204template <class T, class Policy, class Lanczos>
205T 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//
288template <class T>
289struct upper_incomplete_gamma_fract
290{
291private:
292 T z, a;
293 int k;
294public:
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
310template <class T>
311inline 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
320template <class T>
321struct lower_incomplete_gamma_series
322{
323private:
324 T a, z, result;
325public:
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
338template <class T, class Policy>
339inline 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//
356template<class T>
357std::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
367template<class T>
368int 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
379template <class T, class Policy>
380T 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.
436template <class T, class Policy>
437T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign = 0);
438
439template <class T, class Policy>
440T 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
560template <class T, class Policy>
561inline 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
589template <class T, class Policy>
590T 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//
703template <class T, class Policy, class Lanczos>
704T 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
752template <class T, class Policy>
753inline 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//
768template <class T>
769struct 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
784private:
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//
792template <class T, class Policy>
793T 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//
845template <class T, class Policy, class Lanczos>
846T 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//
929template <class T, class Policy>
930T 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//
998template <class T, class Policy>
999inline 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//
1027template <class T, class Policy>
1028inline 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//
1055template <class T, class Policy>
1056T 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//
1092template <class T>
1093struct 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
1108template <class T, class Policy>
1109T 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//
1123template <class T, class Policy>
1124T 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//
1493template <class T, class Policy, class Lanczos>
1494T 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//
1548template <class T, class Policy>
1549T 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
1604template <class T, class Policy>
1605T 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
1662template <class T, class Policy>
1663T 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
1731template <class T, class Policy>
1732T 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
1772template <class T, class Policy>
1773inline 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
1789template <class T, class Policy>
1790struct 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
1829template <class T, class Policy>
1830const typename igamma_initializer<T, Policy>::init igamma_initializer<T, Policy>::initializer;
1831
1832template <class T, class Policy>
1833struct 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
1873template <class T, class Policy>
1874const typename lgamma_initializer<T, Policy>::init lgamma_initializer<T, Policy>::initializer;
1875
1876template <class T1, class T2, class Policy>
1877inline 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
1899template <class T1, class T2>
1900inline 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
1909template <class T>
1910inline typename tools::promote_args<T>::type
1911 tgamma(T z)
1912{
1913 return tgamma(z, policies::policy<>());
1914}
1915
1916template <class T, class Policy>
1917inline 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
1936template <class T>
1937inline typename tools::promote_args<T>::type
1938 lgamma(T z, int* sign)
1939{
1940 return lgamma(z, sign, policies::policy<>());
1941}
1942
1943template <class T, class Policy>
1944inline typename tools::promote_args<T>::type
1945 lgamma(T x, const Policy& pol)
1946{
1947 return ::boost::math::lgamma(x, 0, pol);
1948}
1949
1950template <class T>
1951inline typename tools::promote_args<T>::type
1952 lgamma(T x)
1953{
1954 return ::boost::math::lgamma(x, 0, policies::policy<>());
1955}
1956
1957template <class T, class Policy>
1958inline 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
1975template <class T>
1976inline 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//
1985template <class T1, class T2>
1986inline 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}
1996template <class T1, class T2, class Policy>
1997inline 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//
2005template <class T1, class T2, class Policy>
2006inline 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}
2027template <class T1, class T2>
2028inline 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//
2036template <class T1, class T2, class Policy>
2037inline 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}
2058template <class T1, class T2>
2059inline 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//
2067template <class T1, class T2, class Policy>
2068inline 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}
2089template <class T1, class T2>
2090inline 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:
2097template <class T1, class T2, class Policy>
2098inline 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}
2113template <class T1, class T2>
2114inline 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}
2119template <class T1, class T2, class Policy>
2120inline 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}
2134template <class T1, class T2>
2135inline 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
2141template <class T1, class T2, class Policy>
2142inline 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}
2157template <class T1, class T2>
2158inline 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

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