Chris@16: 
Chris@101: //  Copyright John Maddock 2006-7, 2013-14.
Chris@101: //  Copyright Paul A. Bristow 2007, 2013-14.
Chris@101: //  Copyright Nikhar Agrawal 2013-14
Chris@101: //  Copyright Christopher Kormanyos 2013-14
Chris@16: 
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_MATH_SF_GAMMA_HPP
Chris@16: #define BOOST_MATH_SF_GAMMA_HPP
Chris@16: 
Chris@16: #ifdef _MSC_VER
Chris@16: #pragma once
Chris@16: #endif
Chris@16: 
Chris@16: #include <boost/config.hpp>
Chris@16: #include <boost/math/tools/series.hpp>
Chris@16: #include <boost/math/tools/fraction.hpp>
Chris@16: #include <boost/math/tools/precision.hpp>
Chris@16: #include <boost/math/tools/promotion.hpp>
Chris@16: #include <boost/math/policies/error_handling.hpp>
Chris@16: #include <boost/math/constants/constants.hpp>
Chris@16: #include <boost/math/special_functions/math_fwd.hpp>
Chris@16: #include <boost/math/special_functions/log1p.hpp>
Chris@16: #include <boost/math/special_functions/trunc.hpp>
Chris@16: #include <boost/math/special_functions/powm1.hpp>
Chris@16: #include <boost/math/special_functions/sqrt1pm1.hpp>
Chris@16: #include <boost/math/special_functions/lanczos.hpp>
Chris@16: #include <boost/math/special_functions/fpclassify.hpp>
Chris@16: #include <boost/math/special_functions/detail/igamma_large.hpp>
Chris@16: #include <boost/math/special_functions/detail/unchecked_factorial.hpp>
Chris@16: #include <boost/math/special_functions/detail/lgamma_small.hpp>
Chris@101: #include <boost/math/special_functions/bernoulli.hpp>
Chris@16: #include <boost/type_traits/is_convertible.hpp>
Chris@16: #include <boost/assert.hpp>
Chris@16: #include <boost/mpl/greater.hpp>
Chris@16: #include <boost/mpl/equal_to.hpp>
Chris@16: #include <boost/mpl/greater.hpp>
Chris@16: 
Chris@16: #include <boost/config/no_tr1/cmath.hpp>
Chris@16: #include <algorithm>
Chris@16: 
Chris@16: #ifdef BOOST_MSVC
Chris@16: # pragma warning(push)
Chris@16: # pragma warning(disable: 4702) // unreachable code (return after domain_error throw).
Chris@16: # pragma warning(disable: 4127) // conditional expression is constant.
Chris@16: # pragma warning(disable: 4100) // unreferenced formal parameter.
Chris@16: // Several variables made comments,
Chris@16: // but some difficulty as whether referenced on not may depend on macro values.
Chris@16: // So to be safe, 4100 warnings suppressed.
Chris@16: // TODO - revisit this?
Chris@16: #endif
Chris@16: 
Chris@16: namespace boost{ namespace math{
Chris@16: 
Chris@16: namespace detail{
Chris@16: 
Chris@16: template <class T>
Chris@16: inline bool is_odd(T v, const boost::true_type&)
Chris@16: {
Chris@16:    int i = static_cast<int>(v);
Chris@16:    return i&1;
Chris@16: }
Chris@16: template <class T>
Chris@16: inline bool is_odd(T v, const boost::false_type&)
Chris@16: {
Chris@16:    // Oh dear can't cast T to int!
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    T modulus = v - 2 * floor(v/2);
Chris@16:    return static_cast<bool>(modulus != 0);
Chris@16: }
Chris@16: template <class T>
Chris@16: inline bool is_odd(T v)
Chris@16: {
Chris@16:    return is_odd(v, ::boost::is_convertible<T, int>());
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: T sinpx(T z)
Chris@16: {
Chris@16:    // Ad hoc function calculates x * sin(pi * x),
Chris@16:    // taking extra care near when x is near a whole number.
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    int sign = 1;
Chris@16:    if(z < 0)
Chris@16:    {
Chris@16:       z = -z;
Chris@16:    }
Chris@16:    T fl = floor(z);
Chris@16:    T dist;
Chris@16:    if(is_odd(fl))
Chris@16:    {
Chris@16:       fl += 1;
Chris@16:       dist = fl - z;
Chris@16:       sign = -sign;
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       dist = z - fl;
Chris@16:    }
Chris@16:    BOOST_ASSERT(fl >= 0);
Chris@16:    if(dist > 0.5)
Chris@16:       dist = 1 - dist;
Chris@16:    T result = sin(dist*boost::math::constants::pi<T>());
Chris@16:    return sign*z*result;
Chris@16: } // template <class T> T sinpx(T z)
Chris@16: //
Chris@16: // tgamma(z), with Lanczos support:
Chris@16: //
Chris@16: template <class T, class Policy, class Lanczos>
Chris@16: T gamma_imp(T z, const Policy& pol, const Lanczos& l)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    T result = 1;
Chris@16: 
Chris@16: #ifdef BOOST_MATH_INSTRUMENT
Chris@16:    static bool b = false;
Chris@16:    if(!b)
Chris@16:    {
Chris@16:       std::cout << "tgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;
Chris@16:       b = true;
Chris@16:    }
Chris@16: #endif
Chris@16:    static const char* function = "boost::math::tgamma<%1%>(%1%)";
Chris@16: 
Chris@16:    if(z <= 0)
Chris@16:    {
Chris@16:       if(floor(z) == z)
Chris@16:          return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
Chris@16:       if(z <= -20)
Chris@16:       {
Chris@16:          result = gamma_imp(T(-z), pol, l) * sinpx(z);
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:          if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
Chris@101:             return -boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
Chris@16:          result = -boost::math::constants::pi<T>() / result;
Chris@16:          if(result == 0)
Chris@16:             return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
Chris@16:          if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL)
Chris@16:             return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", result, pol);
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:          return result;
Chris@16:       }
Chris@16: 
Chris@16:       // shift z to > 1:
Chris@16:       while(z < 0)
Chris@16:       {
Chris@16:          result /= z;
Chris@16:          z += 1;
Chris@16:       }
Chris@16:    }
Chris@16:    BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:    if((floor(z) == z) && (z < max_factorial<T>::value))
Chris@16:    {
Chris@16:       result *= unchecked_factorial<T>(itrunc(z, pol) - 1);
Chris@16:       BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:    }
Chris@101:    else if (z < tools::root_epsilon<T>())
Chris@101:    {
Chris@101:       if (z < 1 / tools::max_value<T>())
Chris@101:          result = policies::raise_overflow_error<T>(function, 0, pol);
Chris@101: 	   result *= 1 / z - constants::euler<T>();
Chris@101:    }
Chris@16:    else
Chris@16:    {
Chris@16:       result *= Lanczos::lanczos_sum(z);
Chris@16:       T zgh = (z + static_cast<T>(Lanczos::g()) - boost::math::constants::half<T>());
Chris@16:       T lzgh = log(zgh);
Chris@16:       BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:       BOOST_MATH_INSTRUMENT_VARIABLE(tools::log_max_value<T>());
Chris@16:       if(z * lzgh > tools::log_max_value<T>())
Chris@16:       {
Chris@16:          // we're going to overflow unless this is done with care:
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
Chris@16:          if(lzgh * z / 2 > tools::log_max_value<T>())
Chris@101:             return boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
Chris@16:          T hp = pow(zgh, (z / 2) - T(0.25));
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(hp);
Chris@16:          result *= hp / exp(zgh);
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:          if(tools::max_value<T>() / hp < result)
Chris@101:             return boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
Chris@16:          result *= hp;
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(pow(zgh, z - boost::math::constants::half<T>()));
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(exp(zgh));
Chris@16:          result *= pow(zgh, z - boost::math::constants::half<T>()) / exp(zgh);
Chris@16:          BOOST_MATH_INSTRUMENT_VARIABLE(result);
Chris@16:       }
Chris@16:    }
Chris@16:    return result;
Chris@16: }
Chris@16: //
Chris@16: // lgamma(z) with Lanczos support:
Chris@16: //
Chris@16: template <class T, class Policy, class Lanczos>
Chris@16: T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = 0)
Chris@16: {
Chris@16: #ifdef BOOST_MATH_INSTRUMENT
Chris@16:    static bool b = false;
Chris@16:    if(!b)
Chris@16:    {
Chris@16:       std::cout << "lgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;
Chris@16:       b = true;
Chris@16:    }
Chris@16: #endif
Chris@16: 
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    static const char* function = "boost::math::lgamma<%1%>(%1%)";
Chris@16: 
Chris@16:    T result = 0;
Chris@16:    int sresult = 1;
Chris@101:    if(z <= -tools::root_epsilon<T>())
Chris@16:    {
Chris@16:       // reflection formula:
Chris@16:       if(floor(z) == z)
Chris@16:          return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%.", z, pol);
Chris@16: 
Chris@16:       T t = sinpx(z);
Chris@16:       z = -z;
Chris@16:       if(t < 0)
Chris@16:       {
Chris@16:          t = -t;
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          sresult = -sresult;
Chris@16:       }
Chris@16:       result = log(boost::math::constants::pi<T>()) - lgamma_imp(z, pol, l) - log(t);
Chris@16:    }
Chris@101:    else if (z < tools::root_epsilon<T>())
Chris@101:    {
Chris@101: 	   if (0 == z)
Chris@101: 		   return policies::raise_pole_error<T>(function, "Evaluation of lgamma at %1%.", z, pol);
Chris@101:       if (fabs(z) < 1 / tools::max_value<T>())
Chris@101:          result = -log(fabs(z));
Chris@101:       else
Chris@101: 	      result = log(fabs(1 / z - constants::euler<T>()));
Chris@101: 	   if (z < 0)
Chris@101: 		sresult = -1;
Chris@101:    }
Chris@16:    else if(z < 15)
Chris@16:    {
Chris@16:       typedef typename policies::precision<T, Policy>::type precision_type;
Chris@16:       typedef typename mpl::if_<
Chris@16:          mpl::and_<
Chris@16:             mpl::less_equal<precision_type, mpl::int_<64> >, 
Chris@16:             mpl::greater<precision_type, mpl::int_<0> > 
Chris@16:          >,
Chris@16:          mpl::int_<64>,
Chris@16:          typename mpl::if_<
Chris@16:             mpl::and_<
Chris@16:                mpl::less_equal<precision_type, mpl::int_<113> >,
Chris@16:                mpl::greater<precision_type, mpl::int_<0> > 
Chris@16:             >,
Chris@16:             mpl::int_<113>, mpl::int_<0> >::type
Chris@16:           >::type tag_type;
Chris@16:       result = lgamma_small_imp<T>(z, T(z - 1), T(z - 2), tag_type(), pol, l);
Chris@16:    }
Chris@16:    else if((z >= 3) && (z < 100) && (std::numeric_limits<T>::max_exponent >= 1024))
Chris@16:    {
Chris@16:       // taking the log of tgamma reduces the error, no danger of overflow here:
Chris@16:       result = log(gamma_imp(z, pol, l));
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       // regular evaluation:
Chris@16:       T zgh = static_cast<T>(z + Lanczos::g() - boost::math::constants::half<T>());
Chris@16:       result = log(zgh) - 1;
Chris@16:       result *= z - 0.5f;
Chris@16:       result += log(Lanczos::lanczos_sum_expG_scaled(z));
Chris@16:    }
Chris@16: 
Chris@16:    if(sign)
Chris@16:       *sign = sresult;
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: //
Chris@16: // Incomplete gamma functions follow:
Chris@16: //
Chris@16: template <class T>
Chris@16: struct upper_incomplete_gamma_fract
Chris@16: {
Chris@16: private:
Chris@16:    T z, a;
Chris@16:    int k;
Chris@16: public:
Chris@16:    typedef std::pair<T,T> result_type;
Chris@16: 
Chris@16:    upper_incomplete_gamma_fract(T a1, T z1)
Chris@16:       : z(z1-a1+1), a(a1), k(0)
Chris@16:    {
Chris@16:    }
Chris@16: 
Chris@16:    result_type operator()()
Chris@16:    {
Chris@16:       ++k;
Chris@16:       z += 2;
Chris@16:       return result_type(k * (a - k), z);
Chris@16:    }
Chris@16: };
Chris@16: 
Chris@16: template <class T>
Chris@16: inline T upper_gamma_fraction(T a, T z, T eps)
Chris@16: {
Chris@16:    // Multiply result by z^a * e^-z to get the full
Chris@16:    // upper incomplete integral.  Divide by tgamma(z)
Chris@16:    // to normalise.
Chris@16:    upper_incomplete_gamma_fract<T> f(a, z);
Chris@16:    return 1 / (z - a + 1 + boost::math::tools::continued_fraction_a(f, eps));
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: struct lower_incomplete_gamma_series
Chris@16: {
Chris@16: private:
Chris@16:    T a, z, result;
Chris@16: public:
Chris@16:    typedef T result_type;
Chris@16:    lower_incomplete_gamma_series(T a1, T z1) : a(a1), z(z1), result(1){}
Chris@16: 
Chris@16:    T operator()()
Chris@16:    {
Chris@16:       T r = result;
Chris@16:       a += 1;
Chris@16:       result *= z/a;
Chris@16:       return r;
Chris@16:    }
Chris@16: };
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline T lower_gamma_series(T a, T z, const Policy& pol, T init_value = 0)
Chris@16: {
Chris@16:    // Multiply result by ((z^a) * (e^-z) / a) to get the full
Chris@16:    // lower incomplete integral. Then divide by tgamma(a)
Chris@16:    // to get the normalised value.
Chris@16:    lower_incomplete_gamma_series<T> s(a, z);
Chris@16:    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
Chris@16:    T factor = policies::get_epsilon<T, Policy>();
Chris@16:    T result = boost::math::tools::sum_series(s, factor, max_iter, init_value);
Chris@16:    policies::check_series_iterations<T>("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol);
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: //
Chris@101: // Fully generic tgamma and lgamma use Stirling's approximation
Chris@101: // with Bernoulli numbers.
Chris@16: //
Chris@101: template<class T>
Chris@101: std::size_t highest_bernoulli_index()
Chris@101: {
Chris@101:    const float digits10_of_type = (std::numeric_limits<T>::is_specialized
Chris@101:                                       ? static_cast<float>(std::numeric_limits<T>::digits10)
Chris@101:                                       : static_cast<float>(boost::math::tools::digits<T>() * 0.301F));
Chris@101: 
Chris@101:    // Find the high index n for Bn to produce the desired precision in Stirling's calculation.
Chris@101:    return static_cast<std::size_t>(18.0F + (0.6F * digits10_of_type));
Chris@101: }
Chris@101: 
Chris@101: template<class T>
Chris@101: T minimum_argument_for_bernoulli_recursion()
Chris@101: {
Chris@101:    const float digits10_of_type = (std::numeric_limits<T>::is_specialized
Chris@101:                                       ? static_cast<float>(std::numeric_limits<T>::digits10)
Chris@101:                                       : static_cast<float>(boost::math::tools::digits<T>() * 0.301F));
Chris@101: 
Chris@101:    return T(digits10_of_type * 1.7F);
Chris@101: }
Chris@101: 
Chris@101: // Forward declaration of the lgamma_imp template specialization.
Chris@16: template <class T, class Policy>
Chris@101: T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign = 0);
Chris@101: 
Chris@101: template <class T, class Policy>
Chris@101: T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&)
Chris@16: {
Chris@101:    BOOST_MATH_STD_USING
Chris@101: 
Chris@16:    static const char* function = "boost::math::tgamma<%1%>(%1%)";
Chris@101: 
Chris@101:    // Check if the argument of tgamma is identically zero.
Chris@101:    const bool is_at_zero = (z == 0);
Chris@101: 
Chris@101:    if(is_at_zero)
Chris@101:       return policies::raise_domain_error<T>(function, "Evaluation of tgamma at zero %1%.", z, pol);
Chris@101: 
Chris@101:    const bool b_neg = (z < 0);
Chris@101: 
Chris@101:    const bool floor_of_z_is_equal_to_z = (floor(z) == z);
Chris@101: 
Chris@101:    // Special case handling of small factorials:
Chris@101:    if((!b_neg) && floor_of_z_is_equal_to_z && (z < boost::math::max_factorial<T>::value))
Chris@16:    {
Chris@101:       return boost::math::unchecked_factorial<T>(itrunc(z) - 1);
Chris@16:    }
Chris@101: 
Chris@101:    // Make a local, unsigned copy of the input argument.
Chris@101:    T zz((!b_neg) ? z : -z);
Chris@101: 
Chris@101:    // Special case for ultra-small z:
Chris@101:    if(zz < tools::cbrt_epsilon<T>())
Chris@16:    {
Chris@101:       const T a0(1);
Chris@101:       const T a1(boost::math::constants::euler<T>());
Chris@101:       const T six_euler_squared((boost::math::constants::euler<T>() * boost::math::constants::euler<T>()) * 6);
Chris@101:       const T a2((six_euler_squared -  boost::math::constants::pi_sqr<T>()) / 12);
Chris@101: 
Chris@101:       const T inverse_tgamma_series = z * ((a2 * z + a1) * z + a0);
Chris@101: 
Chris@101:       return 1 / inverse_tgamma_series;
Chris@16:    }
Chris@101: 
Chris@101:    // Scale the argument up for the calculation of lgamma,
Chris@101:    // and use downward recursion later for the final result.
Chris@101:    const T min_arg_for_recursion = minimum_argument_for_bernoulli_recursion<T>();
Chris@101: 
Chris@101:    int n_recur;
Chris@101: 
Chris@101:    if(zz < min_arg_for_recursion)
Chris@16:    {
Chris@101:       n_recur = boost::math::itrunc(min_arg_for_recursion - zz) + 1;
Chris@101: 
Chris@101:       zz += n_recur;
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@101:       n_recur = 0;
Chris@101:    }
Chris@101: 
Chris@101:    const T log_gamma_value = lgamma_imp(zz, pol, lanczos::undefined_lanczos());
Chris@101: 
Chris@101:    if(log_gamma_value > tools::log_max_value<T>())
Chris@101:       return policies::raise_overflow_error<T>(function, 0, pol);
Chris@101: 
Chris@101:    T gamma_value = exp(log_gamma_value);
Chris@101: 
Chris@101:    // Rescale the result using downward recursion if necessary.
Chris@101:    if(n_recur)
Chris@101:    {
Chris@101:       // The order of divides is important, if we keep subtracting 1 from zz
Chris@101:       // we DO NOT get back to z (cancellation error).  Further if z < epsilon
Chris@101:       // we would end up dividing by zero.  Also in order to prevent spurious
Chris@101:       // overflow with the first division, we must save dividing by |z| till last,
Chris@101:       // so the optimal order of divides is z+1, z+2, z+3...z+n_recur-1,z.
Chris@101:       zz = fabs(z) + 1;
Chris@101:       for(int k = 1; k < n_recur; ++k)
Chris@101:       {
Chris@101:          gamma_value /= zz;
Chris@101:          zz += 1;
Chris@101:       }
Chris@101:       gamma_value /= fabs(z);
Chris@101:    }
Chris@101: 
Chris@101:    // Return the result, accounting for possible negative arguments.
Chris@101:    if(b_neg)
Chris@101:    {
Chris@101:       // Provide special error analysis for:
Chris@101:       // * arguments in the neighborhood of a negative integer
Chris@101:       // * arguments exactly equal to a negative integer.
Chris@101: 
Chris@101:       // Check if the argument of tgamma is exactly equal to a negative integer.
Chris@101:       if(floor_of_z_is_equal_to_z)
Chris@101:          return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
Chris@101: 
Chris@101:       gamma_value *= sinpx(z);
Chris@101: 
Chris@101:       BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value);
Chris@101: 
Chris@101:       const bool result_is_too_large_to_represent = (   (abs(gamma_value) < 1)
Chris@101:                                                      && ((tools::max_value<T>() * abs(gamma_value)) < boost::math::constants::pi<T>()));
Chris@101: 
Chris@101:       if(result_is_too_large_to_represent)
Chris@16:          return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
Chris@101: 
Chris@101:       gamma_value = -boost::math::constants::pi<T>() / gamma_value;
Chris@101:       BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value);
Chris@101: 
Chris@101:       if(gamma_value == 0)
Chris@101:          return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
Chris@101: 
Chris@101:       if((boost::math::fpclassify)(gamma_value) == static_cast<int>(FP_SUBNORMAL))
Chris@101:          return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", gamma_value, pol);
Chris@16:    }
Chris@101: 
Chris@101:    return gamma_value;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@101: T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    static const char* function = "boost::math::lgamma<%1%>(%1%)";
Chris@101: 
Chris@101:    // Check if the argument of lgamma is identically zero.
Chris@101:    const bool is_at_zero = (z == 0);
Chris@101: 
Chris@101:    if(is_at_zero)
Chris@101:       return policies::raise_domain_error<T>(function, "Evaluation of lgamma at zero %1%.", z, pol);
Chris@101: 
Chris@101:    const bool b_neg = (z < 0);
Chris@101: 
Chris@101:    const bool floor_of_z_is_equal_to_z = (floor(z) == z);
Chris@101: 
Chris@101:    // Special case handling of small factorials:
Chris@101:    if((!b_neg) && floor_of_z_is_equal_to_z && (z < boost::math::max_factorial<T>::value))
Chris@16:    {
Chris@101:       return log(boost::math::unchecked_factorial<T>(itrunc(z) - 1));
Chris@101:    }
Chris@101: 
Chris@101:    // Make a local, unsigned copy of the input argument.
Chris@101:    T zz((!b_neg) ? z : -z);
Chris@101: 
Chris@101:    const T min_arg_for_recursion = minimum_argument_for_bernoulli_recursion<T>();
Chris@101: 
Chris@101:    T log_gamma_value;
Chris@101: 
Chris@101:    if (zz < min_arg_for_recursion)
Chris@101:    {
Chris@101: 	   // Here we simply take the logarithm of tgamma(). This is somewhat
Chris@101: 	   // inefficient, but simple. The rationale is that the argument here
Chris@101: 	   // is relatively small and overflow is not expected to be likely.
Chris@101:       if (z > -tools::root_epsilon<T>())
Chris@101:       {
Chris@101:          // Reflection formula may fail if z is very close to zero, let the series
Chris@101:          // expansion for tgamma close to zero do the work:
Chris@101:          log_gamma_value = log(abs(gamma_imp(z, pol, lanczos::undefined_lanczos())));
Chris@101:          if (sign)
Chris@101:          {
Chris@101:              *sign = z < 0 ? -1 : 1;
Chris@101:          }
Chris@101:          return log_gamma_value;
Chris@101:       }
Chris@101: 	   else
Chris@101:       {
Chris@101:          // No issue with spurious overflow in reflection formula, 
Chris@101:          // just fall through to regular code:
Chris@101:          log_gamma_value = log(abs(gamma_imp(zz, pol, lanczos::undefined_lanczos())));
Chris@101:       }
Chris@101:    }
Chris@101:    else
Chris@101:    {
Chris@101:       // Perform the Bernoulli series expansion of Stirling's approximation.
Chris@101: 
Chris@101:       const std::size_t number_of_bernoullis_b2n = highest_bernoulli_index<T>();
Chris@101: 
Chris@101:             T one_over_x_pow_two_n_minus_one = 1 / zz;
Chris@101:       const T one_over_x2                    = one_over_x_pow_two_n_minus_one * one_over_x_pow_two_n_minus_one;
Chris@101:             T sum                            = (boost::math::bernoulli_b2n<T>(1) / 2) * one_over_x_pow_two_n_minus_one;
Chris@101:       const T target_epsilon_to_break_loop   = (sum * boost::math::tools::epsilon<T>()) * T(1.0E-10F);
Chris@101: 
Chris@101:       for(std::size_t n = 2U; n < number_of_bernoullis_b2n; ++n)
Chris@101:       {
Chris@101:          one_over_x_pow_two_n_minus_one *= one_over_x2;
Chris@101: 
Chris@101:          const std::size_t n2 = static_cast<std::size_t>(n * 2U);
Chris@101: 
Chris@101:          const T term = (boost::math::bernoulli_b2n<T>(static_cast<int>(n)) * one_over_x_pow_two_n_minus_one) / (n2 * (n2 - 1U));
Chris@101: 
Chris@101:          if((n >= 8U) && (abs(term) < target_epsilon_to_break_loop))
Chris@101:          {
Chris@101:             // We have reached the desired precision in Stirling's expansion.
Chris@101:             // Adding additional terms to the sum of this divergent asymptotic
Chris@101:             // expansion will not improve the result.
Chris@101: 
Chris@101:             // Break from the loop.
Chris@101:             break;
Chris@101:          }
Chris@101: 
Chris@101:          sum += term;
Chris@101:       }
Chris@101: 
Chris@101:       // Complete Stirling's approximation.
Chris@101:       const T half_ln_two_pi = log(boost::math::constants::two_pi<T>()) / 2;
Chris@101: 
Chris@101:       log_gamma_value = ((((zz - boost::math::constants::half<T>()) * log(zz)) - zz) + half_ln_two_pi) + sum;
Chris@101:    }
Chris@101: 
Chris@101:    int sign_of_result = 1;
Chris@101: 
Chris@101:    if(b_neg)
Chris@101:    {
Chris@101:       // Provide special error analysis if the argument is exactly
Chris@101:       // equal to a negative integer.
Chris@101: 
Chris@101:       // Check if the argument of lgamma is exactly equal to a negative integer.
Chris@101:       if(floor_of_z_is_equal_to_z)
Chris@101:          return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%.", z, pol);
Chris@101: 
Chris@101:       T t = sinpx(z);
Chris@101: 
Chris@16:       if(t < 0)
Chris@16:       {
Chris@16:          t = -t;
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@101:          sign_of_result = -sign_of_result;
Chris@16:       }
Chris@101: 
Chris@101:       log_gamma_value = - log_gamma_value
Chris@101:                         + log(boost::math::constants::pi<T>())
Chris@101:                         - log(t);
Chris@16:    }
Chris@101: 
Chris@101:    if(sign != static_cast<int*>(0U)) { *sign = sign_of_result; }
Chris@101: 
Chris@101:    return log_gamma_value;
Chris@16: }
Chris@101: 
Chris@16: //
Chris@16: // This helper calculates tgamma(dz+1)-1 without cancellation errors,
Chris@16: // used by the upper incomplete gamma with z < 1:
Chris@16: //
Chris@16: template <class T, class Policy, class Lanczos>
Chris@16: T tgammap1m1_imp(T dz, Policy const& pol, const Lanczos& l)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    typedef typename policies::precision<T,Policy>::type precision_type;
Chris@16: 
Chris@16:    typedef typename mpl::if_<
Chris@16:       mpl::or_<
Chris@16:          mpl::less_equal<precision_type, mpl::int_<0> >,
Chris@16:          mpl::greater<precision_type, mpl::int_<113> >
Chris@16:       >,
Chris@16:       typename mpl::if_<
Chris@16:          is_same<Lanczos, lanczos::lanczos24m113>,
Chris@16:          mpl::int_<113>,
Chris@16:          mpl::int_<0>
Chris@16:       >::type,
Chris@16:       typename mpl::if_<
Chris@16:          mpl::less_equal<precision_type, mpl::int_<64> >,
Chris@16:          mpl::int_<64>, mpl::int_<113> >::type
Chris@16:        >::type tag_type;
Chris@16: 
Chris@16:    T result;
Chris@16:    if(dz < 0)
Chris@16:    {
Chris@16:       if(dz < -0.5)
Chris@16:       {
Chris@16:          // Best method is simply to subtract 1 from tgamma:
Chris@16:          result = boost::math::tgamma(1+dz, pol) - 1;
Chris@16:          BOOST_MATH_INSTRUMENT_CODE(result);
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          // Use expm1 on lgamma:
Chris@16:          result = boost::math::expm1(-boost::math::log1p(dz, pol) 
Chris@16:             + lgamma_small_imp<T>(dz+2, dz + 1, dz, tag_type(), pol, l));
Chris@16:          BOOST_MATH_INSTRUMENT_CODE(result);
Chris@16:       }
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       if(dz < 2)
Chris@16:       {
Chris@16:          // Use expm1 on lgamma:
Chris@16:          result = boost::math::expm1(lgamma_small_imp<T>(dz+1, dz, dz-1, tag_type(), pol, l), pol);
Chris@16:          BOOST_MATH_INSTRUMENT_CODE(result);
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          // Best method is simply to subtract 1 from tgamma:
Chris@16:          result = boost::math::tgamma(1+dz, pol) - 1;
Chris@16:          BOOST_MATH_INSTRUMENT_CODE(result);
Chris@16:       }
Chris@16:    }
Chris@16: 
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline T tgammap1m1_imp(T dz, Policy const& pol,
Chris@16:                  const ::boost::math::lanczos::undefined_lanczos& l)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING // ADL of std names
Chris@16:    //
Chris@16:    // There should be a better solution than this, but the
Chris@16:    // algebra isn't easy for the general case....
Chris@16:    // Start by subracting 1 from tgamma:
Chris@16:    //
Chris@16:    T result = gamma_imp(T(1 + dz), pol, l) - 1;
Chris@16:    BOOST_MATH_INSTRUMENT_CODE(result);
Chris@16:    //
Chris@16:    // Test the level of cancellation error observed: we loose one bit
Chris@16:    // for each power of 2 the result is less than 1.  If we would get
Chris@16:    // more bits from our most precise lgamma rational approximation, 
Chris@16:    // then use that instead:
Chris@16:    //
Chris@16:    BOOST_MATH_INSTRUMENT_CODE((dz > -0.5));
Chris@16:    BOOST_MATH_INSTRUMENT_CODE((dz < 2));
Chris@16:    BOOST_MATH_INSTRUMENT_CODE((ldexp(1.0, boost::math::policies::digits<T, Policy>()) * fabs(result) < 1e34));
Chris@16:    if((dz > -0.5) && (dz < 2) && (ldexp(1.0, boost::math::policies::digits<T, Policy>()) * fabs(result) < 1e34))
Chris@16:    {
Chris@16:       result = tgammap1m1_imp(dz, pol, boost::math::lanczos::lanczos24m113());
Chris@16:       BOOST_MATH_INSTRUMENT_CODE(result);
Chris@16:    }
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: //
Chris@16: // Series representation for upper fraction when z is small:
Chris@16: //
Chris@16: template <class T>
Chris@16: struct small_gamma2_series
Chris@16: {
Chris@16:    typedef T result_type;
Chris@16: 
Chris@16:    small_gamma2_series(T a_, T x_) : result(-x_), x(-x_), apn(a_+1), n(1){}
Chris@16: 
Chris@16:    T operator()()
Chris@16:    {
Chris@16:       T r = result / (apn);
Chris@16:       result *= x;
Chris@16:       result /= ++n;
Chris@16:       apn += 1;
Chris@16:       return r;
Chris@16:    }
Chris@16: 
Chris@16: private:
Chris@16:    T result, x, apn;
Chris@16:    int n;
Chris@16: };
Chris@16: //
Chris@16: // calculate power term prefix (z^a)(e^-z) used in the non-normalised
Chris@16: // incomplete gammas:
Chris@16: //
Chris@16: template <class T, class Policy>
Chris@16: T full_igamma_prefix(T a, T z, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    T prefix;
Chris@16:    T alz = a * log(z);
Chris@16: 
Chris@16:    if(z >= 1)
Chris@16:    {
Chris@16:       if((alz < tools::log_max_value<T>()) && (-z > tools::log_min_value<T>()))
Chris@16:       {
Chris@16:          prefix = pow(z, a) * exp(-z);
Chris@16:       }
Chris@16:       else if(a >= 1)
Chris@16:       {
Chris@16:          prefix = pow(z / exp(z/a), a);
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          prefix = exp(alz - z);
Chris@16:       }
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       if(alz > tools::log_min_value<T>())
Chris@16:       {
Chris@16:          prefix = pow(z, a) * exp(-z);
Chris@16:       }
Chris@16:       else if(z/a < tools::log_max_value<T>())
Chris@16:       {
Chris@16:          prefix = pow(z / exp(z/a), a);
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          prefix = exp(alz - z);
Chris@16:       }
Chris@16:    }
Chris@16:    //
Chris@16:    // This error handling isn't very good: it happens after the fact
Chris@16:    // rather than before it...
Chris@16:    //
Chris@16:    if((boost::math::fpclassify)(prefix) == (int)FP_INFINITE)
Chris@101:       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);
Chris@16: 
Chris@16:    return prefix;
Chris@16: }
Chris@16: //
Chris@16: // Compute (z^a)(e^-z)/tgamma(a)
Chris@16: // most if the error occurs in this function:
Chris@16: //
Chris@16: template <class T, class Policy, class Lanczos>
Chris@16: T regularised_gamma_prefix(T a, T z, const Policy& pol, const Lanczos& l)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    T agh = a + static_cast<T>(Lanczos::g()) - T(0.5);
Chris@16:    T prefix;
Chris@16:    T d = ((z - a) - static_cast<T>(Lanczos::g()) + T(0.5)) / agh;
Chris@16: 
Chris@16:    if(a < 1)
Chris@16:    {
Chris@16:       //
Chris@16:       // We have to treat a < 1 as a special case because our Lanczos
Chris@16:       // approximations are optimised against the factorials with a > 1,
Chris@16:       // and for high precision types especially (128-bit reals for example)
Chris@16:       // very small values of a can give rather eroneous results for gamma
Chris@16:       // unless we do this:
Chris@16:       //
Chris@16:       // TODO: is this still required?  Lanczos approx should be better now?
Chris@16:       //
Chris@16:       if(z <= tools::log_min_value<T>())
Chris@16:       {
Chris@16:          // Oh dear, have to use logs, should be free of cancellation errors though:
Chris@16:          return exp(a * log(z) - z - lgamma_imp(a, pol, l));
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          // direct calculation, no danger of overflow as gamma(a) < 1/a
Chris@16:          // for small a.
Chris@16:          return pow(z, a) * exp(-z) / gamma_imp(a, pol, l);
Chris@16:       }
Chris@16:    }
Chris@16:    else if((fabs(d*d*a) <= 100) && (a > 150))
Chris@16:    {
Chris@16:       // special case for large a and a ~ z.
Chris@16:       prefix = a * boost::math::log1pmx(d, pol) + z * static_cast<T>(0.5 - Lanczos::g()) / agh;
Chris@16:       prefix = exp(prefix);
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       //
Chris@16:       // general case.
Chris@16:       // direct computation is most accurate, but use various fallbacks
Chris@16:       // for different parts of the problem domain:
Chris@16:       //
Chris@16:       T alz = a * log(z / agh);
Chris@16:       T amz = a - z;
Chris@16:       if(((std::min)(alz, amz) <= tools::log_min_value<T>()) || ((std::max)(alz, amz) >= tools::log_max_value<T>()))
Chris@16:       {
Chris@16:          T amza = amz / a;
Chris@16:          if(((std::min)(alz, amz)/2 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/2 < tools::log_max_value<T>()))
Chris@16:          {
Chris@16:             // compute square root of the result and then square it:
Chris@16:             T sq = pow(z / agh, a / 2) * exp(amz / 2);
Chris@16:             prefix = sq * sq;
Chris@16:          }
Chris@16:          else if(((std::min)(alz, amz)/4 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/4 < tools::log_max_value<T>()) && (z > a))
Chris@16:          {
Chris@16:             // compute the 4th root of the result then square it twice:
Chris@16:             T sq = pow(z / agh, a / 4) * exp(amz / 4);
Chris@16:             prefix = sq * sq;
Chris@16:             prefix *= prefix;
Chris@16:          }
Chris@16:          else if((amza > tools::log_min_value<T>()) && (amza < tools::log_max_value<T>()))
Chris@16:          {
Chris@16:             prefix = pow((z * exp(amza)) / agh, a);
Chris@16:          }
Chris@16:          else
Chris@16:          {
Chris@16:             prefix = exp(alz + amz);
Chris@16:          }
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          prefix = pow(z / agh, a) * exp(amz);
Chris@16:       }
Chris@16:    }
Chris@16:    prefix *= sqrt(agh / boost::math::constants::e<T>()) / Lanczos::lanczos_sum_expG_scaled(a);
Chris@16:    return prefix;
Chris@16: }
Chris@16: //
Chris@16: // And again, without Lanczos support:
Chris@16: //
Chris@16: template <class T, class Policy>
Chris@16: T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined_lanczos&)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    T limit = (std::max)(T(10), a);
Chris@16:    T sum = detail::lower_gamma_series(a, limit, pol) / a;
Chris@16:    sum += detail::upper_gamma_fraction(a, limit, ::boost::math::policies::get_epsilon<T, Policy>());
Chris@16: 
Chris@16:    if(a < 10)
Chris@16:    {
Chris@16:       // special case for small a:
Chris@16:       T prefix = pow(z / 10, a);
Chris@16:       prefix *= exp(10-z);
Chris@16:       if(0 == prefix)
Chris@16:       {
Chris@16:          prefix = pow((z * exp((10-z)/a)) / 10, a);
Chris@16:       }
Chris@16:       prefix /= sum;
Chris@16:       return prefix;
Chris@16:    }
Chris@16: 
Chris@16:    T zoa = z / a;
Chris@16:    T amz = a - z;
Chris@16:    T alzoa = a * log(zoa);
Chris@16:    T prefix;
Chris@16:    if(((std::min)(alzoa, amz) <= tools::log_min_value<T>()) || ((std::max)(alzoa, amz) >= tools::log_max_value<T>()))
Chris@16:    {
Chris@16:       T amza = amz / a;
Chris@16:       if((amza <= tools::log_min_value<T>()) || (amza >= tools::log_max_value<T>()))
Chris@16:       {
Chris@16:          prefix = exp(alzoa + amz);
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          prefix = pow(zoa * exp(amza), a);
Chris@16:       }
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       prefix = pow(zoa, a) * exp(amz);
Chris@16:    }
Chris@16:    prefix /= sum;
Chris@16:    return prefix;
Chris@16: }
Chris@16: //
Chris@16: // Upper gamma fraction for very small a:
Chris@16: //
Chris@16: template <class T, class Policy>
Chris@16: inline T tgamma_small_upper_part(T a, T x, const Policy& pol, T* pgam = 0, bool invert = false, T* pderivative = 0)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING  // ADL of std functions.
Chris@16:    //
Chris@16:    // Compute the full upper fraction (Q) when a is very small:
Chris@16:    //
Chris@16:    T result;
Chris@16:    result = boost::math::tgamma1pm1(a, pol);
Chris@16:    if(pgam)
Chris@16:       *pgam = (result + 1) / a;
Chris@16:    T p = boost::math::powm1(x, a, pol);
Chris@16:    result -= p;
Chris@16:    result /= a;
Chris@16:    detail::small_gamma2_series<T> s(a, x);
Chris@16:    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>() - 10;
Chris@16:    p += 1;
Chris@16:    if(pderivative)
Chris@16:       *pderivative = p / (*pgam * exp(x));
Chris@16:    T init_value = invert ? *pgam : 0;
Chris@16:    result = -p * tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, (init_value - result) / p);
Chris@16:    policies::check_series_iterations<T>("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)", max_iter, pol);
Chris@16:    if(invert)
Chris@16:       result = -result;
Chris@16:    return result;
Chris@16: }
Chris@16: //
Chris@16: // Upper gamma fraction for integer a:
Chris@16: //
Chris@16: template <class T, class Policy>
Chris@16: inline T finite_gamma_q(T a, T x, Policy const& pol, T* pderivative = 0)
Chris@16: {
Chris@16:    //
Chris@16:    // Calculates normalised Q when a is an integer:
Chris@16:    //
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    T e = exp(-x);
Chris@16:    T sum = e;
Chris@16:    if(sum != 0)
Chris@16:    {
Chris@16:       T term = sum;
Chris@16:       for(unsigned n = 1; n < a; ++n)
Chris@16:       {
Chris@16:          term /= n;
Chris@16:          term *= x;
Chris@16:          sum += term;
Chris@16:       }
Chris@16:    }
Chris@16:    if(pderivative)
Chris@16:    {
Chris@16:       *pderivative = e * pow(x, a) / boost::math::unchecked_factorial<T>(itrunc(T(a - 1), pol));
Chris@16:    }
Chris@16:    return sum;
Chris@16: }
Chris@16: //
Chris@16: // Upper gamma fraction for half integer a:
Chris@16: //
Chris@16: template <class T, class Policy>
Chris@16: T finite_half_gamma_q(T a, T x, T* p_derivative, const Policy& pol)
Chris@16: {
Chris@16:    //
Chris@16:    // Calculates normalised Q when a is a half-integer:
Chris@16:    //
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    T e = boost::math::erfc(sqrt(x), pol);
Chris@16:    if((e != 0) && (a > 1))
Chris@16:    {
Chris@16:       T term = exp(-x) / sqrt(constants::pi<T>() * x);
Chris@16:       term *= x;
Chris@16:       static const T half = T(1) / 2;
Chris@16:       term /= half;
Chris@16:       T sum = term;
Chris@16:       for(unsigned n = 2; n < a; ++n)
Chris@16:       {
Chris@16:          term /= n - half;
Chris@16:          term *= x;
Chris@16:          sum += term;
Chris@16:       }
Chris@16:       e += sum;
Chris@16:       if(p_derivative)
Chris@16:       {
Chris@16:          *p_derivative = 0;
Chris@16:       }
Chris@16:    }
Chris@16:    else if(p_derivative)
Chris@16:    {
Chris@16:       // We'll be dividing by x later, so calculate derivative * x:
Chris@16:       *p_derivative = sqrt(x) * exp(-x) / constants::root_pi<T>();
Chris@16:    }
Chris@16:    return e;
Chris@16: }
Chris@16: //
Chris@16: // Main incomplete gamma entry point, handles all four incomplete gamma's:
Chris@16: //
Chris@16: template <class T, class Policy>
Chris@16: T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, 
Chris@16:                        const Policy& pol, T* p_derivative)
Chris@16: {
Chris@16:    static const char* function = "boost::math::gamma_p<%1%>(%1%, %1%)";
Chris@16:    if(a <= 0)
Chris@101:       return policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
Chris@16:    if(x < 0)
Chris@101:       return policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
Chris@16: 
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@16:    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
Chris@16: 
Chris@16:    T result = 0; // Just to avoid warning C4701: potentially uninitialized local variable 'result' used
Chris@16: 
Chris@101:    if(a >= max_factorial<T>::value && !normalised)
Chris@101:    {
Chris@101:       //
Chris@101:       // When we're computing the non-normalized incomplete gamma
Chris@101:       // and a is large the result is rather hard to compute unless
Chris@101:       // we use logs.  There are really two options - if x is a long
Chris@101:       // way from a in value then we can reliably use methods 2 and 4
Chris@101:       // below in logarithmic form and go straight to the result.
Chris@101:       // Otherwise we let the regularized gamma take the strain
Chris@101:       // (the result is unlikely to unerflow in the central region anyway)
Chris@101:       // and combine with lgamma in the hopes that we get a finite result.
Chris@101:       //
Chris@101:       if(invert && (a * 4 < x))
Chris@101:       {
Chris@101:          // This is method 4 below, done in logs:
Chris@101:          result = a * log(x) - x;
Chris@101:          if(p_derivative)
Chris@101:             *p_derivative = exp(result);
Chris@101:          result += log(upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>()));
Chris@101:       }
Chris@101:       else if(!invert && (a > 4 * x))
Chris@101:       {
Chris@101:          // This is method 2 below, done in logs:
Chris@101:          result = a * log(x) - x;
Chris@101:          if(p_derivative)
Chris@101:             *p_derivative = exp(result);
Chris@101:          T init_value = 0;
Chris@101:          result += log(detail::lower_gamma_series(a, x, pol, init_value) / a);
Chris@101:       }
Chris@101:       else
Chris@101:       {
Chris@101:          result = gamma_incomplete_imp(a, x, true, invert, pol, p_derivative);
Chris@101:          if(result == 0)
Chris@101:          {
Chris@101:             if(invert)
Chris@101:             {
Chris@101:                // Try http://functions.wolfram.com/06.06.06.0039.01
Chris@101:                result = 1 + 1 / (12 * a) + 1 / (288 * a * a);
Chris@101:                result = log(result) - a + (a - 0.5f) * log(a) + log(boost::math::constants::root_two_pi<T>());
Chris@101:                if(p_derivative)
Chris@101:                   *p_derivative = exp(a * log(x) - x);
Chris@101:             }
Chris@101:             else
Chris@101:             {
Chris@101:                // This is method 2 below, done in logs, we're really outside the
Chris@101:                // range of this method, but since the result is almost certainly
Chris@101:                // infinite, we should probably be OK:
Chris@101:                result = a * log(x) - x;
Chris@101:                if(p_derivative)
Chris@101:                   *p_derivative = exp(result);
Chris@101:                T init_value = 0;
Chris@101:                result += log(detail::lower_gamma_series(a, x, pol, init_value) / a);
Chris@101:             }
Chris@101:          }
Chris@101:          else
Chris@101:          {
Chris@101:             result = log(result) + boost::math::lgamma(a, pol);
Chris@101:          }
Chris@101:       }
Chris@101:       if(result > tools::log_max_value<T>())
Chris@101:          return policies::raise_overflow_error<T>(function, 0, pol);
Chris@101:       return exp(result);
Chris@101:    }
Chris@101: 
Chris@16:    BOOST_ASSERT((p_derivative == 0) || (normalised == true));
Chris@16: 
Chris@16:    bool is_int, is_half_int;
Chris@16:    bool is_small_a = (a < 30) && (a <= x + 1) && (x < tools::log_max_value<T>());
Chris@16:    if(is_small_a)
Chris@16:    {
Chris@16:       T fa = floor(a);
Chris@16:       is_int = (fa == a);
Chris@16:       is_half_int = is_int ? false : (fabs(fa - a) == 0.5f);
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       is_int = is_half_int = false;
Chris@16:    }
Chris@16: 
Chris@16:    int eval_method;
Chris@16:    
Chris@16:    if(is_int && (x > 0.6))
Chris@16:    {
Chris@16:       // calculate Q via finite sum:
Chris@16:       invert = !invert;
Chris@16:       eval_method = 0;
Chris@16:    }
Chris@16:    else if(is_half_int && (x > 0.2))
Chris@16:    {
Chris@16:       // calculate Q via finite sum for half integer a:
Chris@16:       invert = !invert;
Chris@16:       eval_method = 1;
Chris@16:    }
Chris@101:    else if((x < tools::root_epsilon<T>()) && (a > 1))
Chris@101:    {
Chris@101:       eval_method = 6;
Chris@101:    }
Chris@16:    else if(x < 0.5)
Chris@16:    {
Chris@16:       //
Chris@16:       // Changeover criterion chosen to give a changeover at Q ~ 0.33
Chris@16:       //
Chris@16:       if(-0.4 / log(x) < a)
Chris@16:       {
Chris@16:          eval_method = 2;
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          eval_method = 3;
Chris@16:       }
Chris@16:    }
Chris@16:    else if(x < 1.1)
Chris@16:    {
Chris@16:       //
Chris@16:       // Changover here occurs when P ~ 0.75 or Q ~ 0.25:
Chris@16:       //
Chris@16:       if(x * 0.75f < a)
Chris@16:       {
Chris@16:          eval_method = 2;
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          eval_method = 3;
Chris@16:       }
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       //
Chris@16:       // Begin by testing whether we're in the "bad" zone
Chris@16:       // where the result will be near 0.5 and the usual
Chris@16:       // series and continued fractions are slow to converge:
Chris@16:       //
Chris@16:       bool use_temme = false;
Chris@16:       if(normalised && std::numeric_limits<T>::is_specialized && (a > 20))
Chris@16:       {
Chris@16:          T sigma = fabs((x-a)/a);
Chris@16:          if((a > 200) && (policies::digits<T, Policy>() <= 113))
Chris@16:          {
Chris@16:             //
Chris@16:             // This limit is chosen so that we use Temme's expansion
Chris@16:             // only if the result would be larger than about 10^-6.
Chris@16:             // Below that the regular series and continued fractions
Chris@16:             // converge OK, and if we use Temme's method we get increasing
Chris@16:             // errors from the dominant erfc term as it's (inexact) argument
Chris@16:             // increases in magnitude.
Chris@16:             //
Chris@16:             if(20 / a > sigma * sigma)
Chris@16:                use_temme = true;
Chris@16:          }
Chris@16:          else if(policies::digits<T, Policy>() <= 64)
Chris@16:          {
Chris@16:             // Note in this zone we can't use Temme's expansion for 
Chris@16:             // types longer than an 80-bit real:
Chris@16:             // it would require too many terms in the polynomials.
Chris@16:             if(sigma < 0.4)
Chris@16:                use_temme = true;
Chris@16:          }
Chris@16:       }
Chris@16:       if(use_temme)
Chris@16:       {
Chris@16:          eval_method = 5;
Chris@16:       }
Chris@16:       else
Chris@16:       {
Chris@16:          //
Chris@16:          // Regular case where the result will not be too close to 0.5.
Chris@16:          //
Chris@16:          // Changeover here occurs at P ~ Q ~ 0.5
Chris@16:          // Note that series computation of P is about x2 faster than continued fraction
Chris@16:          // calculation of Q, so try and use the CF only when really necessary, especially
Chris@16:          // for small x.
Chris@16:          //
Chris@16:          if(x - (1 / (3 * x)) < a)
Chris@16:          {
Chris@16:             eval_method = 2;
Chris@16:          }
Chris@16:          else
Chris@16:          {
Chris@16:             eval_method = 4;
Chris@16:             invert = !invert;
Chris@16:          }
Chris@16:       }
Chris@16:    }
Chris@16: 
Chris@16:    switch(eval_method)
Chris@16:    {
Chris@16:    case 0:
Chris@16:       {
Chris@16:          result = finite_gamma_q(a, x, pol, p_derivative);
Chris@16:          if(normalised == false)
Chris@16:             result *= boost::math::tgamma(a, pol);
Chris@16:          break;
Chris@16:       }
Chris@16:    case 1:
Chris@16:       {
Chris@16:          result = finite_half_gamma_q(a, x, p_derivative, pol);
Chris@16:          if(normalised == false)
Chris@16:             result *= boost::math::tgamma(a, pol);
Chris@16:          if(p_derivative && (*p_derivative == 0))
Chris@16:             *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
Chris@16:          break;
Chris@16:       }
Chris@16:    case 2:
Chris@16:       {
Chris@16:          // Compute P:
Chris@16:          result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
Chris@16:          if(p_derivative)
Chris@16:             *p_derivative = result;
Chris@16:          if(result != 0)
Chris@16:          {
Chris@101:             //
Chris@101:             // If we're going to be inverting the result then we can
Chris@101:             // reduce the number of series evaluations by quite
Chris@101:             // a few iterations if we set an initial value for the
Chris@101:             // series sum based on what we'll end up subtracting it from
Chris@101:             // at the end.
Chris@101:             // Have to be careful though that this optimization doesn't 
Chris@101:             // lead to spurious numberic overflow.  Note that the
Chris@101:             // scary/expensive overflow checks below are more often
Chris@101:             // than not bypassed in practice for "sensible" input
Chris@101:             // values:
Chris@101:             //
Chris@16:             T init_value = 0;
Chris@101:             bool optimised_invert = false;
Chris@16:             if(invert)
Chris@16:             {
Chris@101:                init_value = (normalised ? 1 : boost::math::tgamma(a, pol));
Chris@101:                if(normalised || (result >= 1) || (tools::max_value<T>() * result > init_value))
Chris@101:                {
Chris@101:                   init_value /= result;
Chris@101:                   if(normalised || (a < 1) || (tools::max_value<T>() / a > init_value))
Chris@101:                   {
Chris@101:                      init_value *= -a;
Chris@101:                      optimised_invert = true;
Chris@101:                   }
Chris@101:                   else
Chris@101:                      init_value = 0;
Chris@101:                }
Chris@101:                else
Chris@101:                   init_value = 0;
Chris@16:             }
Chris@16:             result *= detail::lower_gamma_series(a, x, pol, init_value) / a;
Chris@101:             if(optimised_invert)
Chris@16:             {
Chris@16:                invert = false;
Chris@16:                result = -result;
Chris@16:             }
Chris@16:          }
Chris@16:          break;
Chris@16:       }
Chris@16:    case 3:
Chris@16:       {
Chris@16:          // Compute Q:
Chris@16:          invert = !invert;
Chris@16:          T g;
Chris@16:          result = tgamma_small_upper_part(a, x, pol, &g, invert, p_derivative);
Chris@16:          invert = false;
Chris@16:          if(normalised)
Chris@16:             result /= g;
Chris@16:          break;
Chris@16:       }
Chris@16:    case 4:
Chris@16:       {
Chris@16:          // Compute Q:
Chris@16:          result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
Chris@16:          if(p_derivative)
Chris@16:             *p_derivative = result;
Chris@16:          if(result != 0)
Chris@16:             result *= upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>());
Chris@16:          break;
Chris@16:       }
Chris@16:    case 5:
Chris@16:       {
Chris@16:          //
Chris@16:          // Use compile time dispatch to the appropriate
Chris@16:          // Temme asymptotic expansion.  This may be dead code
Chris@16:          // if T does not have numeric limits support, or has
Chris@16:          // too many digits for the most precise version of
Chris@16:          // these expansions, in that case we'll be calling
Chris@16:          // an empty function.
Chris@16:          //
Chris@16:          typedef typename policies::precision<T, Policy>::type precision_type;
Chris@16: 
Chris@16:          typedef typename mpl::if_<
Chris@16:             mpl::or_<mpl::equal_to<precision_type, mpl::int_<0> >,
Chris@16:             mpl::greater<precision_type, mpl::int_<113> > >,
Chris@16:             mpl::int_<0>,
Chris@16:             typename mpl::if_<
Chris@16:                mpl::less_equal<precision_type, mpl::int_<53> >,
Chris@16:                mpl::int_<53>,
Chris@16:                typename mpl::if_<
Chris@16:                   mpl::less_equal<precision_type, mpl::int_<64> >,
Chris@16:                   mpl::int_<64>,
Chris@16:                   mpl::int_<113>
Chris@16:                >::type
Chris@16:             >::type
Chris@16:          >::type tag_type;
Chris@16: 
Chris@16:          result = igamma_temme_large(a, x, pol, static_cast<tag_type const*>(0));
Chris@16:          if(x >= a)
Chris@16:             invert = !invert;
Chris@16:          if(p_derivative)
Chris@16:             *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
Chris@16:          break;
Chris@16:       }
Chris@101:    case 6:
Chris@101:       {
Chris@101:          // x is so small that P is necessarily very small too,
Chris@101:          // use http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/
Chris@101:          result = !normalised ? pow(x, a) / (a) : pow(x, a) / boost::math::tgamma(a + 1, pol);
Chris@101:          result *= 1 - a * x / (a + 1);
Chris@101:       }
Chris@16:    }
Chris@16: 
Chris@16:    if(normalised && (result > 1))
Chris@16:       result = 1;
Chris@16:    if(invert)
Chris@16:    {
Chris@16:       T gam = normalised ? 1 : boost::math::tgamma(a, pol);
Chris@16:       result = gam - result;
Chris@16:    }
Chris@16:    if(p_derivative)
Chris@16:    {
Chris@16:       //
Chris@16:       // Need to convert prefix term to derivative:
Chris@16:       //
Chris@16:       if((x < 1) && (tools::max_value<T>() * x < *p_derivative))
Chris@16:       {
Chris@16:          // overflow, just return an arbitrarily large value:
Chris@16:          *p_derivative = tools::max_value<T>() / 2;
Chris@16:       }
Chris@16: 
Chris@16:       *p_derivative /= x;
Chris@16:    }
Chris@16: 
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: //
Chris@16: // Ratios of two gamma functions:
Chris@16: //
Chris@16: template <class T, class Policy, class Lanczos>
Chris@101: T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const Lanczos& l)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@101:    if(z < tools::epsilon<T>())
Chris@101:    {
Chris@101:       //
Chris@101:       // We get spurious numeric overflow unless we're very careful, this
Chris@101:       // can occur either inside Lanczos::lanczos_sum(z) or in the
Chris@101:       // final combination of terms, to avoid this, split the product up
Chris@101:       // into 2 (or 3) parts:
Chris@101:       //
Chris@101:       // G(z) / G(L) = 1 / (z * G(L)) ; z < eps, L = z + delta = delta
Chris@101:       //    z * G(L) = z * G(lim) * (G(L)/G(lim)) ; lim = largest factorial
Chris@101:       //
Chris@101:       if(boost::math::max_factorial<T>::value < delta)
Chris@101:       {
Chris@101:          T ratio = tgamma_delta_ratio_imp_lanczos(delta, T(boost::math::max_factorial<T>::value - delta), pol, l);
Chris@101:          ratio *= z;
Chris@101:          ratio *= boost::math::unchecked_factorial<T>(boost::math::max_factorial<T>::value - 1);
Chris@101:          return 1 / ratio;
Chris@101:       }
Chris@101:       else
Chris@101:       {
Chris@101:          return 1 / (z * boost::math::tgamma(z + delta, pol));
Chris@101:       }
Chris@101:    }
Chris@16:    T zgh = z + Lanczos::g() - constants::half<T>();
Chris@16:    T result;
Chris@16:    if(fabs(delta) < 10)
Chris@16:    {
Chris@16:       result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol));
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       result = pow(zgh / (zgh + delta), z - constants::half<T>());
Chris@16:    }
Chris@101:    // Split the calculation up to avoid spurious overflow:
Chris@101:    result *= Lanczos::lanczos_sum(z) / Lanczos::lanczos_sum(T(z + delta));
Chris@16:    result *= pow(constants::e<T>() / (zgh + delta), delta);
Chris@16:    return result;
Chris@16: }
Chris@16: //
Chris@16: // And again without Lanczos support this time:
Chris@16: //
Chris@16: template <class T, class Policy>
Chris@16: T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const lanczos::undefined_lanczos&)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    //
Chris@16:    // The upper gamma fraction is *very* slow for z < 6, actually it's very
Chris@16:    // slow to converge everywhere but recursing until z > 6 gets rid of the
Chris@16:    // worst of it's behaviour.
Chris@16:    //
Chris@16:    T prefix = 1;
Chris@16:    T zd = z + delta;
Chris@16:    while((zd < 6) && (z < 6))
Chris@16:    {
Chris@16:       prefix /= z;
Chris@16:       prefix *= zd;
Chris@16:       z += 1;
Chris@16:       zd += 1;
Chris@16:    }
Chris@16:    if(delta < 10)
Chris@16:    {
Chris@16:       prefix *= exp(-z * boost::math::log1p(delta / z, pol));
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       prefix *= pow(z / zd, z);
Chris@16:    }
Chris@16:    prefix *= pow(constants::e<T>() / zd, delta);
Chris@16:    T sum = detail::lower_gamma_series(z, z, pol) / z;
Chris@16:    sum += detail::upper_gamma_fraction(z, z, ::boost::math::policies::get_epsilon<T, Policy>());
Chris@16:    T sumd = detail::lower_gamma_series(zd, zd, pol) / zd;
Chris@16:    sumd += detail::upper_gamma_fraction(zd, zd, ::boost::math::policies::get_epsilon<T, Policy>());
Chris@16:    sum /= sumd;
Chris@16:    if(fabs(tools::max_value<T>() / prefix) < fabs(sum))
Chris@16:       return policies::raise_overflow_error<T>("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Result of tgamma is too large to represent.", pol);
Chris@16:    return sum * prefix;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@101:    if((z <= 0) || (z + delta <= 0))
Chris@101:    {
Chris@101:       // This isn't very sofisticated, or accurate, but it does work:
Chris@101:       return boost::math::tgamma(z, pol) / boost::math::tgamma(z + delta, pol);
Chris@101:    }
Chris@16: 
Chris@16:    if(floor(delta) == delta)
Chris@16:    {
Chris@16:       if(floor(z) == z)
Chris@16:       {
Chris@16:          //
Chris@16:          // Both z and delta are integers, see if we can just use table lookup
Chris@16:          // of the factorials to get the result:
Chris@16:          //
Chris@16:          if((z <= max_factorial<T>::value) && (z + delta <= max_factorial<T>::value))
Chris@16:          {
Chris@16:             return unchecked_factorial<T>((unsigned)itrunc(z, pol) - 1) / unchecked_factorial<T>((unsigned)itrunc(T(z + delta), pol) - 1);
Chris@16:          }
Chris@16:       }
Chris@16:       if(fabs(delta) < 20)
Chris@16:       {
Chris@16:          //
Chris@16:          // delta is a small integer, we can use a finite product:
Chris@16:          //
Chris@16:          if(delta == 0)
Chris@16:             return 1;
Chris@16:          if(delta < 0)
Chris@16:          {
Chris@16:             z -= 1;
Chris@16:             T result = z;
Chris@16:             while(0 != (delta += 1))
Chris@16:             {
Chris@16:                z -= 1;
Chris@16:                result *= z;
Chris@16:             }
Chris@16:             return result;
Chris@16:          }
Chris@16:          else
Chris@16:          {
Chris@16:             T result = 1 / z;
Chris@16:             while(0 != (delta -= 1))
Chris@16:             {
Chris@16:                z += 1;
Chris@16:                result /= z;
Chris@16:             }
Chris@16:             return result;
Chris@16:          }
Chris@16:       }
Chris@16:    }
Chris@16:    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
Chris@16:    return tgamma_delta_ratio_imp_lanczos(z, delta, pol, lanczos_type());
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: T tgamma_ratio_imp(T x, T y, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16: 
Chris@101:    if((x <= 0) || (boost::math::isinf)(x))
Chris@101:       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);
Chris@101:    if((y <= 0) || (boost::math::isinf)(y))
Chris@101:       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);
Chris@101: 
Chris@101:    if(x <= tools::min_value<T>())
Chris@101:    {
Chris@101:       // Special case for denorms...Ugh.
Chris@101:       T shift = ldexp(T(1), tools::digits<T>());
Chris@101:       return shift * tgamma_ratio_imp(T(x * shift), y, pol);
Chris@101:    }
Chris@16: 
Chris@16:    if((x < max_factorial<T>::value) && (y < max_factorial<T>::value))
Chris@16:    {
Chris@16:       // Rather than subtracting values, lets just call the gamma functions directly:
Chris@16:       return boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
Chris@16:    }
Chris@16:    T prefix = 1;
Chris@16:    if(x < 1)
Chris@16:    {
Chris@16:       if(y < 2 * max_factorial<T>::value)
Chris@16:       {
Chris@16:          // We need to sidestep on x as well, otherwise we'll underflow
Chris@16:          // before we get to factor in the prefix term:
Chris@16:          prefix /= x;
Chris@16:          x += 1;
Chris@16:          while(y >=  max_factorial<T>::value)
Chris@16:          {
Chris@16:             y -= 1;
Chris@16:             prefix /= y;
Chris@16:          }
Chris@16:          return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
Chris@16:       }
Chris@16:       //
Chris@16:       // result is almost certainly going to underflow to zero, try logs just in case:
Chris@16:       //
Chris@16:       return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
Chris@16:    }
Chris@16:    if(y < 1)
Chris@16:    {
Chris@16:       if(x < 2 * max_factorial<T>::value)
Chris@16:       {
Chris@16:          // We need to sidestep on y as well, otherwise we'll overflow
Chris@16:          // before we get to factor in the prefix term:
Chris@16:          prefix *= y;
Chris@16:          y += 1;
Chris@16:          while(x >= max_factorial<T>::value)
Chris@16:          {
Chris@16:             x -= 1;
Chris@16:             prefix *= x;
Chris@16:          }
Chris@16:          return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
Chris@16:       }
Chris@16:       //
Chris@16:       // Result will almost certainly overflow, try logs just in case:
Chris@16:       //
Chris@16:       return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
Chris@16:    }
Chris@16:    //
Chris@16:    // Regular case, x and y both large and similar in magnitude:
Chris@16:    //
Chris@16:    return boost::math::tgamma_delta_ratio(x, y - x, pol);
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: T gamma_p_derivative_imp(T a, T x, const Policy& pol)
Chris@16: {
Chris@101:    BOOST_MATH_STD_USING
Chris@16:    //
Chris@16:    // Usual error checks first:
Chris@16:    //
Chris@16:    if(a <= 0)
Chris@101:       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);
Chris@16:    if(x < 0)
Chris@101:       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);
Chris@16:    //
Chris@16:    // Now special cases:
Chris@16:    //
Chris@16:    if(x == 0)
Chris@16:    {
Chris@16:       return (a > 1) ? 0 :
Chris@16:          (a == 1) ? 1 : policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
Chris@16:    }
Chris@16:    //
Chris@16:    // Normal case:
Chris@16:    //
Chris@16:    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
Chris@16:    T f1 = detail::regularised_gamma_prefix(a, x, pol, lanczos_type());
Chris@16:    if((x < 1) && (tools::max_value<T>() * x < f1))
Chris@16:    {
Chris@16:       // overflow:
Chris@16:       return policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
Chris@16:    }
Chris@101:    if(f1 == 0)
Chris@101:    {
Chris@101:       // Underflow in calculation, use logs instead:
Chris@101:       f1 = a * log(x) - x - lgamma(a, pol) - log(x);
Chris@101:       f1 = exp(f1);
Chris@101:    }
Chris@101:    else
Chris@101:       f1 /= x;
Chris@16: 
Chris@16:    return f1;
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline typename tools::promote_args<T>::type 
Chris@16:    tgamma(T z, const Policy& /* pol */, const mpl::true_)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename tools::promote_args<T>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16:    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%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: struct igamma_initializer
Chris@16: {
Chris@16:    struct init
Chris@16:    {
Chris@16:       init()
Chris@16:       {
Chris@16:          typedef typename policies::precision<T, Policy>::type precision_type;
Chris@16: 
Chris@16:          typedef typename mpl::if_<
Chris@16:             mpl::or_<mpl::equal_to<precision_type, mpl::int_<0> >,
Chris@16:             mpl::greater<precision_type, mpl::int_<113> > >,
Chris@16:             mpl::int_<0>,
Chris@16:             typename mpl::if_<
Chris@16:                mpl::less_equal<precision_type, mpl::int_<53> >,
Chris@16:                mpl::int_<53>,
Chris@16:                typename mpl::if_<
Chris@16:                   mpl::less_equal<precision_type, mpl::int_<64> >,
Chris@16:                   mpl::int_<64>,
Chris@16:                   mpl::int_<113>
Chris@16:                >::type
Chris@16:             >::type
Chris@16:          >::type tag_type;
Chris@16: 
Chris@16:          do_init(tag_type());
Chris@16:       }
Chris@16:       template <int N>
Chris@16:       static void do_init(const mpl::int_<N>&)
Chris@16:       {
Chris@16:          boost::math::gamma_p(static_cast<T>(400), static_cast<T>(400), Policy());
Chris@16:       }
Chris@16:       static void do_init(const mpl::int_<53>&){}
Chris@16:       void force_instantiate()const{}
Chris@16:    };
Chris@16:    static const init initializer;
Chris@16:    static void force_instantiate()
Chris@16:    {
Chris@16:       initializer.force_instantiate();
Chris@16:    }
Chris@16: };
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: const typename igamma_initializer<T, Policy>::init igamma_initializer<T, Policy>::initializer;
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: struct lgamma_initializer
Chris@16: {
Chris@16:    struct init
Chris@16:    {
Chris@16:       init()
Chris@16:       {
Chris@16:          typedef typename policies::precision<T, Policy>::type precision_type;
Chris@16:          typedef typename mpl::if_<
Chris@16:             mpl::and_<
Chris@16:                mpl::less_equal<precision_type, mpl::int_<64> >, 
Chris@16:                mpl::greater<precision_type, mpl::int_<0> > 
Chris@16:             >,
Chris@16:             mpl::int_<64>,
Chris@16:             typename mpl::if_<
Chris@16:                mpl::and_<
Chris@16:                   mpl::less_equal<precision_type, mpl::int_<113> >,
Chris@16:                   mpl::greater<precision_type, mpl::int_<0> > 
Chris@16:                >,
Chris@16:                mpl::int_<113>, mpl::int_<0> >::type
Chris@16:              >::type tag_type;
Chris@16:          do_init(tag_type());
Chris@16:       }
Chris@16:       static void do_init(const mpl::int_<64>&)
Chris@16:       {
Chris@16:          boost::math::lgamma(static_cast<T>(2.5), Policy());
Chris@16:          boost::math::lgamma(static_cast<T>(1.25), Policy());
Chris@16:          boost::math::lgamma(static_cast<T>(1.75), Policy());
Chris@16:       }
Chris@16:       static void do_init(const mpl::int_<113>&)
Chris@16:       {
Chris@16:          boost::math::lgamma(static_cast<T>(2.5), Policy());
Chris@16:          boost::math::lgamma(static_cast<T>(1.25), Policy());
Chris@16:          boost::math::lgamma(static_cast<T>(1.5), Policy());
Chris@16:          boost::math::lgamma(static_cast<T>(1.75), Policy());
Chris@16:       }
Chris@16:       static void do_init(const mpl::int_<0>&)
Chris@16:       {
Chris@16:       }
Chris@16:       void force_instantiate()const{}
Chris@16:    };
Chris@16:    static const init initializer;
Chris@16:    static void force_instantiate()
Chris@16:    {
Chris@16:       initializer.force_instantiate();
Chris@16:    }
Chris@16: };
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: const typename lgamma_initializer<T, Policy>::init lgamma_initializer<T, Policy>::initializer;
Chris@16: 
Chris@16: template <class T1, class T2, class Policy>
Chris@16: inline typename tools::promote_args<T1, T2>::type
Chris@16:    tgamma(T1 a, T2 z, const Policy&, const mpl::false_)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename tools::promote_args<T1, T2>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    igamma_initializer<value_type, forwarding_policy>::force_instantiate();
Chris@16: 
Chris@16:    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
Chris@16:       detail::gamma_incomplete_imp(static_cast<value_type>(a),
Chris@16:       static_cast<value_type>(z), false, true,
Chris@16:       forwarding_policy(), static_cast<value_type*>(0)), "boost::math::tgamma<%1%>(%1%, %1%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T1, class T2>
Chris@16: inline typename tools::promote_args<T1, T2>::type
Chris@16:    tgamma(T1 a, T2 z, const mpl::false_ tag)
Chris@16: {
Chris@16:    return tgamma(a, z, policies::policy<>(), tag);
Chris@16: }
Chris@16: 
Chris@16: 
Chris@16: } // namespace detail
Chris@16: 
Chris@16: template <class T>
Chris@16: inline typename tools::promote_args<T>::type 
Chris@16:    tgamma(T z)
Chris@16: {
Chris@16:    return tgamma(z, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline typename tools::promote_args<T>::type 
Chris@16:    lgamma(T z, int* sign, const Policy&)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename tools::promote_args<T>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    detail::lgamma_initializer<value_type, forwarding_policy>::force_instantiate();
Chris@16: 
Chris@16:    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%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline typename tools::promote_args<T>::type 
Chris@16:    lgamma(T z, int* sign)
Chris@16: {
Chris@16:    return lgamma(z, sign, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline typename tools::promote_args<T>::type 
Chris@16:    lgamma(T x, const Policy& pol)
Chris@16: {
Chris@16:    return ::boost::math::lgamma(x, 0, pol);
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline typename tools::promote_args<T>::type 
Chris@16:    lgamma(T x)
Chris@16: {
Chris@16:    return ::boost::math::lgamma(x, 0, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: inline typename tools::promote_args<T>::type 
Chris@16:    tgamma1pm1(T z, const Policy& /* pol */)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename tools::promote_args<T>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    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%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T>
Chris@16: inline typename tools::promote_args<T>::type 
Chris@16:    tgamma1pm1(T z)
Chris@16: {
Chris@16:    return tgamma1pm1(z, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: //
Chris@16: // Full upper incomplete gamma:
Chris@16: //
Chris@16: template <class T1, class T2>
Chris@16: inline typename tools::promote_args<T1, T2>::type
Chris@16:    tgamma(T1 a, T2 z)
Chris@16: {
Chris@16:    //
Chris@16:    // Type T2 could be a policy object, or a value, select the 
Chris@16:    // right overload based on T2:
Chris@16:    //
Chris@16:    typedef typename policies::is_policy<T2>::type maybe_policy;
Chris@16:    return detail::tgamma(a, z, maybe_policy());
Chris@16: }
Chris@16: template <class T1, class T2, class Policy>
Chris@16: inline typename tools::promote_args<T1, T2>::type
Chris@16:    tgamma(T1 a, T2 z, const Policy& pol)
Chris@16: {
Chris@16:    return detail::tgamma(a, z, pol, mpl::false_());
Chris@16: }
Chris@16: //
Chris@16: // Full lower incomplete gamma:
Chris@16: //
Chris@16: template <class T1, class T2, class Policy>
Chris@16: inline typename tools::promote_args<T1, T2>::type
Chris@16:    tgamma_lower(T1 a, T2 z, const Policy&)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename tools::promote_args<T1, T2>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
Chris@16: 
Chris@16:    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
Chris@16:       detail::gamma_incomplete_imp(static_cast<value_type>(a),
Chris@16:       static_cast<value_type>(z), false, false,
Chris@16:       forwarding_policy(), static_cast<value_type*>(0)), "tgamma_lower<%1%>(%1%, %1%)");
Chris@16: }
Chris@16: template <class T1, class T2>
Chris@16: inline typename tools::promote_args<T1, T2>::type
Chris@16:    tgamma_lower(T1 a, T2 z)
Chris@16: {
Chris@16:    return tgamma_lower(a, z, policies::policy<>());
Chris@16: }
Chris@16: //
Chris@16: // Regularised upper incomplete gamma:
Chris@16: //
Chris@16: template <class T1, class T2, class Policy>
Chris@16: inline typename tools::promote_args<T1, T2>::type
Chris@16:    gamma_q(T1 a, T2 z, const Policy& /* pol */)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename tools::promote_args<T1, T2>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
Chris@16: 
Chris@16:    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
Chris@16:       detail::gamma_incomplete_imp(static_cast<value_type>(a),
Chris@16:       static_cast<value_type>(z), true, true,
Chris@16:       forwarding_policy(), static_cast<value_type*>(0)), "gamma_q<%1%>(%1%, %1%)");
Chris@16: }
Chris@16: template <class T1, class T2>
Chris@16: inline typename tools::promote_args<T1, T2>::type
Chris@16:    gamma_q(T1 a, T2 z)
Chris@16: {
Chris@16:    return gamma_q(a, z, policies::policy<>());
Chris@16: }
Chris@16: //
Chris@16: // Regularised lower incomplete gamma:
Chris@16: //
Chris@16: template <class T1, class T2, class Policy>
Chris@16: inline typename tools::promote_args<T1, T2>::type
Chris@16:    gamma_p(T1 a, T2 z, const Policy&)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename tools::promote_args<T1, T2>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
Chris@16: 
Chris@16:    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
Chris@16:       detail::gamma_incomplete_imp(static_cast<value_type>(a),
Chris@16:       static_cast<value_type>(z), true, false,
Chris@16:       forwarding_policy(), static_cast<value_type*>(0)), "gamma_p<%1%>(%1%, %1%)");
Chris@16: }
Chris@16: template <class T1, class T2>
Chris@16: inline typename tools::promote_args<T1, T2>::type
Chris@16:    gamma_p(T1 a, T2 z)
Chris@16: {
Chris@16:    return gamma_p(a, z, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: // ratios of gamma functions:
Chris@16: template <class T1, class T2, class Policy>
Chris@16: inline typename tools::promote_args<T1, T2>::type 
Chris@16:    tgamma_delta_ratio(T1 z, T2 delta, const Policy& /* pol */)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename tools::promote_args<T1, T2>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    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%)");
Chris@16: }
Chris@16: template <class T1, class T2>
Chris@16: inline typename tools::promote_args<T1, T2>::type 
Chris@16:    tgamma_delta_ratio(T1 z, T2 delta)
Chris@16: {
Chris@16:    return tgamma_delta_ratio(z, delta, policies::policy<>());
Chris@16: }
Chris@16: template <class T1, class T2, class Policy>
Chris@16: inline typename tools::promote_args<T1, T2>::type 
Chris@16:    tgamma_ratio(T1 a, T2 b, const Policy&)
Chris@16: {
Chris@16:    typedef typename tools::promote_args<T1, T2>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    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%)");
Chris@16: }
Chris@16: template <class T1, class T2>
Chris@16: inline typename tools::promote_args<T1, T2>::type 
Chris@16:    tgamma_ratio(T1 a, T2 b)
Chris@16: {
Chris@16:    return tgamma_ratio(a, b, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T1, class T2, class Policy>
Chris@16: inline typename tools::promote_args<T1, T2>::type 
Chris@16:    gamma_p_derivative(T1 a, T2 x, const Policy&)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    typedef typename tools::promote_args<T1, T2>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    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%)");
Chris@16: }
Chris@16: template <class T1, class T2>
Chris@16: inline typename tools::promote_args<T1, T2>::type 
Chris@16:    gamma_p_derivative(T1 a, T2 x)
Chris@16: {
Chris@16:    return gamma_p_derivative(a, x, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: } // namespace math
Chris@16: } // namespace boost
Chris@16: 
Chris@16: #ifdef BOOST_MSVC
Chris@16: # pragma warning(pop)
Chris@16: #endif
Chris@16: 
Chris@16: #include <boost/math/special_functions/detail/igamma_inverse.hpp>
Chris@16: #include <boost/math/special_functions/detail/gamma_inva.hpp>
Chris@16: #include <boost/math/special_functions/erf.hpp>
Chris@16: 
Chris@16: #endif // BOOST_MATH_SF_GAMMA_HPP
Chris@16: 
Chris@16: 
Chris@16: 
Chris@16: