2 // Copyright John Maddock 2006-7, 2013-14.
3 // Copyright Paul A. Bristow 2007, 2013-14.
4 // Copyright Nikhar Agrawal 2013-14
5 // Copyright Christopher Kormanyos 2013-14
7 // Use, modification and distribution are subject to the
8 // Boost Software License, Version 1.0. (See accompanying file
9 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
11 #ifndef BOOST_MATH_SF_GAMMA_HPP
12 #define BOOST_MATH_SF_GAMMA_HPP
18 #include <boost/config.hpp>
19 #include <boost/math/tools/series.hpp>
20 #include <boost/math/tools/fraction.hpp>
21 #include <boost/math/tools/precision.hpp>
22 #include <boost/math/tools/promotion.hpp>
23 #include <boost/math/policies/error_handling.hpp>
24 #include <boost/math/constants/constants.hpp>
25 #include <boost/math/special_functions/math_fwd.hpp>
26 #include <boost/math/special_functions/log1p.hpp>
27 #include <boost/math/special_functions/trunc.hpp>
28 #include <boost/math/special_functions/powm1.hpp>
29 #include <boost/math/special_functions/sqrt1pm1.hpp>
30 #include <boost/math/special_functions/lanczos.hpp>
31 #include <boost/math/special_functions/fpclassify.hpp>
32 #include <boost/math/special_functions/detail/igamma_large.hpp>
33 #include <boost/math/special_functions/detail/unchecked_factorial.hpp>
34 #include <boost/math/special_functions/detail/lgamma_small.hpp>
35 #include <boost/math/special_functions/bernoulli.hpp>
36 #include <boost/math/special_functions/zeta.hpp>
37 #include <boost/type_traits/is_convertible.hpp>
38 #include <boost/assert.hpp>
39 #include <boost/mpl/greater.hpp>
40 #include <boost/mpl/equal_to.hpp>
41 #include <boost/mpl/greater.hpp>
43 #include <boost/config/no_tr1/cmath.hpp>
47 # pragma warning(push)
48 # pragma warning(disable: 4702) // unreachable code (return after domain_error throw).
49 # pragma warning(disable: 4127) // conditional expression is constant.
50 # pragma warning(disable: 4100) // unreferenced formal parameter.
51 // Several variables made comments,
52 // but some difficulty as whether referenced on not may depend on macro values.
53 // So to be safe, 4100 warnings suppressed.
54 // TODO - revisit this?
57 namespace boost{ namespace math{
62 inline bool is_odd(T v, const boost::true_type&)
64 int i = static_cast<int>(v);
68 inline bool is_odd(T v, const boost::false_type&)
70 // Oh dear can't cast T to int!
72 T modulus = v - 2 * floor(v/2);
73 return static_cast<bool>(modulus != 0);
76 inline bool is_odd(T v)
78 return is_odd(v, ::boost::is_convertible<T, int>());
84 // Ad hoc function calculates x * sin(pi * x),
85 // 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 -boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
140 result = -boost::math::constants::pi<T>() / result;
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);
162 else if (z < tools::root_epsilon<T>())
164 if (z < 1 / tools::max_value<T>())
165 result = policies::raise_overflow_error<T>(function, 0, pol);
166 result *= 1 / z - constants::euler<T>();
170 result *= Lanczos::lanczos_sum(z);
171 T zgh = (z + static_cast<T>(Lanczos::g()) - boost::math::constants::half<T>());
173 BOOST_MATH_INSTRUMENT_VARIABLE(result);
174 BOOST_MATH_INSTRUMENT_VARIABLE(tools::log_max_value<T>());
175 if(z * lzgh > tools::log_max_value<T>())
177 // we're going to overflow unless this is done with care:
178 BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
179 if(lzgh * z / 2 > tools::log_max_value<T>())
180 return boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
181 T hp = pow(zgh, (z / 2) - T(0.25));
182 BOOST_MATH_INSTRUMENT_VARIABLE(hp);
183 result *= hp / exp(zgh);
184 BOOST_MATH_INSTRUMENT_VARIABLE(result);
185 if(tools::max_value<T>() / hp < result)
186 return boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
188 BOOST_MATH_INSTRUMENT_VARIABLE(result);
192 BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
193 BOOST_MATH_INSTRUMENT_VARIABLE(pow(zgh, z - boost::math::constants::half<T>()));
194 BOOST_MATH_INSTRUMENT_VARIABLE(exp(zgh));
195 result *= pow(zgh, z - boost::math::constants::half<T>()) / exp(zgh);
196 BOOST_MATH_INSTRUMENT_VARIABLE(result);
202 // lgamma(z) with Lanczos support:
204 template <class T, class Policy, class Lanczos>
205 T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = 0)
207 #ifdef BOOST_MATH_INSTRUMENT
208 static bool b = false;
211 std::cout << "lgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;
218 static const char* function = "boost::math::lgamma<%1%>(%1%)";
222 if(z <= -tools::root_epsilon<T>())
224 // reflection formula:
226 return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%.", z, pol);
238 result = log(boost::math::constants::pi<T>()) - lgamma_imp(z, pol, l) - log(t);
240 else if (z < tools::root_epsilon<T>())
243 return policies::raise_pole_error<T>(function, "Evaluation of lgamma at %1%.", z, pol);
244 if (fabs(z) < 1 / tools::max_value<T>())
245 result = -log(fabs(z));
247 result = log(fabs(1 / z - constants::euler<T>()));
253 typedef typename policies::precision<T, Policy>::type precision_type;
254 typedef typename mpl::if_<
256 mpl::less_equal<precision_type, mpl::int_<64> >,
257 mpl::greater<precision_type, mpl::int_<0> >
262 mpl::less_equal<precision_type, mpl::int_<113> >,
263 mpl::greater<precision_type, mpl::int_<0> >
265 mpl::int_<113>, mpl::int_<0> >::type
267 result = lgamma_small_imp<T>(z, T(z - 1), T(z - 2), tag_type(), pol, l);
269 else if((z >= 3) && (z < 100) && (std::numeric_limits<T>::max_exponent >= 1024))
271 // taking the log of tgamma reduces the error, no danger of overflow here:
272 result = log(gamma_imp(z, pol, l));
276 // regular evaluation:
277 T zgh = static_cast<T>(z + Lanczos::g() - boost::math::constants::half<T>());
278 result = log(zgh) - 1;
280 result += log(Lanczos::lanczos_sum_expG_scaled(z));
289 // Incomplete gamma functions follow:
292 struct upper_incomplete_gamma_fract
298 typedef std::pair<T,T> result_type;
300 upper_incomplete_gamma_fract(T a1, T z1)
301 : z(z1-a1+1), a(a1), k(0)
305 result_type operator()()
309 return result_type(k * (a - k), z);
314 inline T upper_gamma_fraction(T a, T z, T eps)
316 // Multiply result by z^a * e^-z to get the full
317 // upper incomplete integral. Divide by tgamma(z)
319 upper_incomplete_gamma_fract<T> f(a, z);
320 return 1 / (z - a + 1 + boost::math::tools::continued_fraction_a(f, eps));
324 struct lower_incomplete_gamma_series
329 typedef T result_type;
330 lower_incomplete_gamma_series(T a1, T z1) : a(a1), z(z1), result(1){}
341 template <class T, class Policy>
342 inline T lower_gamma_series(T a, T z, const Policy& pol, T init_value = 0)
344 // Multiply result by ((z^a) * (e^-z) / a) to get the full
345 // lower incomplete integral. Then divide by tgamma(a)
346 // to get the normalised value.
347 lower_incomplete_gamma_series<T> s(a, z);
348 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
349 T factor = policies::get_epsilon<T, Policy>();
350 T result = boost::math::tools::sum_series(s, factor, max_iter, init_value);
351 policies::check_series_iterations<T>("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol);
356 // Fully generic tgamma and lgamma use Stirling's approximation
357 // with Bernoulli numbers.
360 std::size_t highest_bernoulli_index()
362 const float digits10_of_type = (std::numeric_limits<T>::is_specialized
363 ? static_cast<float>(std::numeric_limits<T>::digits10)
364 : static_cast<float>(boost::math::tools::digits<T>() * 0.301F));
366 // Find the high index n for Bn to produce the desired precision in Stirling's calculation.
367 return static_cast<std::size_t>(18.0F + (0.6F * digits10_of_type));
371 T minimum_argument_for_bernoulli_recursion()
373 const float digits10_of_type = (std::numeric_limits<T>::is_specialized
374 ? static_cast<float>(std::numeric_limits<T>::digits10)
375 : static_cast<float>(boost::math::tools::digits<T>() * 0.301F));
377 return T(digits10_of_type * 1.7F);
380 // Forward declaration of the lgamma_imp template specialization.
381 template <class T, class Policy>
382 T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign = 0);
384 template <class T, class Policy>
385 T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&)
389 static const char* function = "boost::math::tgamma<%1%>(%1%)";
391 // Check if the argument of tgamma is identically zero.
392 const bool is_at_zero = (z == 0);
394 if((is_at_zero) || ((boost::math::isinf)(z) && (z < 0)))
395 return policies::raise_domain_error<T>(function, "Evaluation of tgamma at %1%.", z, pol);
397 const bool b_neg = (z < 0);
399 const bool floor_of_z_is_equal_to_z = (floor(z) == z);
401 // Special case handling of small factorials:
402 if((!b_neg) && floor_of_z_is_equal_to_z && (z < boost::math::max_factorial<T>::value))
404 return boost::math::unchecked_factorial<T>(itrunc(z) - 1);
407 // Make a local, unsigned copy of the input argument.
408 T zz((!b_neg) ? z : -z);
410 // Special case for ultra-small z:
411 if(zz < tools::cbrt_epsilon<T>())
414 const T a1(boost::math::constants::euler<T>());
415 const T six_euler_squared((boost::math::constants::euler<T>() * boost::math::constants::euler<T>()) * 6);
416 const T a2((six_euler_squared - boost::math::constants::pi_sqr<T>()) / 12);
418 const T inverse_tgamma_series = z * ((a2 * z + a1) * z + a0);
420 return 1 / inverse_tgamma_series;
423 // Scale the argument up for the calculation of lgamma,
424 // and use downward recursion later for the final result.
425 const T min_arg_for_recursion = minimum_argument_for_bernoulli_recursion<T>();
429 if(zz < min_arg_for_recursion)
431 n_recur = boost::math::itrunc(min_arg_for_recursion - zz) + 1;
440 const T log_gamma_value = lgamma_imp(zz, pol, lanczos::undefined_lanczos());
442 if(log_gamma_value > tools::log_max_value<T>())
443 return policies::raise_overflow_error<T>(function, 0, pol);
445 T gamma_value = exp(log_gamma_value);
447 // Rescale the result using downward recursion if necessary.
450 // The order of divides is important, if we keep subtracting 1 from zz
451 // we DO NOT get back to z (cancellation error). Further if z < epsilon
452 // we would end up dividing by zero. Also in order to prevent spurious
453 // overflow with the first division, we must save dividing by |z| till last,
454 // so the optimal order of divides is z+1, z+2, z+3...z+n_recur-1,z.
456 for(int k = 1; k < n_recur; ++k)
461 gamma_value /= fabs(z);
464 // Return the result, accounting for possible negative arguments.
467 // Provide special error analysis for:
468 // * arguments in the neighborhood of a negative integer
469 // * arguments exactly equal to a negative integer.
471 // Check if the argument of tgamma is exactly equal to a negative integer.
472 if(floor_of_z_is_equal_to_z)
473 return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
475 gamma_value *= sinpx(z);
477 BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value);
479 const bool result_is_too_large_to_represent = ( (abs(gamma_value) < 1)
480 && ((tools::max_value<T>() * abs(gamma_value)) < boost::math::constants::pi<T>()));
482 if(result_is_too_large_to_represent)
483 return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
485 gamma_value = -boost::math::constants::pi<T>() / gamma_value;
486 BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value);
489 return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
491 if((boost::math::fpclassify)(gamma_value) == static_cast<int>(FP_SUBNORMAL))
492 return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", gamma_value, pol);
498 template <class T, class Policy>
499 inline T log_gamma_near_1(const T& z, Policy const& pol)
502 // This is for the multiprecision case where there is
503 // no lanczos support...
505 BOOST_MATH_STD_USING // ADL of std names
507 BOOST_ASSERT(fabs(z) < 1);
509 T result = -constants::euler<T>() * z;
511 T power_term = z * z;
517 term = boost::math::zeta<T>(j + 2, pol) * power_term / (j + 2);
524 } while(fabs(result) * tools::epsilon<T>() < fabs(term));
529 template <class T, class Policy>
530 T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign)
534 static const char* function = "boost::math::lgamma<%1%>(%1%)";
536 // Check if the argument of lgamma is identically zero.
537 const bool is_at_zero = (z == 0);
540 return policies::raise_domain_error<T>(function, "Evaluation of lgamma at zero %1%.", z, pol);
541 if((boost::math::isnan)(z))
542 return policies::raise_domain_error<T>(function, "Evaluation of lgamma at %1%.", z, pol);
543 if((boost::math::isinf)(z))
544 return policies::raise_overflow_error<T>(function, 0, pol);
546 const bool b_neg = (z < 0);
548 const bool floor_of_z_is_equal_to_z = (floor(z) == z);
550 // Special case handling of small factorials:
551 if((!b_neg) && floor_of_z_is_equal_to_z && (z < boost::math::max_factorial<T>::value))
553 return log(boost::math::unchecked_factorial<T>(itrunc(z) - 1));
556 // Make a local, unsigned copy of the input argument.
557 T zz((!b_neg) ? z : -z);
559 const T min_arg_for_recursion = minimum_argument_for_bernoulli_recursion<T>();
563 if (zz < min_arg_for_recursion)
565 // Here we simply take the logarithm of tgamma(). This is somewhat
566 // inefficient, but simple. The rationale is that the argument here
567 // is relatively small and overflow is not expected to be likely.
568 if(fabs(z - 1) < 0.25)
570 return log_gamma_near_1(T(zz - 1), pol);
572 else if(fabs(z - 2) < 0.25)
574 return log_gamma_near_1(T(zz - 2), pol) + log(zz - 1);
576 else if (z > -tools::root_epsilon<T>())
578 // Reflection formula may fail if z is very close to zero, let the series
579 // expansion for tgamma close to zero do the work:
580 log_gamma_value = log(abs(gamma_imp(z, pol, lanczos::undefined_lanczos())));
583 *sign = z < 0 ? -1 : 1;
585 return log_gamma_value;
589 // No issue with spurious overflow in reflection formula,
590 // just fall through to regular code:
591 log_gamma_value = log(abs(gamma_imp(zz, pol, lanczos::undefined_lanczos())));
596 // Perform the Bernoulli series expansion of Stirling's approximation.
598 const std::size_t number_of_bernoullis_b2n = highest_bernoulli_index<T>();
600 T one_over_x_pow_two_n_minus_one = 1 / zz;
601 const T one_over_x2 = one_over_x_pow_two_n_minus_one * one_over_x_pow_two_n_minus_one;
602 T sum = (boost::math::bernoulli_b2n<T>(1) / 2) * one_over_x_pow_two_n_minus_one;
603 const T target_epsilon_to_break_loop = (sum * boost::math::tools::epsilon<T>()) * T(1.0E-10F);
605 for(std::size_t n = 2U; n < number_of_bernoullis_b2n; ++n)
607 one_over_x_pow_two_n_minus_one *= one_over_x2;
609 const std::size_t n2 = static_cast<std::size_t>(n * 2U);
611 const T term = (boost::math::bernoulli_b2n<T>(static_cast<int>(n)) * one_over_x_pow_two_n_minus_one) / (n2 * (n2 - 1U));
613 if((n >= 8U) && (abs(term) < target_epsilon_to_break_loop))
615 // We have reached the desired precision in Stirling's expansion.
616 // Adding additional terms to the sum of this divergent asymptotic
617 // expansion will not improve the result.
619 // Break from the loop.
626 // Complete Stirling's approximation.
627 const T half_ln_two_pi = log(boost::math::constants::two_pi<T>()) / 2;
629 log_gamma_value = ((((zz - boost::math::constants::half<T>()) * log(zz)) - zz) + half_ln_two_pi) + sum;
632 int sign_of_result = 1;
636 // Provide special error analysis if the argument is exactly
637 // equal to a negative integer.
639 // Check if the argument of lgamma is exactly equal to a negative integer.
640 if(floor_of_z_is_equal_to_z)
641 return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%.", z, pol);
651 sign_of_result = -sign_of_result;
654 log_gamma_value = - log_gamma_value
655 + log(boost::math::constants::pi<T>())
659 if(sign != static_cast<int*>(0U)) { *sign = sign_of_result; }
661 return log_gamma_value;
665 // This helper calculates tgamma(dz+1)-1 without cancellation errors,
666 // used by the upper incomplete gamma with z < 1:
668 template <class T, class Policy, class Lanczos>
669 T tgammap1m1_imp(T dz, Policy const& pol, const Lanczos& l)
673 typedef typename policies::precision<T,Policy>::type precision_type;
675 typedef typename mpl::if_<
677 mpl::less_equal<precision_type, mpl::int_<0> >,
678 mpl::greater<precision_type, mpl::int_<113> >
681 mpl::and_<is_same<Lanczos, lanczos::lanczos24m113>, mpl::greater<precision_type, mpl::int_<0> > >,
686 mpl::less_equal<precision_type, mpl::int_<64> >,
687 mpl::int_<64>, mpl::int_<113> >::type
695 // Best method is simply to subtract 1 from tgamma:
696 result = boost::math::tgamma(1+dz, pol) - 1;
697 BOOST_MATH_INSTRUMENT_CODE(result);
701 // Use expm1 on lgamma:
702 result = boost::math::expm1(-boost::math::log1p(dz, pol)
703 + lgamma_small_imp<T>(dz+2, dz + 1, dz, tag_type(), pol, l));
704 BOOST_MATH_INSTRUMENT_CODE(result);
711 // Use expm1 on lgamma:
712 result = boost::math::expm1(lgamma_small_imp<T>(dz+1, dz, dz-1, tag_type(), pol, l), pol);
713 BOOST_MATH_INSTRUMENT_CODE(result);
717 // Best method is simply to subtract 1 from tgamma:
718 result = boost::math::tgamma(1+dz, pol) - 1;
719 BOOST_MATH_INSTRUMENT_CODE(result);
726 template <class T, class Policy>
727 inline T tgammap1m1_imp(T z, Policy const& pol,
728 const ::boost::math::lanczos::undefined_lanczos&)
730 BOOST_MATH_STD_USING // ADL of std names
734 return boost::math::expm1(log_gamma_near_1(z, pol));
736 return boost::math::expm1(boost::math::lgamma(1 + z, pol));
740 // Series representation for upper fraction when z is small:
743 struct small_gamma2_series
745 typedef T result_type;
747 small_gamma2_series(T a_, T x_) : result(-x_), x(-x_), apn(a_+1), n(1){}
751 T r = result / (apn);
763 // calculate power term prefix (z^a)(e^-z) used in the non-normalised
764 // incomplete gammas:
766 template <class T, class Policy>
767 T full_igamma_prefix(T a, T z, const Policy& pol)
776 if((alz < tools::log_max_value<T>()) && (-z > tools::log_min_value<T>()))
778 prefix = pow(z, a) * exp(-z);
782 prefix = pow(z / exp(z/a), a);
786 prefix = exp(alz - z);
791 if(alz > tools::log_min_value<T>())
793 prefix = pow(z, a) * exp(-z);
795 else if(z/a < tools::log_max_value<T>())
797 prefix = pow(z / exp(z/a), a);
801 prefix = exp(alz - z);
805 // This error handling isn't very good: it happens after the fact
806 // rather than before it...
808 if((boost::math::fpclassify)(prefix) == (int)FP_INFINITE)
809 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);
814 // Compute (z^a)(e^-z)/tgamma(a)
815 // most if the error occurs in this function:
817 template <class T, class Policy, class Lanczos>
818 T regularised_gamma_prefix(T a, T z, const Policy& pol, const Lanczos& l)
821 T agh = a + static_cast<T>(Lanczos::g()) - T(0.5);
823 T d = ((z - a) - static_cast<T>(Lanczos::g()) + T(0.5)) / agh;
828 // We have to treat a < 1 as a special case because our Lanczos
829 // approximations are optimised against the factorials with a > 1,
830 // and for high precision types especially (128-bit reals for example)
831 // very small values of a can give rather eroneous results for gamma
832 // unless we do this:
834 // TODO: is this still required? Lanczos approx should be better now?
836 if(z <= tools::log_min_value<T>())
838 // Oh dear, have to use logs, should be free of cancellation errors though:
839 return exp(a * log(z) - z - lgamma_imp(a, pol, l));
843 // direct calculation, no danger of overflow as gamma(a) < 1/a
845 return pow(z, a) * exp(-z) / gamma_imp(a, pol, l);
848 else if((fabs(d*d*a) <= 100) && (a > 150))
850 // special case for large a and a ~ z.
851 prefix = a * boost::math::log1pmx(d, pol) + z * static_cast<T>(0.5 - Lanczos::g()) / agh;
852 prefix = exp(prefix);
858 // direct computation is most accurate, but use various fallbacks
859 // for different parts of the problem domain:
861 T alz = a * log(z / agh);
863 if(((std::min)(alz, amz) <= tools::log_min_value<T>()) || ((std::max)(alz, amz) >= tools::log_max_value<T>()))
866 if(((std::min)(alz, amz)/2 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/2 < tools::log_max_value<T>()))
868 // compute square root of the result and then square it:
869 T sq = pow(z / agh, a / 2) * exp(amz / 2);
872 else if(((std::min)(alz, amz)/4 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/4 < tools::log_max_value<T>()) && (z > a))
874 // compute the 4th root of the result then square it twice:
875 T sq = pow(z / agh, a / 4) * exp(amz / 4);
879 else if((amza > tools::log_min_value<T>()) && (amza < tools::log_max_value<T>()))
881 prefix = pow((z * exp(amza)) / agh, a);
885 prefix = exp(alz + amz);
890 prefix = pow(z / agh, a) * exp(amz);
893 prefix *= sqrt(agh / boost::math::constants::e<T>()) / Lanczos::lanczos_sum_expG_scaled(a);
897 // And again, without Lanczos support:
899 template <class T, class Policy>
900 T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined_lanczos&)
904 T limit = (std::max)(T(10), a);
905 T sum = detail::lower_gamma_series(a, limit, pol) / a;
906 sum += detail::upper_gamma_fraction(a, limit, ::boost::math::policies::get_epsilon<T, Policy>());
910 // special case for small a:
911 T prefix = pow(z / 10, a);
915 prefix = pow((z * exp((10-z)/a)) / 10, a);
923 T alzoa = a * log(zoa);
925 if(((std::min)(alzoa, amz) <= tools::log_min_value<T>()) || ((std::max)(alzoa, amz) >= tools::log_max_value<T>()))
928 if((amza <= tools::log_min_value<T>()) || (amza >= tools::log_max_value<T>()))
930 prefix = exp(alzoa + amz);
934 prefix = pow(zoa * exp(amza), a);
939 prefix = pow(zoa, a) * exp(amz);
945 // Upper gamma fraction for very small a:
947 template <class T, class Policy>
948 inline T tgamma_small_upper_part(T a, T x, const Policy& pol, T* pgam = 0, bool invert = false, T* pderivative = 0)
950 BOOST_MATH_STD_USING // ADL of std functions.
952 // Compute the full upper fraction (Q) when a is very small:
955 result = boost::math::tgamma1pm1(a, pol);
957 *pgam = (result + 1) / a;
958 T p = boost::math::powm1(x, a, pol);
961 detail::small_gamma2_series<T> s(a, x);
962 boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>() - 10;
965 *pderivative = p / (*pgam * exp(x));
966 T init_value = invert ? *pgam : 0;
967 result = -p * tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, (init_value - result) / p);
968 policies::check_series_iterations<T>("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)", max_iter, pol);
974 // Upper gamma fraction for integer a:
976 template <class T, class Policy>
977 inline T finite_gamma_q(T a, T x, Policy const& pol, T* pderivative = 0)
980 // Calculates normalised Q when a is an integer:
988 for(unsigned n = 1; n < a; ++n)
997 *pderivative = e * pow(x, a) / boost::math::unchecked_factorial<T>(itrunc(T(a - 1), pol));
1002 // Upper gamma fraction for half integer a:
1004 template <class T, class Policy>
1005 T finite_half_gamma_q(T a, T x, T* p_derivative, const Policy& pol)
1008 // Calculates normalised Q when a is a half-integer:
1010 BOOST_MATH_STD_USING
1011 T e = boost::math::erfc(sqrt(x), pol);
1012 if((e != 0) && (a > 1))
1014 T term = exp(-x) / sqrt(constants::pi<T>() * x);
1016 static const T half = T(1) / 2;
1019 for(unsigned n = 2; n < a; ++n)
1031 else if(p_derivative)
1033 // We'll be dividing by x later, so calculate derivative * x:
1034 *p_derivative = sqrt(x) * exp(-x) / constants::root_pi<T>();
1039 // Main incomplete gamma entry point, handles all four incomplete gamma's:
1041 template <class T, class Policy>
1042 T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
1043 const Policy& pol, T* p_derivative)
1045 static const char* function = "boost::math::gamma_p<%1%>(%1%, %1%)";
1047 return policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
1049 return policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
1051 BOOST_MATH_STD_USING
1053 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1055 T result = 0; // Just to avoid warning C4701: potentially uninitialized local variable 'result' used
1057 if(a >= max_factorial<T>::value && !normalised)
1060 // When we're computing the non-normalized incomplete gamma
1061 // and a is large the result is rather hard to compute unless
1062 // we use logs. There are really two options - if x is a long
1063 // way from a in value then we can reliably use methods 2 and 4
1064 // below in logarithmic form and go straight to the result.
1065 // Otherwise we let the regularized gamma take the strain
1066 // (the result is unlikely to unerflow in the central region anyway)
1067 // and combine with lgamma in the hopes that we get a finite result.
1069 if(invert && (a * 4 < x))
1071 // This is method 4 below, done in logs:
1072 result = a * log(x) - x;
1074 *p_derivative = exp(result);
1075 result += log(upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>()));
1077 else if(!invert && (a > 4 * x))
1079 // This is method 2 below, done in logs:
1080 result = a * log(x) - x;
1082 *p_derivative = exp(result);
1084 result += log(detail::lower_gamma_series(a, x, pol, init_value) / a);
1088 result = gamma_incomplete_imp(a, x, true, invert, pol, p_derivative);
1093 // Try http://functions.wolfram.com/06.06.06.0039.01
1094 result = 1 + 1 / (12 * a) + 1 / (288 * a * a);
1095 result = log(result) - a + (a - 0.5f) * log(a) + log(boost::math::constants::root_two_pi<T>());
1097 *p_derivative = exp(a * log(x) - x);
1101 // This is method 2 below, done in logs, we're really outside the
1102 // range of this method, but since the result is almost certainly
1103 // infinite, we should probably be OK:
1104 result = a * log(x) - x;
1106 *p_derivative = exp(result);
1108 result += log(detail::lower_gamma_series(a, x, pol, init_value) / a);
1113 result = log(result) + boost::math::lgamma(a, pol);
1116 if(result > tools::log_max_value<T>())
1117 return policies::raise_overflow_error<T>(function, 0, pol);
1121 BOOST_ASSERT((p_derivative == 0) || (normalised == true));
1123 bool is_int, is_half_int;
1124 bool is_small_a = (a < 30) && (a <= x + 1) && (x < tools::log_max_value<T>());
1129 is_half_int = is_int ? false : (fabs(fa - a) == 0.5f);
1133 is_int = is_half_int = false;
1138 if(is_int && (x > 0.6))
1140 // calculate Q via finite sum:
1144 else if(is_half_int && (x > 0.2))
1146 // calculate Q via finite sum for half integer a:
1150 else if((x < tools::root_epsilon<T>()) && (a > 1))
1157 // Changeover criterion chosen to give a changeover at Q ~ 0.33
1159 if(-0.4 / log(x) < a)
1171 // Changover here occurs when P ~ 0.75 or Q ~ 0.25:
1185 // Begin by testing whether we're in the "bad" zone
1186 // where the result will be near 0.5 and the usual
1187 // series and continued fractions are slow to converge:
1189 bool use_temme = false;
1190 if(normalised && std::numeric_limits<T>::is_specialized && (a > 20))
1192 T sigma = fabs((x-a)/a);
1193 if((a > 200) && (policies::digits<T, Policy>() <= 113))
1196 // This limit is chosen so that we use Temme's expansion
1197 // only if the result would be larger than about 10^-6.
1198 // Below that the regular series and continued fractions
1199 // converge OK, and if we use Temme's method we get increasing
1200 // errors from the dominant erfc term as it's (inexact) argument
1201 // increases in magnitude.
1203 if(20 / a > sigma * sigma)
1206 else if(policies::digits<T, Policy>() <= 64)
1208 // Note in this zone we can't use Temme's expansion for
1209 // types longer than an 80-bit real:
1210 // it would require too many terms in the polynomials.
1222 // Regular case where the result will not be too close to 0.5.
1224 // Changeover here occurs at P ~ Q ~ 0.5
1225 // Note that series computation of P is about x2 faster than continued fraction
1226 // calculation of Q, so try and use the CF only when really necessary, especially
1229 if(x - (1 / (3 * x)) < a)
1245 result = finite_gamma_q(a, x, pol, p_derivative);
1246 if(normalised == false)
1247 result *= boost::math::tgamma(a, pol);
1252 result = finite_half_gamma_q(a, x, p_derivative, pol);
1253 if(normalised == false)
1254 result *= boost::math::tgamma(a, pol);
1255 if(p_derivative && (*p_derivative == 0))
1256 *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
1262 result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
1264 *p_derivative = result;
1268 // If we're going to be inverting the result then we can
1269 // reduce the number of series evaluations by quite
1270 // a few iterations if we set an initial value for the
1271 // series sum based on what we'll end up subtracting it from
1273 // Have to be careful though that this optimization doesn't
1274 // lead to spurious numberic overflow. Note that the
1275 // scary/expensive overflow checks below are more often
1276 // than not bypassed in practice for "sensible" input
1280 bool optimised_invert = false;
1283 init_value = (normalised ? 1 : boost::math::tgamma(a, pol));
1284 if(normalised || (result >= 1) || (tools::max_value<T>() * result > init_value))
1286 init_value /= result;
1287 if(normalised || (a < 1) || (tools::max_value<T>() / a > init_value))
1290 optimised_invert = true;
1298 result *= detail::lower_gamma_series(a, x, pol, init_value) / a;
1299 if(optimised_invert)
1312 result = tgamma_small_upper_part(a, x, pol, &g, invert, p_derivative);
1321 result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
1323 *p_derivative = result;
1325 result *= upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>());
1331 // Use compile time dispatch to the appropriate
1332 // Temme asymptotic expansion. This may be dead code
1333 // if T does not have numeric limits support, or has
1334 // too many digits for the most precise version of
1335 // these expansions, in that case we'll be calling
1336 // an empty function.
1338 typedef typename policies::precision<T, Policy>::type precision_type;
1340 typedef typename mpl::if_<
1341 mpl::or_<mpl::equal_to<precision_type, mpl::int_<0> >,
1342 mpl::greater<precision_type, mpl::int_<113> > >,
1345 mpl::less_equal<precision_type, mpl::int_<53> >,
1348 mpl::less_equal<precision_type, mpl::int_<64> >,
1355 result = igamma_temme_large(a, x, pol, static_cast<tag_type const*>(0));
1359 *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
1364 // x is so small that P is necessarily very small too,
1365 // use http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/
1366 result = !normalised ? pow(x, a) / (a) : pow(x, a) / boost::math::tgamma(a + 1, pol);
1367 result *= 1 - a * x / (a + 1);
1371 if(normalised && (result > 1))
1375 T gam = normalised ? 1 : boost::math::tgamma(a, pol);
1376 result = gam - result;
1381 // Need to convert prefix term to derivative:
1383 if((x < 1) && (tools::max_value<T>() * x < *p_derivative))
1385 // overflow, just return an arbitrarily large value:
1386 *p_derivative = tools::max_value<T>() / 2;
1396 // Ratios of two gamma functions:
1398 template <class T, class Policy, class Lanczos>
1399 T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const Lanczos& l)
1401 BOOST_MATH_STD_USING
1402 if(z < tools::epsilon<T>())
1405 // We get spurious numeric overflow unless we're very careful, this
1406 // can occur either inside Lanczos::lanczos_sum(z) or in the
1407 // final combination of terms, to avoid this, split the product up
1408 // into 2 (or 3) parts:
1410 // G(z) / G(L) = 1 / (z * G(L)) ; z < eps, L = z + delta = delta
1411 // z * G(L) = z * G(lim) * (G(L)/G(lim)) ; lim = largest factorial
1413 if(boost::math::max_factorial<T>::value < delta)
1415 T ratio = tgamma_delta_ratio_imp_lanczos(delta, T(boost::math::max_factorial<T>::value - delta), pol, l);
1417 ratio *= boost::math::unchecked_factorial<T>(boost::math::max_factorial<T>::value - 1);
1422 return 1 / (z * boost::math::tgamma(z + delta, pol));
1425 T zgh = static_cast<T>(z + Lanczos::g() - constants::half<T>());
1429 if(fabs(delta) < 10)
1430 result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol));
1436 if(fabs(delta) < 10)
1438 result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol));
1442 result = pow(zgh / (zgh + delta), z - constants::half<T>());
1444 // Split the calculation up to avoid spurious overflow:
1445 result *= Lanczos::lanczos_sum(z) / Lanczos::lanczos_sum(T(z + delta));
1447 result *= pow(constants::e<T>() / (zgh + delta), delta);
1451 // And again without Lanczos support this time:
1453 template <class T, class Policy>
1454 T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const lanczos::undefined_lanczos&)
1456 BOOST_MATH_STD_USING
1458 // The upper gamma fraction is *very* slow for z < 6, actually it's very
1459 // slow to converge everywhere but recursing until z > 6 gets rid of the
1460 // worst of it's behaviour.
1464 while((zd < 6) && (z < 6))
1473 prefix *= exp(-z * boost::math::log1p(delta / z, pol));
1477 prefix *= pow(z / zd, z);
1479 prefix *= pow(constants::e<T>() / zd, delta);
1480 T sum = detail::lower_gamma_series(z, z, pol) / z;
1481 sum += detail::upper_gamma_fraction(z, z, ::boost::math::policies::get_epsilon<T, Policy>());
1482 T sumd = detail::lower_gamma_series(zd, zd, pol) / zd;
1483 sumd += detail::upper_gamma_fraction(zd, zd, ::boost::math::policies::get_epsilon<T, Policy>());
1485 if(fabs(tools::max_value<T>() / prefix) < fabs(sum))
1486 return policies::raise_overflow_error<T>("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Result of tgamma is too large to represent.", pol);
1487 return sum * prefix;
1490 template <class T, class Policy>
1491 T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol)
1493 BOOST_MATH_STD_USING
1495 if((z <= 0) || (z + delta <= 0))
1497 // This isn't very sofisticated, or accurate, but it does work:
1498 return boost::math::tgamma(z, pol) / boost::math::tgamma(z + delta, pol);
1501 if(floor(delta) == delta)
1506 // Both z and delta are integers, see if we can just use table lookup
1507 // of the factorials to get the result:
1509 if((z <= max_factorial<T>::value) && (z + delta <= max_factorial<T>::value))
1511 return unchecked_factorial<T>((unsigned)itrunc(z, pol) - 1) / unchecked_factorial<T>((unsigned)itrunc(T(z + delta), pol) - 1);
1514 if(fabs(delta) < 20)
1517 // delta is a small integer, we can use a finite product:
1525 while(0 != (delta += 1))
1535 while(0 != (delta -= 1))
1544 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1545 return tgamma_delta_ratio_imp_lanczos(z, delta, pol, lanczos_type());
1548 template <class T, class Policy>
1549 T tgamma_ratio_imp(T x, T y, const Policy& pol)
1551 BOOST_MATH_STD_USING
1553 if((x <= 0) || (boost::math::isinf)(x))
1554 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);
1555 if((y <= 0) || (boost::math::isinf)(y))
1556 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);
1558 if(x <= tools::min_value<T>())
1560 // Special case for denorms...Ugh.
1561 T shift = ldexp(T(1), tools::digits<T>());
1562 return shift * tgamma_ratio_imp(T(x * shift), y, pol);
1565 if((x < max_factorial<T>::value) && (y < max_factorial<T>::value))
1567 // Rather than subtracting values, lets just call the gamma functions directly:
1568 return boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
1573 if(y < 2 * max_factorial<T>::value)
1575 // We need to sidestep on x as well, otherwise we'll underflow
1576 // before we get to factor in the prefix term:
1579 while(y >= max_factorial<T>::value)
1584 return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
1587 // result is almost certainly going to underflow to zero, try logs just in case:
1589 return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
1593 if(x < 2 * max_factorial<T>::value)
1595 // We need to sidestep on y as well, otherwise we'll overflow
1596 // before we get to factor in the prefix term:
1599 while(x >= max_factorial<T>::value)
1604 return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
1607 // Result will almost certainly overflow, try logs just in case:
1609 return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
1612 // Regular case, x and y both large and similar in magnitude:
1614 return boost::math::tgamma_delta_ratio(x, y - x, pol);
1617 template <class T, class Policy>
1618 T gamma_p_derivative_imp(T a, T x, const Policy& pol)
1620 BOOST_MATH_STD_USING
1622 // Usual error checks first:
1625 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);
1627 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);
1629 // Now special cases:
1633 return (a > 1) ? 0 :
1634 (a == 1) ? 1 : policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
1639 typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1640 T f1 = detail::regularised_gamma_prefix(a, x, pol, lanczos_type());
1641 if((x < 1) && (tools::max_value<T>() * x < f1))
1644 return policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
1648 // Underflow in calculation, use logs instead:
1649 f1 = a * log(x) - x - lgamma(a, pol) - log(x);
1658 template <class T, class Policy>
1659 inline typename tools::promote_args<T>::type
1660 tgamma(T z, const Policy& /* pol */, const mpl::true_)
1662 BOOST_FPU_EXCEPTION_GUARD
1663 typedef typename tools::promote_args<T>::type result_type;
1664 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1665 typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1666 typedef typename policies::normalise<
1668 policies::promote_float<false>,
1669 policies::promote_double<false>,
1670 policies::discrete_quantile<>,
1671 policies::assert_undefined<> >::type forwarding_policy;
1672 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%)");
1675 template <class T, class Policy>
1676 struct igamma_initializer
1682 typedef typename policies::precision<T, Policy>::type precision_type;
1684 typedef typename mpl::if_<
1685 mpl::or_<mpl::equal_to<precision_type, mpl::int_<0> >,
1686 mpl::greater<precision_type, mpl::int_<113> > >,
1689 mpl::less_equal<precision_type, mpl::int_<53> >,
1692 mpl::less_equal<precision_type, mpl::int_<64> >,
1699 do_init(tag_type());
1702 static void do_init(const mpl::int_<N>&)
1704 // If std::numeric_limits<T>::digits is zero, we must not call
1705 // our inituialization code here as the precision presumably
1706 // varies at runtime, and will not have been set yet. Plus the
1707 // code requiring initialization isn't called when digits == 0.
1708 if(std::numeric_limits<T>::digits)
1710 boost::math::gamma_p(static_cast<T>(400), static_cast<T>(400), Policy());
1713 static void do_init(const mpl::int_<53>&){}
1714 void force_instantiate()const{}
1716 static const init initializer;
1717 static void force_instantiate()
1719 initializer.force_instantiate();
1723 template <class T, class Policy>
1724 const typename igamma_initializer<T, Policy>::init igamma_initializer<T, Policy>::initializer;
1726 template <class T, class Policy>
1727 struct lgamma_initializer
1733 typedef typename policies::precision<T, Policy>::type precision_type;
1734 typedef typename mpl::if_<
1736 mpl::less_equal<precision_type, mpl::int_<64> >,
1737 mpl::greater<precision_type, mpl::int_<0> >
1742 mpl::less_equal<precision_type, mpl::int_<113> >,
1743 mpl::greater<precision_type, mpl::int_<0> >
1745 mpl::int_<113>, mpl::int_<0> >::type
1747 do_init(tag_type());
1749 static void do_init(const mpl::int_<64>&)
1751 boost::math::lgamma(static_cast<T>(2.5), Policy());
1752 boost::math::lgamma(static_cast<T>(1.25), Policy());
1753 boost::math::lgamma(static_cast<T>(1.75), Policy());
1755 static void do_init(const mpl::int_<113>&)
1757 boost::math::lgamma(static_cast<T>(2.5), Policy());
1758 boost::math::lgamma(static_cast<T>(1.25), Policy());
1759 boost::math::lgamma(static_cast<T>(1.5), Policy());
1760 boost::math::lgamma(static_cast<T>(1.75), Policy());
1762 static void do_init(const mpl::int_<0>&)
1765 void force_instantiate()const{}
1767 static const init initializer;
1768 static void force_instantiate()
1770 initializer.force_instantiate();
1774 template <class T, class Policy>
1775 const typename lgamma_initializer<T, Policy>::init lgamma_initializer<T, Policy>::initializer;
1777 template <class T1, class T2, class Policy>
1778 inline typename tools::promote_args<T1, T2>::type
1779 tgamma(T1 a, T2 z, const Policy&, const mpl::false_)
1781 BOOST_FPU_EXCEPTION_GUARD
1782 typedef typename tools::promote_args<T1, T2>::type result_type;
1783 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1784 // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1785 typedef typename policies::normalise<
1787 policies::promote_float<false>,
1788 policies::promote_double<false>,
1789 policies::discrete_quantile<>,
1790 policies::assert_undefined<> >::type forwarding_policy;
1792 igamma_initializer<value_type, forwarding_policy>::force_instantiate();
1794 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
1795 detail::gamma_incomplete_imp(static_cast<value_type>(a),
1796 static_cast<value_type>(z), false, true,
1797 forwarding_policy(), static_cast<value_type*>(0)), "boost::math::tgamma<%1%>(%1%, %1%)");
1800 template <class T1, class T2>
1801 inline typename tools::promote_args<T1, T2>::type
1802 tgamma(T1 a, T2 z, const mpl::false_ tag)
1804 return tgamma(a, z, policies::policy<>(), tag);
1808 } // namespace detail
1811 inline typename tools::promote_args<T>::type
1814 return tgamma(z, policies::policy<>());
1817 template <class T, class Policy>
1818 inline typename tools::promote_args<T>::type
1819 lgamma(T z, int* sign, const Policy&)
1821 BOOST_FPU_EXCEPTION_GUARD
1822 typedef typename tools::promote_args<T>::type result_type;
1823 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1824 typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1825 typedef typename policies::normalise<
1827 policies::promote_float<false>,
1828 policies::promote_double<false>,
1829 policies::discrete_quantile<>,
1830 policies::assert_undefined<> >::type forwarding_policy;
1832 detail::lgamma_initializer<value_type, forwarding_policy>::force_instantiate();
1834 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%)");
1838 inline typename tools::promote_args<T>::type
1839 lgamma(T z, int* sign)
1841 return lgamma(z, sign, policies::policy<>());
1844 template <class T, class Policy>
1845 inline typename tools::promote_args<T>::type
1846 lgamma(T x, const Policy& pol)
1848 return ::boost::math::lgamma(x, 0, pol);
1852 inline typename tools::promote_args<T>::type
1855 return ::boost::math::lgamma(x, 0, policies::policy<>());
1858 template <class T, class Policy>
1859 inline typename tools::promote_args<T>::type
1860 tgamma1pm1(T z, const Policy& /* pol */)
1862 BOOST_FPU_EXCEPTION_GUARD
1863 typedef typename tools::promote_args<T>::type result_type;
1864 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1865 typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1866 typedef typename policies::normalise<
1868 policies::promote_float<false>,
1869 policies::promote_double<false>,
1870 policies::discrete_quantile<>,
1871 policies::assert_undefined<> >::type forwarding_policy;
1873 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%)");
1877 inline typename tools::promote_args<T>::type
1880 return tgamma1pm1(z, policies::policy<>());
1884 // Full upper incomplete gamma:
1886 template <class T1, class T2>
1887 inline typename tools::promote_args<T1, T2>::type
1891 // Type T2 could be a policy object, or a value, select the
1892 // right overload based on T2:
1894 typedef typename policies::is_policy<T2>::type maybe_policy;
1895 return detail::tgamma(a, z, maybe_policy());
1897 template <class T1, class T2, class Policy>
1898 inline typename tools::promote_args<T1, T2>::type
1899 tgamma(T1 a, T2 z, const Policy& pol)
1901 return detail::tgamma(a, z, pol, mpl::false_());
1904 // Full lower incomplete gamma:
1906 template <class T1, class T2, class Policy>
1907 inline typename tools::promote_args<T1, T2>::type
1908 tgamma_lower(T1 a, T2 z, const Policy&)
1910 BOOST_FPU_EXCEPTION_GUARD
1911 typedef typename tools::promote_args<T1, T2>::type result_type;
1912 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1913 // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1914 typedef typename policies::normalise<
1916 policies::promote_float<false>,
1917 policies::promote_double<false>,
1918 policies::discrete_quantile<>,
1919 policies::assert_undefined<> >::type forwarding_policy;
1921 detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
1923 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
1924 detail::gamma_incomplete_imp(static_cast<value_type>(a),
1925 static_cast<value_type>(z), false, false,
1926 forwarding_policy(), static_cast<value_type*>(0)), "tgamma_lower<%1%>(%1%, %1%)");
1928 template <class T1, class T2>
1929 inline typename tools::promote_args<T1, T2>::type
1930 tgamma_lower(T1 a, T2 z)
1932 return tgamma_lower(a, z, policies::policy<>());
1935 // Regularised upper incomplete gamma:
1937 template <class T1, class T2, class Policy>
1938 inline typename tools::promote_args<T1, T2>::type
1939 gamma_q(T1 a, T2 z, const Policy& /* pol */)
1941 BOOST_FPU_EXCEPTION_GUARD
1942 typedef typename tools::promote_args<T1, T2>::type result_type;
1943 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1944 // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1945 typedef typename policies::normalise<
1947 policies::promote_float<false>,
1948 policies::promote_double<false>,
1949 policies::discrete_quantile<>,
1950 policies::assert_undefined<> >::type forwarding_policy;
1952 detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
1954 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
1955 detail::gamma_incomplete_imp(static_cast<value_type>(a),
1956 static_cast<value_type>(z), true, true,
1957 forwarding_policy(), static_cast<value_type*>(0)), "gamma_q<%1%>(%1%, %1%)");
1959 template <class T1, class T2>
1960 inline typename tools::promote_args<T1, T2>::type
1963 return gamma_q(a, z, policies::policy<>());
1966 // Regularised lower incomplete gamma:
1968 template <class T1, class T2, class Policy>
1969 inline typename tools::promote_args<T1, T2>::type
1970 gamma_p(T1 a, T2 z, const Policy&)
1972 BOOST_FPU_EXCEPTION_GUARD
1973 typedef typename tools::promote_args<T1, T2>::type result_type;
1974 typedef typename policies::evaluation<result_type, Policy>::type value_type;
1975 // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1976 typedef typename policies::normalise<
1978 policies::promote_float<false>,
1979 policies::promote_double<false>,
1980 policies::discrete_quantile<>,
1981 policies::assert_undefined<> >::type forwarding_policy;
1983 detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
1985 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
1986 detail::gamma_incomplete_imp(static_cast<value_type>(a),
1987 static_cast<value_type>(z), true, false,
1988 forwarding_policy(), static_cast<value_type*>(0)), "gamma_p<%1%>(%1%, %1%)");
1990 template <class T1, class T2>
1991 inline typename tools::promote_args<T1, T2>::type
1994 return gamma_p(a, z, policies::policy<>());
1997 // ratios of gamma functions:
1998 template <class T1, class T2, class Policy>
1999 inline typename tools::promote_args<T1, T2>::type
2000 tgamma_delta_ratio(T1 z, T2 delta, const Policy& /* pol */)
2002 BOOST_FPU_EXCEPTION_GUARD
2003 typedef typename tools::promote_args<T1, T2>::type result_type;
2004 typedef typename policies::evaluation<result_type, Policy>::type value_type;
2005 typedef typename policies::normalise<
2007 policies::promote_float<false>,
2008 policies::promote_double<false>,
2009 policies::discrete_quantile<>,
2010 policies::assert_undefined<> >::type forwarding_policy;
2012 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%)");
2014 template <class T1, class T2>
2015 inline typename tools::promote_args<T1, T2>::type
2016 tgamma_delta_ratio(T1 z, T2 delta)
2018 return tgamma_delta_ratio(z, delta, policies::policy<>());
2020 template <class T1, class T2, class Policy>
2021 inline typename tools::promote_args<T1, T2>::type
2022 tgamma_ratio(T1 a, T2 b, const Policy&)
2024 typedef typename tools::promote_args<T1, T2>::type result_type;
2025 typedef typename policies::evaluation<result_type, Policy>::type value_type;
2026 typedef typename policies::normalise<
2028 policies::promote_float<false>,
2029 policies::promote_double<false>,
2030 policies::discrete_quantile<>,
2031 policies::assert_undefined<> >::type forwarding_policy;
2033 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%)");
2035 template <class T1, class T2>
2036 inline typename tools::promote_args<T1, T2>::type
2037 tgamma_ratio(T1 a, T2 b)
2039 return tgamma_ratio(a, b, policies::policy<>());
2042 template <class T1, class T2, class Policy>
2043 inline typename tools::promote_args<T1, T2>::type
2044 gamma_p_derivative(T1 a, T2 x, const Policy&)
2046 BOOST_FPU_EXCEPTION_GUARD
2047 typedef typename tools::promote_args<T1, T2>::type result_type;
2048 typedef typename policies::evaluation<result_type, Policy>::type value_type;
2049 typedef typename policies::normalise<
2051 policies::promote_float<false>,
2052 policies::promote_double<false>,
2053 policies::discrete_quantile<>,
2054 policies::assert_undefined<> >::type forwarding_policy;
2056 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%)");
2058 template <class T1, class T2>
2059 inline typename tools::promote_args<T1, T2>::type
2060 gamma_p_derivative(T1 a, T2 x)
2062 return gamma_p_derivative(a, x, policies::policy<>());
2066 } // namespace boost
2069 # pragma warning(pop)
2072 #include <boost/math/special_functions/detail/igamma_inverse.hpp>
2073 #include <boost/math/special_functions/detail/gamma_inva.hpp>
2074 #include <boost/math/special_functions/erf.hpp>
2076 #endif // BOOST_MATH_SF_GAMMA_HPP