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

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