boost/math/special_functions/gamma.hpp

   1 //  Copyright John Maddock 2006-7.
   2 //  Copyright Paul A. Bristow 2007.
   3
   4 //  Use, modification and distribution are subject to the
   5 //  Boost Software License, Version 1.0. (See accompanying file
   6 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   7
   8 #ifndef BOOST_MATH_SF_GAMMA_HPP
   9 #define BOOST_MATH_SF_GAMMA_HPP
  10
  11 #ifdef _MSC_VER
  12 #pragma once
  13 #endif
  14
  15 #include <boost/config.hpp>
  16 #ifdef BOOST_MSVC
  17 # pragma warning(push)
  18 # pragma warning(disable: 4127 4701)
  19 //  // For lexical_cast, until fixed in 1.35?
  20 //  // conditional expression is constant &
  21 //  // Potentially uninitialized local variable 'name' used
  22 #endif
  23 #include <boost/lexical_cast.hpp>
  24 #ifdef BOOST_MSVC
  25 # pragma warning(pop)
  26 #endif
  27 #include <boost/math/tools/series.hpp>
  28 #include <boost/math/tools/fraction.hpp>
  29 #include <boost/math/tools/precision.hpp>
  30 #include <boost/math/tools/promotion.hpp>
  31 #include <boost/math/policies/error_handling.hpp>
  32 #include <boost/math/constants/constants.hpp>
  33 #include <boost/math/special_functions/math_fwd.hpp>
  34 #include <boost/math/special_functions/log1p.hpp>
  35 #include <boost/math/special_functions/trunc.hpp>
  36 #include <boost/math/special_functions/powm1.hpp>
  37 #include <boost/math/special_functions/sqrt1pm1.hpp>
  38 #include <boost/math/special_functions/lanczos.hpp>
  39 #include <boost/math/special_functions/fpclassify.hpp>
  40 #include <boost/math/special_functions/detail/igamma_large.hpp>
  41 #include <boost/math/special_functions/detail/unchecked_factorial.hpp>
  42 #include <boost/math/special_functions/detail/lgamma_small.hpp>
  43 #include <boost/type_traits/is_convertible.hpp>
  44 #include <boost/assert.hpp>
  45 #include <boost/mpl/greater.hpp>
  46 #include <boost/mpl/equal_to.hpp>
  47 #include <boost/mpl/greater.hpp>
  48
  49 #include <boost/config/no_tr1/cmath.hpp>
  50 #include <algorithm>
  51
  52 #ifdef BOOST_MATH_INSTRUMENT
  53 #include <iostream>
  54 #include <iomanip>
  55 #include <typeinfo>
  56 #endif
  57
  58 #ifdef BOOST_MSVC
  59 # pragma warning(push)
  60 # pragma warning(disable: 4702) // unreachable code (return after domain_error throw).
  61 # pragma warning(disable: 4127) // conditional expression is constant.
  62 # pragma warning(disable: 4100) // unreferenced formal parameter.
  63 // Several variables made comments,
  64 // but some difficulty as whether referenced on not may depend on macro values.
  65 // So to be safe, 4100 warnings suppressed.
  66 // TODO - revisit this?
  67 #endif
  68
  69 namespace boost{ namespace math{
  70
  71 namespace detail{
  72
  73 template <class T>
  74 inline bool is_odd(T v, const boost::true_type&)
  75 {
  76    int i = static_cast<int>(v);
  77    return i&1;
  78 }
  79 template <class T>
  80 inline bool is_odd(T v, const boost::false_type&)
  81 {
  82    // Oh dear can't cast T to int!
  83    BOOST_MATH_STD_USING
  84    T modulus = v - 2 * floor(v/2);
  85    return static_cast<bool>(modulus != 0);
  86 }
  87 template <class T>
  88 inline bool is_odd(T v)
  89 {
  90    return is_odd(v, ::boost::is_convertible<T, int>());
  91 }
  92
  93 template <class T>
  94 T sinpx(T z)
  95 {
  96    // Ad hoc function calculates x * sin(pi * x),
  97    // taking extra care near when x is near a whole number.
  98    BOOST_MATH_STD_USING
  99    int sign = 1;
 100    if(z < 0)
 101    {
 102       z = -z;
 103    }
 104    else
 105    {
 106       sign = -sign;
 107    }
 108    T fl = floor(z);
 109    T dist;
 110    if(is_odd(fl))
 111    {
 112       fl += 1;
 113       dist = fl - z;
 114       sign = -sign;
 115    }
 116    else
 117    {
 118       dist = z - fl;
 119    }
 120    BOOST_ASSERT(fl >= 0);
 121    if(dist > 0.5)
 122       dist = 1 - dist;
 123    T result = sin(dist*boost::math::constants::pi<T>());
 124    return sign*z*result;
 125 } // template <class T> T sinpx(T z)
 126 //
 127 // tgamma(z), with Lanczos support:
 128 //
 129 template <class T, class Policy, class L>
 130 T gamma_imp(T z, const Policy& pol, const L& l)
 131 {
 132    BOOST_MATH_STD_USING
 133
 134    T result = 1;
 135
 136 #ifdef BOOST_MATH_INSTRUMENT
 137    static bool b = false;
 138    if(!b)
 139    {
 140       std::cout << "tgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;
 141       b = true;
 142    }
 143 #endif
 144    static const char* function = "boost::math::tgamma<%1%>(%1%)";
 145
 146    if(z <= 0)
 147    {
 148       if(floor(z) == z)
 149          return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
 150       if(z <= -20)
 151       {
 152          result = gamma_imp(T(-z), pol, l) * sinpx(z);
 153          if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
 154             return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
 155          result = -boost::math::constants::pi<T>() / result;
 156          if(result == 0)
 157             return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
 158          if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL)
 159             return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", result, pol);
 160          return result;
 161       }
 162
 163       // shift z to > 1:
 164       while(z < 0)
 165       {
 166          result /= z;
 167          z += 1;
 168       }
 169    }
 170    if((floor(z) == z) && (z < max_factorial<T>::value))
 171    {
 172       result *= unchecked_factorial<T>(itrunc(z, pol) - 1);
 173    }
 174    else
 175    {
 176       result *= L::lanczos_sum(z);
 177       if(z * log(z) > tools::log_max_value<T>())
 178       {
 179          // we're going to overflow unless this is done with care:
 180          T zgh = (z + static_cast<T>(L::g()) - boost::math::constants::half<T>());
 181          if(log(zgh) * z / 2 > tools::log_max_value<T>())
 182             return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
 183          T hp = pow(zgh, (z / 2) - T(0.25));
 184          result *= hp / exp(zgh);
 185          if(tools::max_value<T>() / hp < result)
 186             return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
 187          result *= hp;
 188       }
 189       else
 190       {
 191          T zgh = (z + static_cast<T>(L::g()) - boost::math::constants::half<T>());
 192          result *= pow(zgh, z - boost::math::constants::half<T>()) / exp(zgh);
 193       }
 194    }
 195    return result;
 196 }
 197 //
 198 // lgamma(z) with Lanczos support:
 199 //
 200 template <class T, class Policy, class L>
 201 T lgamma_imp(T z, const Policy& pol, const L& l, int* sign = 0)
 202 {
 203 #ifdef BOOST_MATH_INSTRUMENT
 204    static bool b = false;
 205    if(!b)
 206    {
 207       std::cout << "lgamma_imp called with " << typeid(z).name() << " " << typeid(l).name() << std::endl;
 208       b = true;
 209    }
 210 #endif
 211
 212    BOOST_MATH_STD_USING
 213
 214    static const char* function = "boost::math::lgamma<%1%>(%1%)";
 215
 216    T result = 0;
 217    int sresult = 1;
 218    if(z <= 0)
 219    {
 220       // reflection formula:
 221       if(floor(z) == z)
 222          return policies::raise_pole_error<T>(function, "Evaluation of lgamma at a negative integer %1%.", z, pol);
 223
 224       T t = sinpx(z);
 225       z = -z;
 226       if(t < 0)
 227       {
 228          t = -t;
 229       }
 230       else
 231       {
 232          sresult = -sresult;
 233       }
 234       result = log(boost::math::constants::pi<T>()) - lgamma_imp(z, pol, l) - log(t);
 235    }
 236    else if(z < 15)
 237    {
 238       typedef typename policies::precision<T, Policy>::type precision_type;
 239       typedef typename mpl::if_<
 240          mpl::and_<
 241             mpl::less_equal<precision_type, mpl::int_<64> >,
 242             mpl::greater<precision_type, mpl::int_<0> >
 243          >,
 244          mpl::int_<64>,
 245          typename mpl::if_<
 246             mpl::and_<
 247                mpl::less_equal<precision_type, mpl::int_<113> >,
 248                mpl::greater<precision_type, mpl::int_<0> >
 249             >,
 250             mpl::int_<113>, mpl::int_<0> >::type
 251           >::type tag_type;
 252       result = lgamma_small_imp<T>(z, z - 1, z - 2, tag_type(), pol, l);
 253    }
 254    else if((z >= 3) && (z < 100))
 255    {
 256       // taking the log of tgamma reduces the error, no danger of overflow here:
 257       result = log(gamma_imp(z, pol, l));
 258    }
 259    else
 260    {
 261       // regular evaluation:
 262       T zgh = static_cast<T>(z + L::g() - boost::math::constants::half<T>());
 263       result = log(zgh) - 1;
 264       result *= z - 0.5f;
 265       result += log(L::lanczos_sum_expG_scaled(z));
 266    }
 267
 268    if(sign)
 269       *sign = sresult;
 270    return result;
 271 }
 272
 273 //
 274 // Incomplete gamma functions follow:
 275 //
 276 template <class T>
 277 struct upper_incomplete_gamma_fract
 278 {
 279 private:
 280    T z, a;
 281    int k;
 282 public:
 283    typedef std::pair<T,T> result_type;
 284
 285    upper_incomplete_gamma_fract(T a1, T z1)
 286       : z(z1-a1+1), a(a1), k(0)
 287    {
 288    }
 289
 290    result_type operator()()
 291    {
 292       ++k;
 293       z += 2;
 294       return result_type(k * (a - k), z);
 295    }
 296 };
 297
 298 template <class T>
 299 inline T upper_gamma_fraction(T a, T z, T eps)
 300 {
 301    // Multiply result by z^a * e^-z to get the full
 302    // upper incomplete integral.  Divide by tgamma(z)
 303    // to normalise.
 304    upper_incomplete_gamma_fract<T> f(a, z);
 305    return 1 / (z - a + 1 + boost::math::tools::continued_fraction_a(f, eps));
 306 }
 307
 308 template <class T>
 309 struct lower_incomplete_gamma_series
 310 {
 311 private:
 312    T a, z, result;
 313 public:
 314    typedef T result_type;
 315    lower_incomplete_gamma_series(T a1, T z1) : a(a1), z(z1), result(1){}
 316
 317    T operator()()
 318    {
 319       T r = result;
 320       a += 1;
 321       result *= z/a;
 322       return r;
 323    }
 324 };
 325
 326 template <class T, class Policy>
 327 inline T lower_gamma_series(T a, T z, const Policy& pol, T init_value = 0)
 328 {
 329    // Multiply result by ((z^a) * (e^-z) / a) to get the full
 330    // lower incomplete integral. Then divide by tgamma(a)
 331    // to get the normalised value.
 332    lower_incomplete_gamma_series<T> s(a, z);
 333    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
 334    T factor = policies::get_epsilon<T, Policy>();
 335    T result = boost::math::tools::sum_series(s, factor, max_iter, init_value);
 336    policies::check_series_iterations("boost::math::detail::lower_gamma_series<%1%>(%1%)", max_iter, pol);
 337    return result;
 338 }
 339
 340 //
 341 // Fully generic tgamma and lgamma use the incomplete partial
 342 // sums added together:
 343 //
 344 template <class T, class Policy>
 345 T gamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l)
 346 {
 347    static const char* function = "boost::math::tgamma<%1%>(%1%)";
 348    BOOST_MATH_STD_USING
 349    if((z <= 0) && (floor(z) == z))
 350       return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
 351    if(z <= -20)
 352    {
 353       T result = gamma_imp(-z, pol, l) * sinpx(z);
 354       if((fabs(result) < 1) && (tools::max_value<T>() * fabs(result) < boost::math::constants::pi<T>()))
 355          return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
 356       result = -boost::math::constants::pi<T>() / result;
 357       if(result == 0)
 358          return policies::raise_underflow_error<T>(function, "Result of tgamma is too small to represent.", pol);
 359       if((boost::math::fpclassify)(result) == (int)FP_SUBNORMAL)
 360          return policies::raise_denorm_error<T>(function, "Result of tgamma is denormalized.", result, pol);
 361       return result;
 362    }
 363    //
 364    // The upper gamma fraction is *very* slow for z < 6, actually it's very
 365    // slow to converge everywhere but recursing until z > 6 gets rid of the
 366    // worst of it's behaviour.
 367    //
 368    T prefix = 1;
 369    while(z < 6)
 370    {
 371       prefix /= z;
 372       z += 1;
 373    }
 374    BOOST_MATH_INSTRUMENT_CODE(prefix);
 375    if((floor(z) == z) && (z < max_factorial<T>::value))
 376    {
 377       prefix *= unchecked_factorial<T>(itrunc(z, pol) - 1);
 378    }
 379    else
 380    {
 381       prefix = prefix * pow(z / boost::math::constants::e<T>(), z);
 382       BOOST_MATH_INSTRUMENT_CODE(prefix);
 383       T sum = detail::lower_gamma_series(z, z, pol) / z;
 384       BOOST_MATH_INSTRUMENT_CODE(sum);
 385       sum += detail::upper_gamma_fraction(z, z, ::boost::math::policies::get_epsilon<T, Policy>());
 386       BOOST_MATH_INSTRUMENT_CODE(sum);
 387       if(fabs(tools::max_value<T>() / prefix) < fabs(sum))
 388          return policies::raise_overflow_error<T>(function, "Result of tgamma is too large to represent.", pol);
 389       BOOST_MATH_INSTRUMENT_CODE((sum * prefix));
 390       return sum * prefix;
 391    }
 392    return prefix;
 393 }
 394
 395 template <class T, class Policy>
 396 T lgamma_imp(T z, const Policy& pol, const lanczos::undefined_lanczos& l, int*sign)
 397 {
 398    BOOST_MATH_STD_USING
 399
 400    static const char* function = "boost::math::lgamma<%1%>(%1%)";
 401    T result = 0;
 402    int sresult = 1;
 403    if(z <= 0)
 404    {
 405       if(floor(z) == z)
 406          return policies::raise_pole_error<T>(function, "Evaluation of tgamma at a negative integer %1%.", z, pol);
 407       T t = detail::sinpx(z);
 408       z = -z;
 409       if(t < 0)
 410       {
 411          t = -t;
 412       }
 413       else
 414       {
 415          sresult = -sresult;
 416       }
 417       result = log(boost::math::constants::pi<T>()) - lgamma_imp(z, pol, l, 0) - log(t);
 418    }
 419    else if((z != 1) && (z != 2))
 420    {
 421       T limit = (std::max)(z+1, T(10));
 422       T prefix = z * log(limit) - limit;
 423       T sum = detail::lower_gamma_series(z, limit, pol) / z;
 424       sum += detail::upper_gamma_fraction(z, limit, ::boost::math::policies::get_epsilon<T, Policy>());
 425       result = log(sum) + prefix;
 426    }
 427    if(sign)
 428       *sign = sresult;
 429    return result;
 430 }
 431 //
 432 // This helper calculates tgamma(dz+1)-1 without cancellation errors,
 433 // used by the upper incomplete gamma with z < 1:
 434 //
 435 template <class T, class Policy, class L>
 436 T tgammap1m1_imp(T dz, Policy const& pol, const L& l)
 437 {
 438    BOOST_MATH_STD_USING
 439
 440    typedef typename policies::precision<T,Policy>::type precision_type;
 441
 442    typedef typename mpl::if_<
 443       mpl::or_<
 444          mpl::less_equal<precision_type, mpl::int_<0> >,
 445          mpl::greater<precision_type, mpl::int_<113> >
 446       >,
 447       typename mpl::if_<
 448          is_same<L, lanczos::lanczos24m113>,
 449          mpl::int_<113>,
 450          mpl::int_<0>
 451       >::type,
 452       typename mpl::if_<
 453          mpl::less_equal<precision_type, mpl::int_<64> >,
 454          mpl::int_<64>, mpl::int_<113> >::type
 455        >::type tag_type;
 456
 457    T result;
 458    if(dz < 0)
 459    {
 460       if(dz < -0.5)
 461       {
 462          // Best method is simply to subtract 1 from tgamma:
 463          result = boost::math::tgamma(1+dz, pol) - 1;
 464          BOOST_MATH_INSTRUMENT_CODE(result);
 465       }
 466       else
 467       {
 468          // Use expm1 on lgamma:
 469          result = boost::math::expm1(-boost::math::log1p(dz, pol)
 470             + lgamma_small_imp<T>(dz+2, dz + 1, dz, tag_type(), pol, l));
 471          BOOST_MATH_INSTRUMENT_CODE(result);
 472       }
 473    }
 474    else
 475    {
 476       if(dz < 2)
 477       {
 478          // Use expm1 on lgamma:
 479          result = boost::math::expm1(lgamma_small_imp<T>(dz+1, dz, dz-1, tag_type(), pol, l), pol);
 480          BOOST_MATH_INSTRUMENT_CODE(result);
 481       }
 482       else
 483       {
 484          // Best method is simply to subtract 1 from tgamma:
 485          result = boost::math::tgamma(1+dz, pol) - 1;
 486          BOOST_MATH_INSTRUMENT_CODE(result);
 487       }
 488    }
 489
 490    return result;
 491 }
 492
 493 template <class T, class Policy>
 494 inline T tgammap1m1_imp(T dz, Policy const& pol,
 495                  const ::boost::math::lanczos::undefined_lanczos& l)
 496 {
 497    BOOST_MATH_STD_USING // ADL of std names
 498    //
 499    // There should be a better solution than this, but the
 500    // algebra isn't easy for the general case....
 501    // Start by subracting 1 from tgamma:
 502    //
 503    T result = gamma_imp(1 + dz, pol, l) - 1;
 504    BOOST_MATH_INSTRUMENT_CODE(result);
 505    //
 506    // Test the level of cancellation error observed: we loose one bit
 507    // for each power of 2 the result is less than 1.  If we would get
 508    // more bits from our most precise lgamma rational approximation,
 509    // then use that instead:
 510    //
 511    BOOST_MATH_INSTRUMENT_CODE((dz > -0.5));
 512    BOOST_MATH_INSTRUMENT_CODE((dz < 2));
 513    BOOST_MATH_INSTRUMENT_CODE((ldexp(1.0, boost::math::policies::digits<T, Policy>()) * fabs(result) < 1e34));
 514    if((dz > -0.5) && (dz < 2) && (ldexp(1.0, boost::math::policies::digits<T, Policy>()) * fabs(result) < 1e34))
 515    {
 516       result = tgammap1m1_imp(dz, pol, boost::math::lanczos::lanczos24m113());
 517       BOOST_MATH_INSTRUMENT_CODE(result);
 518    }
 519    return result;
 520 }
 521
 522 //
 523 // Series representation for upper fraction when z is small:
 524 //
 525 template <class T>
 526 struct small_gamma2_series
 527 {
 528    typedef T result_type;
 529
 530    small_gamma2_series(T a_, T x_) : result(-x_), x(-x_), apn(a_+1), n(1){}
 531
 532    T operator()()
 533    {
 534       T r = result / (apn);
 535       result *= x;
 536       result /= ++n;
 537       apn += 1;
 538       return r;
 539    }
 540
 541 private:
 542    T result, x, apn;
 543    int n;
 544 };
 545 //
 546 // calculate power term prefix (z^a)(e^-z) used in the non-normalised
 547 // incomplete gammas:
 548 //
 549 template <class T, class Policy>
 550 T full_igamma_prefix(T a, T z, const Policy& pol)
 551 {
 552    BOOST_MATH_STD_USING
 553
 554    T prefix;
 555    T alz = a * log(z);
 556
 557    if(z >= 1)
 558    {
 559       if((alz < tools::log_max_value<T>()) && (-z > tools::log_min_value<T>()))
 560       {
 561          prefix = pow(z, a) * exp(-z);
 562       }
 563       else if(a >= 1)
 564       {
 565          prefix = pow(z / exp(z/a), a);
 566       }
 567       else
 568       {
 569          prefix = exp(alz - z);
 570       }
 571    }
 572    else
 573    {
 574       if(alz > tools::log_min_value<T>())
 575       {
 576          prefix = pow(z, a) * exp(-z);
 577       }
 578       else if(z/a < tools::log_max_value<T>())
 579       {
 580          prefix = pow(z / exp(z/a), a);
 581       }
 582       else
 583       {
 584          prefix = exp(alz - z);
 585       }
 586    }
 587    //
 588    // This error handling isn't very good: it happens after the fact
 589    // rather than before it...
 590    //
 591    if((boost::math::fpclassify)(prefix) == (int)FP_INFINITE)
 592       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);
 593
 594    return prefix;
 595 }
 596 //
 597 // Compute (z^a)(e^-z)/tgamma(a)
 598 // most if the error occurs in this function:
 599 //
 600 template <class T, class Policy, class L>
 601 T regularised_gamma_prefix(T a, T z, const Policy& pol, const L& l)
 602 {
 603    BOOST_MATH_STD_USING
 604    T agh = a + static_cast<T>(L::g()) - T(0.5);
 605    T prefix;
 606    T d = ((z - a) - static_cast<T>(L::g()) + T(0.5)) / agh;
 607
 608    if(a < 1)
 609    {
 610       //
 611       // We have to treat a < 1 as a special case because our Lanczos
 612       // approximations are optimised against the factorials with a > 1,
 613       // and for high precision types especially (128-bit reals for example)
 614       // very small values of a can give rather eroneous results for gamma
 615       // unless we do this:
 616       //
 617       // TODO: is this still required?  Lanczos approx should be better now?
 618       //
 619       if(z <= tools::log_min_value<T>())
 620       {
 621          // Oh dear, have to use logs, should be free of cancellation errors though:
 622          return exp(a * log(z) - z - lgamma_imp(a, pol, l));
 623       }
 624       else
 625       {
 626          // direct calculation, no danger of overflow as gamma(a) < 1/a
 627          // for small a.
 628          return pow(z, a) * exp(-z) / gamma_imp(a, pol, l);
 629       }
 630    }
 631    else if((fabs(d*d*a) <= 100) && (a > 150))
 632    {
 633       // special case for large a and a ~ z.
 634       prefix = a * boost::math::log1pmx(d, pol) + z * static_cast<T>(0.5 - L::g()) / agh;
 635       prefix = exp(prefix);
 636    }
 637    else
 638    {
 639       //
 640       // general case.
 641       // direct computation is most accurate, but use various fallbacks
 642       // for different parts of the problem domain:
 643       //
 644       T alz = a * log(z / agh);
 645       T amz = a - z;
 646       if(((std::min)(alz, amz) <= tools::log_min_value<T>()) || ((std::max)(alz, amz) >= tools::log_max_value<T>()))
 647       {
 648          T amza = amz / a;
 649          if(((std::min)(alz, amz)/2 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/2 < tools::log_max_value<T>()))
 650          {
 651             // compute square root of the result and then square it:
 652             T sq = pow(z / agh, a / 2) * exp(amz / 2);
 653             prefix = sq * sq;
 654          }
 655          else if(((std::min)(alz, amz)/4 > tools::log_min_value<T>()) && ((std::max)(alz, amz)/4 < tools::log_max_value<T>()) && (z > a))
 656          {
 657             // compute the 4th root of the result then square it twice:
 658             T sq = pow(z / agh, a / 4) * exp(amz / 4);
 659             prefix = sq * sq;
 660             prefix *= prefix;
 661          }
 662          else if((amza > tools::log_min_value<T>()) && (amza < tools::log_max_value<T>()))
 663          {
 664             prefix = pow((z * exp(amza)) / agh, a);
 665          }
 666          else
 667          {
 668             prefix = exp(alz + amz);
 669          }
 670       }
 671       else
 672       {
 673          prefix = pow(z / agh, a) * exp(amz);
 674       }
 675    }
 676    prefix *= sqrt(agh / boost::math::constants::e<T>()) / L::lanczos_sum_expG_scaled(a);
 677    return prefix;
 678 }
 679 //
 680 // And again, without Lanczos support:
 681 //
 682 template <class T, class Policy>
 683 T regularised_gamma_prefix(T a, T z, const Policy& pol, const lanczos::undefined_lanczos&)
 684 {
 685    BOOST_MATH_STD_USING
 686
 687    T limit = (std::max)(T(10), a);
 688    T sum = detail::lower_gamma_series(a, limit, pol) / a;
 689    sum += detail::upper_gamma_fraction(a, limit, ::boost::math::policies::get_epsilon<T, Policy>());
 690
 691    if(a < 10)
 692    {
 693       // special case for small a:
 694       T prefix = pow(z / 10, a);
 695       prefix *= exp(10-z);
 696       if(0 == prefix)
 697       {
 698          prefix = pow((z * exp((10-z)/a)) / 10, a);
 699       }
 700       prefix /= sum;
 701       return prefix;
 702    }
 703
 704    T zoa = z / a;
 705    T amz = a - z;
 706    T alzoa = a * log(zoa);
 707    T prefix;
 708    if(((std::min)(alzoa, amz) <= tools::log_min_value<T>()) || ((std::max)(alzoa, amz) >= tools::log_max_value<T>()))
 709    {
 710       T amza = amz / a;
 711       if((amza <= tools::log_min_value<T>()) || (amza >= tools::log_max_value<T>()))
 712       {
 713          prefix = exp(alzoa + amz);
 714       }
 715       else
 716       {
 717          prefix = pow(zoa * exp(amza), a);
 718       }
 719    }
 720    else
 721    {
 722       prefix = pow(zoa, a) * exp(amz);
 723    }
 724    prefix /= sum;
 725    return prefix;
 726 }
 727 //
 728 // Upper gamma fraction for very small a:
 729 //
 730 template <class T, class Policy>
 731 inline T tgamma_small_upper_part(T a, T x, const Policy& pol, T* pgam = 0, bool invert = false, T* pderivative = 0)
 732 {
 733    BOOST_MATH_STD_USING  // ADL of std functions.
 734    //
 735    // Compute the full upper fraction (Q) when a is very small:
 736    //
 737    T result;
 738    result = boost::math::tgamma1pm1(a, pol);
 739    if(pgam)
 740       *pgam = (result + 1) / a;
 741    T p = boost::math::powm1(x, a, pol);
 742    result -= p;
 743    result /= a;
 744    detail::small_gamma2_series<T> s(a, x);
 745    boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>() - 10;
 746    p += 1;
 747    if(pderivative)
 748       *pderivative = p / (*pgam * exp(x));
 749    T init_value = invert ? *pgam : 0;
 750    result = -p * tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, (init_value - result) / p);
 751    policies::check_series_iterations("boost::math::tgamma_small_upper_part<%1%>(%1%, %1%)", max_iter, pol);
 752    if(invert)
 753       result = -result;
 754    return result;
 755 }
 756 //
 757 // Upper gamma fraction for integer a:
 758 //
 759 template <class T, class Policy>
 760 inline T finite_gamma_q(T a, T x, Policy const& pol, T* pderivative = 0)
 761 {
 762    //
 763    // Calculates normalised Q when a is an integer:
 764    //
 765    BOOST_MATH_STD_USING
 766    T e = exp(-x);
 767    T sum = e;
 768    if(sum != 0)
 769    {
 770       T term = sum;
 771       for(unsigned n = 1; n < a; ++n)
 772       {
 773          term /= n;
 774          term *= x;
 775          sum += term;
 776       }
 777    }
 778    if(pderivative)
 779    {
 780       *pderivative = e * pow(x, a) / boost::math::unchecked_factorial<T>(itrunc(T(a - 1), pol));
 781    }
 782    return sum;
 783 }
 784 //
 785 // Upper gamma fraction for half integer a:
 786 //
 787 template <class T, class Policy>
 788 T finite_half_gamma_q(T a, T x, T* p_derivative, const Policy& pol)
 789 {
 790    //
 791    // Calculates normalised Q when a is a half-integer:
 792    //
 793    BOOST_MATH_STD_USING
 794    T e = boost::math::erfc(sqrt(x), pol);
 795    if((e != 0) && (a > 1))
 796    {
 797       T term = exp(-x) / sqrt(constants::pi<T>() * x);
 798       term *= x;
 799       static const T half = T(1) / 2;
 800       term /= half;
 801       T sum = term;
 802       for(unsigned n = 2; n < a; ++n)
 803       {
 804          term /= n - half;
 805          term *= x;
 806          sum += term;
 807       }
 808       e += sum;
 809       if(p_derivative)
 810       {
 811          *p_derivative = 0;
 812       }
 813    }
 814    else if(p_derivative)
 815    {
 816       // We'll be dividing by x later, so calculate derivative * x:
 817       *p_derivative = sqrt(x) * exp(-x) / constants::root_pi<T>();
 818    }
 819    return e;
 820 }
 821 //
 822 // Main incomplete gamma entry point, handles all four incomplete gamma's:
 823 //
 824 template <class T, class Policy>
 825 T gamma_incomplete_imp(T a, T x, bool normalised, bool invert,
 826                        const Policy& pol, T* p_derivative)
 827 {
 828    static const char* function = "boost::math::gamma_p<%1%>(%1%, %1%)";
 829    if(a <= 0)
 830       policies::raise_domain_error<T>(function, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a, pol);
 831    if(x < 0)
 832       policies::raise_domain_error<T>(function, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x, pol);
 833
 834    BOOST_MATH_STD_USING
 835
 836    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
 837
 838    T result;
 839
 840    BOOST_ASSERT((p_derivative == 0) || (normalised == true));
 841
 842    bool is_int, is_half_int;
 843    bool is_small_a = (a < 30) && (a <= x + 1);
 844    if(is_small_a)
 845    {
 846       T fa = floor(a);
 847       is_int = (fa == a);
 848       is_half_int = is_int ? false : (fabs(fa - a) == 0.5f);
 849    }
 850    else
 851    {
 852       is_int = is_half_int = false;
 853    }
 854
 855    int eval_method;
 856
 857    if(is_int && (x > 0.6))
 858    {
 859       // calculate Q via finite sum:
 860       invert = !invert;
 861       eval_method = 0;
 862    }
 863    else if(is_half_int && (x > 0.2))
 864    {
 865       // calculate Q via finite sum for half integer a:
 866       invert = !invert;
 867       eval_method = 1;
 868    }
 869    else if(x < 0.5)
 870    {
 871       //
 872       // Changeover criterion chosen to give a changeover at Q ~ 0.33
 873       //
 874       if(-0.4 / log(x) < a)
 875       {
 876          eval_method = 2;
 877       }
 878       else
 879       {
 880          eval_method = 3;
 881       }
 882    }
 883    else if(x < 1.1)
 884    {
 885       //
 886       // Changover here occurs when P ~ 0.75 or Q ~ 0.25:
 887       //
 888       if(x * 0.75f < a)
 889       {
 890          eval_method = 2;
 891       }
 892       else
 893       {
 894          eval_method = 3;
 895       }
 896    }
 897    else
 898    {
 899       //
 900       // Begin by testing whether we're in the "bad" zone
 901       // where the result will be near 0.5 and the usual
 902       // series and continued fractions are slow to converge:
 903       //
 904       bool use_temme = false;
 905       if(normalised && std::numeric_limits<T>::is_specialized && (a > 20))
 906       {
 907          T sigma = fabs((x-a)/a);
 908          if((a > 200) && (policies::digits<T, Policy>() <= 113))
 909          {
 910             //
 911             // This limit is chosen so that we use Temme's expansion
 912             // only if the result would be larger than about 10^-6.
 913             // Below that the regular series and continued fractions
 914             // converge OK, and if we use Temme's method we get increasing
 915             // errors from the dominant erfc term as it's (inexact) argument
 916             // increases in magnitude.
 917             //
 918             if(20 / a > sigma * sigma)
 919                use_temme = true;
 920          }
 921          else if(policies::digits<T, Policy>() <= 64)
 922          {
 923             // Note in this zone we can't use Temme's expansion for
 924             // types longer than an 80-bit real:
 925             // it would require too many terms in the polynomials.
 926             if(sigma < 0.4)
 927                use_temme = true;
 928          }
 929       }
 930       if(use_temme)
 931       {
 932          eval_method = 5;
 933       }
 934       else
 935       {
 936          //
 937          // Regular case where the result will not be too close to 0.5.
 938          //
 939          // Changeover here occurs at P ~ Q ~ 0.5
 940          // Note that series computation of P is about x2 faster than continued fraction
 941          // calculation of Q, so try and use the CF only when really necessary, especially
 942          // for small x.
 943          //
 944          if(x - (1 / (3 * x)) < a)
 945          {
 946             eval_method = 2;
 947          }
 948          else
 949          {
 950             eval_method = 4;
 951             invert = !invert;
 952          }
 953       }
 954    }
 955
 956    switch(eval_method)
 957    {
 958    case 0:
 959       {
 960          result = finite_gamma_q(a, x, pol, p_derivative);
 961          if(normalised == false)
 962             result *= boost::math::tgamma(a, pol);
 963          break;
 964       }
 965    case 1:
 966       {
 967          result = finite_half_gamma_q(a, x, p_derivative, pol);
 968          if(normalised == false)
 969             result *= boost::math::tgamma(a, pol);
 970          if(p_derivative && (*p_derivative == 0))
 971             *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
 972          break;
 973       }
 974    case 2:
 975       {
 976          // Compute P:
 977          result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
 978          if(p_derivative)
 979             *p_derivative = result;
 980          if(result != 0)
 981          {
 982             T init_value = 0;
 983             if(invert)
 984             {
 985                init_value = -a * (normalised ? 1 : boost::math::tgamma(a, pol)) / result;
 986             }
 987             result *= detail::lower_gamma_series(a, x, pol, init_value) / a;
 988             if(invert)
 989             {
 990                invert = false;
 991                result = -result;
 992             }
 993          }
 994          break;
 995       }
 996    case 3:
 997       {
 998          // Compute Q:
 999          invert = !invert;
1000          T g;
1001          result = tgamma_small_upper_part(a, x, pol, &g, invert, p_derivative);
1002          invert = false;
1003          if(normalised)
1004             result /= g;
1005          break;
1006       }
1007    case 4:
1008       {
1009          // Compute Q:
1010          result = normalised ? regularised_gamma_prefix(a, x, pol, lanczos_type()) : full_igamma_prefix(a, x, pol);
1011          if(p_derivative)
1012             *p_derivative = result;
1013          if(result != 0)
1014             result *= upper_gamma_fraction(a, x, policies::get_epsilon<T, Policy>());
1015          break;
1016       }
1017    case 5:
1018       {
1019          //
1020          // Use compile time dispatch to the appropriate
1021          // Temme asymptotic expansion.  This may be dead code
1022          // if T does not have numeric limits support, or has
1023          // too many digits for the most precise version of
1024          // these expansions, in that case we'll be calling
1025          // an empty function.
1026          //
1027          typedef typename policies::precision<T, Policy>::type precision_type;
1028
1029          typedef typename mpl::if_<
1030             mpl::or_<mpl::equal_to<precision_type, mpl::int_<0> >,
1031             mpl::greater<precision_type, mpl::int_<113> > >,
1032             mpl::int_<0>,
1033             typename mpl::if_<
1034                mpl::less_equal<precision_type, mpl::int_<53> >,
1035                mpl::int_<53>,
1036                typename mpl::if_<
1037                   mpl::less_equal<precision_type, mpl::int_<64> >,
1038                   mpl::int_<64>,
1039                   mpl::int_<113>
1040                >::type
1041             >::type
1042          >::type tag_type;
1043
1044          result = igamma_temme_large(a, x, pol, static_cast<tag_type const*>(0));
1045          if(x >= a)
1046             invert = !invert;
1047          if(p_derivative)
1048             *p_derivative = regularised_gamma_prefix(a, x, pol, lanczos_type());
1049          break;
1050       }
1051    }
1052
1053    if(normalised && (result > 1))
1054       result = 1;
1055    if(invert)
1056    {
1057       T gam = normalised ? 1 : boost::math::tgamma(a, pol);
1058       result = gam - result;
1059    }
1060    if(p_derivative)
1061    {
1062       //
1063       // Need to convert prefix term to derivative:
1064       //
1065       if((x < 1) && (tools::max_value<T>() * x < *p_derivative))
1066       {
1067          // overflow, just return an arbitrarily large value:
1068          *p_derivative = tools::max_value<T>() / 2;
1069       }
1070
1071       *p_derivative /= x;
1072    }
1073
1074    return result;
1075 }
1076
1077 //
1078 // Ratios of two gamma functions:
1079 //
1080 template <class T, class Policy, class L>
1081 T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const L&)
1082 {
1083    BOOST_MATH_STD_USING
1084    T zgh = z + L::g() - constants::half<T>();
1085    T result;
1086    if(fabs(delta) < 10)
1087    {
1088       result = exp((constants::half<T>() - z) * boost::math::log1p(delta / zgh, pol));
1089    }
1090    else
1091    {
1092       result = pow(zgh / (zgh + delta), z - constants::half<T>());
1093    }
1094    result *= pow(constants::e<T>() / (zgh + delta), delta);
1095    result *= L::lanczos_sum(z) / L::lanczos_sum(z + delta);
1096    return result;
1097 }
1098 //
1099 // And again without Lanczos support this time:
1100 //
1101 template <class T, class Policy>
1102 T tgamma_delta_ratio_imp_lanczos(T z, T delta, const Policy& pol, const lanczos::undefined_lanczos&)
1103 {
1104    BOOST_MATH_STD_USING
1105    //
1106    // The upper gamma fraction is *very* slow for z < 6, actually it's very
1107    // slow to converge everywhere but recursing until z > 6 gets rid of the
1108    // worst of it's behaviour.
1109    //
1110    T prefix = 1;
1111    T zd = z + delta;
1112    while((zd < 6) && (z < 6))
1113    {
1114       prefix /= z;
1115       prefix *= zd;
1116       z += 1;
1117       zd += 1;
1118    }
1119    if(delta < 10)
1120    {
1121       prefix *= exp(-z * boost::math::log1p(delta / z, pol));
1122    }
1123    else
1124    {
1125       prefix *= pow(z / zd, z);
1126    }
1127    prefix *= pow(constants::e<T>() / zd, delta);
1128    T sum = detail::lower_gamma_series(z, z, pol) / z;
1129    sum += detail::upper_gamma_fraction(z, z, ::boost::math::policies::get_epsilon<T, Policy>());
1130    T sumd = detail::lower_gamma_series(zd, zd, pol) / zd;
1131    sumd += detail::upper_gamma_fraction(zd, zd, ::boost::math::policies::get_epsilon<T, Policy>());
1132    sum /= sumd;
1133    if(fabs(tools::max_value<T>() / prefix) < fabs(sum))
1134       return policies::raise_overflow_error<T>("boost::math::tgamma_delta_ratio<%1%>(%1%, %1%)", "Result of tgamma is too large to represent.", pol);
1135    return sum * prefix;
1136 }
1137
1138 template <class T, class Policy>
1139 T tgamma_delta_ratio_imp(T z, T delta, const Policy& pol)
1140 {
1141    BOOST_MATH_STD_USING
1142
1143    if(z <= 0)
1144       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);
1145    if(z+delta <= 0)
1146       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);
1147
1148    if(floor(delta) == delta)
1149    {
1150       if(floor(z) == z)
1151       {
1152          //
1153          // Both z and delta are integers, see if we can just use table lookup
1154          // of the factorials to get the result:
1155          //
1156          if((z <= max_factorial<T>::value) && (z + delta <= max_factorial<T>::value))
1157          {
1158             return unchecked_factorial<T>((unsigned)itrunc(z, pol) - 1) / unchecked_factorial<T>((unsigned)itrunc(T(z + delta), pol) - 1);
1159          }
1160       }
1161       if(fabs(delta) < 20)
1162       {
1163          //
1164          // delta is a small integer, we can use a finite product:
1165          //
1166          if(delta == 0)
1167             return 1;
1168          if(delta < 0)
1169          {
1170             z -= 1;
1171             T result = z;
1172             while(0 != (delta += 1))
1173             {
1174                z -= 1;
1175                result *= z;
1176             }
1177             return result;
1178          }
1179          else
1180          {
1181             T result = 1 / z;
1182             while(0 != (delta -= 1))
1183             {
1184                z += 1;
1185                result /= z;
1186             }
1187             return result;
1188          }
1189       }
1190    }
1191    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1192    return tgamma_delta_ratio_imp_lanczos(z, delta, pol, lanczos_type());
1193 }
1194
1195 template <class T, class Policy>
1196 T gamma_p_derivative_imp(T a, T x, const Policy& pol)
1197 {
1198    //
1199    // Usual error checks first:
1200    //
1201    if(a <= 0)
1202       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);
1203    if(x < 0)
1204       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);
1205    //
1206    // Now special cases:
1207    //
1208    if(x == 0)
1209    {
1210       return (a > 1) ? 0 :
1211          (a == 1) ? 1 : policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
1212    }
1213    //
1214    // Normal case:
1215    //
1216    typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
1217    T f1 = detail::regularised_gamma_prefix(a, x, pol, lanczos_type());
1218    if((x < 1) && (tools::max_value<T>() * x < f1))
1219    {
1220       // overflow:
1221       return policies::raise_overflow_error<T>("boost::math::gamma_p_derivative<%1%>(%1%, %1%)", 0, pol);
1222    }
1223
1224    f1 /= x;
1225
1226    return f1;
1227 }
1228
1229 template <class T, class Policy>
1230 inline typename tools::promote_args<T>::type
1231    tgamma(T z, const Policy& /* pol */, const mpl::true_)
1232 {
1233    BOOST_FPU_EXCEPTION_GUARD
1234    typedef typename tools::promote_args<T>::type result_type;
1235    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1236    typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1237    typedef typename policies::normalise<
1238       Policy,
1239       policies::promote_float<false>,
1240       policies::promote_double<false>,
1241       policies::discrete_quantile<>,
1242       policies::assert_undefined<> >::type forwarding_policy;
1243    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%)");
1244 }
1245
1246 template <class T1, class T2, class Policy>
1247 inline typename tools::promote_args<T1, T2>::type
1248    tgamma(T1 a, T2 z, const Policy&, const mpl::false_)
1249 {
1250    BOOST_FPU_EXCEPTION_GUARD
1251    typedef typename tools::promote_args<T1, T2>::type result_type;
1252    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1253    typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1254    typedef typename policies::normalise<
1255       Policy,
1256       policies::promote_float<false>,
1257       policies::promote_double<false>,
1258       policies::discrete_quantile<>,
1259       policies::assert_undefined<> >::type forwarding_policy;
1260    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
1261       detail::gamma_incomplete_imp(static_cast<value_type>(a),
1262       static_cast<value_type>(z), false, true,
1263       forwarding_policy(), static_cast<value_type*>(0)), "boost::math::tgamma<%1%>(%1%, %1%)");
1264 }
1265
1266 template <class T1, class T2>
1267 inline typename tools::promote_args<T1, T2>::type
1268    tgamma(T1 a, T2 z, const mpl::false_ tag)
1269 {
1270    return tgamma(a, z, policies::policy<>(), tag);
1271 }
1272
1273 } // namespace detail
1274
1275 template <class T>
1276 inline typename tools::promote_args<T>::type
1277    tgamma(T z)
1278 {
1279    return tgamma(z, policies::policy<>());
1280 }
1281
1282 template <class T, class Policy>
1283 inline typename tools::promote_args<T>::type
1284    lgamma(T z, int* sign, const Policy&)
1285 {
1286    BOOST_FPU_EXCEPTION_GUARD
1287    typedef typename tools::promote_args<T>::type result_type;
1288    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1289    typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1290    typedef typename policies::normalise<
1291       Policy,
1292       policies::promote_float<false>,
1293       policies::promote_double<false>,
1294       policies::discrete_quantile<>,
1295       policies::assert_undefined<> >::type forwarding_policy;
1296    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%)");
1297 }
1298
1299 template <class T>
1300 inline typename tools::promote_args<T>::type
1301    lgamma(T z, int* sign)
1302 {
1303    return lgamma(z, sign, policies::policy<>());
1304 }
1305
1306 template <class T, class Policy>
1307 inline typename tools::promote_args<T>::type
1308    lgamma(T x, const Policy& pol)
1309 {
1310    return ::boost::math::lgamma(x, 0, pol);
1311 }
1312
1313 template <class T>
1314 inline typename tools::promote_args<T>::type
1315    lgamma(T x)
1316 {
1317    return ::boost::math::lgamma(x, 0, policies::policy<>());
1318 }
1319
1320 template <class T, class Policy>
1321 inline typename tools::promote_args<T>::type
1322    tgamma1pm1(T z, const Policy& /* pol */)
1323 {
1324    BOOST_FPU_EXCEPTION_GUARD
1325    typedef typename tools::promote_args<T>::type result_type;
1326    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1327    typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1328    typedef typename policies::normalise<
1329       Policy,
1330       policies::promote_float<false>,
1331       policies::promote_double<false>,
1332       policies::discrete_quantile<>,
1333       policies::assert_undefined<> >::type forwarding_policy;
1334
1335    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%)");
1336 }
1337
1338 template <class T>
1339 inline typename tools::promote_args<T>::type
1340    tgamma1pm1(T z)
1341 {
1342    return tgamma1pm1(z, policies::policy<>());
1343 }
1344
1345 //
1346 // Full upper incomplete gamma:
1347 //
1348 template <class T1, class T2>
1349 inline typename tools::promote_args<T1, T2>::type
1350    tgamma(T1 a, T2 z)
1351 {
1352    //
1353    // Type T2 could be a policy object, or a value, select the
1354    // right overload based on T2:
1355    //
1356    typedef typename policies::is_policy<T2>::type maybe_policy;
1357    return detail::tgamma(a, z, maybe_policy());
1358 }
1359 template <class T1, class T2, class Policy>
1360 inline typename tools::promote_args<T1, T2>::type
1361    tgamma(T1 a, T2 z, const Policy& pol)
1362 {
1363    return detail::tgamma(a, z, pol, mpl::false_());
1364 }
1365 //
1366 // Full lower incomplete gamma:
1367 //
1368 template <class T1, class T2, class Policy>
1369 inline typename tools::promote_args<T1, T2>::type
1370    tgamma_lower(T1 a, T2 z, const Policy&)
1371 {
1372    BOOST_FPU_EXCEPTION_GUARD
1373    typedef typename tools::promote_args<T1, T2>::type result_type;
1374    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1375    typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1376    typedef typename policies::normalise<
1377       Policy,
1378       policies::promote_float<false>,
1379       policies::promote_double<false>,
1380       policies::discrete_quantile<>,
1381       policies::assert_undefined<> >::type forwarding_policy;
1382
1383    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
1384       detail::gamma_incomplete_imp(static_cast<value_type>(a),
1385       static_cast<value_type>(z), false, false,
1386       forwarding_policy(), static_cast<value_type*>(0)), "tgamma_lower<%1%>(%1%, %1%)");
1387 }
1388 template <class T1, class T2>
1389 inline typename tools::promote_args<T1, T2>::type
1390    tgamma_lower(T1 a, T2 z)
1391 {
1392    return tgamma_lower(a, z, policies::policy<>());
1393 }
1394 //
1395 // Regularised upper incomplete gamma:
1396 //
1397 template <class T1, class T2, class Policy>
1398 inline typename tools::promote_args<T1, T2>::type
1399    gamma_q(T1 a, T2 z, const Policy& /* pol */)
1400 {
1401    BOOST_FPU_EXCEPTION_GUARD
1402    typedef typename tools::promote_args<T1, T2>::type result_type;
1403    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1404    typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1405    typedef typename policies::normalise<
1406       Policy,
1407       policies::promote_float<false>,
1408       policies::promote_double<false>,
1409       policies::discrete_quantile<>,
1410       policies::assert_undefined<> >::type forwarding_policy;
1411
1412    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
1413       detail::gamma_incomplete_imp(static_cast<value_type>(a),
1414       static_cast<value_type>(z), true, true,
1415       forwarding_policy(), static_cast<value_type*>(0)), "gamma_q<%1%>(%1%, %1%)");
1416 }
1417 template <class T1, class T2>
1418 inline typename tools::promote_args<T1, T2>::type
1419    gamma_q(T1 a, T2 z)
1420 {
1421    return gamma_q(a, z, policies::policy<>());
1422 }
1423 //
1424 // Regularised lower incomplete gamma:
1425 //
1426 template <class T1, class T2, class Policy>
1427 inline typename tools::promote_args<T1, T2>::type
1428    gamma_p(T1 a, T2 z, const Policy&)
1429 {
1430    BOOST_FPU_EXCEPTION_GUARD
1431    typedef typename tools::promote_args<T1, T2>::type result_type;
1432    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1433    typedef typename lanczos::lanczos<value_type, Policy>::type evaluation_type;
1434    typedef typename policies::normalise<
1435       Policy,
1436       policies::promote_float<false>,
1437       policies::promote_double<false>,
1438       policies::discrete_quantile<>,
1439       policies::assert_undefined<> >::type forwarding_policy;
1440
1441    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
1442       detail::gamma_incomplete_imp(static_cast<value_type>(a),
1443       static_cast<value_type>(z), true, false,
1444       forwarding_policy(), static_cast<value_type*>(0)), "gamma_p<%1%>(%1%, %1%)");
1445 }
1446 template <class T1, class T2>
1447 inline typename tools::promote_args<T1, T2>::type
1448    gamma_p(T1 a, T2 z)
1449 {
1450    return gamma_p(a, z, policies::policy<>());
1451 }
1452
1453 // ratios of gamma functions:
1454 template <class T1, class T2, class Policy>
1455 inline typename tools::promote_args<T1, T2>::type
1456    tgamma_delta_ratio(T1 z, T2 delta, const Policy& /* pol */)
1457 {
1458    BOOST_FPU_EXCEPTION_GUARD
1459    typedef typename tools::promote_args<T1, T2>::type result_type;
1460    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1461    typedef typename policies::normalise<
1462       Policy,
1463       policies::promote_float<false>,
1464       policies::promote_double<false>,
1465       policies::discrete_quantile<>,
1466       policies::assert_undefined<> >::type forwarding_policy;
1467
1468    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%)");
1469 }
1470 template <class T1, class T2>
1471 inline typename tools::promote_args<T1, T2>::type
1472    tgamma_delta_ratio(T1 z, T2 delta)
1473 {
1474    return tgamma_delta_ratio(z, delta, policies::policy<>());
1475 }
1476 template <class T1, class T2, class Policy>
1477 inline typename tools::promote_args<T1, T2>::type
1478    tgamma_ratio(T1 a, T2 b, const Policy&)
1479 {
1480    typedef typename tools::promote_args<T1, T2>::type result_type;
1481    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1482    typedef typename policies::normalise<
1483       Policy,
1484       policies::promote_float<false>,
1485       policies::promote_double<false>,
1486       policies::discrete_quantile<>,
1487       policies::assert_undefined<> >::type forwarding_policy;
1488
1489    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%)");
1490 }
1491 template <class T1, class T2>
1492 inline typename tools::promote_args<T1, T2>::type
1493    tgamma_ratio(T1 a, T2 b)
1494 {
1495    return tgamma_ratio(a, b, policies::policy<>());
1496 }
1497
1498 template <class T1, class T2, class Policy>
1499 inline typename tools::promote_args<T1, T2>::type
1500    gamma_p_derivative(T1 a, T2 x, const Policy&)
1501 {
1502    BOOST_FPU_EXCEPTION_GUARD
1503    typedef typename tools::promote_args<T1, T2>::type result_type;
1504    typedef typename policies::evaluation<result_type, Policy>::type value_type;
1505    typedef typename policies::normalise<
1506       Policy,
1507       policies::promote_float<false>,
1508       policies::promote_double<false>,
1509       policies::discrete_quantile<>,
1510       policies::assert_undefined<> >::type forwarding_policy;
1511
1512    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%)");
1513 }
1514 template <class T1, class T2>
1515 inline typename tools::promote_args<T1, T2>::type
1516    gamma_p_derivative(T1 a, T2 x)
1517 {
1518    return gamma_p_derivative(a, x, policies::policy<>());
1519 }
1520
1521 } // namespace math
1522 } // namespace boost
1523
1524 #ifdef BOOST_MSVC
1525 # pragma warning(pop)
1526 #endif
1527
1528 #include <boost/math/special_functions/detail/igamma_inverse.hpp>
1529 #include <boost/math/special_functions/detail/gamma_inva.hpp>
1530 #include <boost/math/special_functions/erf.hpp>
1531
1532 #endif // BOOST_MATH_SF_GAMMA_HPP
1533
1534
1535
1536