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
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)
10 #ifndef BOOST_MATH_SF_GAMMA_HPP
11 #define BOOST_MATH_SF_GAMMA_HPP
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>
40 #include <type_traits>
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 // Several variables made comments,
48 // but some difficulty as whether referenced on not may depend on macro values.
49 // So to be safe, 4100 warnings suppressed.
50 // TODO - revisit this?
53 namespace boost{ namespace math{
58 inline bool is_odd(T v, const std::true_type&)
60 int i = static_cast<int>(v);
64 inline bool is_odd(T v, const std::false_type&)
66 // Oh dear can't cast T to int!
68 T modulus = v - 2 * floor(v/2);
69 return static_cast<bool>(modulus != 0);
72 inline bool is_odd(T v)
74 return is_odd(v, ::std::is_convertible<T, int>());
80 // Ad hoc function calculates x * sin(pi * x),
81 // taking extra care near when x is near a whole number.
100 BOOST_MATH_ASSERT(fl >= 0);
103 T result = sin(dist*boost::math::constants::pi<T>());
104 return sign*z*result;
105 } // template <class T> T sinpx(T z)
107 // tgamma(z), with Lanczos support:
109 template <class T, class Policy, class Lanczos>
110 T gamma_imp(T z, const Policy& pol, const Lanczos& l)
116 #ifdef BOOST_MATH_INSTRUMENT
117 static bool b = false;
120 std::cout << "tgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;
124 static const char* function = "boost::math::tgamma<%1%>(%1%)";
129 return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
132 result = gamma_imp(T(-z), pol, l) * sinpx(z);
133 BOOST_MATH_INSTRUMENT_VARIABLE(result);
134 if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
135 return -boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
136 result = -boost::math::constants::pi<T>() / result;
138 return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
139 if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL)
140 return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", result, pol);
141 BOOST_MATH_INSTRUMENT_VARIABLE(result);
152 BOOST_MATH_INSTRUMENT_VARIABLE(result);
153 if((floor(z) == z) && (z < max_factorial<T>::value))
155 result *= unchecked_factorial<T>(itrunc(z, pol) - 1);
156 BOOST_MATH_INSTRUMENT_VARIABLE(result);
158 else if (z < tools::root_epsilon<T>())
160 if (z < 1 / tools::max_value<T>())
161 result = policies::raise_overflow_error<T>(function, 0, pol);
162 result *= 1 / z - constants::euler<T>();
166 result *= Lanczos::lanczos_sum(z);
167 T zgh = (z + static_cast<T>(Lanczos::g()) - boost::math::constants::half<T>());
169 BOOST_MATH_INSTRUMENT_VARIABLE(result);
170 BOOST_MATH_INSTRUMENT_VARIABLE(tools::log_max_value<T>());
171 if(z * lzgh > tools::log_max_value<T>())
173 // we're going to overflow unless this is done with care:
174 BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
175 if(lzgh * z / 2 > tools::log_max_value<T>())
176 return boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
177 T hp = pow(zgh, (z / 2) - T(0.25));
178 BOOST_MATH_INSTRUMENT_VARIABLE(hp);
179 result *= hp / exp(zgh);
180 BOOST_MATH_INSTRUMENT_VARIABLE(result);
181 if(tools::max_value<T>() / hp < result)
182 return boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
184 BOOST_MATH_INSTRUMENT_VARIABLE(result);
188 BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
189 BOOST_MATH_INSTRUMENT_VARIABLE(pow(zgh, z - boost::math::constants::half<T>()));
190 BOOST_MATH_INSTRUMENT_VARIABLE(exp(zgh));
191 result *= pow(zgh, z - boost::math::constants::half<T>()) / exp(zgh);
192 BOOST_MATH_INSTRUMENT_VARIABLE(result);
198 // lgamma(z) with Lanczos support:
200 template <class T, class Policy, class Lanczos>
201 T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = 0)
203 #ifdef BOOST_MATH_INSTRUMENT
204 static bool b = false;
207 std::cout << "lgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;
214 static const char* function = "boost::math::lgamma<%1%>(%1%)";
218 if(z <= -tools::root_epsilon<T>())
220 // reflection formula:
222 return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%.", z, pol);
234 result = log(boost::math::constants::pi<T>()) - lgamma_imp(z, pol, l) - log(t);
236 else if (z < tools::root_epsilon<T>())
239 return policies::raise_pole_error<T>(function, "Evaluation of lgamma at %1%.", z, pol);
240 if (4 * fabs(z) < tools::epsilon<T>())
241 result = -log(fabs(z));
243 result = log(fabs(1 / z - constants::euler<T>()));
249 typedef typename policies::precision<T, Policy>::type precision_type;
250 typedef std::integral_constant<int,
251 precision_type::value <= 0 ? 0 :
252 precision_type::value <= 64 ? 64 :
253 precision_type::value <= 113 ? 113 : 0
256 result = lgamma_small_imp<T>(z, T(z - 1), T(z - 2), tag_type(), pol, l);
258 else if((z >= 3) && (z < 100) && (std::numeric_limits<T>::max_exponent >= 1024))
260 // taking the log of tgamma reduces the error, no danger of overflow here:
261 result = log(gamma_imp(z, pol, l));
265 // regular evaluation:
266 T zgh = static_cast<T>(z + Lanczos::g() - boost::math::constants::half<T>());
267 result = log(zgh) - 1;
270 // Only add on the lanczos sum part if we're going to need it:
272 if(result * tools::epsilon<T>() < 20)
273 result += log(Lanczos::lanczos_sum_expG_scaled(z));
282 // Incomplete gamma functions follow:
285 struct upper_incomplete_gamma_fract
291 typedef std::pair<T,T> result_type;
293 upper_incomplete_gamma_fract(T a1, T z1)
294 : z(z1-a1+1), a(a1), k(0)
298 result_type operator()()
302 return result_type(k * (a - k), z);
307 inline T upper_gamma_fraction(T a, T z, T eps)
309 // Multiply result by z^a * e^-z to get the full
310 // upper incomplete integral. Divide by tgamma(z)
312 upper_incomplete_gamma_fract<T> f(a, z);
313 return 1 / (z - a + 1 + boost::math::tools::continued_fraction_a(f, eps));
317 struct lower_incomplete_gamma_series
322 typedef T result_type;
323 lower_incomplete_gamma_series(T a1, T z1) : a(a1), z(z1), result(1){}
334 template <class T, class Policy>
335 inline T lower_gamma_series(T a, T z, const Policy& pol, T init_value = 0)
337 // Multiply result by ((z^a) * (e^-z) / a) to get the full
338 // lower incomplete integral. Then divide by tgamma(a)
339 // to get the normalised value.
340 lower_incomplete_gamma_series<T> s(a, z);
341 std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
342 T factor = policies::get_epsilon<T, Policy>();
343 T result = boost::math::tools::sum_series(s, factor, max_iter, init_value);
344 policies::check_series_iterations<T>("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol);
349 // Fully generic tgamma and lgamma use Stirling's approximation
350 // with Bernoulli numbers.
353 std::size_t highest_bernoulli_index()
355 const float digits10_of_type = (std::numeric_limits<T>::is_specialized
356 ? static_cast<float>(std::numeric_limits<T>::digits10)
357 : static_cast<float>(boost::math::tools::digits<T>() * 0.301F));
359 // Find the high index n for Bn to produce the desired precision in Stirling's calculation.
360 return static_cast<std::size_t>(18.0F + (0.6F * digits10_of_type));
364 int minimum_argument_for_bernoulli_recursion()
368 const float digits10_of_type = (std::numeric_limits<T>::is_specialized
369 ? (float) std::numeric_limits<T>::digits10
370 : (float) (boost::math::tools::digits<T>() * 0.301F));
372 int min_arg = (int) (digits10_of_type * 1.7F);
374 if(digits10_of_type < 50.0F)
376 // The following code sequence has been modified
377 // within the context of issue 396.
379 // The calculation of the test-variable limit has now
380 // been protected against overflow/underflow dangers.
382 // The previous line looked like this and did, in fact,
383 // underflow ldexp when using certain multiprecision types.
385 // const float limit = std::ceil(std::pow(1.0f / std::ldexp(1.0f, 1-boost::math::tools::digits<T>()), 1.0f / 20.0f));
387 // The new safe version of the limit check is now here.
388 const float d2_minus_one = ((digits10_of_type / 0.301F) - 1.0F);
389 const float limit = ceil(exp((d2_minus_one * log(2.0F)) / 20.0F));
391 min_arg = (int) ((std::min)(digits10_of_type * 1.7F, limit));
397 template <class T, class Policy>
398 T scaled_tgamma_no_lanczos(const T& z, const Policy& pol, bool islog = false)
402 // Calculates tgamma(z) / (z/e)^z
403 // Requires that our argument is large enough for Sterling's approximation to hold.
404 // Used internally when combining gamma's of similar magnitude without logarithms.
406 BOOST_MATH_ASSERT(minimum_argument_for_bernoulli_recursion<T>() <= z);
408 // Perform the Bernoulli series expansion of Stirling's approximation.
410 const std::size_t number_of_bernoullis_b2n = policies::get_max_series_iterations<Policy>();
412 T one_over_x_pow_two_n_minus_one = 1 / z;
413 const T one_over_x2 = one_over_x_pow_two_n_minus_one * one_over_x_pow_two_n_minus_one;
414 T sum = (boost::math::bernoulli_b2n<T>(1) / 2) * one_over_x_pow_two_n_minus_one;
415 const T target_epsilon_to_break_loop = sum * boost::math::tools::epsilon<T>();
416 const T half_ln_two_pi_over_z = sqrt(boost::math::constants::two_pi<T>() / z);
417 T last_term = 2 * sum;
419 for (std::size_t n = 2U;; ++n)
421 one_over_x_pow_two_n_minus_one *= one_over_x2;
423 const std::size_t n2 = static_cast<std::size_t>(n * 2U);
425 const T term = (boost::math::bernoulli_b2n<T>(static_cast<int>(n)) * one_over_x_pow_two_n_minus_one) / (n2 * (n2 - 1U));
427 if ((n >= 3U) && (abs(term) < target_epsilon_to_break_loop))
429 // We have reached the desired precision in Stirling's expansion.
430 // Adding additional terms to the sum of this divergent asymptotic
431 // expansion will not improve the result.
433 // Break from the loop.
436 if (n > number_of_bernoullis_b2n)
437 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);
441 // Sanity check for divergence:
442 T fterm = fabs(term);
443 if(fterm > last_term)
444 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);
448 // Complete Stirling's approximation.
449 T scaled_gamma_value = islog ? T(sum + log(half_ln_two_pi_over_z)) : T(exp(sum) * half_ln_two_pi_over_z);
450 return scaled_gamma_value;
453 // Forward declaration of the lgamma_imp template specialization.
454 template <class T, class Policy>
455 T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign = 0);
457 template <class T, class Policy>
458 T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&)
462 static const char* function = "boost::math::tgamma<%1%>(%1%)";
464 // Check if the argument of tgamma is identically zero.
465 const bool is_at_zero = (z == 0);
467 if((boost::math::isnan)(z) || (is_at_zero) || ((boost::math::isinf)(z) && (z < 0)))
468 return policies::raise_domain_error<T>(function, "Evaluation of tgamma at %1%.", z, pol);
470 const bool b_neg = (z < 0);
472 const bool floor_of_z_is_equal_to_z = (floor(z) == z);
474 // Special case handling of small factorials:
475 if((!b_neg) && floor_of_z_is_equal_to_z && (z < boost::math::max_factorial<T>::value))
477 return boost::math::unchecked_factorial<T>(itrunc(z) - 1);
480 // Make a local, unsigned copy of the input argument.
481 T zz((!b_neg) ? z : -z);
483 // Special case for ultra-small z:
484 if(zz < tools::cbrt_epsilon<T>())
487 const T a1(boost::math::constants::euler<T>());
488 const T six_euler_squared((boost::math::constants::euler<T>() * boost::math::constants::euler<T>()) * 6);
489 const T a2((six_euler_squared - boost::math::constants::pi_sqr<T>()) / 12);
491 const T inverse_tgamma_series = z * ((a2 * z + a1) * z + a0);
493 return 1 / inverse_tgamma_series;
496 // Scale the argument up for the calculation of lgamma,
497 // and use downward recursion later for the final result.
498 const int min_arg_for_recursion = minimum_argument_for_bernoulli_recursion<T>();
502 if(zz < min_arg_for_recursion)
504 n_recur = boost::math::itrunc(min_arg_for_recursion - zz) + 1;
514 if (zz > tools::log_max_value<T>())
515 return policies::raise_overflow_error<T>(function, 0, pol);
516 if (log(zz) * zz / 2 > tools::log_max_value<T>())
517 return policies::raise_overflow_error<T>(function, 0, pol);
519 T gamma_value = scaled_tgamma_no_lanczos(zz, pol);
520 T power_term = pow(zz, zz / 2);
521 T exp_term = exp(-zz);
522 gamma_value *= (power_term * exp_term);
523 if(!n_recur && (tools::max_value<T>() / power_term < gamma_value))
524 return policies::raise_overflow_error<T>(function, 0, pol);
525 gamma_value *= power_term;
527 // Rescale the result using downward recursion if necessary.
530 // The order of divides is important, if we keep subtracting 1 from zz
531 // we DO NOT get back to z (cancellation error). Further if z < epsilon
532 // we would end up dividing by zero. Also in order to prevent spurious
533 // overflow with the first division, we must save dividing by |z| till last,
534 // so the optimal order of divides is z+1, z+2, z+3...z+n_recur-1,z.
536 for(int k = 1; k < n_recur; ++k)
541 gamma_value /= fabs(z);
544 // Return the result, accounting for possible negative arguments.
547 // Provide special error analysis for:
548 // * arguments in the neighborhood of a negative integer
549 // * arguments exactly equal to a negative integer.
551 // Check if the argument of tgamma is exactly equal to a negative integer.
552 if(floor_of_z_is_equal_to_z)
553 return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
555 gamma_value *= sinpx(z);
557 BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value);
559 const bool result_is_too_large_to_represent = ( (abs(gamma_value) < 1)
560 && ((tools::max_value<T>() * abs(gamma_value)) < boost::math::constants::pi<T>()));
562 if(result_is_too_large_to_represent)
563 return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
565 gamma_value = -boost::math::constants::pi<T>() / gamma_value;
566 BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value);
569 return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
571 if((boost::math::fpclassify)(gamma_value) == static_cast<int>(FP_SUBNORMAL))
572 return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", gamma_value, pol);
578 template <class T, class Policy>
579 inline T log_gamma_near_1(const T& z, Policy const& pol)
582 // This is for the multiprecision case where there is
583 // no lanczos support, use a taylor series at z = 1,
584 // see https://www.wolframalpha.com/input/?i=taylor+series+lgamma(x)+at+x+%3D+1
586 BOOST_MATH_STD_USING // ADL of std names
588 BOOST_MATH_ASSERT(fabs(z) < 1);
590 T result = -constants::euler<T>() * z;
592 T power_term = z * z / 2;
598 term = power_term * boost::math::polygamma(n - 1, T(1), pol);
602 } while (fabs(result) * tools::epsilon<T>() < fabs(term));
607 template <class T, class Policy>
608 T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign)
612 static const char* function = "boost::math::lgamma<%1%>(%1%)";
614 // Check if the argument of lgamma is identically zero.
615 const bool is_at_zero = (z == 0);
618 return policies::raise_domain_error<T>(function, "Evaluation of lgamma at zero %1%.", z, pol);
619 if((boost::math::isnan)(z))
620 return policies::raise_domain_error<T>(function, "Evaluation of lgamma at %1%.", z, pol);
621 if((boost::math::isinf)(z))
622 return policies::raise_overflow_error<T>(function, 0, pol);
624 const bool b_neg = (z < 0);
626 const bool floor_of_z_is_equal_to_z = (floor(z) == z);
628 // Special case handling of small factorials:
629 if((!b_neg) && floor_of_z_is_equal_to_z && (z < boost::math::max_factorial<T>::value))
633 return log(boost::math::unchecked_factorial<T>(itrunc(z) - 1));
636 // Make a local, unsigned copy of the input argument.
637 T zz((!b_neg) ? z : -z);
639 const int min_arg_for_recursion = minimum_argument_for_bernoulli_recursion<T>();
643 if (zz < min_arg_for_recursion)
645 // Here we simply take the logarithm of tgamma(). This is somewhat
646 // inefficient, but simple. The rationale is that the argument here
647 // is relatively small and overflow is not expected to be likely.
650 if(fabs(z - 1) < 0.25)
652 log_gamma_value = log_gamma_near_1(T(zz - 1), pol);
654 else if(fabs(z - 2) < 0.25)
656 log_gamma_value = log_gamma_near_1(T(zz - 2), pol) + log(zz - 1);
658 else if (z > -tools::root_epsilon<T>())
660 // Reflection formula may fail if z is very close to zero, let the series
661 // expansion for tgamma close to zero do the work:
663 *sign = z < 0 ? -1 : 1;
664 return log(abs(gamma_imp(z, pol, lanczos::undefined_lanczos())));
668 // No issue with spurious overflow in reflection formula,
669 // just fall through to regular code:
670 T g = gamma_imp(zz, pol, lanczos::undefined_lanczos());
673 *sign = g < 0 ? -1 : 1;
675 log_gamma_value = log(abs(g));
680 // Perform the Bernoulli series expansion of Stirling's approximation.
681 T sum = scaled_tgamma_no_lanczos(zz, pol, true);
682 log_gamma_value = zz * (log(zz) - 1) + sum;
685 int sign_of_result = 1;
689 // Provide special error analysis if the argument is exactly
690 // equal to a negative integer.
692 // Check if the argument of lgamma is exactly equal to a negative integer.
693 if(floor_of_z_is_equal_to_z)
694 return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%.", z, pol);
704 sign_of_result = -sign_of_result;
707 log_gamma_value = - log_gamma_value
708 + log(boost::math::constants::pi<T>())
712 if(sign != static_cast<int*>(0U)) { *sign = sign_of_result; }
714 return log_gamma_value;
718 // This helper calculates tgamma(dz+1)-1 without cancellation errors,
719 // used by the upper incomplete gamma with z < 1:
721 template <class T, class Policy, class Lanczos>
722 T tgammap1m1_imp(T dz, Policy const& pol, const Lanczos& l)
726 typedef typename policies::precision<T,Policy>::type precision_type;
728 typedef std::integral_constant<int,
729 precision_type::value <= 0 ? 0 :
730 precision_type::value <= 64 ? 64 :
731 precision_type::value <= 113 ? 113 : 0
739 // Best method is simply to subtract 1 from tgamma:
740 result = boost::math::tgamma(1+dz, pol) - 1;
741 BOOST_MATH_INSTRUMENT_CODE(result);
745 // Use expm1 on lgamma:
746 result = boost::math::expm1(-boost::math::log1p(dz, pol)
747 + lgamma_small_imp<T>(dz+2, dz + 1, dz, tag_type(), pol, l), pol);
748 BOOST_MATH_INSTRUMENT_CODE(result);
755 // Use expm1 on lgamma:
756 result = boost::math::expm1(lgamma_small_imp<T>(dz+1, dz, dz-1, tag_type(), pol, l), pol);
757 BOOST_MATH_INSTRUMENT_CODE(result);
761 // Best method is simply to subtract 1 from tgamma:
762 result = boost::math::tgamma(1+dz, pol) - 1;
763 BOOST_MATH_INSTRUMENT_CODE(result);
770 template <class T, class Policy>
771 inline T tgammap1m1_imp(T z, Policy const& pol,
772 const ::boost::math::lanczos::undefined_lanczos&)
774 BOOST_MATH_STD_USING // ADL of std names
778 return boost::math::expm1(log_gamma_near_1(z, pol));
780 return boost::math::expm1(boost::math::lgamma(1 + z, pol));
784 // Series representation for upper fraction when z is small:
787 struct small_gamma2_series
789 typedef T result_type;
791 small_gamma2_series(T a_, T x_) : result(-x_), x(-x_), apn(a_+1), n(1){}
795 T r = result / (apn);
807 // calculate power term prefix (z^a)(e^-z) used in the non-normalised
808 // incomplete gammas:
810 template <class T, class Policy>
811 T full_igamma_prefix(T a, T z, const Policy& pol)
816 if (z > tools::max_value<T>())
822 if((alz < tools::log_max_value<T>()) && (-z > tools::log_min_value<T>()))
824 prefix = pow(z, a) * exp(-z);
828 prefix = pow(z / exp(z/a), a);
832 prefix = exp(alz - z);
837 if(alz > tools::log_min_value<T>())
839 prefix = pow(z, a) * exp(-z);
841 else if(z/a < tools::log_max_value<T>())
843 prefix = pow(z / exp(z/a), a);
847 prefix = exp(alz - z);
851 // This error handling isn't very good: it happens after the fact
852 // rather than before it...
854 if((boost::math::fpclassify)(prefix) == (int)FP_INFINITE)
855 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);
860 // Compute (z^a)(e^-z)/tgamma(a)
861 // most if the error occurs in this function:
863 template <class T, class Policy, class Lanczos>
864 T regularised_gamma_prefix(T a, T z, const Policy& pol, const Lanczos& l)
867 if (z >= tools::max_value<T>())
869 T agh = a + static_cast<T>(Lanczos::g()) - T(0.5);
871 T d = ((z - a) - static_cast<T>(Lanczos::g()) + T(0.5)) / agh;
876 // We have to treat a < 1 as a special case because our Lanczos
877 // approximations are optimised against the factorials with a > 1,
878 // and for high precision types especially (128-bit reals for example)
879 // very small values of a can give rather erroneous results for gamma
880 // unless we do this:
882 // TODO: is this still required? Lanczos approx should be better now?
884 if(z <= tools::log_min_value<T>())
886 // Oh dear, have to use logs, should be free of cancellation errors though:
887 return exp(a * log(z) - z - lgamma_imp(a, pol, l));
891 // direct calculation, no danger of overflow as gamma(a) < 1/a
893 return pow(z, a) * exp(-z) / gamma_imp(a, pol, l);
896 else if((fabs(d*d*a) <= 100) && (a > 150))
898 // special case for large a and a ~ z.
899 prefix = a * boost::math::log1pmx(d, pol) + z * static_cast<T>(0.5 - Lanczos::g()) / agh;
900 prefix = exp(prefix);
906 // direct computation is most accurate, but use various fallbacks
907 // for different parts of the problem domain:
909 T alz = a * log(z / agh);
911 if(((std::min)(alz, amz) <= tools::log_min_value<T>()) || ((std::max)(alz, amz) >= tools::log_max_value<T>()))
914 if(((std::min)(alz, amz)/2 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/2 < tools::log_max_value<T>()))
916 // compute square root of the result and then square it:
917 T sq = pow(z / agh, a / 2) * exp(amz / 2);
920 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 // compute the 4th root of the result then square it twice:
923 T sq = pow(z / agh, a / 4) * exp(amz / 4);
927 else if((amza > tools::log_min_value<T>()) && (amza < tools::log_max_value<T>()))
929 prefix = pow((z * exp(amza)) / agh, a);
933 prefix = exp(alz + amz);
938 prefix = pow(z / agh, a) * exp(amz);
941 prefix *= sqrt(agh / boost::math::constants::e<T>()) / Lanczos::lanczos_sum_expG_scaled(a);
945 // And again, without Lanczos support:
947 template <class T, class Policy>
948 T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined_lanczos& l)
952 if((a < 1) && (z < 1))
954 // No overflow possible since the power terms tend to unity as a,z -> 0
955 return pow(z, a) * exp(-z) / boost::math::tgamma(a, pol);
957 else if(a > minimum_argument_for_bernoulli_recursion<T>())
959 T scaled_gamma = scaled_tgamma_no_lanczos(a, pol);
960 T power_term = pow(z / a, a / 2);
962 if ((0 == power_term) || (fabs(a_minus_z) > tools::log_max_value<T>()))
964 // The result is probably zero, but we need to be sure:
965 return exp(a * log(z / a) + a_minus_z - log(scaled_gamma));
967 return (power_term * exp(a_minus_z)) * (power_term / scaled_gamma);
972 // Usual case is to calculate the prefix at a+shift and recurse down
973 // to the value we want:
975 const int min_z = minimum_argument_for_bernoulli_recursion<T>();
976 long shift = 1 + ltrunc(min_z - a);
977 T result = regularised_gamma_prefix(T(a + shift), z, pol, l);
980 for (long i = 0; i < shift; ++i)
990 // We failed, most probably we have z << 1, try again, this time
991 // we calculate z^a e^-z / tgamma(a+shift), combining power terms
992 // as we go. And again recurse down to the result.
994 T scaled_gamma = scaled_tgamma_no_lanczos(T(a + shift), pol);
995 T power_term_1 = pow(z / (a + shift), a);
996 T power_term_2 = pow(a + shift, -shift);
997 T power_term_3 = exp(a + shift - z);
998 if ((0 == power_term_1) || (0 == power_term_2) || (0 == power_term_3) || (fabs(a + shift - z) > tools::log_max_value<T>()))
1000 // We have no test case that gets here, most likely the type T
1001 // has a high precision but low exponent range:
1002 return exp(a * log(z) - z - boost::math::lgamma(a, pol));
1004 result = power_term_1 * power_term_2 * power_term_3 / scaled_gamma;
1005 for (long i = 0; i < shift; ++i)
1014 // Upper gamma fraction for very small a:
1016 template <class T, class Policy>
1017 inline T tgamma_small_upper_part(T a, T x, const Policy& pol, T* pgam = 0, bool invert = false, T* pderivative = 0)
1019 BOOST_MATH_STD_USING // ADL of std functions.
1021 // Compute the full upper fraction (Q) when a is very small:
1024 result = boost::math::tgamma1pm1(a, pol);
1026 *pgam = (result + 1) / a;
1027 T p = boost::math::powm1(x, a, pol);
1030 detail::small_gamma2_series<T> s(a, x);
1031 std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>() - 10;
1034 *pderivative = p / (*pgam * exp(x));
1035 T init_value = invert ? *pgam : 0;
1036 result = -p * tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, (init_value - result) / p);
1037 policies::check_series_iterations<T>("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)", max_iter, pol);
1043 // Upper gamma fraction for integer a:
1045 template <class T, class Policy>
1046 inline T finite_gamma_q(T a, T x, Policy const& pol, T* pderivative = 0)
1049 // Calculates normalised Q when a is an integer:
1051 BOOST_MATH_STD_USING
1057 for(unsigned n = 1; n < a; ++n)
1066 *pderivative = e * pow(x, a) / boost::math::unchecked_factorial<T>(itrunc(T(a - 1), pol));
1071 // Upper gamma fraction for half integer a:
1073 template <class T, class Policy>
1074 T finite_half_gamma_q(T a, T x, T* p_derivative, const Policy& pol)
1077 // Calculates normalised Q when a is a half-integer:
1079 BOOST_MATH_STD_USING
1080 T e = boost::math::erfc(sqrt(x), pol);
1081 if((e != 0) && (a > 1))
1083 T term = exp(-x) / sqrt(constants::pi<T>() * x);
1085 static const T half = T(1) / 2;
1088 for(unsigned n = 2; n < a; ++n)
1100 else if(p_derivative)
1102 // We'll be dividing by x later, so calculate derivative * x:
1103 *p_derivative = sqrt(x) * exp(-x) / constants::root_pi<T>();
1108 // Asymptotic approximation for large argument, see: https://dlmf.nist.gov/8.11#E2
1111 struct incomplete_tgamma_large_x_series
1113 typedef T result_type;
1114 incomplete_tgamma_large_x_series(const T& a, const T& x)
1115 : a_poch(a - 1), z(x), term(1) {}
1126 template <class T, class Policy>
1127 T incomplete_tgamma_large_x(const T& a, const T& x, const Policy& pol)
1129 BOOST_MATH_STD_USING
1130 incomplete_tgamma_large_x_series<T> s(a, x);
1131 std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
1132 T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
1133 boost::math::policies::check_series_iterations<T>("boost::math::tgamma<%1%>(%1%,%1%)", max_iter, pol);
1139 // Main incomplete gamma entry point, handles all four incomplete gamma's:
1141 template <class T, class Policy>
1142 T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
1143 const Policy& pol, T* p_derivative)
1145 static const char* function = "boost::math::gamma_p<%1%>(%1%, %1%)";
1147 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 return policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
1151 BOOST_MATH_STD_USING
1153 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1155 T result = 0; // Just to avoid warning C4701: potentially uninitialized local variable 'result' used
1157 if(a >= max_factorial<T>::value && !normalised)
1160 // When we're computing the non-normalized incomplete gamma
1161 // and a is large the result is rather hard to compute unless
1162 // we use logs. There are really two options - if x is a long
1163 // way from a in value then we can reliably use methods 2 and 4
1164 // below in logarithmic form and go straight to the result.
1165 // Otherwise we let the regularized gamma take the strain
1166 // (the result is unlikely to underflow in the central region anyway)
1167 // and combine with lgamma in the hopes that we get a finite result.
1169 if(invert && (a * 4 < x))
1171 // This is method 4 below, done in logs:
1172 result = a * log(x) - x;
1174 *p_derivative = exp(result);
1175 result += log(upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>()));
1177 else if(!invert && (a > 4 * x))
1179 // This is method 2 below, done in logs:
1180 result = a * log(x) - x;
1182 *p_derivative = exp(result);
1184 result += log(detail::lower_gamma_series(a, x, pol, init_value) / a);
1188 result = gamma_incomplete_imp(a, x, true, invert, pol, p_derivative);
1193 // Try http://functions.wolfram.com/06.06.06.0039.01
1194 result = 1 + 1 / (12 * a) + 1 / (288 * a * a);
1195 result = log(result) - a + (a - 0.5f) * log(a) + log(boost::math::constants::root_two_pi<T>());
1197 *p_derivative = exp(a * log(x) - x);
1201 // This is method 2 below, done in logs, we're really outside the
1202 // range of this method, but since the result is almost certainly
1203 // infinite, we should probably be OK:
1204 result = a * log(x) - x;
1206 *p_derivative = exp(result);
1208 result += log(detail::lower_gamma_series(a, x, pol, init_value) / a);
1213 result = log(result) + boost::math::lgamma(a, pol);
1216 if(result > tools::log_max_value<T>())
1217 return policies::raise_overflow_error<T>(function, 0, pol);
1221 BOOST_MATH_ASSERT((p_derivative == 0) || normalised);
1223 bool is_int, is_half_int;
1224 bool is_small_a = (a < 30) && (a <= x + 1) && (x < tools::log_max_value<T>());
1229 is_half_int = is_int ? false : (fabs(fa - a) == 0.5f);
1233 is_int = is_half_int = false;
1238 if(is_int && (x > 0.6))
1240 // calculate Q via finite sum:
1244 else if(is_half_int && (x > 0.2))
1246 // calculate Q via finite sum for half integer a:
1250 else if((x < tools::root_epsilon<T>()) && (a > 1))
1254 else if ((x > 1000) && ((a < x) || (fabs(a - 50) / x < 1)))
1256 // calculate Q via asymptotic approximation:
1263 // Changeover criterion chosen to give a changeover at Q ~ 0.33
1265 if(-0.4 / log(x) < a)
1277 // Changover here occurs when P ~ 0.75 or Q ~ 0.25:
1291 // Begin by testing whether we're in the "bad" zone
1292 // where the result will be near 0.5 and the usual
1293 // series and continued fractions are slow to converge:
1295 bool use_temme = false;
1296 if(normalised && std::numeric_limits<T>::is_specialized && (a > 20))
1298 T sigma = fabs((x-a)/a);
1299 if((a > 200) && (policies::digits<T, Policy>() <= 113))
1302 // This limit is chosen so that we use Temme's expansion
1303 // only if the result would be larger than about 10^-6.
1304 // Below that the regular series and continued fractions
1305 // converge OK, and if we use Temme's method we get increasing
1306 // errors from the dominant erfc term as it's (inexact) argument
1307 // increases in magnitude.
1309 if(20 / a > sigma * sigma)
1312 else if(policies::digits<T, Policy>() <= 64)
1314 // Note in this zone we can't use Temme's expansion for
1315 // types longer than an 80-bit real:
1316 // it would require too many terms in the polynomials.
1328 // Regular case where the result will not be too close to 0.5.
1330 // Changeover here occurs at P ~ Q ~ 0.5
1331 // Note that series computation of P is about x2 faster than continued fraction
1332 // calculation of Q, so try and use the CF only when really necessary, especially
1335 if(x - (1 / (3 * x)) < a)
1351 result = finite_gamma_q(a, x, pol, p_derivative);
1353 result *= boost::math::tgamma(a, pol);
1358 result = finite_half_gamma_q(a, x, p_derivative, pol);
1360 result *= boost::math::tgamma(a, pol);
1361 if(p_derivative && (*p_derivative == 0))
1362 *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
1368 result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
1370 *p_derivative = result;
1374 // If we're going to be inverting the result then we can
1375 // reduce the number of series evaluations by quite
1376 // a few iterations if we set an initial value for the
1377 // series sum based on what we'll end up subtracting it from
1379 // Have to be careful though that this optimization doesn't
1380 // lead to spurious numeric overflow. Note that the
1381 // scary/expensive overflow checks below are more often
1382 // than not bypassed in practice for "sensible" input
1386 bool optimised_invert = false;
1389 init_value = (normalised ? 1 : boost::math::tgamma(a, pol));
1390 if(normalised || (result >= 1) || (tools::max_value<T>() * result > init_value))
1392 init_value /= result;
1393 if(normalised || (a < 1) || (tools::max_value<T>() / a > init_value))
1396 optimised_invert = true;
1404 result *= detail::lower_gamma_series(a, x, pol, init_value) / a;
1405 if(optimised_invert)
1418 result = tgamma_small_upper_part(a, x, pol, &g, invert, p_derivative);
1427 result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
1429 *p_derivative = result;
1431 result *= upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>());
1437 // Use compile time dispatch to the appropriate
1438 // Temme asymptotic expansion. This may be dead code
1439 // if T does not have numeric limits support, or has
1440 // too many digits for the most precise version of
1441 // these expansions, in that case we'll be calling
1442 // an empty function.
1444 typedef typename policies::precision<T, Policy>::type precision_type;
1446 typedef std::integral_constant<int,
1447 precision_type::value <= 0 ? 0 :
1448 precision_type::value <= 53 ? 53 :
1449 precision_type::value <= 64 ? 64 :
1450 precision_type::value <= 113 ? 113 : 0
1453 result = igamma_temme_large(a, x, pol, static_cast<tag_type const*>(0));
1457 *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
1462 // x is so small that P is necessarily very small too,
1463 // use http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/
1465 result = pow(x, a) / (a);
1470 result = pow(x, a) / boost::math::tgamma(a + 1, pol);
1472 catch (const std::overflow_error&)
1477 result *= 1 - a * x / (a + 1);
1479 *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
1486 result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
1488 *p_derivative = result;
1491 result *= incomplete_tgamma_large_x(a, x, pol);
1496 if(normalised && (result > 1))
1500 T gam = normalised ? 1 : boost::math::tgamma(a, pol);
1501 result = gam - result;
1506 // Need to convert prefix term to derivative:
1508 if((x < 1) && (tools::max_value<T>() * x < *p_derivative))
1510 // overflow, just return an arbitrarily large value:
1511 *p_derivative = tools::max_value<T>() / 2;
1521 // Ratios of two gamma functions:
1523 template <class T, class Policy, class Lanczos>
1524 T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const Lanczos& l)
1526 BOOST_MATH_STD_USING
1527 if(z < tools::epsilon<T>())
1530 // We get spurious numeric overflow unless we're very careful, this
1531 // can occur either inside Lanczos::lanczos_sum(z) or in the
1532 // final combination of terms, to avoid this, split the product up
1533 // into 2 (or 3) parts:
1535 // G(z) / G(L) = 1 / (z * G(L)) ; z < eps, L = z + delta = delta
1536 // z * G(L) = z * G(lim) * (G(L)/G(lim)) ; lim = largest factorial
1538 if(boost::math::max_factorial<T>::value < delta)
1540 T ratio = tgamma_delta_ratio_imp_lanczos(delta, T(boost::math::max_factorial<T>::value - delta), pol, l);
1542 ratio *= boost::math::unchecked_factorial<T>(boost::math::max_factorial<T>::value - 1);
1547 return 1 / (z * boost::math::tgamma(z + delta, pol));
1550 T zgh = static_cast<T>(z + Lanczos::g() - constants::half<T>());
1554 if (fabs(delta / zgh) < boost::math::tools::epsilon<T>())
1557 // result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol));
1559 // log1p(delta / zgh) = delta / zgh = delta / z
1560 // multiplying we get -delta.
1561 result = exp(-delta);
1564 // from the pow formula below... but this may actually be wrong, we just can't really calculate it :(
1569 if(fabs(delta) < 10)
1571 result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol));
1575 result = pow(zgh / (zgh + delta), z - constants::half<T>());
1577 // Split the calculation up to avoid spurious overflow:
1578 result *= Lanczos::lanczos_sum(z) / Lanczos::lanczos_sum(T(z + delta));
1580 result *= pow(constants::e<T>() / (zgh + delta), delta);
1584 // And again without Lanczos support this time:
1586 template <class T, class Policy>
1587 T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const lanczos::undefined_lanczos& l)
1589 BOOST_MATH_STD_USING
1592 // We adjust z and delta so that both z and z+delta are large enough for
1593 // Sterling's approximation to hold. We can then calculate the ratio
1594 // for the adjusted values, and rescale back down to z and z+delta.
1596 // Get the required shifts first:
1598 long numerator_shift = 0;
1599 long denominator_shift = 0;
1600 const int min_z = minimum_argument_for_bernoulli_recursion<T>();
1603 numerator_shift = 1 + ltrunc(min_z - z);
1604 if (min_z > z + delta)
1605 denominator_shift = 1 + ltrunc(min_z - z - delta);
1607 // If the shifts are zero, then we can just combine scaled tgamma's
1608 // and combine the remaining terms:
1610 if (numerator_shift == 0 && denominator_shift == 0)
1612 T scaled_tgamma_num = scaled_tgamma_no_lanczos(z, pol);
1613 T scaled_tgamma_denom = scaled_tgamma_no_lanczos(T(z + delta), pol);
1614 T result = scaled_tgamma_num / scaled_tgamma_denom;
1615 result *= exp(z * boost::math::log1p(-delta / (z + delta), pol)) * pow((delta + z) / constants::e<T>(), -delta);
1619 // We're going to have to rescale first, get the adjusted z and delta values,
1620 // plus the ratio for the adjusted values:
1622 T zz = z + numerator_shift;
1623 T dd = delta - (numerator_shift - denominator_shift);
1624 T ratio = tgamma_delta_ratio_imp_lanczos(zz, dd, pol, l);
1626 // Use gamma recurrence relations to get back to the original
1629 for (long long i = 0; i < numerator_shift; ++i)
1632 if (i < denominator_shift)
1633 ratio *= (z + delta + i);
1635 for (long long i = numerator_shift; i < denominator_shift; ++i)
1637 ratio *= (z + delta + i);
1642 template <class T, class Policy>
1643 T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol)
1645 BOOST_MATH_STD_USING
1647 if((z <= 0) || (z + delta <= 0))
1649 // This isn't very sophisticated, or accurate, but it does work:
1650 return boost::math::tgamma(z, pol) / boost::math::tgamma(z + delta, pol);
1653 if(floor(delta) == delta)
1658 // Both z and delta are integers, see if we can just use table lookup
1659 // of the factorials to get the result:
1661 if((z <= max_factorial<T>::value) && (z + delta <= max_factorial<T>::value))
1663 return unchecked_factorial<T>((unsigned)itrunc(z, pol) - 1) / unchecked_factorial<T>((unsigned)itrunc(T(z + delta), pol) - 1);
1666 if(fabs(delta) < 20)
1669 // delta is a small integer, we can use a finite product:
1677 while(0 != (delta += 1))
1687 while(0 != (delta -= 1))
1696 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1697 return tgamma_delta_ratio_imp_lanczos(z, delta, pol, lanczos_type());
1700 template <class T, class Policy>
1701 T tgamma_ratio_imp(T x, T y, const Policy& pol)
1703 BOOST_MATH_STD_USING
1705 if((x <= 0) || (boost::math::isinf)(x))
1706 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);
1707 if((y <= 0) || (boost::math::isinf)(y))
1708 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);
1710 if(x <= tools::min_value<T>())
1712 // Special case for denorms...Ugh.
1713 T shift = ldexp(T(1), tools::digits<T>());
1714 return shift * tgamma_ratio_imp(T(x * shift), y, pol);
1717 if((x < max_factorial<T>::value) && (y < max_factorial<T>::value))
1719 // Rather than subtracting values, lets just call the gamma functions directly:
1720 return boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
1725 if(y < 2 * max_factorial<T>::value)
1727 // We need to sidestep on x as well, otherwise we'll underflow
1728 // before we get to factor in the prefix term:
1731 while(y >= max_factorial<T>::value)
1736 return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
1739 // result is almost certainly going to underflow to zero, try logs just in case:
1741 return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
1745 if(x < 2 * max_factorial<T>::value)
1747 // We need to sidestep on y as well, otherwise we'll overflow
1748 // before we get to factor in the prefix term:
1751 while(x >= max_factorial<T>::value)
1756 return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
1759 // Result will almost certainly overflow, try logs just in case:
1761 return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
1764 // Regular case, x and y both large and similar in magnitude:
1766 return boost::math::tgamma_delta_ratio(x, y - x, pol);
1769 template <class T, class Policy>
1770 T gamma_p_derivative_imp(T a, T x, const Policy& pol)
1772 BOOST_MATH_STD_USING
1774 // Usual error checks first:
1777 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);
1779 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);
1781 // Now special cases:
1785 return (a > 1) ? 0 :
1786 (a == 1) ? 1 : policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
1791 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1792 T f1 = detail::regularised_gamma_prefix(a, x, pol, lanczos_type());
1793 if((x < 1) && (tools::max_value<T>() * x < f1))
1796 return policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
1800 // Underflow in calculation, use logs instead:
1801 f1 = a * log(x) - x - lgamma(a, pol) - log(x);
1810 template <class T, class Policy>
1811 inline typename tools::promote_args<T>::type
1812 tgamma(T z, const Policy& /* pol */, const std::true_type)
1814 BOOST_FPU_EXCEPTION_GUARD
1815 typedef typename tools::promote_args<T>::type result_type;
1816 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1817 typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1818 typedef typename policies::normalise<
1820 policies::promote_float<false>,
1821 policies::promote_double<false>,
1822 policies::discrete_quantile<>,
1823 policies::assert_undefined<> >::type forwarding_policy;
1824 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%)");
1827 template <class T, class Policy>
1828 struct igamma_initializer
1834 typedef typename policies::precision<T, Policy>::type precision_type;
1836 typedef std::integral_constant<int,
1837 precision_type::value <= 0 ? 0 :
1838 precision_type::value <= 53 ? 53 :
1839 precision_type::value <= 64 ? 64 :
1840 precision_type::value <= 113 ? 113 : 0
1843 do_init(tag_type());
1846 static void do_init(const std::integral_constant<int, N>&)
1848 // If std::numeric_limits<T>::digits is zero, we must not call
1849 // our initialization code here as the precision presumably
1850 // varies at runtime, and will not have been set yet. Plus the
1851 // code requiring initialization isn't called when digits == 0.
1852 if(std::numeric_limits<T>::digits)
1854 boost::math::gamma_p(static_cast<T>(400), static_cast<T>(400), Policy());
1857 static void do_init(const std::integral_constant<int, 53>&){}
1858 void force_instantiate()const{}
1860 static const init initializer;
1861 static void force_instantiate()
1863 initializer.force_instantiate();
1867 template <class T, class Policy>
1868 const typename igamma_initializer<T, Policy>::init igamma_initializer<T, Policy>::initializer;
1870 template <class T, class Policy>
1871 struct lgamma_initializer
1877 typedef typename policies::precision<T, Policy>::type precision_type;
1878 typedef std::integral_constant<int,
1879 precision_type::value <= 0 ? 0 :
1880 precision_type::value <= 64 ? 64 :
1881 precision_type::value <= 113 ? 113 : 0
1884 do_init(tag_type());
1886 static void do_init(const std::integral_constant<int, 64>&)
1888 boost::math::lgamma(static_cast<T>(2.5), Policy());
1889 boost::math::lgamma(static_cast<T>(1.25), Policy());
1890 boost::math::lgamma(static_cast<T>(1.75), Policy());
1892 static void do_init(const std::integral_constant<int, 113>&)
1894 boost::math::lgamma(static_cast<T>(2.5), Policy());
1895 boost::math::lgamma(static_cast<T>(1.25), Policy());
1896 boost::math::lgamma(static_cast<T>(1.5), Policy());
1897 boost::math::lgamma(static_cast<T>(1.75), Policy());
1899 static void do_init(const std::integral_constant<int, 0>&)
1902 void force_instantiate()const{}
1904 static const init initializer;
1905 static void force_instantiate()
1907 initializer.force_instantiate();
1911 template <class T, class Policy>
1912 const typename lgamma_initializer<T, Policy>::init lgamma_initializer<T, Policy>::initializer;
1914 template <class T1, class T2, class Policy>
1915 inline typename tools::promote_args<T1, T2>::type
1916 tgamma(T1 a, T2 z, const Policy&, const std::false_type)
1918 BOOST_FPU_EXCEPTION_GUARD
1919 typedef typename tools::promote_args<T1, T2>::type result_type;
1920 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1921 // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1922 typedef typename policies::normalise<
1924 policies::promote_float<false>,
1925 policies::promote_double<false>,
1926 policies::discrete_quantile<>,
1927 policies::assert_undefined<> >::type forwarding_policy;
1929 igamma_initializer<value_type, forwarding_policy>::force_instantiate();
1931 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
1932 detail::gamma_incomplete_imp(static_cast<value_type>(a),
1933 static_cast<value_type>(z), false, true,
1934 forwarding_policy(), static_cast<value_type*>(0)), "boost::math::tgamma<%1%>(%1%, %1%)");
1937 template <class T1, class T2>
1938 inline typename tools::promote_args<T1, T2>::type
1939 tgamma(T1 a, T2 z, const std::false_type& tag)
1941 return tgamma(a, z, policies::policy<>(), tag);
1945 } // namespace detail
1948 inline typename tools::promote_args<T>::type
1951 return tgamma(z, policies::policy<>());
1954 template <class T, class Policy>
1955 inline typename tools::promote_args<T>::type
1956 lgamma(T z, int* sign, const Policy&)
1958 BOOST_FPU_EXCEPTION_GUARD
1959 typedef typename tools::promote_args<T>::type result_type;
1960 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1961 typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1962 typedef typename policies::normalise<
1964 policies::promote_float<false>,
1965 policies::promote_double<false>,
1966 policies::discrete_quantile<>,
1967 policies::assert_undefined<> >::type forwarding_policy;
1969 detail::lgamma_initializer<value_type, forwarding_policy>::force_instantiate();
1971 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%)");
1975 inline typename tools::promote_args<T>::type
1976 lgamma(T z, int* sign)
1978 return lgamma(z, sign, policies::policy<>());
1981 template <class T, class Policy>
1982 inline typename tools::promote_args<T>::type
1983 lgamma(T x, const Policy& pol)
1985 return ::boost::math::lgamma(x, 0, pol);
1989 inline typename tools::promote_args<T>::type
1992 return ::boost::math::lgamma(x, 0, policies::policy<>());
1995 template <class T, class Policy>
1996 inline typename tools::promote_args<T>::type
1997 tgamma1pm1(T z, const Policy& /* pol */)
1999 BOOST_FPU_EXCEPTION_GUARD
2000 typedef typename tools::promote_args<T>::type result_type;
2001 typedef typename policies::evaluation<result_type, Policy>::type value_type;
2002 typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
2003 typedef typename policies::normalise<
2005 policies::promote_float<false>,
2006 policies::promote_double<false>,
2007 policies::discrete_quantile<>,
2008 policies::assert_undefined<> >::type forwarding_policy;
2010 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%)");
2014 inline typename tools::promote_args<T>::type
2017 return tgamma1pm1(z, policies::policy<>());
2021 // Full upper incomplete gamma:
2023 template <class T1, class T2>
2024 inline typename tools::promote_args<T1, T2>::type
2028 // Type T2 could be a policy object, or a value, select the
2029 // right overload based on T2:
2031 typedef typename policies::is_policy<T2>::type maybe_policy;
2032 return detail::tgamma(a, z, maybe_policy());
2034 template <class T1, class T2, class Policy>
2035 inline typename tools::promote_args<T1, T2>::type
2036 tgamma(T1 a, T2 z, const Policy& pol)
2038 return detail::tgamma(a, z, pol, std::false_type());
2041 // Full lower incomplete gamma:
2043 template <class T1, class T2, class Policy>
2044 inline typename tools::promote_args<T1, T2>::type
2045 tgamma_lower(T1 a, T2 z, const Policy&)
2047 BOOST_FPU_EXCEPTION_GUARD
2048 typedef typename tools::promote_args<T1, T2>::type result_type;
2049 typedef typename policies::evaluation<result_type, Policy>::type value_type;
2050 // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
2051 typedef typename policies::normalise<
2053 policies::promote_float<false>,
2054 policies::promote_double<false>,
2055 policies::discrete_quantile<>,
2056 policies::assert_undefined<> >::type forwarding_policy;
2058 detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
2060 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
2061 detail::gamma_incomplete_imp(static_cast<value_type>(a),
2062 static_cast<value_type>(z), false, false,
2063 forwarding_policy(), static_cast<value_type*>(0)), "tgamma_lower<%1%>(%1%, %1%)");
2065 template <class T1, class T2>
2066 inline typename tools::promote_args<T1, T2>::type
2067 tgamma_lower(T1 a, T2 z)
2069 return tgamma_lower(a, z, policies::policy<>());
2072 // Regularised upper incomplete gamma:
2074 template <class T1, class T2, class Policy>
2075 inline typename tools::promote_args<T1, T2>::type
2076 gamma_q(T1 a, T2 z, const Policy& /* pol */)
2078 BOOST_FPU_EXCEPTION_GUARD
2079 typedef typename tools::promote_args<T1, T2>::type result_type;
2080 typedef typename policies::evaluation<result_type, Policy>::type value_type;
2081 // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
2082 typedef typename policies::normalise<
2084 policies::promote_float<false>,
2085 policies::promote_double<false>,
2086 policies::discrete_quantile<>,
2087 policies::assert_undefined<> >::type forwarding_policy;
2089 detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
2091 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
2092 detail::gamma_incomplete_imp(static_cast<value_type>(a),
2093 static_cast<value_type>(z), true, true,
2094 forwarding_policy(), static_cast<value_type*>(0)), "gamma_q<%1%>(%1%, %1%)");
2096 template <class T1, class T2>
2097 inline typename tools::promote_args<T1, T2>::type
2100 return gamma_q(a, z, policies::policy<>());
2103 // Regularised lower incomplete gamma:
2105 template <class T1, class T2, class Policy>
2106 inline typename tools::promote_args<T1, T2>::type
2107 gamma_p(T1 a, T2 z, const Policy&)
2109 BOOST_FPU_EXCEPTION_GUARD
2110 typedef typename tools::promote_args<T1, T2>::type result_type;
2111 typedef typename policies::evaluation<result_type, Policy>::type value_type;
2112 // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
2113 typedef typename policies::normalise<
2115 policies::promote_float<false>,
2116 policies::promote_double<false>,
2117 policies::discrete_quantile<>,
2118 policies::assert_undefined<> >::type forwarding_policy;
2120 detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
2122 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
2123 detail::gamma_incomplete_imp(static_cast<value_type>(a),
2124 static_cast<value_type>(z), true, false,
2125 forwarding_policy(), static_cast<value_type*>(0)), "gamma_p<%1%>(%1%, %1%)");
2127 template <class T1, class T2>
2128 inline typename tools::promote_args<T1, T2>::type
2131 return gamma_p(a, z, policies::policy<>());
2134 // ratios of gamma functions:
2135 template <class T1, class T2, class Policy>
2136 inline typename tools::promote_args<T1, T2>::type
2137 tgamma_delta_ratio(T1 z, T2 delta, const Policy& /* pol */)
2139 BOOST_FPU_EXCEPTION_GUARD
2140 typedef typename tools::promote_args<T1, T2>::type result_type;
2141 typedef typename policies::evaluation<result_type, Policy>::type value_type;
2142 typedef typename policies::normalise<
2144 policies::promote_float<false>,
2145 policies::promote_double<false>,
2146 policies::discrete_quantile<>,
2147 policies::assert_undefined<> >::type forwarding_policy;
2149 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%)");
2151 template <class T1, class T2>
2152 inline typename tools::promote_args<T1, T2>::type
2153 tgamma_delta_ratio(T1 z, T2 delta)
2155 return tgamma_delta_ratio(z, delta, policies::policy<>());
2157 template <class T1, class T2, class Policy>
2158 inline typename tools::promote_args<T1, T2>::type
2159 tgamma_ratio(T1 a, T2 b, const Policy&)
2161 typedef typename tools::promote_args<T1, T2>::type result_type;
2162 typedef typename policies::evaluation<result_type, Policy>::type value_type;
2163 typedef typename policies::normalise<
2165 policies::promote_float<false>,
2166 policies::promote_double<false>,
2167 policies::discrete_quantile<>,
2168 policies::assert_undefined<> >::type forwarding_policy;
2170 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%)");
2172 template <class T1, class T2>
2173 inline typename tools::promote_args<T1, T2>::type
2174 tgamma_ratio(T1 a, T2 b)
2176 return tgamma_ratio(a, b, policies::policy<>());
2179 template <class T1, class T2, class Policy>
2180 inline typename tools::promote_args<T1, T2>::type
2181 gamma_p_derivative(T1 a, T2 x, const Policy&)
2183 BOOST_FPU_EXCEPTION_GUARD
2184 typedef typename tools::promote_args<T1, T2>::type result_type;
2185 typedef typename policies::evaluation<result_type, Policy>::type value_type;
2186 typedef typename policies::normalise<
2188 policies::promote_float<false>,
2189 policies::promote_double<false>,
2190 policies::discrete_quantile<>,
2191 policies::assert_undefined<> >::type forwarding_policy;
2193 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%)");
2195 template <class T1, class T2>
2196 inline typename tools::promote_args<T1, T2>::type
2197 gamma_p_derivative(T1 a, T2 x)
2199 return gamma_p_derivative(a, x, policies::policy<>());
2203 } // namespace boost
2206 # pragma warning(pop)
2209 #include <boost/math/special_functions/detail/igamma_inverse.hpp>
2210 #include <boost/math/special_functions/detail/gamma_inva.hpp>
2211 #include <boost/math/special_functions/erf.hpp>
2213 #endif // BOOST_MATH_SF_GAMMA_HPP