igamma_large.hpp

来自「Boost provides free peer-reviewed portab」· HPP 代码 · 共 770 行 · 第 1/2 页

//  Copyright John Maddock 2006.
//  Use, modification and distribution are subject to the
//  Boost Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
//
// This file implements the asymptotic expansions of the incomplete
// gamma functions P(a, x) and Q(a, x), used when a is large and
// x ~ a.
//
// The primary reference is:
//
// "The Asymptotic Expansion of the Incomplete Gamma Functions"
// N. M. Temme.
// Siam J. Math Anal. Vol 10 No 4, July 1979, p757.
//
// A different way of evaluating these expansions,
// plus a lot of very useful background information is in:
// 
// "A Set of Algorithms For the Incomplete Gamma Functions."
// N. M. Temme.
// Probability in the Engineering and Informational Sciences,
// 8, 1994, 291.
//
// An alternative implementation is in:
//
// "Computation of the Incomplete Gamma Function Ratios and their Inverse."
// A. R. Didonato and A. H. Morris.
// ACM TOMS, Vol 12, No 4, Dec 1986, p377.
//
// There are various versions of the same code below, each accurate
// to a different precision.  To understand the code, refer to Didonato
// and Morris, from Eq 17 and 18 onwards.
//
// The coefficients used here are not taken from Didonato and Morris:
// the domain over which these expansions are used is slightly different
// to theirs, and their constants are not quite accurate enough for
// 128-bit long double's.  Instead the coefficients were calculated
// using the methods described by Temme p762 from Eq 3.8 onwards.
// The values obtained agree with those obtained by Didonato and Morris
// (at least to the first 30 digits that they provide).
// At double precision the degrees of polynomial required for full
// machine precision are close to those recomended to Didonato and Morris,
// but of course many more terms are needed for larger types.
//
#ifndef BOOST_MATH_DETAIL_IGAMMA_LARGE
#define BOOST_MATH_DETAIL_IGAMMA_LARGE

#ifdef _MSC_VER
#pragma once
#endif

namespace boost{ namespace math{ namespace detail{

// This version will never be called (at runtime), it's a stub used
// when T is unsuitable to be passed to these routines:
//
template <class T, class Policy>
inline T igamma_temme_large(T, T, const Policy& /* pol */, mpl::int_<0> const *)
{
   // stub function, should never actually be called
   BOOST_ASSERT(0);
   return 0;
}
//
// This version is accurate for up to 64-bit mantissa's, 
// (80-bit long double, or 10^-20).
//
template <class T, class Policy>
T igamma_temme_large(T a, T x, const Policy& pol, mpl::int_<64> const *)
{
   BOOST_MATH_STD_USING // ADL of std functions
   T sigma = (x - a) / a;
   T phi = -boost::math::log1pmx(sigma, pol);
   T y = a * phi;
   T z = sqrt(2 * phi);
   if(x < a)
      z = -z;

   T workspace[13];

   static const T C0[] = {
      -0.333333333333333333333L,
      0.0833333333333333333333L,
      -0.0148148148148148148148L,
      0.00115740740740740740741L,
      0.000352733686067019400353L,
      -0.0001787551440329218107L,
      0.39192631785224377817e-4L,
      -0.218544851067999216147e-5L,
      -0.18540622107151599607e-5L,
      0.829671134095308600502e-6L,
      -0.176659527368260793044e-6L,
      0.670785354340149858037e-8L,
      0.102618097842403080426e-7L,
      -0.438203601845335318655e-8L,
      0.914769958223679023418e-9L,
      -0.255141939949462497669e-10L,
      -0.583077213255042506746e-10L,
      0.243619480206674162437e-10L,
      -0.502766928011417558909e-11L,
   };
   workspace[0] = tools::evaluate_polynomial(C0, z);

   static const T C1[] = {
      -0.00185185185185185185185L,
      -0.00347222222222222222222L,
      0.00264550264550264550265L,
      -0.000990226337448559670782L,
      0.000205761316872427983539L,
      -0.40187757201646090535e-6L,
      -0.18098550334489977837e-4L,
      0.764916091608111008464e-5L,
      -0.161209008945634460038e-5L,
      0.464712780280743434226e-8L,
      0.137863344691572095931e-6L,
      -0.575254560351770496402e-7L,
      0.119516285997781473243e-7L,
      -0.175432417197476476238e-10L,
      -0.100915437106004126275e-8L,
      0.416279299184258263623e-9L,
      -0.856390702649298063807e-10L,
   };
   workspace[1] = tools::evaluate_polynomial(C1, z);

   static const T C2[] = {
      0.00413359788359788359788L,
      -0.00268132716049382716049L,
      0.000771604938271604938272L,
      0.200938786008230452675e-5L,
      -0.000107366532263651605215L,
      0.529234488291201254164e-4L,
      -0.127606351886187277134e-4L,
      0.342357873409613807419e-7L,
      0.137219573090629332056e-5L,
      -0.629899213838005502291e-6L,
      0.142806142060642417916e-6L,
      -0.204770984219908660149e-9L,
      -0.140925299108675210533e-7L,
      0.622897408492202203356e-8L,
      -0.136704883966171134993e-8L,
   };
   workspace[2] = tools::evaluate_polynomial(C2, z);

   static const T C3[] = {
      0.000649434156378600823045L,
      0.000229472093621399176955L,
      -0.000469189494395255712128L,
      0.000267720632062838852962L,
      -0.756180167188397641073e-4L,
      -0.239650511386729665193e-6L,
      0.110826541153473023615e-4L,
      -0.56749528269915965675e-5L,
      0.142309007324358839146e-5L,
      -0.278610802915281422406e-10L,
      -0.169584040919302772899e-6L,
      0.809946490538808236335e-7L,
      -0.191111684859736540607e-7L,
   };
   workspace[3] = tools::evaluate_polynomial(C3, z);

   static const T C4[] = {
      -0.000861888290916711698605L,
      0.000784039221720066627474L,
      -0.000299072480303190179733L,
      -0.146384525788434181781e-5L,
      0.664149821546512218666e-4L,
      -0.396836504717943466443e-4L,
      0.113757269706784190981e-4L,
      0.250749722623753280165e-9L,
      -0.169541495365583060147e-5L,
      0.890750753220530968883e-6L,
      -0.229293483400080487057e-6L,
   };
   workspace[4] = tools::evaluate_polynomial(C4, z);

   static const T C5[] = {
      -0.000336798553366358150309L,
      -0.697281375836585777429e-4L,
      0.000277275324495939207873L,
      -0.000199325705161888477003L,
      0.679778047793720783882e-4L,
      0.141906292064396701483e-6L,
      -0.135940481897686932785e-4L,
      0.801847025633420153972e-5L,
      -0.229148117650809517038e-5L,
   };
   workspace[5] = tools::evaluate_polynomial(C5, z);

   static const T C6[] = {
      0.000531307936463992223166L,
      -0.000592166437353693882865L,
      0.000270878209671804482771L,
      0.790235323266032787212e-6L,
      -0.815396936756196875093e-4L,
      0.561168275310624965004e-4L,
      -0.183291165828433755673e-4L,
      -0.307961345060330478256e-8L,
      0.346515536880360908674e-5L,
      -0.20291327396058603727e-5L,
      0.57887928631490037089e-6L,
   };
   workspace[6] = tools::evaluate_polynomial(C6, z);

   static const T C7[] = {
      0.000344367606892377671254L,
      0.517179090826059219337e-4L,
      -0.000334931610811422363117L,
      0.000281269515476323702274L,
      -0.000109765822446847310235L,
      -0.127410090954844853795e-6L,
      0.277444515115636441571e-4L,
      -0.182634888057113326614e-4L,
      0.578769494973505239894e-5L,
   };
   workspace[7] = tools::evaluate_polynomial(C7, z);

   static const T C8[] = {
      -0.000652623918595309418922L,
      0.000839498720672087279993L,
      -0.000438297098541721005061L,
      -0.696909145842055197137e-6L,
      0.000166448466420675478374L,
      -0.000127835176797692185853L,
      0.462995326369130429061e-4L,
   };
   workspace[8] = tools::evaluate_polynomial(C8, z);

   static const T C9[] = {
      -0.000596761290192746250124L,
      -0.720489541602001055909e-4L,
      0.000678230883766732836162L,
      -0.0006401475260262758451L,
      0.000277501076343287044992L,
   };
   workspace[9] = tools::evaluate_polynomial(C9, z);

   static const T C10[] = {
      0.00133244544948006563713L,
      -0.0019144384985654775265L,
      0.00110893691345966373396L,
   };
   workspace[10] = tools::evaluate_polynomial(C10, z);

   static const T C11[] = {
      0.00157972766073083495909L,
      0.000162516262783915816899L,
      -0.00206334210355432762645L,
      0.00213896861856890981541L,
      -0.00101085593912630031708L,
   };
   workspace[11] = tools::evaluate_polynomial(C11, z);

   static const T C12[] = {
      -0.00407251211951401664727L,
      0.00640336283380806979482L,
      -0.00404101610816766177474L,
   };
   workspace[12] = tools::evaluate_polynomial(C12, z);

   T result = tools::evaluate_polynomial(workspace, 1/a);
   result *= exp(-y) / sqrt(2 * constants::pi<T>() * a);
   if(x < a)
      result = -result;

   result += boost::math::erfc(sqrt(y), pol) / 2;

   return result;
}
//
// This one is accurate for 53-bit mantissa's
// (IEEE double precision or 10^-17).
//
template <class T, class Policy>
T igamma_temme_large(T a, T x, const Policy& pol, mpl::int_<53> const *)
{
   BOOST_MATH_STD_USING // ADL of std functions
   T sigma = (x - a) / a;
   T phi = -boost::math::log1pmx(sigma, pol);
   T y = a * phi;
   T z = sqrt(2 * phi);
   if(x < a)
      z = -z;

   T workspace[10];

   static const T C0[] = {
      static_cast<T>(-0.33333333333333333L),
      static_cast<T>(0.083333333333333333L),
      static_cast<T>(-0.014814814814814815L),
      static_cast<T>(0.0011574074074074074L),
      static_cast<T>(0.0003527336860670194L),
      static_cast<T>(-0.00017875514403292181L),
      static_cast<T>(0.39192631785224378e-4L),
      static_cast<T>(-0.21854485106799922e-5L),
      static_cast<T>(-0.185406221071516e-5L),
      static_cast<T>(0.8296711340953086e-6L),
      static_cast<T>(-0.17665952736826079e-6L),
      static_cast<T>(0.67078535434014986e-8L),
      static_cast<T>(0.10261809784240308e-7L),
      static_cast<T>(-0.43820360184533532e-8L),
      static_cast<T>(0.91476995822367902e-9L),
   };
   workspace[0] = tools::evaluate_polynomial(C0, z);

   static const T C1[] = {
      static_cast<T>(-0.0018518518518518519L),
      static_cast<T>(-0.0034722222222222222L),
      static_cast<T>(0.0026455026455026455L),
      static_cast<T>(-0.00099022633744855967L),
      static_cast<T>(0.00020576131687242798L),
      static_cast<T>(-0.40187757201646091e-6L),
      static_cast<T>(-0.18098550334489978e-4L),
      static_cast<T>(0.76491609160811101e-5L),
      static_cast<T>(-0.16120900894563446e-5L),
      static_cast<T>(0.46471278028074343e-8L),
      static_cast<T>(0.1378633446915721e-6L),
      static_cast<T>(-0.5752545603517705e-7L),
      static_cast<T>(0.11951628599778147e-7L),
   };
   workspace[1] = tools::evaluate_polynomial(C1, z);

   static const T C2[] = {
      static_cast<T>(0.0041335978835978836L),
      static_cast<T>(-0.0026813271604938272L),
      static_cast<T>(0.00077160493827160494L),
      static_cast<T>(0.20093878600823045e-5L),
      static_cast<T>(-0.00010736653226365161L),
      static_cast<T>(0.52923448829120125e-4L),
      static_cast<T>(-0.12760635188618728e-4L),
      static_cast<T>(0.34235787340961381e-7L),
      static_cast<T>(0.13721957309062933e-5L),
      static_cast<T>(-0.6298992138380055e-6L),
      static_cast<T>(0.14280614206064242e-6L),
   };
   workspace[2] = tools::evaluate_polynomial(C2, z);

   static const T C3[] = {
      static_cast<T>(0.00064943415637860082L),
      static_cast<T>(0.00022947209362139918L),
      static_cast<T>(-0.00046918949439525571L),
      static_cast<T>(0.00026772063206283885L),
      static_cast<T>(-0.75618016718839764e-4L),
      static_cast<T>(-0.23965051138672967e-6L),
      static_cast<T>(0.11082654115347302e-4L),
      static_cast<T>(-0.56749528269915966e-5L),
      static_cast<T>(0.14230900732435884e-5L),
   };
   workspace[3] = tools::evaluate_polynomial(C3, z);

   static const T C4[] = {
      static_cast<T>(-0.0008618882909167117L),
      static_cast<T>(0.00078403922172006663L),
      static_cast<T>(-0.00029907248030319018L),
      static_cast<T>(-0.14638452578843418e-5L),
      static_cast<T>(0.66414982154651222e-4L),
      static_cast<T>(-0.39683650471794347e-4L),
      static_cast<T>(0.11375726970678419e-4L),
   };
   workspace[4] = tools::evaluate_polynomial(C4, z);

   static const T C5[] = {
      static_cast<T>(-0.00033679855336635815L),
      static_cast<T>(-0.69728137583658578e-4L),
      static_cast<T>(0.00027727532449593921L),
      static_cast<T>(-0.00019932570516188848L),
      static_cast<T>(0.67977804779372078e-4L),
      static_cast<T>(0.1419062920643967e-6L),
      static_cast<T>(-0.13594048189768693e-4L),
      static_cast<T>(0.80184702563342015e-5L),
      static_cast<T>(-0.22914811765080952e-5L),
   };
   workspace[5] = tools::evaluate_polynomial(C5, z);

   static const T C6[] = {
      static_cast<T>(0.00053130793646399222L),
      static_cast<T>(-0.00059216643735369388L),
      static_cast<T>(0.00027087820967180448L),
      static_cast<T>(0.79023532326603279e-6L),
      static_cast<T>(-0.81539693675619688e-4L),
      static_cast<T>(0.56116827531062497e-4L),
      static_cast<T>(-0.18329116582843376e-4L),
   };
   workspace[6] = tools::evaluate_polynomial(C6, z);

   static const T C7[] = {