2 // Copyright John Maddock 2006-7.
3 // Copyright Paul A. Bristow 2007.
5 // Use, modification and distribution are subject to the
6 // Boost Software License, Version 1.0. (See accompanying file
7 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
9 #ifndef BOOST_MATH_SF_GAMMA_HPP
10 #define BOOST_MATH_SF_GAMMA_HPP
16 #include <boost/config.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/policies/error_handling.hpp>
22 #include <boost/math/constants/constants.hpp>
23 #include <boost/math/special_functions/math_fwd.hpp>
24 #include <boost/math/special_functions/log1p.hpp>
25 #include <boost/math/special_functions/trunc.hpp>
26 #include <boost/math/special_functions/powm1.hpp>
27 #include <boost/math/special_functions/sqrt1pm1.hpp>
28 #include <boost/math/special_functions/lanczos.hpp>
29 #include <boost/math/special_functions/fpclassify.hpp>
30 #include <boost/math/special_functions/detail/igamma_large.hpp>
31 #include <boost/math/special_functions/detail/unchecked_factorial.hpp>
32 #include <boost/math/special_functions/detail/lgamma_small.hpp>
33 #include <boost/type_traits/is_convertible.hpp>
34 #include <boost/assert.hpp>
35 #include <boost/mpl/greater.hpp>
36 #include <boost/mpl/equal_to.hpp>
37 #include <boost/mpl/greater.hpp>
39 #include <boost/config/no_tr1/cmath.hpp>
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 boost::true_type&)
60 int i = static_cast<int>(v);
64 inline bool is_odd(T v, const boost::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, ::boost::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.
104 BOOST_ASSERT(fl >= 0);
107 T result = sin(dist*boost::math::constants::pi<T>());
108 return sign*z*result;
109 } // template <class T> T sinpx(T z)
111 // tgamma(z), with Lanczos support:
113 template <class T, class Policy, class Lanczos>
114 T gamma_imp(T z, const Policy& pol, const Lanczos& l)
120 #ifdef BOOST_MATH_INSTRUMENT
121 static bool b = false;
124 std::cout << "tgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;
128 static const char* function = "boost::math::tgamma<%1%>(%1%)";
133 return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
136 result = gamma_imp(T(-z), pol, l) * sinpx(z);
137 BOOST_MATH_INSTRUMENT_VARIABLE(result);
138 if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
139 return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
140 result = -boost::math::constants::pi<T>() / result;
142 return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
143 if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL)
144 return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", result, pol);
145 BOOST_MATH_INSTRUMENT_VARIABLE(result);
156 BOOST_MATH_INSTRUMENT_VARIABLE(result);
157 if((floor(z) == z) && (z < max_factorial<T>::value))
159 result *= unchecked_factorial<T>(itrunc(z, pol) - 1);
160 BOOST_MATH_INSTRUMENT_VARIABLE(result);
164 result *= Lanczos::lanczos_sum(z);
165 T zgh = (z + static_cast<T>(Lanczos::g()) - boost::math::constants::half<T>());
167 BOOST_MATH_INSTRUMENT_VARIABLE(result);
168 BOOST_MATH_INSTRUMENT_VARIABLE(tools::log_max_value<T>());
169 if(z * lzgh > tools::log_max_value<T>())
171 // we're going to overflow unless this is done with care:
172 BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
173 if(lzgh * z / 2 > tools::log_max_value<T>())
174 return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
175 T hp = pow(zgh, (z / 2) - T(0.25));
176 BOOST_MATH_INSTRUMENT_VARIABLE(hp);
177 result *= hp / exp(zgh);
178 BOOST_MATH_INSTRUMENT_VARIABLE(result);
179 if(tools::max_value<T>() / hp < result)
180 return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
182 BOOST_MATH_INSTRUMENT_VARIABLE(result);
186 BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
187 BOOST_MATH_INSTRUMENT_VARIABLE(pow(zgh, z - boost::math::constants::half<T>()));
188 BOOST_MATH_INSTRUMENT_VARIABLE(exp(zgh));
189 result *= pow(zgh, z - boost::math::constants::half<T>()) / exp(zgh);
190 BOOST_MATH_INSTRUMENT_VARIABLE(result);
196 // lgamma(z) with Lanczos support:
198 template <class T, class Policy, class Lanczos>
199 T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = 0)
201 #ifdef BOOST_MATH_INSTRUMENT
202 static bool b = false;
205 std::cout << "lgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;
212 static const char* function = "boost::math::lgamma<%1%>(%1%)";
218 // reflection formula:
220 return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%.", z, pol);
232 result = log(boost::math::constants::pi<T>()) - lgamma_imp(z, pol, l) - log(t);
236 typedef typename policies::precision<T, Policy>::type precision_type;
237 typedef typename mpl::if_<
239 mpl::less_equal<precision_type, mpl::int_<64> >,
240 mpl::greater<precision_type, mpl::int_<0> >
245 mpl::less_equal<precision_type, mpl::int_<113> >,
246 mpl::greater<precision_type, mpl::int_<0> >
248 mpl::int_<113>, mpl::int_<0> >::type
250 result = lgamma_small_imp<T>(z, T(z - 1), T(z - 2), tag_type(), pol, l);
252 else if((z >= 3) && (z < 100) && (std::numeric_limits<T>::max_exponent >= 1024))
254 // taking the log of tgamma reduces the error, no danger of overflow here:
255 result = log(gamma_imp(z, pol, l));
259 // regular evaluation:
260 T zgh = static_cast<T>(z + Lanczos::g() - boost::math::constants::half<T>());
261 result = log(zgh) - 1;
263 result += log(Lanczos::lanczos_sum_expG_scaled(z));
272 // Incomplete gamma functions follow:
275 struct upper_incomplete_gamma_fract
281 typedef std::pair<T,T> result_type;
283 upper_incomplete_gamma_fract(T a1, T z1)
284 : z(z1-a1+1), a(a1), k(0)
288 result_type operator()()
292 return result_type(k * (a - k), z);
297 inline T upper_gamma_fraction(T a, T z, T eps)
299 // Multiply result by z^a * e^-z to get the full
300 // upper incomplete integral. Divide by tgamma(z)
302 upper_incomplete_gamma_fract<T> f(a, z);
303 return 1 / (z - a + 1 + boost::math::tools::continued_fraction_a(f, eps));
307 struct lower_incomplete_gamma_series
312 typedef T result_type;
313 lower_incomplete_gamma_series(T a1, T z1) : a(a1), z(z1), result(1){}
324 template <class T, class Policy>
325 inline T lower_gamma_series(T a, T z, const Policy& pol, T init_value = 0)
327 // Multiply result by ((z^a) * (e^-z) / a) to get the full
328 // lower incomplete integral. Then divide by tgamma(a)
329 // to get the normalised value.
330 lower_incomplete_gamma_series<T> s(a, z);
331 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
332 T factor = policies::get_epsilon<T, Policy>();
333 T result = boost::math::tools::sum_series(s, factor, max_iter, init_value);
334 policies::check_series_iterations<T>("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol);
339 // Fully generic tgamma and lgamma use the incomplete partial
340 // sums added together:
342 template <class T, class Policy>
343 T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l)
345 static const char* function = "boost::math::tgamma<%1%>(%1%)";
347 if((z <= 0) && (floor(z) == z))
348 return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
351 T result = gamma_imp(T(-z), pol, l) * sinpx(z);
352 if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
353 return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
354 result = -boost::math::constants::pi<T>() / result;
356 return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
357 if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL)
358 return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", result, pol);
362 // The upper gamma fraction is *very* slow for z < 6, actually it's very
363 // slow to converge everywhere but recursing until z > 6 gets rid of the
364 // worst of it's behaviour.
372 BOOST_MATH_INSTRUMENT_CODE(prefix);
373 if((floor(z) == z) && (z < max_factorial<T>::value))
375 prefix *= unchecked_factorial<T>(itrunc(z, pol) - 1);
379 prefix = prefix * pow(z / boost::math::constants::e<T>(), z);
380 BOOST_MATH_INSTRUMENT_CODE(prefix);
381 T sum = detail::lower_gamma_series(z, z, pol) / z;
382 BOOST_MATH_INSTRUMENT_CODE(sum);
383 sum += detail::upper_gamma_fraction(z, z, ::boost::math::policies::get_epsilon<T, Policy>());
384 BOOST_MATH_INSTRUMENT_CODE(sum);
385 if(fabs(tools::max_value<T>() / prefix) < fabs(sum))
386 return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
387 BOOST_MATH_INSTRUMENT_CODE((sum * prefix));
393 template <class T, class Policy>
394 T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l, int*sign)
398 static const char* function = "boost::math::lgamma<%1%>(%1%)";
404 return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
405 T t = detail::sinpx(z);
415 result = log(boost::math::constants::pi<T>()) - lgamma_imp(z, pol, l, 0) - log(t);
417 else if((z != 1) && (z != 2))
419 T limit = (std::max)(T(z+1), T(10));
420 T prefix = z * log(limit) - limit;
421 T sum = detail::lower_gamma_series(z, limit, pol) / z;
422 sum += detail::upper_gamma_fraction(z, limit, ::boost::math::policies::get_epsilon<T, Policy>());
423 result = log(sum) + prefix;
430 // This helper calculates tgamma(dz+1)-1 without cancellation errors,
431 // used by the upper incomplete gamma with z < 1:
433 template <class T, class Policy, class Lanczos>
434 T tgammap1m1_imp(T dz, Policy const& pol, const Lanczos& l)
438 typedef typename policies::precision<T,Policy>::type precision_type;
440 typedef typename mpl::if_<
442 mpl::less_equal<precision_type, mpl::int_<0> >,
443 mpl::greater<precision_type, mpl::int_<113> >
446 is_same<Lanczos, lanczos::lanczos24m113>,
451 mpl::less_equal<precision_type, mpl::int_<64> >,
452 mpl::int_<64>, mpl::int_<113> >::type
460 // Best method is simply to subtract 1 from tgamma:
461 result = boost::math::tgamma(1+dz, pol) - 1;
462 BOOST_MATH_INSTRUMENT_CODE(result);
466 // Use expm1 on lgamma:
467 result = boost::math::expm1(-boost::math::log1p(dz, pol)
468 + lgamma_small_imp<T>(dz+2, dz + 1, dz, tag_type(), pol, l));
469 BOOST_MATH_INSTRUMENT_CODE(result);
476 // Use expm1 on lgamma:
477 result = boost::math::expm1(lgamma_small_imp<T>(dz+1, dz, dz-1, tag_type(), pol, l), pol);
478 BOOST_MATH_INSTRUMENT_CODE(result);
482 // Best method is simply to subtract 1 from tgamma:
483 result = boost::math::tgamma(1+dz, pol) - 1;
484 BOOST_MATH_INSTRUMENT_CODE(result);
491 template <class T, class Policy>
492 inline T tgammap1m1_imp(T dz, Policy const& pol,
493 const ::boost::math::lanczos::undefined_lanczos& l)
495 BOOST_MATH_STD_USING // ADL of std names
497 // There should be a better solution than this, but the
498 // algebra isn't easy for the general case....
499 // Start by subracting 1 from tgamma:
501 T result = gamma_imp(T(1 + dz), pol, l) - 1;
502 BOOST_MATH_INSTRUMENT_CODE(result);
504 // Test the level of cancellation error observed: we loose one bit
505 // for each power of 2 the result is less than 1. If we would get
506 // more bits from our most precise lgamma rational approximation,
507 // then use that instead:
509 BOOST_MATH_INSTRUMENT_CODE((dz > -0.5));
510 BOOST_MATH_INSTRUMENT_CODE((dz < 2));
511 BOOST_MATH_INSTRUMENT_CODE((ldexp(1.0, boost::math::policies::digits<T, Policy>()) * fabs(result) < 1e34));
512 if((dz > -0.5) && (dz < 2) && (ldexp(1.0, boost::math::policies::digits<T, Policy>()) * fabs(result) < 1e34))
514 result = tgammap1m1_imp(dz, pol, boost::math::lanczos::lanczos24m113());
515 BOOST_MATH_INSTRUMENT_CODE(result);
521 // Series representation for upper fraction when z is small:
524 struct small_gamma2_series
526 typedef T result_type;
528 small_gamma2_series(T a_, T x_) : result(-x_), x(-x_), apn(a_+1), n(1){}
532 T r = result / (apn);
544 // calculate power term prefix (z^a)(e^-z) used in the non-normalised
545 // incomplete gammas:
547 template <class T, class Policy>
548 T full_igamma_prefix(T a, T z, const Policy& pol)
557 if((alz < tools::log_max_value<T>()) && (-z > tools::log_min_value<T>()))
559 prefix = pow(z, a) * exp(-z);
563 prefix = pow(z / exp(z/a), a);
567 prefix = exp(alz - z);
572 if(alz > tools::log_min_value<T>())
574 prefix = pow(z, a) * exp(-z);
576 else if(z/a < tools::log_max_value<T>())
578 prefix = pow(z / exp(z/a), a);
582 prefix = exp(alz - z);
586 // This error handling isn't very good: it happens after the fact
587 // rather than before it...
589 if((boost::math::fpclassify)(prefix) == (int)FP_INFINITE)
590 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);
595 // Compute (z^a)(e^-z)/tgamma(a)
596 // most if the error occurs in this function:
598 template <class T, class Policy, class Lanczos>
599 T regularised_gamma_prefix(T a, T z, const Policy& pol, const Lanczos& l)
602 T agh = a + static_cast<T>(Lanczos::g()) - T(0.5);
604 T d = ((z - a) - static_cast<T>(Lanczos::g()) + T(0.5)) / agh;
609 // We have to treat a < 1 as a special case because our Lanczos
610 // approximations are optimised against the factorials with a > 1,
611 // and for high precision types especially (128-bit reals for example)
612 // very small values of a can give rather eroneous results for gamma
613 // unless we do this:
615 // TODO: is this still required? Lanczos approx should be better now?
617 if(z <= tools::log_min_value<T>())
619 // Oh dear, have to use logs, should be free of cancellation errors though:
620 return exp(a * log(z) - z - lgamma_imp(a, pol, l));
624 // direct calculation, no danger of overflow as gamma(a) < 1/a
626 return pow(z, a) * exp(-z) / gamma_imp(a, pol, l);
629 else if((fabs(d*d*a) <= 100) && (a > 150))
631 // special case for large a and a ~ z.
632 prefix = a * boost::math::log1pmx(d, pol) + z * static_cast<T>(0.5 - Lanczos::g()) / agh;
633 prefix = exp(prefix);
639 // direct computation is most accurate, but use various fallbacks
640 // for different parts of the problem domain:
642 T alz = a * log(z / agh);
644 if(((std::min)(alz, amz) <= tools::log_min_value<T>()) || ((std::max)(alz, amz) >= tools::log_max_value<T>()))
647 if(((std::min)(alz, amz)/2 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/2 < tools::log_max_value<T>()))
649 // compute square root of the result and then square it:
650 T sq = pow(z / agh, a / 2) * exp(amz / 2);
653 else if(((std::min)(alz, amz)/4 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/4 < tools::log_max_value<T>()) && (z > a))
655 // compute the 4th root of the result then square it twice:
656 T sq = pow(z / agh, a / 4) * exp(amz / 4);
660 else if((amza > tools::log_min_value<T>()) && (amza < tools::log_max_value<T>()))
662 prefix = pow((z * exp(amza)) / agh, a);
666 prefix = exp(alz + amz);
671 prefix = pow(z / agh, a) * exp(amz);
674 prefix *= sqrt(agh / boost::math::constants::e<T>()) / Lanczos::lanczos_sum_expG_scaled(a);
678 // And again, without Lanczos support:
680 template <class T, class Policy>
681 T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined_lanczos&)
685 T limit = (std::max)(T(10), a);
686 T sum = detail::lower_gamma_series(a, limit, pol) / a;
687 sum += detail::upper_gamma_fraction(a, limit, ::boost::math::policies::get_epsilon<T, Policy>());
691 // special case for small a:
692 T prefix = pow(z / 10, a);
696 prefix = pow((z * exp((10-z)/a)) / 10, a);
704 T alzoa = a * log(zoa);
706 if(((std::min)(alzoa, amz) <= tools::log_min_value<T>()) || ((std::max)(alzoa, amz) >= tools::log_max_value<T>()))
709 if((amza <= tools::log_min_value<T>()) || (amza >= tools::log_max_value<T>()))
711 prefix = exp(alzoa + amz);
715 prefix = pow(zoa * exp(amza), a);
720 prefix = pow(zoa, a) * exp(amz);
726 // Upper gamma fraction for very small a:
728 template <class T, class Policy>
729 inline T tgamma_small_upper_part(T a, T x, const Policy& pol, T* pgam = 0, bool invert = false, T* pderivative = 0)
731 BOOST_MATH_STD_USING // ADL of std functions.
733 // Compute the full upper fraction (Q) when a is very small:
736 result = boost::math::tgamma1pm1(a, pol);
738 *pgam = (result + 1) / a;
739 T p = boost::math::powm1(x, a, pol);
742 detail::small_gamma2_series<T> s(a, x);
743 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>() - 10;
746 *pderivative = p / (*pgam * exp(x));
747 T init_value = invert ? *pgam : 0;
748 result = -p * tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, (init_value - result) / p);
749 policies::check_series_iterations<T>("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)", max_iter, pol);
755 // Upper gamma fraction for integer a:
757 template <class T, class Policy>
758 inline T finite_gamma_q(T a, T x, Policy const& pol, T* pderivative = 0)
761 // Calculates normalised Q when a is an integer:
769 for(unsigned n = 1; n < a; ++n)
778 *pderivative = e * pow(x, a) / boost::math::unchecked_factorial<T>(itrunc(T(a - 1), pol));
783 // Upper gamma fraction for half integer a:
785 template <class T, class Policy>
786 T finite_half_gamma_q(T a, T x, T* p_derivative, const Policy& pol)
789 // Calculates normalised Q when a is a half-integer:
792 T e = boost::math::erfc(sqrt(x), pol);
793 if((e != 0) && (a > 1))
795 T term = exp(-x) / sqrt(constants::pi<T>() * x);
797 static const T half = T(1) / 2;
800 for(unsigned n = 2; n < a; ++n)
812 else if(p_derivative)
814 // We'll be dividing by x later, so calculate derivative * x:
815 *p_derivative = sqrt(x) * exp(-x) / constants::root_pi<T>();
820 // Main incomplete gamma entry point, handles all four incomplete gamma's:
822 template <class T, class Policy>
823 T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
824 const Policy& pol, T* p_derivative)
826 static const char* function = "boost::math::gamma_p<%1%>(%1%, %1%)";
828 policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
830 policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
834 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
836 T result = 0; // Just to avoid warning C4701: potentially uninitialized local variable 'result' used
838 BOOST_ASSERT((p_derivative == 0) || (normalised == true));
840 bool is_int, is_half_int;
841 bool is_small_a = (a < 30) && (a <= x + 1) && (x < tools::log_max_value<T>());
846 is_half_int = is_int ? false : (fabs(fa - a) == 0.5f);
850 is_int = is_half_int = false;
855 if(is_int && (x > 0.6))
857 // calculate Q via finite sum:
861 else if(is_half_int && (x > 0.2))
863 // calculate Q via finite sum for half integer a:
870 // Changeover criterion chosen to give a changeover at Q ~ 0.33
872 if(-0.4 / log(x) < a)
884 // Changover here occurs when P ~ 0.75 or Q ~ 0.25:
898 // Begin by testing whether we're in the "bad" zone
899 // where the result will be near 0.5 and the usual
900 // series and continued fractions are slow to converge:
902 bool use_temme = false;
903 if(normalised && std::numeric_limits<T>::is_specialized && (a > 20))
905 T sigma = fabs((x-a)/a);
906 if((a > 200) && (policies::digits<T, Policy>() <= 113))
909 // This limit is chosen so that we use Temme's expansion
910 // only if the result would be larger than about 10^-6.
911 // Below that the regular series and continued fractions
912 // converge OK, and if we use Temme's method we get increasing
913 // errors from the dominant erfc term as it's (inexact) argument
914 // increases in magnitude.
916 if(20 / a > sigma * sigma)
919 else if(policies::digits<T, Policy>() <= 64)
921 // Note in this zone we can't use Temme's expansion for
922 // types longer than an 80-bit real:
923 // it would require too many terms in the polynomials.
935 // Regular case where the result will not be too close to 0.5.
937 // Changeover here occurs at P ~ Q ~ 0.5
938 // Note that series computation of P is about x2 faster than continued fraction
939 // calculation of Q, so try and use the CF only when really necessary, especially
942 if(x - (1 / (3 * x)) < a)
958 result = finite_gamma_q(a, x, pol, p_derivative);
959 if(normalised == false)
960 result *= boost::math::tgamma(a, pol);
965 result = finite_half_gamma_q(a, x, p_derivative, pol);
966 if(normalised == false)
967 result *= boost::math::tgamma(a, pol);
968 if(p_derivative && (*p_derivative == 0))
969 *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
975 result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
977 *p_derivative = result;
983 init_value = -a * (normalised ? 1 : boost::math::tgamma(a, pol)) / result;
985 result *= detail::lower_gamma_series(a, x, pol, init_value) / a;
999 result = tgamma_small_upper_part(a, x, pol, &g, invert, p_derivative);
1008 result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
1010 *p_derivative = result;
1012 result *= upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>());
1018 // Use compile time dispatch to the appropriate
1019 // Temme asymptotic expansion. This may be dead code
1020 // if T does not have numeric limits support, or has
1021 // too many digits for the most precise version of
1022 // these expansions, in that case we'll be calling
1023 // an empty function.
1025 typedef typename policies::precision<T, Policy>::type precision_type;
1027 typedef typename mpl::if_<
1028 mpl::or_<mpl::equal_to<precision_type, mpl::int_<0> >,
1029 mpl::greater<precision_type, mpl::int_<113> > >,
1032 mpl::less_equal<precision_type, mpl::int_<53> >,
1035 mpl::less_equal<precision_type, mpl::int_<64> >,
1042 result = igamma_temme_large(a, x, pol, static_cast<tag_type const*>(0));
1046 *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
1051 if(normalised && (result > 1))
1055 T gam = normalised ? 1 : boost::math::tgamma(a, pol);
1056 result = gam - result;
1061 // Need to convert prefix term to derivative:
1063 if((x < 1) && (tools::max_value<T>() * x < *p_derivative))
1065 // overflow, just return an arbitrarily large value:
1066 *p_derivative = tools::max_value<T>() / 2;
1076 // Ratios of two gamma functions:
1078 template <class T, class Policy, class Lanczos>
1079 T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const Lanczos&)
1081 BOOST_MATH_STD_USING
1082 T zgh = z + Lanczos::g() - constants::half<T>();
1084 if(fabs(delta) < 10)
1086 result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol));
1090 result = pow(zgh / (zgh + delta), z - constants::half<T>());
1092 result *= pow(constants::e<T>() / (zgh + delta), delta);
1093 result *= Lanczos::lanczos_sum(z) / Lanczos::lanczos_sum(T(z + delta));
1097 // And again without Lanczos support this time:
1099 template <class T, class Policy>
1100 T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const lanczos::undefined_lanczos&)
1102 BOOST_MATH_STD_USING
1104 // The upper gamma fraction is *very* slow for z < 6, actually it's very
1105 // slow to converge everywhere but recursing until z > 6 gets rid of the
1106 // worst of it's behaviour.
1110 while((zd < 6) && (z < 6))
1119 prefix *= exp(-z * boost::math::log1p(delta / z, pol));
1123 prefix *= pow(z / zd, z);
1125 prefix *= pow(constants::e<T>() / zd, delta);
1126 T sum = detail::lower_gamma_series(z, z, pol) / z;
1127 sum += detail::upper_gamma_fraction(z, z, ::boost::math::policies::get_epsilon<T, Policy>());
1128 T sumd = detail::lower_gamma_series(zd, zd, pol) / zd;
1129 sumd += detail::upper_gamma_fraction(zd, zd, ::boost::math::policies::get_epsilon<T, Policy>());
1131 if(fabs(tools::max_value<T>() / prefix) < fabs(sum))
1132 return policies::raise_overflow_error<T>("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Result of tgamma is too large to represent.", pol);
1133 return sum * prefix;
1136 template <class T, class Policy>
1137 T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol)
1139 BOOST_MATH_STD_USING
1142 policies::raise_domain_error<T>("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got a=%1%).", z, pol);
1144 policies::raise_domain_error<T>("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got b=%1%).", z+delta, pol);
1146 if(floor(delta) == delta)
1151 // Both z and delta are integers, see if we can just use table lookup
1152 // of the factorials to get the result:
1154 if((z <= max_factorial<T>::value) && (z + delta <= max_factorial<T>::value))
1156 return unchecked_factorial<T>((unsigned)itrunc(z, pol) - 1) / unchecked_factorial<T>((unsigned)itrunc(T(z + delta), pol) - 1);
1159 if(fabs(delta) < 20)
1162 // delta is a small integer, we can use a finite product:
1170 while(0 != (delta += 1))
1180 while(0 != (delta -= 1))
1189 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1190 return tgamma_delta_ratio_imp_lanczos(z, delta, pol, lanczos_type());
1193 template <class T, class Policy>
1194 T tgamma_ratio_imp(T x, T y, const Policy& pol)
1196 BOOST_MATH_STD_USING
1198 if((x <= tools::min_value<T>()) || (boost::math::isinf)(x))
1199 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);
1200 if((y <= tools::min_value<T>()) || (boost::math::isinf)(y))
1201 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);
1203 if((x < max_factorial<T>::value) && (y < max_factorial<T>::value))
1205 // Rather than subtracting values, lets just call the gamma functions directly:
1206 return boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
1211 if(y < 2 * max_factorial<T>::value)
1213 // We need to sidestep on x as well, otherwise we'll underflow
1214 // before we get to factor in the prefix term:
1217 while(y >= max_factorial<T>::value)
1222 return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
1225 // result is almost certainly going to underflow to zero, try logs just in case:
1227 return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
1231 if(x < 2 * max_factorial<T>::value)
1233 // We need to sidestep on y as well, otherwise we'll overflow
1234 // before we get to factor in the prefix term:
1237 while(x >= max_factorial<T>::value)
1242 return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
1245 // Result will almost certainly overflow, try logs just in case:
1247 return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
1250 // Regular case, x and y both large and similar in magnitude:
1252 return boost::math::tgamma_delta_ratio(x, y - x, pol);
1255 template <class T, class Policy>
1256 T gamma_p_derivative_imp(T a, T x, const Policy& pol)
1259 // Usual error checks first:
1262 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);
1264 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);
1266 // Now special cases:
1270 return (a > 1) ? 0 :
1271 (a == 1) ? 1 : policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
1276 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1277 T f1 = detail::regularised_gamma_prefix(a, x, pol, lanczos_type());
1278 if((x < 1) && (tools::max_value<T>() * x < f1))
1281 return policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
1289 template <class T, class Policy>
1290 inline typename tools::promote_args<T>::type
1291 tgamma(T z, const Policy& /* pol */, const mpl::true_)
1293 BOOST_FPU_EXCEPTION_GUARD
1294 typedef typename tools::promote_args<T>::type result_type;
1295 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1296 typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1297 typedef typename policies::normalise<
1299 policies::promote_float<false>,
1300 policies::promote_double<false>,
1301 policies::discrete_quantile<>,
1302 policies::assert_undefined<> >::type forwarding_policy;
1303 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%)");
1306 template <class T, class Policy>
1307 struct igamma_initializer
1313 typedef typename policies::precision<T, Policy>::type precision_type;
1315 typedef typename mpl::if_<
1316 mpl::or_<mpl::equal_to<precision_type, mpl::int_<0> >,
1317 mpl::greater<precision_type, mpl::int_<113> > >,
1320 mpl::less_equal<precision_type, mpl::int_<53> >,
1323 mpl::less_equal<precision_type, mpl::int_<64> >,
1330 do_init(tag_type());
1333 static void do_init(const mpl::int_<N>&)
1335 boost::math::gamma_p(static_cast<T>(400), static_cast<T>(400), Policy());
1337 static void do_init(const mpl::int_<53>&){}
1338 void force_instantiate()const{}
1340 static const init initializer;
1341 static void force_instantiate()
1343 initializer.force_instantiate();
1347 template <class T, class Policy>
1348 const typename igamma_initializer<T, Policy>::init igamma_initializer<T, Policy>::initializer;
1350 template <class T, class Policy>
1351 struct lgamma_initializer
1357 typedef typename policies::precision<T, Policy>::type precision_type;
1358 typedef typename mpl::if_<
1360 mpl::less_equal<precision_type, mpl::int_<64> >,
1361 mpl::greater<precision_type, mpl::int_<0> >
1366 mpl::less_equal<precision_type, mpl::int_<113> >,
1367 mpl::greater<precision_type, mpl::int_<0> >
1369 mpl::int_<113>, mpl::int_<0> >::type
1371 do_init(tag_type());
1373 static void do_init(const mpl::int_<64>&)
1375 boost::math::lgamma(static_cast<T>(2.5), Policy());
1376 boost::math::lgamma(static_cast<T>(1.25), Policy());
1377 boost::math::lgamma(static_cast<T>(1.75), Policy());
1379 static void do_init(const mpl::int_<113>&)
1381 boost::math::lgamma(static_cast<T>(2.5), Policy());
1382 boost::math::lgamma(static_cast<T>(1.25), Policy());
1383 boost::math::lgamma(static_cast<T>(1.5), Policy());
1384 boost::math::lgamma(static_cast<T>(1.75), Policy());
1386 static void do_init(const mpl::int_<0>&)
1389 void force_instantiate()const{}
1391 static const init initializer;
1392 static void force_instantiate()
1394 initializer.force_instantiate();
1398 template <class T, class Policy>
1399 const typename lgamma_initializer<T, Policy>::init lgamma_initializer<T, Policy>::initializer;
1401 template <class T1, class T2, class Policy>
1402 inline typename tools::promote_args<T1, T2>::type
1403 tgamma(T1 a, T2 z, const Policy&, const mpl::false_)
1405 BOOST_FPU_EXCEPTION_GUARD
1406 typedef typename tools::promote_args<T1, T2>::type result_type;
1407 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1408 // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1409 typedef typename policies::normalise<
1411 policies::promote_float<false>,
1412 policies::promote_double<false>,
1413 policies::discrete_quantile<>,
1414 policies::assert_undefined<> >::type forwarding_policy;
1416 igamma_initializer<value_type, forwarding_policy>::force_instantiate();
1418 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
1419 detail::gamma_incomplete_imp(static_cast<value_type>(a),
1420 static_cast<value_type>(z), false, true,
1421 forwarding_policy(), static_cast<value_type*>(0)), "boost::math::tgamma<%1%>(%1%, %1%)");
1424 template <class T1, class T2>
1425 inline typename tools::promote_args<T1, T2>::type
1426 tgamma(T1 a, T2 z, const mpl::false_ tag)
1428 return tgamma(a, z, policies::policy<>(), tag);
1432 } // namespace detail
1435 inline typename tools::promote_args<T>::type
1438 return tgamma(z, policies::policy<>());
1441 template <class T, class Policy>
1442 inline typename tools::promote_args<T>::type
1443 lgamma(T z, int* sign, const Policy&)
1445 BOOST_FPU_EXCEPTION_GUARD
1446 typedef typename tools::promote_args<T>::type result_type;
1447 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1448 typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1449 typedef typename policies::normalise<
1451 policies::promote_float<false>,
1452 policies::promote_double<false>,
1453 policies::discrete_quantile<>,
1454 policies::assert_undefined<> >::type forwarding_policy;
1456 detail::lgamma_initializer<value_type, forwarding_policy>::force_instantiate();
1458 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%)");
1462 inline typename tools::promote_args<T>::type
1463 lgamma(T z, int* sign)
1465 return lgamma(z, sign, policies::policy<>());
1468 template <class T, class Policy>
1469 inline typename tools::promote_args<T>::type
1470 lgamma(T x, const Policy& pol)
1472 return ::boost::math::lgamma(x, 0, pol);
1476 inline typename tools::promote_args<T>::type
1479 return ::boost::math::lgamma(x, 0, policies::policy<>());
1482 template <class T, class Policy>
1483 inline typename tools::promote_args<T>::type
1484 tgamma1pm1(T z, const Policy& /* pol */)
1486 BOOST_FPU_EXCEPTION_GUARD
1487 typedef typename tools::promote_args<T>::type result_type;
1488 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1489 typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1490 typedef typename policies::normalise<
1492 policies::promote_float<false>,
1493 policies::promote_double<false>,
1494 policies::discrete_quantile<>,
1495 policies::assert_undefined<> >::type forwarding_policy;
1497 return policies::checked_narrowing_cast<typename remove_cv<result_type>::type, forwarding_policy>(detail::tgammap1m1_imp(static_cast<value_type>(z), forwarding_policy(), evaluation_type()), "boost::math::tgamma1pm1<%!%>(%1%)");
1501 inline typename tools::promote_args<T>::type
1504 return tgamma1pm1(z, policies::policy<>());
1508 // Full upper incomplete gamma:
1510 template <class T1, class T2>
1511 inline typename tools::promote_args<T1, T2>::type
1515 // Type T2 could be a policy object, or a value, select the
1516 // right overload based on T2:
1518 typedef typename policies::is_policy<T2>::type maybe_policy;
1519 return detail::tgamma(a, z, maybe_policy());
1521 template <class T1, class T2, class Policy>
1522 inline typename tools::promote_args<T1, T2>::type
1523 tgamma(T1 a, T2 z, const Policy& pol)
1525 return detail::tgamma(a, z, pol, mpl::false_());
1528 // Full lower incomplete gamma:
1530 template <class T1, class T2, class Policy>
1531 inline typename tools::promote_args<T1, T2>::type
1532 tgamma_lower(T1 a, T2 z, const Policy&)
1534 BOOST_FPU_EXCEPTION_GUARD
1535 typedef typename tools::promote_args<T1, T2>::type result_type;
1536 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1537 // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1538 typedef typename policies::normalise<
1540 policies::promote_float<false>,
1541 policies::promote_double<false>,
1542 policies::discrete_quantile<>,
1543 policies::assert_undefined<> >::type forwarding_policy;
1545 detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
1547 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
1548 detail::gamma_incomplete_imp(static_cast<value_type>(a),
1549 static_cast<value_type>(z), false, false,
1550 forwarding_policy(), static_cast<value_type*>(0)), "tgamma_lower<%1%>(%1%, %1%)");
1552 template <class T1, class T2>
1553 inline typename tools::promote_args<T1, T2>::type
1554 tgamma_lower(T1 a, T2 z)
1556 return tgamma_lower(a, z, policies::policy<>());
1559 // Regularised upper incomplete gamma:
1561 template <class T1, class T2, class Policy>
1562 inline typename tools::promote_args<T1, T2>::type
1563 gamma_q(T1 a, T2 z, const Policy& /* pol */)
1565 BOOST_FPU_EXCEPTION_GUARD
1566 typedef typename tools::promote_args<T1, T2>::type result_type;
1567 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1568 // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1569 typedef typename policies::normalise<
1571 policies::promote_float<false>,
1572 policies::promote_double<false>,
1573 policies::discrete_quantile<>,
1574 policies::assert_undefined<> >::type forwarding_policy;
1576 detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
1578 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
1579 detail::gamma_incomplete_imp(static_cast<value_type>(a),
1580 static_cast<value_type>(z), true, true,
1581 forwarding_policy(), static_cast<value_type*>(0)), "gamma_q<%1%>(%1%, %1%)");
1583 template <class T1, class T2>
1584 inline typename tools::promote_args<T1, T2>::type
1587 return gamma_q(a, z, policies::policy<>());
1590 // Regularised lower incomplete gamma:
1592 template <class T1, class T2, class Policy>
1593 inline typename tools::promote_args<T1, T2>::type
1594 gamma_p(T1 a, T2 z, const Policy&)
1596 BOOST_FPU_EXCEPTION_GUARD
1597 typedef typename tools::promote_args<T1, T2>::type result_type;
1598 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1599 // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1600 typedef typename policies::normalise<
1602 policies::promote_float<false>,
1603 policies::promote_double<false>,
1604 policies::discrete_quantile<>,
1605 policies::assert_undefined<> >::type forwarding_policy;
1607 detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
1609 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
1610 detail::gamma_incomplete_imp(static_cast<value_type>(a),
1611 static_cast<value_type>(z), true, false,
1612 forwarding_policy(), static_cast<value_type*>(0)), "gamma_p<%1%>(%1%, %1%)");
1614 template <class T1, class T2>
1615 inline typename tools::promote_args<T1, T2>::type
1618 return gamma_p(a, z, policies::policy<>());
1621 // ratios of gamma functions:
1622 template <class T1, class T2, class Policy>
1623 inline typename tools::promote_args<T1, T2>::type
1624 tgamma_delta_ratio(T1 z, T2 delta, const Policy& /* pol */)
1626 BOOST_FPU_EXCEPTION_GUARD
1627 typedef typename tools::promote_args<T1, T2>::type result_type;
1628 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1629 typedef typename policies::normalise<
1631 policies::promote_float<false>,
1632 policies::promote_double<false>,
1633 policies::discrete_quantile<>,
1634 policies::assert_undefined<> >::type forwarding_policy;
1636 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%)");
1638 template <class T1, class T2>
1639 inline typename tools::promote_args<T1, T2>::type
1640 tgamma_delta_ratio(T1 z, T2 delta)
1642 return tgamma_delta_ratio(z, delta, policies::policy<>());
1644 template <class T1, class T2, class Policy>
1645 inline typename tools::promote_args<T1, T2>::type
1646 tgamma_ratio(T1 a, T2 b, const Policy&)
1648 typedef typename tools::promote_args<T1, T2>::type result_type;
1649 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1650 typedef typename policies::normalise<
1652 policies::promote_float<false>,
1653 policies::promote_double<false>,
1654 policies::discrete_quantile<>,
1655 policies::assert_undefined<> >::type forwarding_policy;
1657 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%)");
1659 template <class T1, class T2>
1660 inline typename tools::promote_args<T1, T2>::type
1661 tgamma_ratio(T1 a, T2 b)
1663 return tgamma_ratio(a, b, policies::policy<>());
1666 template <class T1, class T2, class Policy>
1667 inline typename tools::promote_args<T1, T2>::type
1668 gamma_p_derivative(T1 a, T2 x, const Policy&)
1670 BOOST_FPU_EXCEPTION_GUARD
1671 typedef typename tools::promote_args<T1, T2>::type result_type;
1672 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1673 typedef typename policies::normalise<
1675 policies::promote_float<false>,
1676 policies::promote_double<false>,
1677 policies::discrete_quantile<>,
1678 policies::assert_undefined<> >::type forwarding_policy;
1680 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%)");
1682 template <class T1, class T2>
1683 inline typename tools::promote_args<T1, T2>::type
1684 gamma_p_derivative(T1 a, T2 x)
1686 return gamma_p_derivative(a, x, policies::policy<>());
1690 } // namespace boost
1693 # pragma warning(pop)
1696 #include <boost/math/special_functions/detail/igamma_inverse.hpp>
1697 #include <boost/math/special_functions/detail/gamma_inva.hpp>
1698 #include <boost/math/special_functions/erf.hpp>
1700 #endif // BOOST_MATH_SF_GAMMA_HPP