text: //  Copyright John Maddock 2006-7.
  //  Copyright Paul A. Bristow 2007.
  
  //  Use, modification and distribution are subject to the
  //  Boost Software License, Version 1.0. (See accompanying file
  //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  
  #ifndef BOOST_MATH_SF_GAMMA_HPP
  #define BOOST_MATH_SF_GAMMA_HPP
  
  #ifdef _MSC_VER
  #pragma once
  #endif
  
  #include <boost/config.hpp>
  #ifdef BOOST_MSVC
  # pragma warning(push)
  # pragma warning(disable: 4127 4701)
  //  // For lexical_cast, until fixed in 1.35?
  //  // conditional expression is constant &
  //  // Potentially uninitialized local variable 'name' used
  #endif
  #include <boost/lexical_cast.hpp>
  #ifdef BOOST_MSVC
  # pragma warning(pop)
  #endif
  #include <boost/math/tools/series.hpp>
  #include <boost/math/tools/fraction.hpp>
  #include <boost/math/tools/precision.hpp>
  #include <boost/math/tools/promotion.hpp>
  #include <boost/math/policies/error_handling.hpp>
  #include <boost/math/constants/constants.hpp>
  #include <boost/math/special_functions/math_fwd.hpp>
  #include <boost/math/special_functions/log1p.hpp>
  #include <boost/math/special_functions/trunc.hpp>
  #include <boost/math/special_functions/powm1.hpp>
  #include <boost/math/special_functions/sqrt1pm1.hpp>
  #include <boost/math/special_functions/lanczos.hpp>
  #include <boost/math/special_functions/fpclassify.hpp>
  #include <boost/math/special_functions/detail/igamma_large.hpp>
  #include <boost/math/special_functions/detail/unchecked_factorial.hpp>
  #include <boost/math/special_functions/detail/lgamma_small.hpp>
  #include <boost/type_traits/is_convertible.hpp>
  #include <boost/assert.hpp>
  #include <boost/mpl/greater.hpp>
  #include <boost/mpl/equal_to.hpp>
  #include <boost/mpl/greater.hpp>
  
  #include <boost/config/no_tr1/cmath.hpp>
  #include <algorithm>
  
  #ifdef BOOST_MATH_INSTRUMENT
  #include <iostream>
  #include <iomanip>
  #include <typeinfo>
  #endif
  
  #ifdef BOOST_MSVC
  # pragma warning(push)
  # pragma warning(disable: 4702) // unreachable code (return after domain_error throw).
  # pragma warning(disable: 4127) // conditional expression is constant.
  # pragma warning(disable: 4100) // unreferenced formal parameter.
  // Several variables made comments,
  // but some difficulty as whether referenced on not may depend on macro values.
  // So to be safe, 4100 warnings suppressed.
  // TODO - revisit this?
  #endif
  
  namespace boost{ namespace math{
  
  namespace detail{
  
  template <class T>
  inline bool is_odd(T v, const boost::true_type&)
  {
     int i = static_cast<int>(v);
     return i&1;
  }
  template <class T>
  inline bool is_odd(T v, const boost::false_type&)
  {
     // Oh dear can't cast T to int!
     BOOST_MATH_STD_USING
     T modulus = v - 2 * floor(v/2);
     return static_cast<bool>(modulus != 0);
  }
  template <class T>
  inline bool is_odd(T v)
  {
     return is_odd(v, ::boost::is_convertible<T, int>());
  }
  
  template <class T>
  T sinpx(T z)
  {
     // Ad hoc function calculates x * sin(pi * x),
     // taking extra care near when x is near a whole number.
     BOOST_MATH_STD_USING
     int sign = 1;
     if(z < 0)
     {
        z = -z;
     }
     else
     {
        sign = -sign;
     }
     T fl = floor(z);
     T dist;
     if(is_odd(fl))
     {
        fl += 1;
        dist = fl - z;
        sign = -sign;
     }
     else
     {
        dist = z - fl;
     }
     BOOST_ASSERT(fl >= 0);
     if(dist > 0.5)
        dist = 1 - dist;
     T result = sin(dist*boost::math::constants::pi<T>());
     return sign*z*result;
  } // template <class T> T sinpx(T z)
  //
  // tgamma(z), with Lanczos support:
  //
  template <class T, class Policy, class L>
  T gamma_imp(T z, const Policy& pol, const L& l)
  {
     BOOST_MATH_STD_USING
  
     T result = 1;
  
  #ifdef BOOST_MATH_INSTRUMENT
     static bool b = false;
     if(!b)
     {
        std::cout << "tgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;
        b = true;
     }
  #endif
     static const char* function = "boost::math::tgamma<%1%>(%1%)";
  
     if(z <= 0)
     {
        if(floor(z) == z)
           return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
        if(z <= -20)
        {
           result = gamma_imp(T(-z), pol, l) * sinpx(z);
           if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
              return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
           result = -boost::math::constants::pi<T>() / result;
           if(result == 0)
              return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
           if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL)
              return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", result, pol);
           return result;
        }
  
        // shift z to > 1:
        while(z < 0)
        {
           result /= z;
           z += 1;
        }
     }
     if((floor(z) == z) && (z < max_factorial<T>::value))
     {
        result *= unchecked_factorial<T>(itrunc(z, pol) - 1);
     }
     else
     {
        result *= L::lanczos_sum(z);
        if(z * log(z) > tools::log_max_value<T>())
        {
           // we're going to overflow unless this is done with care:
           T zgh = (z + static_cast<T>(L::g()) - boost::math::constants::half<T>());
           if(log(zgh) * z / 2 > tools::log_max_value<T>())
              return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
           T hp = pow(zgh, (z / 2) - T(0.25));
           result *= hp / exp(zgh);
           if(tools::max_value<T>() / hp < result)
              return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
           result *= hp;
        }
        else
        {
           T zgh = (z + static_cast<T>(L::g()) - boost::math::constants::half<T>());
           result *= pow(zgh, z - boost::math::constants::half<T>()) / exp(zgh);
        }
     }
     return result;
  }
  //
  // lgamma(z) with Lanczos support:
  //
  template <class T, class Policy, class L>
  T lgamma_imp(T z, const Policy& pol, const L& l, int* sign = 0)
  {
  #ifdef BOOST_MATH_INSTRUMENT
     static bool b = false;
     if(!b)
     {
        std::cout << "lgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;
        b = true;
     }
  #endif
  
     BOOST_MATH_STD_USING
  
     static const char* function = "boost::math::lgamma<%1%>(%1%)";
  
     T result = 0;
     int sresult = 1;
     if(z <= 0)
     {
        // reflection formula:
        if(floor(z) == z)
           return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%.", z, pol);
  
        T t = sinpx(z);
        z = -z;
        if(t < 0)
        {
           t = -t;
        }
        else
        {
           sresult = -sresult;
        }
        result = log(boost::math::constants::pi<T>()) - lgamma_imp(z, pol, l) - log(t);
     }
     else if(z < 15)
     {
        typedef typename policies::precision<T, Policy>::type precision_type;
        typedef typename mpl::if_<
           mpl::and_<
              mpl::less_equal<precision_type, mpl::int_<64> >, 
              mpl::greater<precision_type, mpl::int_<0> > 
           >,
           mpl::int_<64>,
           typename mpl::if_<
              mpl::and_<
                 mpl::less_equal<precision_type, mpl::int_<113> >,
                 mpl::greater<precision_type, mpl::int_<0> > 
              >,
              mpl::int_<113>, mpl::int_<0> >::type
            >::type tag_type;
        result = lgamma_small_imp<T>(z, z - 1, z - 2, tag_type(), pol, l);
     }
     else if((z >= 3) && (z < 100))
     {
        // taking the log of tgamma reduces the error, no danger of overflow here:
        result = log(gamma_imp(z, pol, l));
     }
     else
     {
        // regular evaluation:
        T zgh = static_cast<T>(z + L::g() - boost::math::constants::half<T>());
        result = log(zgh) - 1;
        result *= z - 0.5f;
        result += log(L::lanczos_sum_expG_scaled(z));
     }
  
     if(sign)
        *sign = sresult;
     return result;
  }
  
  //
  // Incomplete gamma functions follow:
  //
  template <class T>
  struct upper_incomplete_gamma_fract
  {
  private:
     T z, a;
     int k;
  public:
     typedef std::pair<T,T> result_type;
  
     upper_incomplete_gamma_fract(T a1, T z1)
        : z(z1-a1+1), a(a1), k(0)
     {
     }
  
     result_type operator()()
     {
        ++k;
        z += 2;
        return result_type(k * (a - k), z);
     }
  };
  
  template <class T>
  inline T upper_gamma_fraction(T a, T z, T eps)
  {
     // Multiply result by z^a * e^-z to get the full
     // upper incomplete integral.  Divide by tgamma(z)
     // to normalise.
     upper_incomplete_gamma_fract<T> f(a, z);
     return 1 / (z - a + 1 + boost::math::tools::continued_fraction_a(f, eps));
  }
  
  template <class T>
  struct lower_incomplete_gamma_series
  {
  private:
     T a, z, result;
  public:
     typedef T result_type;
     lower_incomplete_gamma_series(T a1, T z1) : a(a1), z(z1), result(1){}
  
     T operator()()
     {
        T r = result;
        a += 1;
        result *= z/a;
        return r;
     }
  };
  
  template <class T, class Policy>
  inline T lower_gamma_series(T a, T z, const Policy& pol, T init_value = 0)
  {
     // Multiply result by ((z^a) * (e^-z) / a) to get the full
     // lower incomplete integral. Then divide by tgamma(a)
     // to get the normalised value.
     lower_incomplete_gamma_series<T> s(a, z);
     boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
     T factor = policies::get_epsilon<T, Policy>();
     T result = boost::math::tools::sum_series(s, factor, max_iter, init_value);
     policies::check_series_iterations("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol);
     return result;
  }
  
  //
  // Fully generic tgamma and lgamma use the incomplete partial
  // sums added together:
  //
  template <class T, class Policy>
  T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l)
  {
     static const char* function = "boost::math::tgamma<%1%>(%1%)";
     BOOST_MATH_STD_USING
     if((z <= 0) && (floor(z) == z))
        return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
     if(z <= -20)
     {
        T result = gamma_imp(-z, pol, l) * sinpx(z);
        if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
           return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
        result = -boost::math::constants::pi<T>() / result;
        if(result == 0)
           return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
        if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL)
           return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", result, pol);
        return result;
     }
     //
     // The upper gamma fraction is *very* slow for z < 6, actually it's very
     // slow to converge everywhere but recursing until z > 6 gets rid of the
     // worst of it's behaviour.
     //
     T prefix = 1;
     while(z < 6)
     {
        prefix /= z;
        z += 1;
     }
     BOOST_MATH_INSTRUMENT_CODE(prefix);
     if((floor(z) == z) && (z < max_factorial<T>::value))
     {
        prefix *= unchecked_factorial<T>(itrunc(z, pol) - 1);
     }
     else
     {
        prefix = prefix * pow(z / boost::math::constants::e<T>(), z);
        BOOST_MATH_INSTRUMENT_CODE(prefix);
        T sum = detail::lower_gamma_series(z, z, pol) / z;
        BOOST_MATH_INSTRUMENT_CODE(sum);
        sum += detail::upper_gamma_fraction(z, z, ::boost::math::policies::get_epsilon<T, Policy>());
        BOOST_MATH_INSTRUMENT_CODE(sum);
        if(fabs(tools::max_value<T>() / prefix) < fabs(sum))
           return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
        BOOST_MATH_INSTRUMENT_CODE((sum * prefix));
        return sum * prefix;
     }
     return prefix;
  }
  
  template <class T, class Policy>
  T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l, int*sign)
  {
     BOOST_MATH_STD_USING
  
     static const char* function = "boost::math::lgamma<%1%>(%1%)";
     T result = 0;
     int sresult = 1;
     if(z <= 0)
     {
        if(floor(z) == z)
           return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
        T t = detail::sinpx(z);
        z = -z;
        if(t < 0)
        {
           t = -t;
        }
        else
        {
           sresult = -sresult;
        }
        result = log(boost::math::constants::pi<T>()) - lgamma_imp(z, pol, l, 0) - log(t);
     }
     else if((z != 1) && (z != 2))
     {
        T limit = (std::max)(z+1, T(10));
        T prefix = z * log(limit) - limit;
        T sum = detail::lower_gamma_series(z, limit, pol) / z;
        sum += detail::upper_gamma_fraction(z, limit, ::boost::math::policies::get_epsilon<T, Policy>());
        result = log(sum) + prefix;
     }
     if(sign)
        *sign = sresult;
     return result;
  }
  //
  // This helper calculates tgamma(dz+1)-1 without cancellation errors,
  // used by the upper incomplete gamma with z < 1:
  //
  template <class T, class Policy, class L>
  T tgammap1m1_imp(T dz, Policy const& pol, const L& l)
  {
     BOOST_MATH_STD_USING
  
     typedef typename policies::precision<T,Policy>::type precision_type;
  
     typedef typename mpl::if_<
        mpl::or_<
           mpl::less_equal<precision_type, mpl::int_<0> >,
           mpl::greater<precision_type, mpl::int_<113> >
        >,
        typename mpl::if_<
           is_same<L, lanczos::lanczos24m113>,
           mpl::int_<113>,
           mpl::int_<0>
        >::type,
        typename mpl::if_<
           mpl::less_equal<precision_type, mpl::int_<64> >,
           mpl::int_<64>, mpl::int_<113> >::type
         >::type tag_type;
  
     T result;
     if(dz < 0)
     {
        if(dz < -0.5)
        {
           // Best method is simply to subtract 1 from tgamma:
           result = boost::math::tgamma(1+dz, pol) - 1;
           BOOST_MATH_INSTRUMENT_CODE(result);
        }
        else
        {
           // Use expm1 on lgamma:
           result = boost::math::expm1(-boost::math::log1p(dz, pol) 
              + lgamma_small_imp<T>(dz+2, dz + 1, dz, tag_type(), pol, l));
           BOOST_MATH_INSTRUMENT_CODE(result);
        }
     }
     else
     {
        if(dz < 2)
        {
           // Use expm1 on lgamma:
           result = boost::math::expm1(lgamma_small_imp<T>(dz+1, dz, dz-1, tag_type(), pol, l), pol);
           BOOST_MATH_INSTRUMENT_CODE(result);
        }
        else
        {
           // Best method is simply to subtract 1 from tgamma:
           result = boost::math::tgamma(1+dz, pol) - 1;
           BOOST_MATH_INSTRUMENT_CODE(result);
        }
     }
  
     return result;
  }
  
  template <class T, class Policy>
  inline T tgammap1m1_imp(T dz, Policy const& pol,
                   const ::boost::math::lanczos::undefined_lanczos& l)
  {
     BOOST_MATH_STD_USING // ADL of std names
     //
     // There should be a better solution than this, but the
     // algebra isn't easy for the general case....
     // Start by subracting 1 from tgamma:
     //
     T result = gamma_imp(1 + dz, pol, l) - 1;
     BOOST_MATH_INSTRUMENT_CODE(result);
     //
     // Test the level of cancellation error observed: we loose one bit
     // for each power of 2 the result is less than 1.  If we would get
     // more bits from our most precise lgamma rational approximation, 
     // then use that instead:
     //
     BOOST_MATH_INSTRUMENT_CODE((dz > -0.5));
     BOOST_MATH_INSTRUMENT_CODE((dz < 2));
     BOOST_MATH_INSTRUMENT_CODE((ldexp(1.0, boost::math::policies::digits<T, Policy>()) * fabs(result) < 1e34));
     if((dz > -0.5) && (dz < 2) && (ldexp(1.0, boost::math::policies::digits<T, Policy>()) * fabs(result) < 1e34))
     {
        result = tgammap1m1_imp(dz, pol, boost::math::lanczos::lanczos24m113());
        BOOST_MATH_INSTRUMENT_CODE(result);
     }
     return result;
  }
  
  //
  // Series representation for upper fraction when z is small:
  //
  template <class T>
  struct small_gamma2_series
  {
     typedef T result_type;
  
     small_gamma2_series(T a_, T x_) : result(-x_), x(-x_), apn(a_+1), n(1){}
  
     T operator()()
     {
        T r = result / (apn);
        result *= x;
        result /= ++n;
        apn += 1;
        return r;
     }
  
  private:
     T result, x, apn;
     int n;
  };
  //
  // calculate power term prefix (z^a)(e^-z) used in the non-normalised
  // incomplete gammas:
  //
  template <class T, class Policy>
  T full_igamma_prefix(T a, T z, const Policy& pol)
  {
     BOOST_MATH_STD_USING
  
     T prefix;
     T alz = a * log(z);
  
     if(z >= 1)
     {
        if((alz < tools::log_max_value<T>()) && (-z > tools::log_min_value<T>()))
        {
           prefix = pow(z, a) * exp(-z);
        }
        else if(a >= 1)
        {
           prefix = pow(z / exp(z/a), a);
        }
        else
        {
           prefix = exp(alz - z);
        }
     }
     else
     {
        if(alz > tools::log_min_value<T>())
        {
           prefix = pow(z, a) * exp(-z);
        }
        else if(z/a < tools::log_max_value<T>())
        {
           prefix = pow(z / exp(z/a), a);
        }
        else
        {
           prefix = exp(alz - z);
        }
     }
     //
     // This error handling isn't very good: it happens after the fact
     // rather than before it...
     //
     if((boost::math::fpclassify)(prefix) == (int)FP_INFINITE)
        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);
  
     return prefix;
  }
  //
  // Compute (z^a)(e^-z)/tgamma(a)
  // most if the error occurs in this function:
  //
  template <class T, class Policy, class L>
  T regularised_gamma_prefix(T a, T z, const Policy& pol, const L& l)
  {
     BOOST_MATH_STD_USING
     T agh = a + static_cast<T>(L::g()) - T(0.5);
     T prefix;
     T d = ((z - a) - static_cast<T>(L::g()) + T(0.5)) / agh;
  
     if(a < 1)
     {
        //
        // We have to treat a < 1 as a special case because our Lanczos
        // approximations are optimised against the factorials with a > 1,
        // and for high precision types especially (128-bit reals for example)
        // very small values of a can give rather eroneous results for gamma
        // unless we do this:
        //
        // TODO: is this still required?  Lanczos approx should be better now?
        //
        if(z <= tools::log_min_value<T>())
        {
           // Oh dear, have to use logs, should be free of cancellation errors though:
           return exp(a * log(z) - z - lgamma_imp(a, pol, l));
        }
        else
        {
           // direct calculation, no danger of overflow as gamma(a) < 1/a
           // for small a.
           return pow(z, a) * exp(-z) / gamma_imp(a, pol, l);
        }
     }
     else if((fabs(d*d*a) <= 100) && (a > 150))
     {
        // special case for large a and a ~ z.
        prefix = a * boost::math::log1pmx(d, pol) + z * static_cast<T>(0.5 - L::g()) / agh;
        prefix = exp(prefix);
     }
     else
     {
        //
        // general case.
        // direct computation is most accurate, but use various fallbacks
        // for different parts of the problem domain:
        //
        T alz = a * log(z / agh);
        T amz = a - z;
        if(((std::min)(alz, amz) <= tools::log_min_value<T>()) || ((std::max)(alz, amz) >= tools::log_max_value<T>()))
        {
           T amza = amz / a;
           if(((std::min)(alz, amz)/2 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/2 < tools::log_max_value<T>()))
           {
              // compute square root of the result and then square it:
              T sq = pow(z / agh, a / 2) * exp(amz / 2);
              prefix = sq * sq;
           }
           else if(((std::min)(alz, amz)/4 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/4 < tools::log_max_value<T>()) && (z > a))
           {
              // compute the 4th root of the result then square it twice:
              T sq = pow(z / agh, a / 4) * exp(amz / 4);
              prefix = sq * sq;
              prefix *= prefix;
           }
           else if((amza > tools::log_min_value<T>()) && (amza < tools::log_max_value<T>()))
           {
              prefix = pow((z * exp(amza)) / agh, a);
           }
           else
           {
              prefix = exp(alz + amz);
           }
        }
        else
        {
           prefix = pow(z / agh, a) * exp(amz);
        }
     }
     prefix *= sqrt(agh / boost::math::constants::e<T>()) / L::lanczos_sum_expG_scaled(a);
     return prefix;
  }
  //
  // And again, without Lanczos support:
  //
  template <class T, class Policy>
  T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined_lanczos&)
  {
     BOOST_MATH_STD_USING
  
     T limit = (std::max)(T(10), a);
     T sum = detail::lower_gamma_series(a, limit, pol) / a;
     sum += detail::upper_gamma_fraction(a, limit, ::boost::math::policies::get_epsilon<T, Policy>());
  
     if(a < 10)
     {
        // special case for small a:
        T prefix = pow(z / 10, a);
        prefix *= exp(10-z);
        if(0 == prefix)
        {
           prefix = pow((z * exp((10-z)/a)) / 10, a);
        }
        prefix /= sum;
        return prefix;
     }
  
     T zoa = z / a;
     T amz = a - z;
     T alzoa = a * log(zoa);
     T prefix;
     if(((std::min)(alzoa, amz) <= tools::log_min_value<T>()) || ((std::max)(alzoa, amz) >= tools::log_max_value<T>()))
     {
        T amza = amz / a;
        if((amza <= tools::log_min_value<T>()) || (amza >= tools::log_max_value<T>()))
        {
           prefix = exp(alzoa + amz);
        }
        else
        {
           prefix = pow(zoa * exp(amza), a);
        }
     }
     else
     {
        prefix = pow(zoa, a) * exp(amz);
     }
     prefix /= sum;
     return prefix;
  }
  //
  // Upper gamma fraction for very small a:
  //
  template <class T, class Policy>
  inline T tgamma_small_upper_part(T a, T x, const Policy& pol, T* pgam = 0, bool invert = false, T* pderivative = 0)
  {
     BOOST_MATH_STD_USING  // ADL of std functions.
     //
     // Compute the full upper fraction (Q) when a is very small:
     //
     T result;
     result = boost::math::tgamma1pm1(a, pol);
     if(pgam)
        *pgam = (result + 1) / a;
     T p = boost::math::powm1(x, a, pol);
     result -= p;
     result /= a;
     detail::small_gamma2_series<T> s(a, x);
     boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>() - 10;
     p += 1;
     if(pderivative)
        *pderivative = p / (*pgam * exp(x));
     T init_value = invert ? *pgam : 0;
     result = -p * tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, (init_value - result) / p);
     policies::check_series_iterations("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)", max_iter, pol);
     if(invert)
        result = -result;
     return result;
  }
  //
  // Upper gamma fraction for integer a:
  //
  template <class T, class Policy>
  inline T finite_gamma_q(T a, T x, Policy const& pol, T* pderivative = 0)
  {
     //
     // Calculates normalised Q when a is an integer:
     //
     BOOST_MATH_STD_USING
     T e = exp(-x);
     T sum = e;
     if(sum != 0)
     {
        T term = sum;
        for(unsigned n = 1; n < a; ++n)
        {
           term /= n;
           term *= x;
           sum += term;
        }
     }
     if(pderivative)
     {
        *pderivative = e * pow(x, a) / boost::math::unchecked_factorial<T>(itrunc(T(a - 1), pol));
     }
     return sum;
  }
  //
  // Upper gamma fraction for half integer a:
  //
  template <class T, class Policy>
  T finite_half_gamma_q(T a, T x, T* p_derivative, const Policy& pol)
  {
     //
     // Calculates normalised Q when a is a half-integer:
     //
     BOOST_MATH_STD_USING
     T e = boost::math::erfc(sqrt(x), pol);
     if((e != 0) && (a > 1))
     {
        T term = exp(-x) / sqrt(constants::pi<T>() * x);
        term *= x;
        static const T half = T(1) / 2;
        term /= half;
        T sum = term;
        for(unsigned n = 2; n < a; ++n)
        {
           term /= n - half;
           term *= x;
           sum += term;
        }
        e += sum;
        if(p_derivative)
        {
           *p_derivative = 0;
        }
     }
     else if(p_derivative)
     {
        // We'll be dividing by x later, so calculate derivative * x:
        *p_derivative = sqrt(x) * exp(-x) / constants::root_pi<T>();
     }
     return e;
  }
  //
  // Main incomplete gamma entry point, handles all four incomplete gamma's:
  //
  template <class T, class Policy>
  T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, 
                         const Policy& pol, T* p_derivative)
  {
     static const char* function = "boost::math::gamma_p<%1%>(%1%, %1%)";
     if(a <= 0)
        policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
     if(x < 0)
        policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
  
     BOOST_MATH_STD_USING
  
     typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
  
     T result = 0; // Just to avoid warning C4701: potentially uninitialized local variable 'result' used
  
     BOOST_ASSERT((p_derivative == 0) || (normalised == true));
  
     bool is_int, is_half_int;
     bool is_small_a = (a < 30) && (a <= x + 1);
     if(is_small_a)
     {
        T fa = floor(a);
        is_int = (fa == a);
        is_half_int = is_int ? false : (fabs(fa - a) == 0.5f);
     }
     else
     {
        is_int = is_half_int = false;
     }
  
     int eval_method;
     
     if(is_int && (x > 0.6))
     {
        // calculate Q via finite sum:
        invert = !invert;
        eval_method = 0;
     }
     else if(is_half_int && (x > 0.2))
     {
        // calculate Q via finite sum for half integer a:
        invert = !invert;
        eval_method = 1;
     }
     else if(x < 0.5)
     {
        //
        // Changeover criterion chosen to give a changeover at Q ~ 0.33
        //
        if(-0.4 / log(x) < a)
        {
           eval_method = 2;
        }
        else
        {
           eval_method = 3;
        }
     }
     else if(x < 1.1)
     {
        //
        // Changover here occurs when P ~ 0.75 or Q ~ 0.25:
        //
        if(x * 0.75f < a)
        {
           eval_method = 2;
        }
        else
        {
           eval_method = 3;
        }
     }
     else
     {
        //
        // Begin by testing whether we're in the "bad" zone
        // where the result will be near 0.5 and the usual
        // series and continued fractions are slow to converge:
        //
        bool use_temme = false;
        if(normalised && std::numeric_limits<T>::is_specialized && (a > 20))
        {
           T sigma = fabs((x-a)/a);
           if((a > 200) && (policies::digits<T, Policy>() <= 113))
           {
              //
              // This limit is chosen so that we use Temme's expansion
              // only if the result would be larger than about 10^-6.
              // Below that the regular series and continued fractions
              // converge OK, and if we use Temme's method we get increasing
              // errors from the dominant erfc term as it's (inexact) argument
              // increases in magnitude.
              //
              if(20 / a > sigma * sigma)
                 use_temme = true;
           }
           else if(policies::digits<T, Policy>() <= 64)
           {
              // Note in this zone we can't use Temme's expansion for 
              // types longer than an 80-bit real:
              // it would require too many terms in the polynomials.
              if(sigma < 0.4)
                 use_temme = true;
           }
        }
        if(use_temme)
        {
           eval_method = 5;
        }
        else
        {
           //
           // Regular case where the result will not be too close to 0.5.
           //
           // Changeover here occurs at P ~ Q ~ 0.5
           // Note that series computation of P is about x2 faster than continued fraction
           // calculation of Q, so try and use the CF only when really necessary, especially
           // for small x.
           //
           if(x - (1 / (3 * x)) < a)
           {
              eval_method = 2;
           }
           else
           {
              eval_method = 4;
              invert = !invert;
           }
        }
     }
  
     switch(eval_method)
     {
     case 0:
        {
           result = finite_gamma_q(a, x, pol, p_derivative);
           if(normalised == false)
              result *= boost::math::tgamma(a, pol);
           break;
        }
     case 1:
        {
           result = finite_half_gamma_q(a, x, p_derivative, pol);
           if(normalised == false)
              result *= boost::math::tgamma(a, pol);
           if(p_derivative && (*p_derivative == 0))
              *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
           break;
        }
     case 2:
        {
           // Compute P:
           result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
           if(p_derivative)
              *p_derivative = result;
           if(result != 0)
           {
              T init_value = 0;
              if(invert)
              {
                 init_value = -a * (normalised ? 1 : boost::math::tgamma(a, pol)) / result;
              }
              result *= detail::lower_gamma_series(a, x, pol, init_value) / a;
              if(invert)
              {
                 invert = false;
                 result = -result;
              }
           }
           break;
        }
     case 3:
        {
           // Compute Q:
           invert = !invert;
           T g;
           result = tgamma_small_upper_part(a, x, pol, &g, invert, p_derivative);
           invert = false;
           if(normalised)
              result /= g;
           break;
        }
     case 4:
        {
           // Compute Q:
           result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
           if(p_derivative)
              *p_derivative = result;
           if(result != 0)
              result *= upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>());
           break;
        }
     case 5:
        {
           //
           // Use compile time dispatch to the appropriate
           // Temme asymptotic expansion.  This may be dead code
           // if T does not have numeric limits support, or has
           // too many digits for the most precise version of
           // these expansions, in that case we'll be calling
           // an empty function.
           //
           typedef typename policies::precision<T, Policy>::type precision_type;
  
           typedef typename mpl::if_<
              mpl::or_<mpl::equal_to<precision_type, mpl::int_<0> >,
              mpl::greater<precision_type, mpl::int_<113> > >,
              mpl::int_<0>,
              typename mpl::if_<
                 mpl::less_equal<precision_type, mpl::int_<53> >,
                 mpl::int_<53>,
                 typename mpl::if_<
                    mpl::less_equal<precision_type, mpl::int_<64> >,
                    mpl::int_<64>,
                    mpl::int_<113>
                 >::type
              >::type
           >::type tag_type;
  
           result = igamma_temme_large(a, x, pol, static_cast<tag_type const*>(0));
           if(x >= a)
              invert = !invert;
           if(p_derivative)
              *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
           break;
        }
     }
  
     if(normalised && (result > 1))
        result = 1;
     if(invert)
     {
        T gam = normalised ? 1 : boost::math::tgamma(a, pol);
        result = gam - result;
     }
     if(p_derivative)
     {
        //
        // Need to convert prefix term to derivative:
        //
        if((x < 1) && (tools::max_value<T>() * x < *p_derivative))
        {
           // overflow, just return an arbitrarily large value:
           *p_derivative = tools::max_value<T>() / 2;
        }
  
        *p_derivative /= x;
     }
  
     return result;
  }
  
  //
  // Ratios of two gamma functions:
  //
  template <class T, class Policy, class L>
  T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const L&)
  {
     BOOST_MATH_STD_USING
     T zgh = z + L::g() - constants::half<T>();
     T result;
     if(fabs(delta) < 10)
     {
        result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol));
     }
     else
     {
        result = pow(zgh / (zgh + delta), z - constants::half<T>());
     }
     result *= pow(constants::e<T>() / (zgh + delta), delta);
     result *= L::lanczos_sum(z) / L::lanczos_sum(z + delta);
     return result;
  }
  //
  // And again without Lanczos support this time:
  //
  template <class T, class Policy>
  T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const lanczos::undefined_lanczos&)
  {
     BOOST_MATH_STD_USING
     //
     // The upper gamma fraction is *very* slow for z < 6, actually it's very
     // slow to converge everywhere but recursing until z > 6 gets rid of the
     // worst of it's behaviour.
     //
     T prefix = 1;
     T zd = z + delta;
     while((zd < 6) && (z < 6))
     {
        prefix /= z;
        prefix *= zd;
        z += 1;
        zd += 1;
     }
     if(delta < 10)
     {
        prefix *= exp(-z * boost::math::log1p(delta / z, pol));
     }
     else
     {
        prefix *= pow(z / zd, z);
     }
     prefix *= pow(constants::e<T>() / zd, delta);
     T sum = detail::lower_gamma_series(z, z, pol) / z;
     sum += detail::upper_gamma_fraction(z, z, ::boost::math::policies::get_epsilon<T, Policy>());
     T sumd = detail::lower_gamma_series(zd, zd, pol) / zd;
     sumd += detail::upper_gamma_fraction(zd, zd, ::boost::math::policies::get_epsilon<T, Policy>());
     sum /= sumd;
     if(fabs(tools::max_value<T>() / prefix) < fabs(sum))
        return policies::raise_overflow_error<T>("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Result of tgamma is too large to represent.", pol);
     return sum * prefix;
  }
  
  template <class T, class Policy>
  T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol)
  {
     BOOST_MATH_STD_USING
  
     if(z <= 0)
        policies::raise_domain_error<T>("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got a=%1%).", z, pol);
     if(z+delta <= 0)
        policies::raise_domain_error<T>("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Gamma function ratios only implemented for positive arguments (got b=%1%).", z+delta, pol);
  
     if(floor(delta) == delta)
     {
        if(floor(z) == z)
        {
           //
           // Both z and delta are integers, see if we can just use table lookup
           // of the factorials to get the result:
           //
           if((z <= max_factorial<T>::value) && (z + delta <= max_factorial<T>::value))
           {
              return unchecked_factorial<T>((unsigned)itrunc(z, pol) - 1) / unchecked_factorial<T>((unsigned)itrunc(T(z + delta), pol) - 1);
           }
        }
        if(fabs(delta) < 20)
        {
           //
           // delta is a small integer, we can use a finite product:
           //
           if(delta == 0)
              return 1;
           if(delta < 0)
           {
              z -= 1;
              T result = z;
              while(0 != (delta += 1))
              {
                 z -= 1;
                 result *= z;
              }
              return result;
           }
           else
           {
              T result = 1 / z;
              while(0 != (delta -= 1))
              {
                 z += 1;
                 result /= z;
              }
              return result;
           }
        }
     }
     typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
     return tgamma_delta_ratio_imp_lanczos(z, delta, pol, lanczos_type());
  }
  
  template <class T, class Policy>
  T gamma_p_derivative_imp(T a, T x, const Policy& pol)
  {
     //
     // Usual error checks first:
     //
     if(a <= 0)
        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);
     if(x < 0)
        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);
     //
     // Now special cases:
     //
     if(x == 0)
     {
        return (a > 1) ? 0 :
           (a == 1) ? 1 : policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
     }
     //
     // Normal case:
     //
     typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
     T f1 = detail::regularised_gamma_prefix(a, x, pol, lanczos_type());
     if((x < 1) && (tools::max_value<T>() * x < f1))
     {
        // overflow:
        return policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
     }
  
     f1 /= x;
  
     return f1;
  }
  
  template <class T, class Policy>
  inline typename tools::promote_args<T>::type 
     tgamma(T z, const Policy& /* pol */, const mpl::true_)
  {
     BOOST_FPU_EXCEPTION_GUARD
     typedef typename tools::promote_args<T>::type result_type;
     typedef typename policies::evaluation<result_type, Policy>::type value_type;
     typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
     typedef typename policies::normalise<
        Policy, 
        policies::promote_float<false>, 
        policies::promote_double<false>, 
        policies::discrete_quantile<>,
        policies::assert_undefined<> >::type forwarding_policy;
     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%)");
  }
  
  template <class T1, class T2, class Policy>
  inline typename tools::promote_args<T1, T2>::type
     tgamma(T1 a, T2 z, const Policy&, const mpl::false_)
  {
     BOOST_FPU_EXCEPTION_GUARD
     typedef typename tools::promote_args<T1, T2>::type result_type;
     typedef typename policies::evaluation<result_type, Policy>::type value_type;
     typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
     typedef typename policies::normalise<
        Policy, 
        policies::promote_float<false>, 
        policies::promote_double<false>, 
        policies::discrete_quantile<>,
        policies::assert_undefined<> >::type forwarding_policy;
     return policies::checked_narrowing_cast<result_type, forwarding_policy>(
        detail::gamma_incomplete_imp(static_cast<value_type>(a),
        static_cast<value_type>(z), false, true,
        forwarding_policy(), static_cast<value_type*>(0)), "boost::math::tgamma<%1%>(%1%, %1%)");
  }
  
  template <class T1, class T2>
  inline typename tools::promote_args<T1, T2>::type
     tgamma(T1 a, T2 z, const mpl::false_ tag)
  {
     return tgamma(a, z, policies::policy<>(), tag);
  }
  
  } // namespace detail
  
  template <class T>
  inline typename tools::promote_args<T>::type 
     tgamma(T z)
  {
     return tgamma(z, policies::policy<>());
  }
  
  template <class T, class Policy>
  inline typename tools::promote_args<T>::type 
     lgamma(T z, int* sign, const Policy&)
  {
     BOOST_FPU_EXCEPTION_GUARD
     typedef typename tools::promote_args<T>::type result_type;
     typedef typename policies::evaluation<result_type, Policy>::type value_type;
     typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
     typedef typename policies::normalise<
        Policy, 
        policies::promote_float<false>, 
        policies::promote_double<false>, 
        policies::discrete_quantile<>,
        policies::assert_undefined<> >::type forwarding_policy;
     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%)");
  }
  
  template <class T>
  inline typename tools::promote_args<T>::type 
     lgamma(T z, int* sign)
  {
     return lgamma(z, sign, policies::policy<>());
  }
  
  template <class T, class Policy>
  inline typename tools::promote_args<T>::type 
     lgamma(T x, const Policy& pol)
  {
     return ::boost::math::lgamma(x, 0, pol);
  }
  
  template <class T>
  inline typename tools::promote_args<T>::type 
     lgamma(T x)
  {
     return ::boost::math::lgamma(x, 0, policies::policy<>());
  }
  
  template <class T, class Policy>
  inline typename tools::promote_args<T>::type 
     tgamma1pm1(T z, const Policy& /* pol */)
  {
     BOOST_FPU_EXCEPTION_GUARD
     typedef typename tools::promote_args<T>::type result_type;
     typedef typename policies::evaluation<result_type, Policy>::type value_type;
     typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
     typedef typename policies::normalise<
        Policy, 
        policies::promote_float<false>, 
        policies::promote_double<false>, 
        policies::discrete_quantile<>,
        policies::assert_undefined<> >::type forwarding_policy;
  
     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%)");
  }
  
  template <class T>
  inline typename tools::promote_args<T>::type 
     tgamma1pm1(T z)
  {
     return tgamma1pm1(z, policies::policy<>());
  }
  
  //
  // Full upper incomplete gamma:
  //
  template <class T1, class T2>
  inline typename tools::promote_args<T1, T2>::type
     tgamma(T1 a, T2 z)
  {
     //
     // Type T2 could be a policy object, or a value, select the 
     // right overload based on T2:
     //
     typedef typename policies::is_policy<T2>::type maybe_policy;
     return detail::tgamma(a, z, maybe_policy());
  }
  template <class T1, class T2, class Policy>
  inline typename tools::promote_args<T1, T2>::type
     tgamma(T1 a, T2 z, const Policy& pol)
  {
     return detail::tgamma(a, z, pol, mpl::false_());
  }
  //
  // Full lower incomplete gamma:
  //
  template <class T1, class T2, class Policy>
  inline typename tools::promote_args<T1, T2>::type
     tgamma_lower(T1 a, T2 z, const Policy&)
  {
     BOOST_FPU_EXCEPTION_GUARD
     typedef typename tools::promote_args<T1, T2>::type result_type;
     typedef typename policies::evaluation<result_type, Policy>::type value_type;
     typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
     typedef typename policies::normalise<
        Policy, 
        policies::promote_float<false>, 
        policies::promote_double<false>, 
        policies::discrete_quantile<>,
        policies::assert_undefined<> >::type forwarding_policy;
  
     return policies::checked_narrowing_cast<result_type, forwarding_policy>(
        detail::gamma_incomplete_imp(static_cast<value_type>(a),
        static_cast<value_type>(z), false, false,
        forwarding_policy(), static_cast<value_type*>(0)), "tgamma_lower<%1%>(%1%, %1%)");
  }
  template <class T1, class T2>
  inline typename tools::promote_args<T1, T2>::type
     tgamma_lower(T1 a, T2 z)
  {
     return tgamma_lower(a, z, policies::policy<>());
  }
  //
  // Regularised upper incomplete gamma:
  //
  template <class T1, class T2, class Policy>
  inline typename tools::promote_args<T1, T2>::type
     gamma_q(T1 a, T2 z, const Policy& /* pol */)
  {
     BOOST_FPU_EXCEPTION_GUARD
     typedef typename tools::promote_args<T1, T2>::type result_type;
     typedef typename policies::evaluation<result_type, Policy>::type value_type;
     typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
     typedef typename policies::normalise<
        Policy, 
        policies::promote_float<false>, 
        policies::promote_double<false>, 
        policies::discrete_quantile<>,
        policies::assert_undefined<> >::type forwarding_policy;
  
     return policies::checked_narrowing_cast<result_type, forwarding_policy>(
        detail::gamma_incomplete_imp(static_cast<value_type>(a),
        static_cast<value_type>(z), true, true,
        forwarding_policy(), static_cast<value_type*>(0)), "gamma_q<%1%>(%1%, %1%)");
  }
  template <class T1, class T2>
  inline typename tools::promote_args<T1, T2>::type
     gamma_q(T1 a, T2 z)
  {
     return gamma_q(a, z, policies::policy<>());
  }
  //
  // Regularised lower incomplete gamma:
  //
  template <class T1, class T2, class Policy>
  inline typename tools::promote_args<T1, T2>::type
     gamma_p(T1 a, T2 z, const Policy&)
  {
     BOOST_FPU_EXCEPTION_GUARD
     typedef typename tools::promote_args<T1, T2>::type result_type;
     typedef typename policies::evaluation<result_type, Policy>::type value_type;
     typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
     typedef typename policies::normalise<
        Policy, 
        policies::promote_float<false>, 
        policies::promote_double<false>, 
        policies::discrete_quantile<>,
        policies::assert_undefined<> >::type forwarding_policy;
  
     return policies::checked_narrowing_cast<result_type, forwarding_policy>(
        detail::gamma_incomplete_imp(static_cast<value_type>(a),
        static_cast<value_type>(z), true, false,
        forwarding_policy(), static_cast<value_type*>(0)), "gamma_p<%1%>(%1%, %1%)");
  }
  template <class T1, class T2>
  inline typename tools::promote_args<T1, T2>::type
     gamma_p(T1 a, T2 z)
  {
     return gamma_p(a, z, policies::policy<>());
  }
  
  // ratios of gamma functions:
  template <class T1, class T2, class Policy>
  inline typename tools::promote_args<T1, T2>::type 
     tgamma_delta_ratio(T1 z, T2 delta, const Policy& /* pol */)
  {
     BOOST_FPU_EXCEPTION_GUARD
     typedef typename tools::promote_args<T1, T2>::type result_type;
     typedef typename policies::evaluation<result_type, Policy>::type value_type;
     typedef typename policies::normalise<
        Policy, 
        policies::promote_float<false>, 
        policies::promote_double<false>, 
        policies::discrete_quantile<>,
        policies::assert_undefined<> >::type forwarding_policy;
  
     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%)");
  }
  template <class T1, class T2>
  inline typename tools::promote_args<T1, T2>::type 
     tgamma_delta_ratio(T1 z, T2 delta)
  {
     return tgamma_delta_ratio(z, delta, policies::policy<>());
  }
  template <class T1, class T2, class Policy>
  inline typename tools::promote_args<T1, T2>::type 
     tgamma_ratio(T1 a, T2 b, const Policy&)
  {
     typedef typename tools::promote_args<T1, T2>::type result_type;
     typedef typename policies::evaluation<result_type, Policy>::type value_type;
     typedef typename policies::normalise<
        Policy, 
        policies::promote_float<false>, 
        policies::promote_double<false>, 
        policies::discrete_quantile<>,
        policies::assert_undefined<> >::type forwarding_policy;
  
     return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_delta_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(static_cast<value_type>(b) - static_cast<value_type>(a)), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
  }
  template <class T1, class T2>
  inline typename tools::promote_args<T1, T2>::type 
     tgamma_ratio(T1 a, T2 b)
  {
     return tgamma_ratio(a, b, policies::policy<>());
  }
  
  template <class T1, class T2, class Policy>
  inline typename tools::promote_args<T1, T2>::type 
     gamma_p_derivative(T1 a, T2 x, const Policy&)
  {
     BOOST_FPU_EXCEPTION_GUARD
     typedef typename tools::promote_args<T1, T2>::type result_type;
     typedef typename policies::evaluation<result_type, Policy>::type value_type;
     typedef typename policies::normalise<
        Policy, 
        policies::promote_float<false>, 
        policies::promote_double<false>, 
        policies::discrete_quantile<>,
        policies::assert_undefined<> >::type forwarding_policy;
  
     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%)");
  }
  template <class T1, class T2>
  inline typename tools::promote_args<T1, T2>::type 
     gamma_p_derivative(T1 a, T2 x)
  {
     return gamma_p_derivative(a, x, policies::policy<>());
  }
  
  } // namespace math
  } // namespace boost
  
  #ifdef BOOST_MSVC
  # pragma warning(pop)
  #endif
  
  #include <boost/math/special_functions/detail/igamma_inverse.hpp>
  #include <boost/math/special_functions/detail/gamma_inva.hpp>
  #include <boost/math/special_functions/erf.hpp>
  
  #endif // BOOST_MATH_SF_GAMMA_HPP
  
  
  
  
  