boost/math/special_functions/gamma.hpp

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