boost/math/special_functions/gamma.hpp

   1
   2 //  Copyright John Maddock 2006-7, 2013-14.
   3 //  Copyright Paul A. Bristow 2007, 2013-14.
   4 //  Copyright Nikhar Agrawal 2013-14
   5 //  Copyright Christopher Kormanyos 2013-14
   6
   7 //  Use, modification and distribution are subject to the
   8 //  Boost Software License, Version 1.0. (See accompanying file
   9 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  10
  11 #ifndef BOOST_MATH_SF_GAMMA_HPP
  12 #define BOOST_MATH_SF_GAMMA_HPP
  13
  14 #ifdef _MSC_VER
  15 #pragma once
  16 #endif
  17
  18 #include <boost/config.hpp>
  19 #include <boost/math/tools/series.hpp>
  20 #include <boost/math/tools/fraction.hpp>
  21 #include <boost/math/tools/precision.hpp>
  22 #include <boost/math/tools/promotion.hpp>
  23 #include <boost/math/policies/error_handling.hpp>
  24 #include <boost/math/constants/constants.hpp>
  25 #include <boost/math/special_functions/math_fwd.hpp>
  26 #include <boost/math/special_functions/log1p.hpp>
  27 #include <boost/math/special_functions/trunc.hpp>
  28 #include <boost/math/special_functions/powm1.hpp>
  29 #include <boost/math/special_functions/sqrt1pm1.hpp>
  30 #include <boost/math/special_functions/lanczos.hpp>
  31 #include <boost/math/special_functions/fpclassify.hpp>
  32 #include <boost/math/special_functions/detail/igamma_large.hpp>
  33 #include <boost/math/special_functions/detail/unchecked_factorial.hpp>
  34 #include <boost/math/special_functions/detail/lgamma_small.hpp>
  35 #include <boost/math/special_functions/bernoulli.hpp>
  36 #include <boost/math/special_functions/polygamma.hpp>
  37 #include <boost/type_traits/is_convertible.hpp>
  38 #include <boost/assert.hpp>
  39 #include <boost/mpl/greater.hpp>
  40 #include <boost/mpl/equal_to.hpp>
  41 #include <boost/mpl/greater.hpp>
  42
  43 #include <boost/config/no_tr1/cmath.hpp>
  44 #include <algorithm>
  45
  46 #ifdef BOOST_MSVC
  47 # pragma warning(push)
  48 # pragma warning(disable: 4702) // unreachable code (return after domain_error throw).
  49 # pragma warning(disable: 4127) // conditional expression is constant.
  50 # pragma warning(disable: 4100) // unreferenced formal parameter.
  51 // Several variables made comments,
  52 // but some difficulty as whether referenced on not may depend on macro values.
  53 // So to be safe, 4100 warnings suppressed.
  54 // TODO - revisit this?
  55 #endif
  56
  57 namespace boost{ namespace math{
  58
  59 namespace detail{
  60
  61 template <class T>
  62 inline bool is_odd(T v, const boost::true_type&)
  63 {
  64    int i = static_cast<int>(v);
  65    return i&1;
  66 }
  67 template <class T>
  68 inline bool is_odd(T v, const boost::false_type&)
  69 {
  70    // Oh dear can't cast T to int!
  71    BOOST_MATH_STD_USING
  72    T modulus = v - 2 * floor(v/2);
  73    return static_cast<bool>(modulus != 0);
  74 }
  75 template <class T>
  76 inline bool is_odd(T v)
  77 {
  78    return is_odd(v, ::boost::is_convertible<T, int>());
  79 }
  80
  81 template <class T>
  82 T sinpx(T z)
  83 {
  84    // Ad hoc function calculates x * sin(pi * x),
  85    // taking extra care near when x is near a whole number.
  86    BOOST_MATH_STD_USING
  87    int sign = 1;
  88    if(z < 0)
  89    {
  90       z = -z;
  91    }
  92    T fl = floor(z);
  93    T dist;
  94    if(is_odd(fl))
  95    {
  96       fl += 1;
  97       dist = fl - z;
  98       sign = -sign;
  99    }
 100    else
 101    {
 102       dist = z - fl;
 103    }
 104    BOOST_ASSERT(fl >= 0);
 105    if(dist > 0.5)
 106       dist = 1 - dist;
 107    T result = sin(dist*boost::math::constants::pi<T>());
 108    return sign*z*result;
 109 } // template <class T> T sinpx(T z)
 110 //
 111 // tgamma(z), with Lanczos support:
 112 //
 113 template <class T, class Policy, class Lanczos>
 114 T gamma_imp(T z, const Policy& pol, const Lanczos& l)
 115 {
 116    BOOST_MATH_STD_USING
 117
 118    T result = 1;
 119
 120 #ifdef BOOST_MATH_INSTRUMENT
 121    static bool b = false;
 122    if(!b)
 123    {
 124       std::cout << "tgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;
 125       b = true;
 126    }
 127 #endif
 128    static const char* function = "boost::math::tgamma<%1%>(%1%)";
 129
 130    if(z <= 0)
 131    {
 132       if(floor(z) == z)
 133          return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
 134       if(z <= -20)
 135       {
 136          result = gamma_imp(T(-z), pol, l) * sinpx(z);
 137          BOOST_MATH_INSTRUMENT_VARIABLE(result);
 138          if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
 139             return -boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
 140          result = -boost::math::constants::pi<T>() / result;
 141          if(result == 0)
 142             return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
 143          if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL)
 144             return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", result, pol);
 145          BOOST_MATH_INSTRUMENT_VARIABLE(result);
 146          return result;
 147       }
 148
 149       // shift z to > 1:
 150       while(z < 0)
 151       {
 152          result /= z;
 153          z += 1;
 154       }
 155    }
 156    BOOST_MATH_INSTRUMENT_VARIABLE(result);
 157    if((floor(z) == z) && (z < max_factorial<T>::value))
 158    {
 159       result *= unchecked_factorial<T>(itrunc(z, pol) - 1);
 160       BOOST_MATH_INSTRUMENT_VARIABLE(result);
 161    }
 162    else if (z < tools::root_epsilon<T>())
 163    {
 164       if (z < 1 / tools::max_value<T>())
 165          result = policies::raise_overflow_error<T>(function, 0, pol);
 166       result *= 1 / z - constants::euler<T>();
 167    }
 168    else
 169    {
 170       result *= Lanczos::lanczos_sum(z);
 171       T zgh = (z + static_cast<T>(Lanczos::g()) - boost::math::constants::half<T>());
 172       T lzgh = log(zgh);
 173       BOOST_MATH_INSTRUMENT_VARIABLE(result);
 174       BOOST_MATH_INSTRUMENT_VARIABLE(tools::log_max_value<T>());
 175       if(z * lzgh > tools::log_max_value<T>())
 176       {
 177          // we're going to overflow unless this is done with care:
 178          BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
 179          if(lzgh * z / 2 > tools::log_max_value<T>())
 180             return boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
 181          T hp = pow(zgh, (z / 2) - T(0.25));
 182          BOOST_MATH_INSTRUMENT_VARIABLE(hp);
 183          result *= hp / exp(zgh);
 184          BOOST_MATH_INSTRUMENT_VARIABLE(result);
 185          if(tools::max_value<T>() / hp < result)
 186             return boost::math::sign(result) * policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
 187          result *= hp;
 188          BOOST_MATH_INSTRUMENT_VARIABLE(result);
 189       }
 190       else
 191       {
 192          BOOST_MATH_INSTRUMENT_VARIABLE(zgh);
 193          BOOST_MATH_INSTRUMENT_VARIABLE(pow(zgh, z - boost::math::constants::half<T>()));
 194          BOOST_MATH_INSTRUMENT_VARIABLE(exp(zgh));
 195          result *= pow(zgh, z - boost::math::constants::half<T>()) / exp(zgh);
 196          BOOST_MATH_INSTRUMENT_VARIABLE(result);
 197       }
 198    }
 199    return result;
 200 }
 201 //
 202 // lgamma(z) with Lanczos support:
 203 //
 204 template <class T, class Policy, class Lanczos>
 205 T lgamma_imp(T z, const Policy& pol, const Lanczos& l, int* sign = 0)
 206 {
 207 #ifdef BOOST_MATH_INSTRUMENT
 208    static bool b = false;
 209    if(!b)
 210    {
 211       std::cout << "lgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;
 212       b = true;
 213    }
 214 #endif
 215
 216    BOOST_MATH_STD_USING
 217
 218    static const char* function = "boost::math::lgamma<%1%>(%1%)";
 219
 220    T result = 0;
 221    int sresult = 1;
 222    if(z <= -tools::root_epsilon<T>())
 223    {
 224       // reflection formula:
 225       if(floor(z) == z)
 226          return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%.", z, pol);
 227
 228       T t = sinpx(z);
 229       z = -z;
 230       if(t < 0)
 231       {
 232          t = -t;
 233       }
 234       else
 235       {
 236          sresult = -sresult;
 237       }
 238       result = log(boost::math::constants::pi<T>()) - lgamma_imp(z, pol, l) - log(t);
 239    }
 240    else if (z < tools::root_epsilon<T>())
 241    {
 242       if (0 == z)
 243          return policies::raise_pole_error<T>(function, "Evaluation of lgamma at %1%.", z, pol);
 244       if (fabs(z) < 1 / tools::max_value<T>())
 245          result = -log(fabs(z));
 246       else
 247          result = log(fabs(1 / z - constants::euler<T>()));
 248       if (z < 0)
 249          sresult = -1;
 250    }
 251    else if(z < 15)
 252    {
 253       typedef typename policies::precision<T, Policy>::type precision_type;
 254       typedef typename mpl::if_<
 255          mpl::and_<
 256             mpl::less_equal<precision_type, mpl::int_<64> >,
 257             mpl::greater<precision_type, mpl::int_<0> >
 258          >,
 259          mpl::int_<64>,
 260          typename mpl::if_<
 261             mpl::and_<
 262                mpl::less_equal<precision_type, mpl::int_<113> >,
 263                mpl::greater<precision_type, mpl::int_<0> >
 264             >,
 265             mpl::int_<113>, mpl::int_<0> >::type
 266           >::type tag_type;
 267       result = lgamma_small_imp<T>(z, T(z - 1), T(z - 2), tag_type(), pol, l);
 268    }
 269    else if((z >= 3) && (z < 100) && (std::numeric_limits<T>::max_exponent >= 1024))
 270    {
 271       // taking the log of tgamma reduces the error, no danger of overflow here:
 272       result = log(gamma_imp(z, pol, l));
 273    }
 274    else
 275    {
 276       // regular evaluation:
 277       T zgh = static_cast<T>(z + Lanczos::g() - boost::math::constants::half<T>());
 278       result = log(zgh) - 1;
 279       result *= z - 0.5f;
 280       //
 281       // Only add on the lanczos sum part if we're going to need it:
 282       //
 283       if(result * tools::epsilon<T>() < 20)
 284          result += log(Lanczos::lanczos_sum_expG_scaled(z));
 285    }
 286
 287    if(sign)
 288       *sign = sresult;
 289    return result;
 290 }
 291
 292 //
 293 // Incomplete gamma functions follow:
 294 //
 295 template <class T>
 296 struct upper_incomplete_gamma_fract
 297 {
 298 private:
 299    T z, a;
 300    int k;
 301 public:
 302    typedef std::pair<T,T> result_type;
 303
 304    upper_incomplete_gamma_fract(T a1, T z1)
 305       : z(z1-a1+1), a(a1), k(0)
 306    {
 307    }
 308
 309    result_type operator()()
 310    {
 311       ++k;
 312       z += 2;
 313       return result_type(k * (a - k), z);
 314    }
 315 };
 316
 317 template <class T>
 318 inline T upper_gamma_fraction(T a, T z, T eps)
 319 {
 320    // Multiply result by z^a * e^-z to get the full
 321    // upper incomplete integral.  Divide by tgamma(z)
 322    // to normalise.
 323    upper_incomplete_gamma_fract<T> f(a, z);
 324    return 1 / (z - a + 1 + boost::math::tools::continued_fraction_a(f, eps));
 325 }
 326
 327 template <class T>
 328 struct lower_incomplete_gamma_series
 329 {
 330 private:
 331    T a, z, result;
 332 public:
 333    typedef T result_type;
 334    lower_incomplete_gamma_series(T a1, T z1) : a(a1), z(z1), result(1){}
 335
 336    T operator()()
 337    {
 338       T r = result;
 339       a += 1;
 340       result *= z/a;
 341       return r;
 342    }
 343 };
 344
 345 template <class T, class Policy>
 346 inline T lower_gamma_series(T a, T z, const Policy& pol, T init_value = 0)
 347 {
 348    // Multiply result by ((z^a) * (e^-z) / a) to get the full
 349    // lower incomplete integral. Then divide by tgamma(a)
 350    // to get the normalised value.
 351    lower_incomplete_gamma_series<T> s(a, z);
 352    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
 353    T factor = policies::get_epsilon<T, Policy>();
 354    T result = boost::math::tools::sum_series(s, factor, max_iter, init_value);
 355    policies::check_series_iterations<T>("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol);
 356    return result;
 357 }
 358
 359 //
 360 // Fully generic tgamma and lgamma use Stirling's approximation
 361 // with Bernoulli numbers.
 362 //
 363 template<class T>
 364 std::size_t highest_bernoulli_index()
 365 {
 366    const float digits10_of_type = (std::numeric_limits<T>::is_specialized
 367                                       ? static_cast<float>(std::numeric_limits<T>::digits10)
 368                                       : static_cast<float>(boost::math::tools::digits<T>() * 0.301F));
 369
 370    // Find the high index n for Bn to produce the desired precision in Stirling's calculation.
 371    return static_cast<std::size_t>(18.0F + (0.6F * digits10_of_type));
 372 }
 373
 374 template<class T>
 375 int minimum_argument_for_bernoulli_recursion()
 376 {
 377    const float digits10_of_type = (std::numeric_limits<T>::is_specialized
 378                                       ? static_cast<float>(std::numeric_limits<T>::digits10)
 379                                       : static_cast<float>(boost::math::tools::digits<T>() * 0.301F));
 380
 381    const float limit = std::ceil(std::pow(1.0f / std::ldexp(1.0f, 1-boost::math::tools::digits<T>()), 1.0f / 20.0f));
 382
 383    return (int)((std::min)(digits10_of_type * 1.7F, limit));
 384 }
 385
 386 template <class T, class Policy>
 387 T scaled_tgamma_no_lanczos(const T& z, const Policy& pol, bool islog = false)
 388 {
 389    BOOST_MATH_STD_USING
 390    //
 391    // Calculates tgamma(z) / (z/e)^z
 392    // Requires that our argument is large enough for Sterling's approximation to hold.
 393    // Used internally when combining gamma's of similar magnitude without logarithms.
 394    //
 395    BOOST_ASSERT(minimum_argument_for_bernoulli_recursion<T>() <= z);
 396
 397    // Perform the Bernoulli series expansion of Stirling's approximation.
 398
 399    const std::size_t number_of_bernoullis_b2n = policies::get_max_series_iterations<Policy>();
 400
 401    T one_over_x_pow_two_n_minus_one = 1 / z;
 402    const T one_over_x2 = one_over_x_pow_two_n_minus_one * one_over_x_pow_two_n_minus_one;
 403    T sum = (boost::math::bernoulli_b2n<T>(1) / 2) * one_over_x_pow_two_n_minus_one;
 404    const T target_epsilon_to_break_loop = sum * boost::math::tools::epsilon<T>();
 405    const T half_ln_two_pi_over_z = sqrt(boost::math::constants::two_pi<T>() / z);
 406    T last_term = 2 * sum;
 407
 408    for (std::size_t n = 2U;; ++n)
 409    {
 410       one_over_x_pow_two_n_minus_one *= one_over_x2;
 411
 412       const std::size_t n2 = static_cast<std::size_t>(n * 2U);
 413
 414       const T term = (boost::math::bernoulli_b2n<T>(static_cast<int>(n)) * one_over_x_pow_two_n_minus_one) / (n2 * (n2 - 1U));
 415
 416       if ((n >= 3U) && (abs(term) < target_epsilon_to_break_loop))
 417       {
 418          // We have reached the desired precision in Stirling's expansion.
 419          // Adding additional terms to the sum of this divergent asymptotic
 420          // expansion will not improve the result.
 421
 422          // Break from the loop.
 423          break;
 424       }
 425       if (n > number_of_bernoullis_b2n)
 426          return policies::raise_evaluation_error("scaled_tgamma_no_lanczos<%1%>()", "Exceeded maximum series iterations without reaching convergence, best approximation was %1%", T(exp(sum) * half_ln_two_pi_over_z), pol);
 427
 428       sum += term;
 429
 430       // Sanity check for divergence:
 431       T fterm = fabs(term);
 432       if(fterm > last_term)
 433          return policies::raise_evaluation_error("scaled_tgamma_no_lanczos<%1%>()", "Series became divergent without reaching convergence, best approximation was %1%", T(exp(sum) * half_ln_two_pi_over_z), pol);
 434       last_term = fterm;
 435    }
 436
 437    // Complete Stirling's approximation.
 438    T scaled_gamma_value = islog ? T(sum + log(half_ln_two_pi_over_z)) : T(exp(sum) * half_ln_two_pi_over_z);
 439    return scaled_gamma_value;
 440 }
 441
 442 // Forward declaration of the lgamma_imp template specialization.
 443 template <class T, class Policy>
 444 T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign = 0);
 445
 446 template <class T, class Policy>
 447 T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&)
 448 {
 449    BOOST_MATH_STD_USING
 450
 451    static const char* function = "boost::math::tgamma<%1%>(%1%)";
 452
 453    // Check if the argument of tgamma is identically zero.
 454    const bool is_at_zero = (z == 0);
 455
 456    if((boost::math::isnan)(z) || (is_at_zero) || ((boost::math::isinf)(z) && (z < 0)))
 457       return policies::raise_domain_error<T>(function, "Evaluation of tgamma at %1%.", z, pol);
 458
 459    const bool b_neg = (z < 0);
 460
 461    const bool floor_of_z_is_equal_to_z = (floor(z) == z);
 462
 463    // Special case handling of small factorials:
 464    if((!b_neg) && floor_of_z_is_equal_to_z && (z < boost::math::max_factorial<T>::value))
 465    {
 466       return boost::math::unchecked_factorial<T>(itrunc(z) - 1);
 467    }
 468
 469    // Make a local, unsigned copy of the input argument.
 470    T zz((!b_neg) ? z : -z);
 471
 472    // Special case for ultra-small z:
 473    if(zz < tools::cbrt_epsilon<T>())
 474    {
 475       const T a0(1);
 476       const T a1(boost::math::constants::euler<T>());
 477       const T six_euler_squared((boost::math::constants::euler<T>() * boost::math::constants::euler<T>()) * 6);
 478       const T a2((six_euler_squared -  boost::math::constants::pi_sqr<T>()) / 12);
 479
 480       const T inverse_tgamma_series = z * ((a2 * z + a1) * z + a0);
 481
 482       return 1 / inverse_tgamma_series;
 483    }
 484
 485    // Scale the argument up for the calculation of lgamma,
 486    // and use downward recursion later for the final result.
 487    const int min_arg_for_recursion = minimum_argument_for_bernoulli_recursion<T>();
 488
 489    int n_recur;
 490
 491    if(zz < min_arg_for_recursion)
 492    {
 493       n_recur = boost::math::itrunc(min_arg_for_recursion - zz) + 1;
 494
 495       zz += n_recur;
 496    }
 497    else
 498    {
 499       n_recur = 0;
 500    }
 501    if (!n_recur)
 502    {
 503       if (zz > tools::log_max_value<T>())
 504          return policies::raise_overflow_error<T>(function, 0, pol);
 505       if (log(zz) * zz / 2 > tools::log_max_value<T>())
 506          return policies::raise_overflow_error<T>(function, 0, pol);
 507    }
 508    T gamma_value = scaled_tgamma_no_lanczos(zz, pol);
 509    T power_term = pow(zz, zz / 2);
 510    T exp_term = exp(-zz);
 511    gamma_value *= (power_term * exp_term);
 512    if(!n_recur && (tools::max_value<T>() / power_term < gamma_value))
 513       return policies::raise_overflow_error<T>(function, 0, pol);
 514    gamma_value *= power_term;
 515
 516    // Rescale the result using downward recursion if necessary.
 517    if(n_recur)
 518    {
 519       // The order of divides is important, if we keep subtracting 1 from zz
 520       // we DO NOT get back to z (cancellation error).  Further if z < epsilon
 521       // we would end up dividing by zero.  Also in order to prevent spurious
 522       // overflow with the first division, we must save dividing by |z| till last,
 523       // so the optimal order of divides is z+1, z+2, z+3...z+n_recur-1,z.
 524       zz = fabs(z) + 1;
 525       for(int k = 1; k < n_recur; ++k)
 526       {
 527          gamma_value /= zz;
 528          zz += 1;
 529       }
 530       gamma_value /= fabs(z);
 531    }
 532
 533    // Return the result, accounting for possible negative arguments.
 534    if(b_neg)
 535    {
 536       // Provide special error analysis for:
 537       // * arguments in the neighborhood of a negative integer
 538       // * arguments exactly equal to a negative integer.
 539
 540       // Check if the argument of tgamma is exactly equal to a negative integer.
 541       if(floor_of_z_is_equal_to_z)
 542          return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
 543
 544       gamma_value *= sinpx(z);
 545
 546       BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value);
 547
 548       const bool result_is_too_large_to_represent = (   (abs(gamma_value) < 1)
 549                                                      && ((tools::max_value<T>() * abs(gamma_value)) < boost::math::constants::pi<T>()));
 550
 551       if(result_is_too_large_to_represent)
 552          return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
 553
 554       gamma_value = -boost::math::constants::pi<T>() / gamma_value;
 555       BOOST_MATH_INSTRUMENT_VARIABLE(gamma_value);
 556
 557       if(gamma_value == 0)
 558          return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
 559
 560       if((boost::math::fpclassify)(gamma_value) == static_cast<int>(FP_SUBNORMAL))
 561          return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", gamma_value, pol);
 562    }
 563
 564    return gamma_value;
 565 }
 566
 567 template <class T, class Policy>
 568 inline T log_gamma_near_1(const T& z, Policy const& pol)
 569 {
 570    //
 571    // This is for the multiprecision case where there is
 572    // no lanczos support, use a taylor series at z = 1,
 573    // see https://www.wolframalpha.com/input/?i=taylor+series+lgamma(x)+at+x+%3D+1
 574    //
 575    BOOST_MATH_STD_USING // ADL of std names
 576
 577    BOOST_ASSERT(fabs(z) < 1);
 578
 579    T result = -constants::euler<T>() * z;
 580
 581    T power_term = z * z / 2;
 582    int n = 2;
 583    T term = 0;
 584
 585    do
 586    {
 587       term = power_term * boost::math::polygamma(n - 1, T(1));
 588       result += term;
 589       ++n;
 590       power_term *= z / n;
 591    } while (fabs(result) * tools::epsilon<T>() < fabs(term));
 592
 593    return result;
 594 }
 595
 596 template <class T, class Policy>
 597 T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos&, int* sign)
 598 {
 599    BOOST_MATH_STD_USING
 600
 601    static const char* function = "boost::math::lgamma<%1%>(%1%)";
 602
 603    // Check if the argument of lgamma is identically zero.
 604    const bool is_at_zero = (z == 0);
 605
 606    if(is_at_zero)
 607       return policies::raise_domain_error<T>(function, "Evaluation of lgamma at zero %1%.", z, pol);
 608    if((boost::math::isnan)(z))
 609       return policies::raise_domain_error<T>(function, "Evaluation of lgamma at %1%.", z, pol);
 610    if((boost::math::isinf)(z))
 611       return policies::raise_overflow_error<T>(function, 0, pol);
 612
 613    const bool b_neg = (z < 0);
 614
 615    const bool floor_of_z_is_equal_to_z = (floor(z) == z);
 616
 617    // Special case handling of small factorials:
 618    if((!b_neg) && floor_of_z_is_equal_to_z && (z < boost::math::max_factorial<T>::value))
 619    {
 620       if (sign)
 621          *sign = 1;
 622       return log(boost::math::unchecked_factorial<T>(itrunc(z) - 1));
 623    }
 624
 625    // Make a local, unsigned copy of the input argument.
 626    T zz((!b_neg) ? z : -z);
 627
 628    const int min_arg_for_recursion = minimum_argument_for_bernoulli_recursion<T>();
 629
 630    T log_gamma_value;
 631
 632    if (zz < min_arg_for_recursion)
 633    {
 634       // Here we simply take the logarithm of tgamma(). This is somewhat
 635       // inefficient, but simple. The rationale is that the argument here
 636       // is relatively small and overflow is not expected to be likely.
 637       if (sign)
 638          * sign = 1;
 639       if(fabs(z - 1) < 0.25)
 640       {
 641          log_gamma_value = log_gamma_near_1(T(zz - 1), pol);
 642       }
 643       else if(fabs(z - 2) < 0.25)
 644       {
 645          log_gamma_value = log_gamma_near_1(T(zz - 2), pol) + log(zz - 1);
 646       }
 647       else if (z > -tools::root_epsilon<T>())
 648       {
 649          // Reflection formula may fail if z is very close to zero, let the series
 650          // expansion for tgamma close to zero do the work:
 651          if (sign)
 652             *sign = z < 0 ? -1 : 1;
 653          return log(abs(gamma_imp(z, pol, lanczos::undefined_lanczos())));
 654       }
 655       else
 656       {
 657          // No issue with spurious overflow in reflection formula,
 658          // just fall through to regular code:
 659          T g = gamma_imp(zz, pol, lanczos::undefined_lanczos());
 660          if (sign)
 661          {
 662             *sign = g < 0 ? -1 : 1;
 663          }
 664          log_gamma_value = log(abs(g));
 665       }
 666    }
 667    else
 668    {
 669       // Perform the Bernoulli series expansion of Stirling's approximation.
 670       T sum = scaled_tgamma_no_lanczos(zz, pol, true);
 671       log_gamma_value = zz * (log(zz) - 1) + sum;
 672    }
 673
 674    int sign_of_result = 1;
 675
 676    if(b_neg)
 677    {
 678       // Provide special error analysis if the argument is exactly
 679       // equal to a negative integer.
 680
 681       // Check if the argument of lgamma is exactly equal to a negative integer.
 682       if(floor_of_z_is_equal_to_z)
 683          return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%.", z, pol);
 684
 685       T t = sinpx(z);
 686
 687       if(t < 0)
 688       {
 689          t = -t;
 690       }
 691       else
 692       {
 693          sign_of_result = -sign_of_result;
 694       }
 695
 696       log_gamma_value = - log_gamma_value
 697                         + log(boost::math::constants::pi<T>())
 698                         - log(t);
 699    }
 700
 701    if(sign != static_cast<int*>(0U)) { *sign = sign_of_result; }
 702
 703    return log_gamma_value;
 704 }
 705
 706 //
 707 // This helper calculates tgamma(dz+1)-1 without cancellation errors,
 708 // used by the upper incomplete gamma with z < 1:
 709 //
 710 template <class T, class Policy, class Lanczos>
 711 T tgammap1m1_imp(T dz, Policy const& pol, const Lanczos& l)
 712 {
 713    BOOST_MATH_STD_USING
 714
 715    typedef typename policies::precision<T,Policy>::type precision_type;
 716
 717    typedef typename mpl::if_<
 718       mpl::or_<
 719          mpl::less_equal<precision_type, mpl::int_<0> >,
 720          mpl::greater<precision_type, mpl::int_<113> >
 721       >,
 722       typename mpl::if_<
 723          mpl::and_<is_same<Lanczos, lanczos::lanczos24m113>, mpl::greater<precision_type, mpl::int_<0> > >,
 724          mpl::int_<113>,
 725          mpl::int_<0>
 726       >::type,
 727       typename mpl::if_<
 728          mpl::less_equal<precision_type, mpl::int_<64> >,
 729          mpl::int_<64>, mpl::int_<113> >::type
 730        >::type tag_type;
 731
 732    T result;
 733    if(dz < 0)
 734    {
 735       if(dz < -0.5)
 736       {
 737          // Best method is simply to subtract 1 from tgamma:
 738          result = boost::math::tgamma(1+dz, pol) - 1;
 739          BOOST_MATH_INSTRUMENT_CODE(result);
 740       }
 741       else
 742       {
 743          // Use expm1 on lgamma:
 744          result = boost::math::expm1(-boost::math::log1p(dz, pol)
 745             + lgamma_small_imp<T>(dz+2, dz + 1, dz, tag_type(), pol, l));
 746          BOOST_MATH_INSTRUMENT_CODE(result);
 747       }
 748    }
 749    else
 750    {
 751       if(dz < 2)
 752       {
 753          // Use expm1 on lgamma:
 754          result = boost::math::expm1(lgamma_small_imp<T>(dz+1, dz, dz-1, tag_type(), pol, l), pol);
 755          BOOST_MATH_INSTRUMENT_CODE(result);
 756       }
 757       else
 758       {
 759          // Best method is simply to subtract 1 from tgamma:
 760          result = boost::math::tgamma(1+dz, pol) - 1;
 761          BOOST_MATH_INSTRUMENT_CODE(result);
 762       }
 763    }
 764
 765    return result;
 766 }
 767
 768 template <class T, class Policy>
 769 inline T tgammap1m1_imp(T z, Policy const& pol,
 770                  const ::boost::math::lanczos::undefined_lanczos&)
 771 {
 772    BOOST_MATH_STD_USING // ADL of std names
 773
 774    if(fabs(z) < 0.55)
 775    {
 776       return boost::math::expm1(log_gamma_near_1(z, pol));
 777    }
 778    return boost::math::expm1(boost::math::lgamma(1 + z, pol));
 779 }
 780
 781 //
 782 // Series representation for upper fraction when z is small:
 783 //
 784 template <class T>
 785 struct small_gamma2_series
 786 {
 787    typedef T result_type;
 788
 789    small_gamma2_series(T a_, T x_) : result(-x_), x(-x_), apn(a_+1), n(1){}
 790
 791    T operator()()
 792    {
 793       T r = result / (apn);
 794       result *= x;
 795       result /= ++n;
 796       apn += 1;
 797       return r;
 798    }
 799
 800 private:
 801    T result, x, apn;
 802    int n;
 803 };
 804 //
 805 // calculate power term prefix (z^a)(e^-z) used in the non-normalised
 806 // incomplete gammas:
 807 //
 808 template <class T, class Policy>
 809 T full_igamma_prefix(T a, T z, const Policy& pol)
 810 {
 811    BOOST_MATH_STD_USING
 812
 813    T prefix;
 814    if (z > tools::max_value<T>())
 815       return 0;
 816    T alz = a * log(z);
 817
 818    if(z >= 1)
 819    {
 820       if((alz < tools::log_max_value<T>()) && (-z > tools::log_min_value<T>()))
 821       {
 822          prefix = pow(z, a) * exp(-z);
 823       }
 824       else if(a >= 1)
 825       {
 826          prefix = pow(z / exp(z/a), a);
 827       }
 828       else
 829       {
 830          prefix = exp(alz - z);
 831       }
 832    }
 833    else
 834    {
 835       if(alz > tools::log_min_value<T>())
 836       {
 837          prefix = pow(z, a) * exp(-z);
 838       }
 839       else if(z/a < tools::log_max_value<T>())
 840       {
 841          prefix = pow(z / exp(z/a), a);
 842       }
 843       else
 844       {
 845          prefix = exp(alz - z);
 846       }
 847    }
 848    //
 849    // This error handling isn't very good: it happens after the fact
 850    // rather than before it...
 851    //
 852    if((boost::math::fpclassify)(prefix) == (int)FP_INFINITE)
 853       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);
 854
 855    return prefix;
 856 }
 857 //
 858 // Compute (z^a)(e^-z)/tgamma(a)
 859 // most if the error occurs in this function:
 860 //
 861 template <class T, class Policy, class Lanczos>
 862 T regularised_gamma_prefix(T a, T z, const Policy& pol, const Lanczos& l)
 863 {
 864    BOOST_MATH_STD_USING
 865    if (z >= tools::max_value<T>())
 866       return 0;
 867    T agh = a + static_cast<T>(Lanczos::g()) - T(0.5);
 868    T prefix;
 869    T d = ((z - a) - static_cast<T>(Lanczos::g()) + T(0.5)) / agh;
 870
 871    if(a < 1)
 872    {
 873       //
 874       // We have to treat a < 1 as a special case because our Lanczos
 875       // approximations are optimised against the factorials with a > 1,
 876       // and for high precision types especially (128-bit reals for example)
 877       // very small values of a can give rather eroneous results for gamma
 878       // unless we do this:
 879       //
 880       // TODO: is this still required?  Lanczos approx should be better now?
 881       //
 882       if(z <= tools::log_min_value<T>())
 883       {
 884          // Oh dear, have to use logs, should be free of cancellation errors though:
 885          return exp(a * log(z) - z - lgamma_imp(a, pol, l));
 886       }
 887       else
 888       {
 889          // direct calculation, no danger of overflow as gamma(a) < 1/a
 890          // for small a.
 891          return pow(z, a) * exp(-z) / gamma_imp(a, pol, l);
 892       }
 893    }
 894    else if((fabs(d*d*a) <= 100) && (a > 150))
 895    {
 896       // special case for large a and a ~ z.
 897       prefix = a * boost::math::log1pmx(d, pol) + z * static_cast<T>(0.5 - Lanczos::g()) / agh;
 898       prefix = exp(prefix);
 899    }
 900    else
 901    {
 902       //
 903       // general case.
 904       // direct computation is most accurate, but use various fallbacks
 905       // for different parts of the problem domain:
 906       //
 907       T alz = a * log(z / agh);
 908       T amz = a - z;
 909       if(((std::min)(alz, amz) <= tools::log_min_value<T>()) || ((std::max)(alz, amz) >= tools::log_max_value<T>()))
 910       {
 911          T amza = amz / a;
 912          if(((std::min)(alz, amz)/2 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/2 < tools::log_max_value<T>()))
 913          {
 914             // compute square root of the result and then square it:
 915             T sq = pow(z / agh, a / 2) * exp(amz / 2);
 916             prefix = sq * sq;
 917          }
 918          else if(((std::min)(alz, amz)/4 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/4 < tools::log_max_value<T>()) && (z > a))
 919          {
 920             // compute the 4th root of the result then square it twice:
 921             T sq = pow(z / agh, a / 4) * exp(amz / 4);
 922             prefix = sq * sq;
 923             prefix *= prefix;
 924          }
 925          else if((amza > tools::log_min_value<T>()) && (amza < tools::log_max_value<T>()))
 926          {
 927             prefix = pow((z * exp(amza)) / agh, a);
 928          }
 929          else
 930          {
 931             prefix = exp(alz + amz);
 932          }
 933       }
 934       else
 935       {
 936          prefix = pow(z / agh, a) * exp(amz);
 937       }
 938    }
 939    prefix *= sqrt(agh / boost::math::constants::e<T>()) / Lanczos::lanczos_sum_expG_scaled(a);
 940    return prefix;
 941 }
 942 //
 943 // And again, without Lanczos support:
 944 //
 945 template <class T, class Policy>
 946 T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined_lanczos& l)
 947 {
 948    BOOST_MATH_STD_USING
 949
 950    if((a < 1) && (z < 1))
 951    {
 952       // No overflow possible since the power terms tend to unity as a,z -> 0
 953       return pow(z, a) * exp(-z) / boost::math::tgamma(a, pol);
 954    }
 955    else if(a > minimum_argument_for_bernoulli_recursion<T>())
 956    {
 957       T scaled_gamma = scaled_tgamma_no_lanczos(a, pol);
 958       T power_term = pow(z / a, a / 2);
 959       T a_minus_z = a - z;
 960       if ((0 == power_term) || (fabs(a_minus_z) > tools::log_max_value<T>()))
 961       {
 962          // The result is probably zero, but we need to be sure:
 963          return exp(a * log(z / a) + a_minus_z - log(scaled_gamma));
 964       }
 965       return (power_term * exp(a_minus_z)) * (power_term / scaled_gamma);
 966    }
 967    else
 968    {
 969       //
 970       // Usual case is to calculate the prefix at a+shift and recurse down
 971       // to the value we want:
 972       //
 973       const int min_z = minimum_argument_for_bernoulli_recursion<T>();
 974       long shift = 1 + ltrunc(min_z - a);
 975       T result = regularised_gamma_prefix(T(a + shift), z, pol, l);
 976       if (result != 0)
 977       {
 978          for (long i = 0; i < shift; ++i)
 979          {
 980             result /= z;
 981             result *= a + i;
 982          }
 983          return result;
 984       }
 985       else
 986       {
 987          //
 988          // We failed, most probably we have z << 1, try again, this time
 989          // we calculate z^a e^-z / tgamma(a+shift), combining power terms
 990          // as we go.  And again recurse down to the result.
 991          //
 992          T scaled_gamma = scaled_tgamma_no_lanczos(T(a + shift), pol);
 993          T power_term_1 = pow(z / (a + shift), a);
 994          T power_term_2 = pow(a + shift, -shift);
 995          T power_term_3 = exp(a + shift - z);
 996          if ((0 == power_term_1) || (0 == power_term_2) || (0 == power_term_3) || (fabs(a + shift - z) > tools::log_max_value<T>()))
 997          {
 998             // We have no test case that gets here, most likely the type T
 999             // has a high precision but low exponent range:
1000             return exp(a * log(z) - z - boost::math::lgamma(a, pol));
1001          }
1002          result = power_term_1 * power_term_2 * power_term_3 / scaled_gamma;
1003          for (long i = 0; i < shift; ++i)
1004          {
1005             result *= a + i;
1006          }
1007          return result;
1008       }
1009    }
1010 }
1011 //
1012 // Upper gamma fraction for very small a:
1013 //
1014 template <class T, class Policy>
1015 inline T tgamma_small_upper_part(T a, T x, const Policy& pol, T* pgam = 0, bool invert = false, T* pderivative = 0)
1016 {
1017    BOOST_MATH_STD_USING  // ADL of std functions.
1018    //
1019    // Compute the full upper fraction (Q) when a is very small:
1020    //
1021    T result;
1022    result = boost::math::tgamma1pm1(a, pol);
1023    if(pgam)
1024       *pgam = (result + 1) / a;
1025    T p = boost::math::powm1(x, a, pol);
1026    result -= p;
1027    result /= a;
1028    detail::small_gamma2_series<T> s(a, x);
1029    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>() - 10;
1030    p += 1;
1031    if(pderivative)
1032       *pderivative = p / (*pgam * exp(x));
1033    T init_value = invert ? *pgam : 0;
1034    result = -p * tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, (init_value - result) / p);
1035    policies::check_series_iterations<T>("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)", max_iter, pol);
1036    if(invert)
1037       result = -result;
1038    return result;
1039 }
1040 //
1041 // Upper gamma fraction for integer a:
1042 //
1043 template <class T, class Policy>
1044 inline T finite_gamma_q(T a, T x, Policy const& pol, T* pderivative = 0)
1045 {
1046    //
1047    // Calculates normalised Q when a is an integer:
1048    //
1049    BOOST_MATH_STD_USING
1050    T e = exp(-x);
1051    T sum = e;
1052    if(sum != 0)
1053    {
1054       T term = sum;
1055       for(unsigned n = 1; n < a; ++n)
1056       {
1057          term /= n;
1058          term *= x;
1059          sum += term;
1060       }
1061    }
1062    if(pderivative)
1063    {
1064       *pderivative = e * pow(x, a) / boost::math::unchecked_factorial<T>(itrunc(T(a - 1), pol));
1065    }
1066    return sum;
1067 }
1068 //
1069 // Upper gamma fraction for half integer a:
1070 //
1071 template <class T, class Policy>
1072 T finite_half_gamma_q(T a, T x, T* p_derivative, const Policy& pol)
1073 {
1074    //
1075    // Calculates normalised Q when a is a half-integer:
1076    //
1077    BOOST_MATH_STD_USING
1078    T e = boost::math::erfc(sqrt(x), pol);
1079    if((e != 0) && (a > 1))
1080    {
1081       T term = exp(-x) / sqrt(constants::pi<T>() * x);
1082       term *= x;
1083       static const T half = T(1) / 2;
1084       term /= half;
1085       T sum = term;
1086       for(unsigned n = 2; n < a; ++n)
1087       {
1088          term /= n - half;
1089          term *= x;
1090          sum += term;
1091       }
1092       e += sum;
1093       if(p_derivative)
1094       {
1095          *p_derivative = 0;
1096       }
1097    }
1098    else if(p_derivative)
1099    {
1100       // We'll be dividing by x later, so calculate derivative * x:
1101       *p_derivative = sqrt(x) * exp(-x) / constants::root_pi<T>();
1102    }
1103    return e;
1104 }
1105 //
1106 // Asymptotic approximation for large argument, see: https://dlmf.nist.gov/8.11#E2
1107 //
1108 template <class T>
1109 struct incomplete_tgamma_large_x_series
1110 {
1111    typedef T result_type;
1112    incomplete_tgamma_large_x_series(const T& a, const T& x)
1113       : a_poch(a - 1), z(x), term(1) {}
1114    T operator()()
1115    {
1116       T result = term;
1117       term *= a_poch / z;
1118       a_poch -= 1;
1119       return result;
1120    }
1121    T a_poch, z, term;
1122 };
1123
1124 template <class T, class Policy>
1125 T incomplete_tgamma_large_x(const T& a, const T& x, const Policy& pol)
1126 {
1127    BOOST_MATH_STD_USING
1128    incomplete_tgamma_large_x_series<T> s(a, x);
1129    boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
1130    T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
1131    boost::math::policies::check_series_iterations<T>("boost::math::tgamma<%1%>(%1%,%1%)", max_iter, pol);
1132    return result;
1133 }
1134
1135
1136 //
1137 // Main incomplete gamma entry point, handles all four incomplete gamma's:
1138 //
1139 template <class T, class Policy>
1140 T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
1141                        const Policy& pol, T* p_derivative)
1142 {
1143    static const char* function = "boost::math::gamma_p<%1%>(%1%, %1%)";
1144    if(a <= 0)
1145       return policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
1146    if(x < 0)
1147       return policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
1148
1149    BOOST_MATH_STD_USING
1150
1151    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1152
1153    T result = 0; // Just to avoid warning C4701: potentially uninitialized local variable 'result' used
1154
1155    if(a >= max_factorial<T>::value && !normalised)
1156    {
1157       //
1158       // When we're computing the non-normalized incomplete gamma
1159       // and a is large the result is rather hard to compute unless
1160       // we use logs.  There are really two options - if x is a long
1161       // way from a in value then we can reliably use methods 2 and 4
1162       // below in logarithmic form and go straight to the result.
1163       // Otherwise we let the regularized gamma take the strain
1164       // (the result is unlikely to unerflow in the central region anyway)
1165       // and combine with lgamma in the hopes that we get a finite result.
1166       //
1167       if(invert && (a * 4 < x))
1168       {
1169          // This is method 4 below, done in logs:
1170          result = a * log(x) - x;
1171          if(p_derivative)
1172             *p_derivative = exp(result);
1173          result += log(upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>()));
1174       }
1175       else if(!invert && (a > 4 * x))
1176       {
1177          // This is method 2 below, done in logs:
1178          result = a * log(x) - x;
1179          if(p_derivative)
1180             *p_derivative = exp(result);
1181          T init_value = 0;
1182          result += log(detail::lower_gamma_series(a, x, pol, init_value) / a);
1183       }
1184       else
1185       {
1186          result = gamma_incomplete_imp(a, x, true, invert, pol, p_derivative);
1187          if(result == 0)
1188          {
1189             if(invert)
1190             {
1191                // Try http://functions.wolfram.com/06.06.06.0039.01
1192                result = 1 + 1 / (12 * a) + 1 / (288 * a * a);
1193                result = log(result) - a + (a - 0.5f) * log(a) + log(boost::math::constants::root_two_pi<T>());
1194                if(p_derivative)
1195                   *p_derivative = exp(a * log(x) - x);
1196             }
1197             else
1198             {
1199                // This is method 2 below, done in logs, we're really outside the
1200                // range of this method, but since the result is almost certainly
1201                // infinite, we should probably be OK:
1202                result = a * log(x) - x;
1203                if(p_derivative)
1204                   *p_derivative = exp(result);
1205                T init_value = 0;
1206                result += log(detail::lower_gamma_series(a, x, pol, init_value) / a);
1207             }
1208          }
1209          else
1210          {
1211             result = log(result) + boost::math::lgamma(a, pol);
1212          }
1213       }
1214       if(result > tools::log_max_value<T>())
1215          return policies::raise_overflow_error<T>(function, 0, pol);
1216       return exp(result);
1217    }
1218
1219    BOOST_ASSERT((p_derivative == 0) || (normalised == true));
1220
1221    bool is_int, is_half_int;
1222    bool is_small_a = (a < 30) && (a <= x + 1) && (x < tools::log_max_value<T>());
1223    if(is_small_a)
1224    {
1225       T fa = floor(a);
1226       is_int = (fa == a);
1227       is_half_int = is_int ? false : (fabs(fa - a) == 0.5f);
1228    }
1229    else
1230    {
1231       is_int = is_half_int = false;
1232    }
1233
1234    int eval_method;
1235
1236    if(is_int && (x > 0.6))
1237    {
1238       // calculate Q via finite sum:
1239       invert = !invert;
1240       eval_method = 0;
1241    }
1242    else if(is_half_int && (x > 0.2))
1243    {
1244       // calculate Q via finite sum for half integer a:
1245       invert = !invert;
1246       eval_method = 1;
1247    }
1248    else if((x < tools::root_epsilon<T>()) && (a > 1))
1249    {
1250       eval_method = 6;
1251    }
1252    else if ((x > 1000) && ((a < x) || (fabs(a - 50) / x < 1)))
1253    {
1254       // calculate Q via asymptotic approximation:
1255       invert = !invert;
1256       eval_method = 7;
1257    }
1258    else if(x < 0.5)
1259    {
1260       //
1261       // Changeover criterion chosen to give a changeover at Q ~ 0.33
1262       //
1263       if(-0.4 / log(x) < a)
1264       {
1265          eval_method = 2;
1266       }
1267       else
1268       {
1269          eval_method = 3;
1270       }
1271    }
1272    else if(x < 1.1)
1273    {
1274       //
1275       // Changover here occurs when P ~ 0.75 or Q ~ 0.25:
1276       //
1277       if(x * 0.75f < a)
1278       {
1279          eval_method = 2;
1280       }
1281       else
1282       {
1283          eval_method = 3;
1284       }
1285    }
1286    else
1287    {
1288       //
1289       // Begin by testing whether we're in the "bad" zone
1290       // where the result will be near 0.5 and the usual
1291       // series and continued fractions are slow to converge:
1292       //
1293       bool use_temme = false;
1294       if(normalised && std::numeric_limits<T>::is_specialized && (a > 20))
1295       {
1296          T sigma = fabs((x-a)/a);
1297          if((a > 200) && (policies::digits<T, Policy>() <= 113))
1298          {
1299             //
1300             // This limit is chosen so that we use Temme's expansion
1301             // only if the result would be larger than about 10^-6.
1302             // Below that the regular series and continued fractions
1303             // converge OK, and if we use Temme's method we get increasing
1304             // errors from the dominant erfc term as it's (inexact) argument
1305             // increases in magnitude.
1306             //
1307             if(20 / a > sigma * sigma)
1308                use_temme = true;
1309          }
1310          else if(policies::digits<T, Policy>() <= 64)
1311          {
1312             // Note in this zone we can't use Temme's expansion for
1313             // types longer than an 80-bit real:
1314             // it would require too many terms in the polynomials.
1315             if(sigma < 0.4)
1316                use_temme = true;
1317          }
1318       }
1319       if(use_temme)
1320       {
1321          eval_method = 5;
1322       }
1323       else
1324       {
1325          //
1326          // Regular case where the result will not be too close to 0.5.
1327          //
1328          // Changeover here occurs at P ~ Q ~ 0.5
1329          // Note that series computation of P is about x2 faster than continued fraction
1330          // calculation of Q, so try and use the CF only when really necessary, especially
1331          // for small x.
1332          //
1333          if(x - (1 / (3 * x)) < a)
1334          {
1335             eval_method = 2;
1336          }
1337          else
1338          {
1339             eval_method = 4;
1340             invert = !invert;
1341          }
1342       }
1343    }
1344
1345    switch(eval_method)
1346    {
1347    case 0:
1348       {
1349          result = finite_gamma_q(a, x, pol, p_derivative);
1350          if(normalised == false)
1351             result *= boost::math::tgamma(a, pol);
1352          break;
1353       }
1354    case 1:
1355       {
1356          result = finite_half_gamma_q(a, x, p_derivative, pol);
1357          if(normalised == false)
1358             result *= boost::math::tgamma(a, pol);
1359          if(p_derivative && (*p_derivative == 0))
1360             *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
1361          break;
1362       }
1363    case 2:
1364       {
1365          // Compute P:
1366          result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
1367          if(p_derivative)
1368             *p_derivative = result;
1369          if(result != 0)
1370          {
1371             //
1372             // If we're going to be inverting the result then we can
1373             // reduce the number of series evaluations by quite
1374             // a few iterations if we set an initial value for the
1375             // series sum based on what we'll end up subtracting it from
1376             // at the end.
1377             // Have to be careful though that this optimization doesn't
1378             // lead to spurious numberic overflow.  Note that the
1379             // scary/expensive overflow checks below are more often
1380             // than not bypassed in practice for "sensible" input
1381             // values:
1382             //
1383             T init_value = 0;
1384             bool optimised_invert = false;
1385             if(invert)
1386             {
1387                init_value = (normalised ? 1 : boost::math::tgamma(a, pol));
1388                if(normalised || (result >= 1) || (tools::max_value<T>() * result > init_value))
1389                {
1390                   init_value /= result;
1391                   if(normalised || (a < 1) || (tools::max_value<T>() / a > init_value))
1392                   {
1393                      init_value *= -a;
1394                      optimised_invert = true;
1395                   }
1396                   else
1397                      init_value = 0;
1398                }
1399                else
1400                   init_value = 0;
1401             }
1402             result *= detail::lower_gamma_series(a, x, pol, init_value) / a;
1403             if(optimised_invert)
1404             {
1405                invert = false;
1406                result = -result;
1407             }
1408          }
1409          break;
1410       }
1411    case 3:
1412       {
1413          // Compute Q:
1414          invert = !invert;
1415          T g;
1416          result = tgamma_small_upper_part(a, x, pol, &g, invert, p_derivative);
1417          invert = false;
1418          if(normalised)
1419             result /= g;
1420          break;
1421       }
1422    case 4:
1423       {
1424          // Compute Q:
1425          result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
1426          if(p_derivative)
1427             *p_derivative = result;
1428          if(result != 0)
1429             result *= upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>());
1430          break;
1431       }
1432    case 5:
1433       {
1434          //
1435          // Use compile time dispatch to the appropriate
1436          // Temme asymptotic expansion.  This may be dead code
1437          // if T does not have numeric limits support, or has
1438          // too many digits for the most precise version of
1439          // these expansions, in that case we'll be calling
1440          // an empty function.
1441          //
1442          typedef typename policies::precision<T, Policy>::type precision_type;
1443
1444          typedef typename mpl::if_<
1445             mpl::or_<mpl::equal_to<precision_type, mpl::int_<0> >,
1446             mpl::greater<precision_type, mpl::int_<113> > >,
1447             mpl::int_<0>,
1448             typename mpl::if_<
1449                mpl::less_equal<precision_type, mpl::int_<53> >,
1450                mpl::int_<53>,
1451                typename mpl::if_<
1452                   mpl::less_equal<precision_type, mpl::int_<64> >,
1453                   mpl::int_<64>,
1454                   mpl::int_<113>
1455                >::type
1456             >::type
1457          >::type tag_type;
1458
1459          result = igamma_temme_large(a, x, pol, static_cast<tag_type const*>(0));
1460          if(x >= a)
1461             invert = !invert;
1462          if(p_derivative)
1463             *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
1464          break;
1465       }
1466    case 6:
1467       {
1468          // x is so small that P is necessarily very small too,
1469          // use http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/
1470          result = !normalised ? pow(x, a) / (a) : pow(x, a) / boost::math::tgamma(a + 1, pol);
1471          result *= 1 - a * x / (a + 1);
1472          if (p_derivative)
1473             *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
1474          break;
1475       }
1476    case 7:
1477    {
1478       // x is large,
1479       // Compute Q:
1480       result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
1481       if (p_derivative)
1482          *p_derivative = result;
1483       result /= x;
1484       if (result != 0)
1485          result *= incomplete_tgamma_large_x(a, x, pol);
1486       break;
1487    }
1488    }
1489
1490    if(normalised && (result > 1))
1491       result = 1;
1492    if(invert)
1493    {
1494       T gam = normalised ? 1 : boost::math::tgamma(a, pol);
1495       result = gam - result;
1496    }
1497    if(p_derivative)
1498    {
1499       //
1500       // Need to convert prefix term to derivative:
1501       //
1502       if((x < 1) && (tools::max_value<T>() * x < *p_derivative))
1503       {
1504          // overflow, just return an arbitrarily large value:
1505          *p_derivative = tools::max_value<T>() / 2;
1506       }
1507
1508       *p_derivative /= x;
1509    }
1510
1511    return result;
1512 }
1513
1514 //
1515 // Ratios of two gamma functions:
1516 //
1517 template <class T, class Policy, class Lanczos>
1518 T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const Lanczos& l)
1519 {
1520    BOOST_MATH_STD_USING
1521    if(z < tools::epsilon<T>())
1522    {
1523       //
1524       // We get spurious numeric overflow unless we're very careful, this
1525       // can occur either inside Lanczos::lanczos_sum(z) or in the
1526       // final combination of terms, to avoid this, split the product up
1527       // into 2 (or 3) parts:
1528       //
1529       // G(z) / G(L) = 1 / (z * G(L)) ; z < eps, L = z + delta = delta
1530       //    z * G(L) = z * G(lim) * (G(L)/G(lim)) ; lim = largest factorial
1531       //
1532       if(boost::math::max_factorial<T>::value < delta)
1533       {
1534          T ratio = tgamma_delta_ratio_imp_lanczos(delta, T(boost::math::max_factorial<T>::value - delta), pol, l);
1535          ratio *= z;
1536          ratio *= boost::math::unchecked_factorial<T>(boost::math::max_factorial<T>::value - 1);
1537          return 1 / ratio;
1538       }
1539       else
1540       {
1541          return 1 / (z * boost::math::tgamma(z + delta, pol));
1542       }
1543    }
1544    T zgh = static_cast<T>(z + Lanczos::g() - constants::half<T>());
1545    T result;
1546    if(z + delta == z)
1547    {
1548       if(fabs(delta) < 10)
1549          result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol));
1550       else
1551          result = 1;
1552    }
1553    else
1554    {
1555       if(fabs(delta) < 10)
1556       {
1557          result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol));
1558       }
1559       else
1560       {
1561          result = pow(zgh / (zgh + delta), z - constants::half<T>());
1562       }
1563       // Split the calculation up to avoid spurious overflow:
1564       result *= Lanczos::lanczos_sum(z) / Lanczos::lanczos_sum(T(z + delta));
1565    }
1566    result *= pow(constants::e<T>() / (zgh + delta), delta);
1567    return result;
1568 }
1569 //
1570 // And again without Lanczos support this time:
1571 //
1572 template <class T, class Policy>
1573 T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const lanczos::undefined_lanczos& l)
1574 {
1575    BOOST_MATH_STD_USING
1576
1577    //
1578    // We adjust z and delta so that both z and z+delta are large enough for
1579    // Sterling's approximation to hold.  We can then calculate the ratio
1580    // for the adjusted values, and rescale back down to z and z+delta.
1581    //
1582    // Get the required shifts first:
1583    //
1584    long numerator_shift = 0;
1585    long denominator_shift = 0;
1586    const int min_z = minimum_argument_for_bernoulli_recursion<T>();
1587
1588    if (min_z > z)
1589       numerator_shift = 1 + ltrunc(min_z - z);
1590    if (min_z > z + delta)
1591       denominator_shift = 1 + ltrunc(min_z - z - delta);
1592    //
1593    // If the shifts are zero, then we can just combine scaled tgamma's
1594    // and combine the remaining terms:
1595    //
1596    if (numerator_shift == 0 && denominator_shift == 0)
1597    {
1598       T scaled_tgamma_num = scaled_tgamma_no_lanczos(z, pol);
1599       T scaled_tgamma_denom = scaled_tgamma_no_lanczos(T(z + delta), pol);
1600       T result = scaled_tgamma_num / scaled_tgamma_denom;
1601       result *= exp(z * boost::math::log1p(-delta / (z + delta), pol)) * pow((delta + z) / constants::e<T>(), -delta);
1602       return result;
1603    }
1604    //
1605    // We're going to have to rescale first, get the adjusted z and delta values,
1606    // plus the ratio for the adjusted values:
1607    //
1608    T zz = z + numerator_shift;
1609    T dd = delta - (numerator_shift - denominator_shift);
1610    T ratio = tgamma_delta_ratio_imp_lanczos(zz, dd, pol, l);
1611    //
1612    // Use gamma recurrence relations to get back to the original
1613    // z and z+delta:
1614    //
1615    for (long long i = 0; i < numerator_shift; ++i)
1616    {
1617       ratio /= (z + i);
1618       if (i < denominator_shift)
1619          ratio *= (z + delta + i);
1620    }
1621    for (long long i = numerator_shift; i < denominator_shift; ++i)
1622    {
1623       ratio *= (z + delta + i);
1624    }
1625    return ratio;
1626 }
1627
1628 template <class T, class Policy>
1629 T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol)
1630 {
1631    BOOST_MATH_STD_USING
1632
1633    if((z <= 0) || (z + delta <= 0))
1634    {
1635       // This isn't very sofisticated, or accurate, but it does work:
1636       return boost::math::tgamma(z, pol) / boost::math::tgamma(z + delta, pol);
1637    }
1638
1639    if(floor(delta) == delta)
1640    {
1641       if(floor(z) == z)
1642       {
1643          //
1644          // Both z and delta are integers, see if we can just use table lookup
1645          // of the factorials to get the result:
1646          //
1647          if((z <= max_factorial<T>::value) && (z + delta <= max_factorial<T>::value))
1648          {
1649             return unchecked_factorial<T>((unsigned)itrunc(z, pol) - 1) / unchecked_factorial<T>((unsigned)itrunc(T(z + delta), pol) - 1);
1650          }
1651       }
1652       if(fabs(delta) < 20)
1653       {
1654          //
1655          // delta is a small integer, we can use a finite product:
1656          //
1657          if(delta == 0)
1658             return 1;
1659          if(delta < 0)
1660          {
1661             z -= 1;
1662             T result = z;
1663             while(0 != (delta += 1))
1664             {
1665                z -= 1;
1666                result *= z;
1667             }
1668             return result;
1669          }
1670          else
1671          {
1672             T result = 1 / z;
1673             while(0 != (delta -= 1))
1674             {
1675                z += 1;
1676                result /= z;
1677             }
1678             return result;
1679          }
1680       }
1681    }
1682    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1683    return tgamma_delta_ratio_imp_lanczos(z, delta, pol, lanczos_type());
1684 }
1685
1686 template <class T, class Policy>
1687 T tgamma_ratio_imp(T x, T y, const Policy& pol)
1688 {
1689    BOOST_MATH_STD_USING
1690
1691    if((x <= 0) || (boost::math::isinf)(x))
1692       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);
1693    if((y <= 0) || (boost::math::isinf)(y))
1694       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);
1695
1696    if(x <= tools::min_value<T>())
1697    {
1698       // Special case for denorms...Ugh.
1699       T shift = ldexp(T(1), tools::digits<T>());
1700       return shift * tgamma_ratio_imp(T(x * shift), y, pol);
1701    }
1702
1703    if((x < max_factorial<T>::value) && (y < max_factorial<T>::value))
1704    {
1705       // Rather than subtracting values, lets just call the gamma functions directly:
1706       return boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
1707    }
1708    T prefix = 1;
1709    if(x < 1)
1710    {
1711       if(y < 2 * max_factorial<T>::value)
1712       {
1713          // We need to sidestep on x as well, otherwise we'll underflow
1714          // before we get to factor in the prefix term:
1715          prefix /= x;
1716          x += 1;
1717          while(y >=  max_factorial<T>::value)
1718          {
1719             y -= 1;
1720             prefix /= y;
1721          }
1722          return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
1723       }
1724       //
1725       // result is almost certainly going to underflow to zero, try logs just in case:
1726       //
1727       return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
1728    }
1729    if(y < 1)
1730    {
1731       if(x < 2 * max_factorial<T>::value)
1732       {
1733          // We need to sidestep on y as well, otherwise we'll overflow
1734          // before we get to factor in the prefix term:
1735          prefix *= y;
1736          y += 1;
1737          while(x >= max_factorial<T>::value)
1738          {
1739             x -= 1;
1740             prefix *= x;
1741          }
1742          return prefix * boost::math::tgamma(x, pol) / boost::math::tgamma(y, pol);
1743       }
1744       //
1745       // Result will almost certainly overflow, try logs just in case:
1746       //
1747       return exp(boost::math::lgamma(x, pol) - boost::math::lgamma(y, pol));
1748    }
1749    //
1750    // Regular case, x and y both large and similar in magnitude:
1751    //
1752    return boost::math::tgamma_delta_ratio(x, y - x, pol);
1753 }
1754
1755 template <class T, class Policy>
1756 T gamma_p_derivative_imp(T a, T x, const Policy& pol)
1757 {
1758    BOOST_MATH_STD_USING
1759    //
1760    // Usual error checks first:
1761    //
1762    if(a <= 0)
1763       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);
1764    if(x < 0)
1765       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);
1766    //
1767    // Now special cases:
1768    //
1769    if(x == 0)
1770    {
1771       return (a > 1) ? 0 :
1772          (a == 1) ? 1 : policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
1773    }
1774    //
1775    // Normal case:
1776    //
1777    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1778    T f1 = detail::regularised_gamma_prefix(a, x, pol, lanczos_type());
1779    if((x < 1) && (tools::max_value<T>() * x < f1))
1780    {
1781       // overflow:
1782       return policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
1783    }
1784    if(f1 == 0)
1785    {
1786       // Underflow in calculation, use logs instead:
1787       f1 = a * log(x) - x - lgamma(a, pol) - log(x);
1788       f1 = exp(f1);
1789    }
1790    else
1791       f1 /= x;
1792
1793    return f1;
1794 }
1795
1796 template <class T, class Policy>
1797 inline typename tools::promote_args<T>::type
1798    tgamma(T z, const Policy& /* pol */, const mpl::true_)
1799 {
1800    BOOST_FPU_EXCEPTION_GUARD
1801    typedef typename tools::promote_args<T>::type result_type;
1802    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1803    typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1804    typedef typename policies::normalise<
1805       Policy,
1806       policies::promote_float<false>,
1807       policies::promote_double<false>,
1808       policies::discrete_quantile<>,
1809       policies::assert_undefined<> >::type forwarding_policy;
1810    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%)");
1811 }
1812
1813 template <class T, class Policy>
1814 struct igamma_initializer
1815 {
1816    struct init
1817    {
1818       init()
1819       {
1820          typedef typename policies::precision<T, Policy>::type precision_type;
1821
1822          typedef typename mpl::if_<
1823             mpl::or_<mpl::equal_to<precision_type, mpl::int_<0> >,
1824             mpl::greater<precision_type, mpl::int_<113> > >,
1825             mpl::int_<0>,
1826             typename mpl::if_<
1827                mpl::less_equal<precision_type, mpl::int_<53> >,
1828                mpl::int_<53>,
1829                typename mpl::if_<
1830                   mpl::less_equal<precision_type, mpl::int_<64> >,
1831                   mpl::int_<64>,
1832                   mpl::int_<113>
1833                >::type
1834             >::type
1835          >::type tag_type;
1836
1837          do_init(tag_type());
1838       }
1839       template <int N>
1840       static void do_init(const mpl::int_<N>&)
1841       {
1842          // If std::numeric_limits<T>::digits is zero, we must not call
1843          // our inituialization code here as the precision presumably
1844          // varies at runtime, and will not have been set yet.  Plus the
1845          // code requiring initialization isn't called when digits == 0.
1846          if(std::numeric_limits<T>::digits)
1847          {
1848             boost::math::gamma_p(static_cast<T>(400), static_cast<T>(400), Policy());
1849          }
1850       }
1851       static void do_init(const mpl::int_<53>&){}
1852       void force_instantiate()const{}
1853    };
1854    static const init initializer;
1855    static void force_instantiate()
1856    {
1857       initializer.force_instantiate();
1858    }
1859 };
1860
1861 template <class T, class Policy>
1862 const typename igamma_initializer<T, Policy>::init igamma_initializer<T, Policy>::initializer;
1863
1864 template <class T, class Policy>
1865 struct lgamma_initializer
1866 {
1867    struct init
1868    {
1869       init()
1870       {
1871          typedef typename policies::precision<T, Policy>::type precision_type;
1872          typedef typename mpl::if_<
1873             mpl::and_<
1874                mpl::less_equal<precision_type, mpl::int_<64> >,
1875                mpl::greater<precision_type, mpl::int_<0> >
1876             >,
1877             mpl::int_<64>,
1878             typename mpl::if_<
1879                mpl::and_<
1880                   mpl::less_equal<precision_type, mpl::int_<113> >,
1881                   mpl::greater<precision_type, mpl::int_<0> >
1882                >,
1883                mpl::int_<113>, mpl::int_<0> >::type
1884              >::type tag_type;
1885          do_init(tag_type());
1886       }
1887       static void do_init(const mpl::int_<64>&)
1888       {
1889          boost::math::lgamma(static_cast<T>(2.5), Policy());
1890          boost::math::lgamma(static_cast<T>(1.25), Policy());
1891          boost::math::lgamma(static_cast<T>(1.75), Policy());
1892       }
1893       static void do_init(const mpl::int_<113>&)
1894       {
1895          boost::math::lgamma(static_cast<T>(2.5), Policy());
1896          boost::math::lgamma(static_cast<T>(1.25), Policy());
1897          boost::math::lgamma(static_cast<T>(1.5), Policy());
1898          boost::math::lgamma(static_cast<T>(1.75), Policy());
1899       }
1900       static void do_init(const mpl::int_<0>&)
1901       {
1902       }
1903       void force_instantiate()const{}
1904    };
1905    static const init initializer;
1906    static void force_instantiate()
1907    {
1908       initializer.force_instantiate();
1909    }
1910 };
1911
1912 template <class T, class Policy>
1913 const typename lgamma_initializer<T, Policy>::init lgamma_initializer<T, Policy>::initializer;
1914
1915 template <class T1, class T2, class Policy>
1916 inline typename tools::promote_args<T1, T2>::type
1917    tgamma(T1 a, T2 z, const Policy&, const mpl::false_)
1918 {
1919    BOOST_FPU_EXCEPTION_GUARD
1920    typedef typename tools::promote_args<T1, T2>::type result_type;
1921    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1922    // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1923    typedef typename policies::normalise<
1924       Policy,
1925       policies::promote_float<false>,
1926       policies::promote_double<false>,
1927       policies::discrete_quantile<>,
1928       policies::assert_undefined<> >::type forwarding_policy;
1929
1930    igamma_initializer<value_type, forwarding_policy>::force_instantiate();
1931
1932    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
1933       detail::gamma_incomplete_imp(static_cast<value_type>(a),
1934       static_cast<value_type>(z), false, true,
1935       forwarding_policy(), static_cast<value_type*>(0)), "boost::math::tgamma<%1%>(%1%, %1%)");
1936 }
1937
1938 template <class T1, class T2>
1939 inline typename tools::promote_args<T1, T2>::type
1940    tgamma(T1 a, T2 z, const mpl::false_ tag)
1941 {
1942    return tgamma(a, z, policies::policy<>(), tag);
1943 }
1944
1945
1946 } // namespace detail
1947
1948 template <class T>
1949 inline typename tools::promote_args<T>::type
1950    tgamma(T z)
1951 {
1952    return tgamma(z, policies::policy<>());
1953 }
1954
1955 template <class T, class Policy>
1956 inline typename tools::promote_args<T>::type
1957    lgamma(T z, int* sign, const Policy&)
1958 {
1959    BOOST_FPU_EXCEPTION_GUARD
1960    typedef typename tools::promote_args<T>::type result_type;
1961    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1962    typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1963    typedef typename policies::normalise<
1964       Policy,
1965       policies::promote_float<false>,
1966       policies::promote_double<false>,
1967       policies::discrete_quantile<>,
1968       policies::assert_undefined<> >::type forwarding_policy;
1969
1970    detail::lgamma_initializer<value_type, forwarding_policy>::force_instantiate();
1971
1972    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%)");
1973 }
1974
1975 template <class T>
1976 inline typename tools::promote_args<T>::type
1977    lgamma(T z, int* sign)
1978 {
1979    return lgamma(z, sign, policies::policy<>());
1980 }
1981
1982 template <class T, class Policy>
1983 inline typename tools::promote_args<T>::type
1984    lgamma(T x, const Policy& pol)
1985 {
1986    return ::boost::math::lgamma(x, 0, pol);
1987 }
1988
1989 template <class T>
1990 inline typename tools::promote_args<T>::type
1991    lgamma(T x)
1992 {
1993    return ::boost::math::lgamma(x, 0, policies::policy<>());
1994 }
1995
1996 template <class T, class Policy>
1997 inline typename tools::promote_args<T>::type
1998    tgamma1pm1(T z, const Policy& /* pol */)
1999 {
2000    BOOST_FPU_EXCEPTION_GUARD
2001    typedef typename tools::promote_args<T>::type result_type;
2002    typedef typename policies::evaluation<result_type, Policy>::type value_type;
2003    typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
2004    typedef typename policies::normalise<
2005       Policy,
2006       policies::promote_float<false>,
2007       policies::promote_double<false>,
2008       policies::discrete_quantile<>,
2009       policies::assert_undefined<> >::type forwarding_policy;
2010
2011    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%)");
2012 }
2013
2014 template <class T>
2015 inline typename tools::promote_args<T>::type
2016    tgamma1pm1(T z)
2017 {
2018    return tgamma1pm1(z, policies::policy<>());
2019 }
2020
2021 //
2022 // Full upper incomplete gamma:
2023 //
2024 template <class T1, class T2>
2025 inline typename tools::promote_args<T1, T2>::type
2026    tgamma(T1 a, T2 z)
2027 {
2028    //
2029    // Type T2 could be a policy object, or a value, select the
2030    // right overload based on T2:
2031    //
2032    typedef typename policies::is_policy<T2>::type maybe_policy;
2033    return detail::tgamma(a, z, maybe_policy());
2034 }
2035 template <class T1, class T2, class Policy>
2036 inline typename tools::promote_args<T1, T2>::type
2037    tgamma(T1 a, T2 z, const Policy& pol)
2038 {
2039    return detail::tgamma(a, z, pol, mpl::false_());
2040 }
2041 //
2042 // Full lower incomplete gamma:
2043 //
2044 template <class T1, class T2, class Policy>
2045 inline typename tools::promote_args<T1, T2>::type
2046    tgamma_lower(T1 a, T2 z, const Policy&)
2047 {
2048    BOOST_FPU_EXCEPTION_GUARD
2049    typedef typename tools::promote_args<T1, T2>::type result_type;
2050    typedef typename policies::evaluation<result_type, Policy>::type value_type;
2051    // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
2052    typedef typename policies::normalise<
2053       Policy,
2054       policies::promote_float<false>,
2055       policies::promote_double<false>,
2056       policies::discrete_quantile<>,
2057       policies::assert_undefined<> >::type forwarding_policy;
2058
2059    detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
2060
2061    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
2062       detail::gamma_incomplete_imp(static_cast<value_type>(a),
2063       static_cast<value_type>(z), false, false,
2064       forwarding_policy(), static_cast<value_type*>(0)), "tgamma_lower<%1%>(%1%, %1%)");
2065 }
2066 template <class T1, class T2>
2067 inline typename tools::promote_args<T1, T2>::type
2068    tgamma_lower(T1 a, T2 z)
2069 {
2070    return tgamma_lower(a, z, policies::policy<>());
2071 }
2072 //
2073 // Regularised upper incomplete gamma:
2074 //
2075 template <class T1, class T2, class Policy>
2076 inline typename tools::promote_args<T1, T2>::type
2077    gamma_q(T1 a, T2 z, const Policy& /* pol */)
2078 {
2079    BOOST_FPU_EXCEPTION_GUARD
2080    typedef typename tools::promote_args<T1, T2>::type result_type;
2081    typedef typename policies::evaluation<result_type, Policy>::type value_type;
2082    // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
2083    typedef typename policies::normalise<
2084       Policy,
2085       policies::promote_float<false>,
2086       policies::promote_double<false>,
2087       policies::discrete_quantile<>,
2088       policies::assert_undefined<> >::type forwarding_policy;
2089
2090    detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
2091
2092    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
2093       detail::gamma_incomplete_imp(static_cast<value_type>(a),
2094       static_cast<value_type>(z), true, true,
2095       forwarding_policy(), static_cast<value_type*>(0)), "gamma_q<%1%>(%1%, %1%)");
2096 }
2097 template <class T1, class T2>
2098 inline typename tools::promote_args<T1, T2>::type
2099    gamma_q(T1 a, T2 z)
2100 {
2101    return gamma_q(a, z, policies::policy<>());
2102 }
2103 //
2104 // Regularised lower incomplete gamma:
2105 //
2106 template <class T1, class T2, class Policy>
2107 inline typename tools::promote_args<T1, T2>::type
2108    gamma_p(T1 a, T2 z, const Policy&)
2109 {
2110    BOOST_FPU_EXCEPTION_GUARD
2111    typedef typename tools::promote_args<T1, T2>::type result_type;
2112    typedef typename policies::evaluation<result_type, Policy>::type value_type;
2113    // typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
2114    typedef typename policies::normalise<
2115       Policy,
2116       policies::promote_float<false>,
2117       policies::promote_double<false>,
2118       policies::discrete_quantile<>,
2119       policies::assert_undefined<> >::type forwarding_policy;
2120
2121    detail::igamma_initializer<value_type, forwarding_policy>::force_instantiate();
2122
2123    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
2124       detail::gamma_incomplete_imp(static_cast<value_type>(a),
2125       static_cast<value_type>(z), true, false,
2126       forwarding_policy(), static_cast<value_type*>(0)), "gamma_p<%1%>(%1%, %1%)");
2127 }
2128 template <class T1, class T2>
2129 inline typename tools::promote_args<T1, T2>::type
2130    gamma_p(T1 a, T2 z)
2131 {
2132    return gamma_p(a, z, policies::policy<>());
2133 }
2134
2135 // ratios of gamma functions:
2136 template <class T1, class T2, class Policy>
2137 inline typename tools::promote_args<T1, T2>::type
2138    tgamma_delta_ratio(T1 z, T2 delta, const Policy& /* pol */)
2139 {
2140    BOOST_FPU_EXCEPTION_GUARD
2141    typedef typename tools::promote_args<T1, T2>::type result_type;
2142    typedef typename policies::evaluation<result_type, Policy>::type value_type;
2143    typedef typename policies::normalise<
2144       Policy,
2145       policies::promote_float<false>,
2146       policies::promote_double<false>,
2147       policies::discrete_quantile<>,
2148       policies::assert_undefined<> >::type forwarding_policy;
2149
2150    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_delta_ratio_imp(static_cast<value_type>(z), static_cast<value_type>(delta), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
2151 }
2152 template <class T1, class T2>
2153 inline typename tools::promote_args<T1, T2>::type
2154    tgamma_delta_ratio(T1 z, T2 delta)
2155 {
2156    return tgamma_delta_ratio(z, delta, policies::policy<>());
2157 }
2158 template <class T1, class T2, class Policy>
2159 inline typename tools::promote_args<T1, T2>::type
2160    tgamma_ratio(T1 a, T2 b, const Policy&)
2161 {
2162    typedef typename tools::promote_args<T1, T2>::type result_type;
2163    typedef typename policies::evaluation<result_type, Policy>::type value_type;
2164    typedef typename policies::normalise<
2165       Policy,
2166       policies::promote_float<false>,
2167       policies::promote_double<false>,
2168       policies::discrete_quantile<>,
2169       policies::assert_undefined<> >::type forwarding_policy;
2170
2171    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::tgamma_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(b), forwarding_policy()), "boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)");
2172 }
2173 template <class T1, class T2>
2174 inline typename tools::promote_args<T1, T2>::type
2175    tgamma_ratio(T1 a, T2 b)
2176 {
2177    return tgamma_ratio(a, b, policies::policy<>());
2178 }
2179
2180 template <class T1, class T2, class Policy>
2181 inline typename tools::promote_args<T1, T2>::type
2182    gamma_p_derivative(T1 a, T2 x, const Policy&)
2183 {
2184    BOOST_FPU_EXCEPTION_GUARD
2185    typedef typename tools::promote_args<T1, T2>::type result_type;
2186    typedef typename policies::evaluation<result_type, Policy>::type value_type;
2187    typedef typename policies::normalise<
2188       Policy,
2189       policies::promote_float<false>,
2190       policies::promote_double<false>,
2191       policies::discrete_quantile<>,
2192       policies::assert_undefined<> >::type forwarding_policy;
2193
2194    return policies::checked_narrowing_cast<result_type, forwarding_policy>(detail::gamma_p_derivative_imp(static_cast<value_type>(a), static_cast<value_type>(x), forwarding_policy()), "boost::math::gamma_p_derivative<%1%>(%1%, %1%)");
2195 }
2196 template <class T1, class T2>
2197 inline typename tools::promote_args<T1, T2>::type
2198    gamma_p_derivative(T1 a, T2 x)
2199 {
2200    return gamma_p_derivative(a, x, policies::policy<>());
2201 }
2202
2203 } // namespace math
2204 } // namespace boost
2205
2206 #ifdef BOOST_MSVC
2207 # pragma warning(pop)
2208 #endif
2209
2210 #include <boost/math/special_functions/detail/igamma_inverse.hpp>
2211 #include <boost/math/special_functions/detail/gamma_inva.hpp>
2212 #include <boost/math/special_functions/erf.hpp>
2213
2214 #endif // BOOST_MATH_SF_GAMMA_HPP
2215
2216
2217
2218