special_functions/detail/igamma_large.hpp

0001 //  Copyright John Maddock 2006.
0002 //  Copyright Matt Borland 2024.
0003 //  Use, modification and distribution are subject to the
0004 //  Boost Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006 //
0007 // This file implements the asymptotic expansions of the incomplete
0008 // gamma functions P(a, x) and Q(a, x), used when a is large and
0009 // x ~ a.
0010 //
0011 // The primary reference is:
0012 //
0013 // "The Asymptotic Expansion of the Incomplete Gamma Functions"
0014 // N. M. Temme.
0015 // Siam J. Math Anal. Vol 10 No 4, July 1979, p757.
0016 //
0017 // A different way of evaluating these expansions,
0018 // plus a lot of very useful background information is in:
0019 //
0020 // "A Set of Algorithms For the Incomplete Gamma Functions."
0021 // N. M. Temme.
0022 // Probability in the Engineering and Informational Sciences,
0023 // 8, 1994, 291.
0024 //
0025 // An alternative implementation is in:
0026 //
0027 // "Computation of the Incomplete Gamma Function Ratios and their Inverse."
0028 // A. R. Didonato and A. H. Morris.
0029 // ACM TOMS, Vol 12, No 4, Dec 1986, p377.
0030 //
0031 // There are various versions of the same code below, each accurate
0032 // to a different precision.  To understand the code, refer to Didonato
0033 // and Morris, from Eq 17 and 18 onwards.
0034 //
0035 // The coefficients used here are not taken from Didonato and Morris:
0036 // the domain over which these expansions are used is slightly different
0037 // to theirs, and their constants are not quite accurate enough for
0038 // 128-bit long double's.  Instead the coefficients were calculated
0039 // using the methods described by Temme p762 from Eq 3.8 onwards.
0040 // The values obtained agree with those obtained by Didonato and Morris
0041 // (at least to the first 30 digits that they provide).
0042 // At double precision the degrees of polynomial required for full
0043 // machine precision are close to those recommended to Didonato and Morris,
0044 // but of course many more terms are needed for larger types.
0045 //
0046 #ifndef BOOST_MATH_DETAIL_IGAMMA_LARGE
0047 #define BOOST_MATH_DETAIL_IGAMMA_LARGE
0048
0049 #ifdef _MSC_VER
0050 #pragma once
0051 #endif
0052
0053 #if defined(__GNUC__) && defined(BOOST_MATH_USE_FLOAT128)
0054 //
0055 // This is the only way we can avoid
0056 // warning: non-standard suffix on floating constant [-Wpedantic]
0057 // when building with -Wall -pedantic.  Neither __extension__
0058 // nor #pragma diagnostic ignored work :(
0059 //
0060 #pragma GCC system_header
0061 #endif
0062
0063 #include <boost/math/tools/config.hpp>
0064 #include <boost/math/tools/type_traits.hpp>
0065
0066 namespace boost{ namespace math{ namespace detail{
0067
0068 // This version will never be called (at runtime), it's a stub used
0069 // when T is unsuitable to be passed to these routines:
0070 //
0071 template <class T, class Policy>
0072 BOOST_MATH_GPU_ENABLED inline T igamma_temme_large(T, T, const Policy& /* pol */, const boost::math::integral_constant<int, 0>&)
0073 {
0074    // stub function, should never actually be called
0075    BOOST_MATH_ASSERT(0);
0076    return 0;
0077 }
0078 //
0079 // This version is accurate for up to 64-bit mantissa's,
0080 // (80-bit long double, or 10^-20).
0081 //
0082
0083 #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0084
0085 template <class T, class Policy>
0086 BOOST_MATH_GPU_ENABLED T igamma_temme_large(T a, T x, const Policy& pol, const boost::math::integral_constant<int, 64>&)
0087 {
0088    BOOST_MATH_STD_USING // ADL of std functions
0089    T sigma = (x - a) / a;
0090    T phi = -boost::math::log1pmx(sigma, pol);
0091    T y = a * phi;
0092    T z = sqrt(2 * phi);
0093    if(x < a)
0094       z = -z;
0095
0096    T workspace[13];
0097
0098    BOOST_MATH_STATIC const T C0[] = {
0099       BOOST_MATH_BIG_CONSTANT(T, 64, -0.333333333333333333333),
0100       BOOST_MATH_BIG_CONSTANT(T, 64, 0.0833333333333333333333),
0101       BOOST_MATH_BIG_CONSTANT(T, 64, -0.0148148148148148148148),
0102       BOOST_MATH_BIG_CONSTANT(T, 64, 0.00115740740740740740741),
0103       BOOST_MATH_BIG_CONSTANT(T, 64, 0.000352733686067019400353),
0104       BOOST_MATH_BIG_CONSTANT(T, 64, -0.0001787551440329218107),
0105       BOOST_MATH_BIG_CONSTANT(T, 64, 0.39192631785224377817e-4),
0106       BOOST_MATH_BIG_CONSTANT(T, 64, -0.218544851067999216147e-5),
0107       BOOST_MATH_BIG_CONSTANT(T, 64, -0.18540622107151599607e-5),
0108       BOOST_MATH_BIG_CONSTANT(T, 64, 0.829671134095308600502e-6),
0109       BOOST_MATH_BIG_CONSTANT(T, 64, -0.176659527368260793044e-6),
0110       BOOST_MATH_BIG_CONSTANT(T, 64, 0.670785354340149858037e-8),
0111       BOOST_MATH_BIG_CONSTANT(T, 64, 0.102618097842403080426e-7),
0112       BOOST_MATH_BIG_CONSTANT(T, 64, -0.438203601845335318655e-8),
0113       BOOST_MATH_BIG_CONSTANT(T, 64, 0.914769958223679023418e-9),
0114       BOOST_MATH_BIG_CONSTANT(T, 64, -0.255141939949462497669e-10),
0115       BOOST_MATH_BIG_CONSTANT(T, 64, -0.583077213255042506746e-10),
0116       BOOST_MATH_BIG_CONSTANT(T, 64, 0.243619480206674162437e-10),
0117       BOOST_MATH_BIG_CONSTANT(T, 64, -0.502766928011417558909e-11),
0118    };
0119    workspace[0] = tools::evaluate_polynomial(C0, z);
0120
0121    BOOST_MATH_STATIC const T C1[] = {
0122       BOOST_MATH_BIG_CONSTANT(T, 64, -0.00185185185185185185185),
0123       BOOST_MATH_BIG_CONSTANT(T, 64, -0.00347222222222222222222),
0124       BOOST_MATH_BIG_CONSTANT(T, 64, 0.00264550264550264550265),
0125       BOOST_MATH_BIG_CONSTANT(T, 64, -0.000990226337448559670782),
0126       BOOST_MATH_BIG_CONSTANT(T, 64, 0.000205761316872427983539),
0127       BOOST_MATH_BIG_CONSTANT(T, 64, -0.40187757201646090535e-6),
0128       BOOST_MATH_BIG_CONSTANT(T, 64, -0.18098550334489977837e-4),
0129       BOOST_MATH_BIG_CONSTANT(T, 64, 0.764916091608111008464e-5),
0130       BOOST_MATH_BIG_CONSTANT(T, 64, -0.161209008945634460038e-5),
0131       BOOST_MATH_BIG_CONSTANT(T, 64, 0.464712780280743434226e-8),
0132       BOOST_MATH_BIG_CONSTANT(T, 64, 0.137863344691572095931e-6),
0133       BOOST_MATH_BIG_CONSTANT(T, 64, -0.575254560351770496402e-7),
0134       BOOST_MATH_BIG_CONSTANT(T, 64, 0.119516285997781473243e-7),
0135       BOOST_MATH_BIG_CONSTANT(T, 64, -0.175432417197476476238e-10),
0136       BOOST_MATH_BIG_CONSTANT(T, 64, -0.100915437106004126275e-8),
0137       BOOST_MATH_BIG_CONSTANT(T, 64, 0.416279299184258263623e-9),
0138       BOOST_MATH_BIG_CONSTANT(T, 64, -0.856390702649298063807e-10),
0139    };
0140    workspace[1] = tools::evaluate_polynomial(C1, z);
0141
0142    BOOST_MATH_STATIC const T C2[] = {
0143       BOOST_MATH_BIG_CONSTANT(T, 64, 0.00413359788359788359788),
0144       BOOST_MATH_BIG_CONSTANT(T, 64, -0.00268132716049382716049),
0145       BOOST_MATH_BIG_CONSTANT(T, 64, 0.000771604938271604938272),
0146       BOOST_MATH_BIG_CONSTANT(T, 64, 0.200938786008230452675e-5),
0147       BOOST_MATH_BIG_CONSTANT(T, 64, -0.000107366532263651605215),
0148       BOOST_MATH_BIG_CONSTANT(T, 64, 0.529234488291201254164e-4),
0149       BOOST_MATH_BIG_CONSTANT(T, 64, -0.127606351886187277134e-4),
0150       BOOST_MATH_BIG_CONSTANT(T, 64, 0.342357873409613807419e-7),
0151       BOOST_MATH_BIG_CONSTANT(T, 64, 0.137219573090629332056e-5),
0152       BOOST_MATH_BIG_CONSTANT(T, 64, -0.629899213838005502291e-6),
0153       BOOST_MATH_BIG_CONSTANT(T, 64, 0.142806142060642417916e-6),
0154       BOOST_MATH_BIG_CONSTANT(T, 64, -0.204770984219908660149e-9),
0155       BOOST_MATH_BIG_CONSTANT(T, 64, -0.140925299108675210533e-7),
0156       BOOST_MATH_BIG_CONSTANT(T, 64, 0.622897408492202203356e-8),
0157       BOOST_MATH_BIG_CONSTANT(T, 64, -0.136704883966171134993e-8),
0158    };
0159    workspace[2] = tools::evaluate_polynomial(C2, z);
0160
0161    BOOST_MATH_STATIC const T C3[] = {
0162       BOOST_MATH_BIG_CONSTANT(T, 64, 0.000649434156378600823045),
0163       BOOST_MATH_BIG_CONSTANT(T, 64, 0.000229472093621399176955),
0164       BOOST_MATH_BIG_CONSTANT(T, 64, -0.000469189494395255712128),
0165       BOOST_MATH_BIG_CONSTANT(T, 64, 0.000267720632062838852962),
0166       BOOST_MATH_BIG_CONSTANT(T, 64, -0.756180167188397641073e-4),
0167       BOOST_MATH_BIG_CONSTANT(T, 64, -0.239650511386729665193e-6),
0168       BOOST_MATH_BIG_CONSTANT(T, 64, 0.110826541153473023615e-4),
0169       BOOST_MATH_BIG_CONSTANT(T, 64, -0.56749528269915965675e-5),
0170       BOOST_MATH_BIG_CONSTANT(T, 64, 0.142309007324358839146e-5),
0171       BOOST_MATH_BIG_CONSTANT(T, 64, -0.278610802915281422406e-10),
0172       BOOST_MATH_BIG_CONSTANT(T, 64, -0.169584040919302772899e-6),
0173       BOOST_MATH_BIG_CONSTANT(T, 64, 0.809946490538808236335e-7),
0174       BOOST_MATH_BIG_CONSTANT(T, 64, -0.191111684859736540607e-7),
0175    };
0176    workspace[3] = tools::evaluate_polynomial(C3, z);
0177
0178    BOOST_MATH_STATIC const T C4[] = {
0179       BOOST_MATH_BIG_CONSTANT(T, 64, -0.000861888290916711698605),
0180       BOOST_MATH_BIG_CONSTANT(T, 64, 0.000784039221720066627474),
0181       BOOST_MATH_BIG_CONSTANT(T, 64, -0.000299072480303190179733),
0182       BOOST_MATH_BIG_CONSTANT(T, 64, -0.146384525788434181781e-5),
0183       BOOST_MATH_BIG_CONSTANT(T, 64, 0.664149821546512218666e-4),
0184       BOOST_MATH_BIG_CONSTANT(T, 64, -0.396836504717943466443e-4),
0185       BOOST_MATH_BIG_CONSTANT(T, 64, 0.113757269706784190981e-4),
0186       BOOST_MATH_BIG_CONSTANT(T, 64, 0.250749722623753280165e-9),
0187       BOOST_MATH_BIG_CONSTANT(T, 64, -0.169541495365583060147e-5),
0188       BOOST_MATH_BIG_CONSTANT(T, 64, 0.890750753220530968883e-6),
0189       BOOST_MATH_BIG_CONSTANT(T, 64, -0.229293483400080487057e-6),
0190    };
0191    workspace[4] = tools::evaluate_polynomial(C4, z);
0192
0193    BOOST_MATH_STATIC const T C5[] = {
0194       BOOST_MATH_BIG_CONSTANT(T, 64, -0.000336798553366358150309),
0195       BOOST_MATH_BIG_CONSTANT(T, 64, -0.697281375836585777429e-4),
0196       BOOST_MATH_BIG_CONSTANT(T, 64, 0.000277275324495939207873),
0197       BOOST_MATH_BIG_CONSTANT(T, 64, -0.000199325705161888477003),
0198       BOOST_MATH_BIG_CONSTANT(T, 64, 0.679778047793720783882e-4),
0199       BOOST_MATH_BIG_CONSTANT(T, 64, 0.141906292064396701483e-6),
0200       BOOST_MATH_BIG_CONSTANT(T, 64, -0.135940481897686932785e-4),
0201       BOOST_MATH_BIG_CONSTANT(T, 64, 0.801847025633420153972e-5),
0202       BOOST_MATH_BIG_CONSTANT(T, 64, -0.229148117650809517038e-5),
0203    };
0204    workspace[5] = tools::evaluate_polynomial(C5, z);
0205
0206    BOOST_MATH_STATIC const T C6[] = {
0207       BOOST_MATH_BIG_CONSTANT(T, 64, 0.000531307936463992223166),
0208       BOOST_MATH_BIG_CONSTANT(T, 64, -0.000592166437353693882865),
0209       BOOST_MATH_BIG_CONSTANT(T, 64, 0.000270878209671804482771),
0210       BOOST_MATH_BIG_CONSTANT(T, 64, 0.790235323266032787212e-6),
0211       BOOST_MATH_BIG_CONSTANT(T, 64, -0.815396936756196875093e-4),
0212       BOOST_MATH_BIG_CONSTANT(T, 64, 0.561168275310624965004e-4),
0213       BOOST_MATH_BIG_CONSTANT(T, 64, -0.183291165828433755673e-4),
0214       BOOST_MATH_BIG_CONSTANT(T, 64, -0.307961345060330478256e-8),
0215       BOOST_MATH_BIG_CONSTANT(T, 64, 0.346515536880360908674e-5),
0216       BOOST_MATH_BIG_CONSTANT(T, 64, -0.20291327396058603727e-5),
0217       BOOST_MATH_BIG_CONSTANT(T, 64, 0.57887928631490037089e-6),
0218    };
0219    workspace[6] = tools::evaluate_polynomial(C6, z);
0220
0221    BOOST_MATH_STATIC const T C7[] = {
0222       BOOST_MATH_BIG_CONSTANT(T, 64, 0.000344367606892377671254),
0223       BOOST_MATH_BIG_CONSTANT(T, 64, 0.517179090826059219337e-4),
0224       BOOST_MATH_BIG_CONSTANT(T, 64, -0.000334931610811422363117),
0225       BOOST_MATH_BIG_CONSTANT(T, 64, 0.000281269515476323702274),
0226       BOOST_MATH_BIG_CONSTANT(T, 64, -0.000109765822446847310235),
0227       BOOST_MATH_BIG_CONSTANT(T, 64, -0.127410090954844853795e-6),
0228       BOOST_MATH_BIG_CONSTANT(T, 64, 0.277444515115636441571e-4),
0229       BOOST_MATH_BIG_CONSTANT(T, 64, -0.182634888057113326614e-4),
0230       BOOST_MATH_BIG_CONSTANT(T, 64, 0.578769494973505239894e-5),
0231    };
0232    workspace[7] = tools::evaluate_polynomial(C7, z);
0233
0234    BOOST_MATH_STATIC const T C8[] = {
0235       BOOST_MATH_BIG_CONSTANT(T, 64, -0.000652623918595309418922),
0236       BOOST_MATH_BIG_CONSTANT(T, 64, 0.000839498720672087279993),
0237       BOOST_MATH_BIG_CONSTANT(T, 64, -0.000438297098541721005061),
0238       BOOST_MATH_BIG_CONSTANT(T, 64, -0.696909145842055197137e-6),
0239       BOOST_MATH_BIG_CONSTANT(T, 64, 0.000166448466420675478374),
0240       BOOST_MATH_BIG_CONSTANT(T, 64, -0.000127835176797692185853),
0241       BOOST_MATH_BIG_CONSTANT(T, 64, 0.462995326369130429061e-4),
0242    };
0243    workspace[8] = tools::evaluate_polynomial(C8, z);
0244
0245    BOOST_MATH_STATIC const T C9[] = {
0246       BOOST_MATH_BIG_CONSTANT(T, 64, -0.000596761290192746250124),
0247       BOOST_MATH_BIG_CONSTANT(T, 64, -0.720489541602001055909e-4),
0248       BOOST_MATH_BIG_CONSTANT(T, 64, 0.000678230883766732836162),
0249       BOOST_MATH_BIG_CONSTANT(T, 64, -0.0006401475260262758451),
0250       BOOST_MATH_BIG_CONSTANT(T, 64, 0.000277501076343287044992),
0251    };
0252    workspace[9] = tools::evaluate_polynomial(C9, z);
0253
0254    BOOST_MATH_STATIC const T C10[] = {
0255       BOOST_MATH_BIG_CONSTANT(T, 64, 0.00133244544948006563713),
0256       BOOST_MATH_BIG_CONSTANT(T, 64, -0.0019144384985654775265),
0257       BOOST_MATH_BIG_CONSTANT(T, 64, 0.00110893691345966373396),
0258    };
0259    workspace[10] = tools::evaluate_polynomial(C10, z);
0260
0261    BOOST_MATH_STATIC const T C11[] = {
0262       BOOST_MATH_BIG_CONSTANT(T, 64, 0.00157972766073083495909),
0263       BOOST_MATH_BIG_CONSTANT(T, 64, 0.000162516262783915816899),
0264       BOOST_MATH_BIG_CONSTANT(T, 64, -0.00206334210355432762645),
0265       BOOST_MATH_BIG_CONSTANT(T, 64, 0.00213896861856890981541),
0266       BOOST_MATH_BIG_CONSTANT(T, 64, -0.00101085593912630031708),
0267    };
0268    workspace[11] = tools::evaluate_polynomial(C11, z);
0269
0270    BOOST_MATH_STATIC const T C12[] = {
0271       BOOST_MATH_BIG_CONSTANT(T, 64, -0.00407251211951401664727),
0272       BOOST_MATH_BIG_CONSTANT(T, 64, 0.00640336283380806979482),
0273       BOOST_MATH_BIG_CONSTANT(T, 64, -0.00404101610816766177474),
0274    };
0275    workspace[12] = tools::evaluate_polynomial(C12, z);
0276
0277    T result = tools::evaluate_polynomial<13, T, T>(workspace, 1/a);
0278    result *= exp(-y) / sqrt(2 * constants::pi<T>() * a);
0279    if(x < a)
0280       result = -result;
0281
0282    result += boost::math::erfc(sqrt(y), pol) / 2;
0283
0284    return result;
0285 }
0286
0287 #endif
0288
0289 //
0290 // This one is accurate for 53-bit mantissa's
0291 // (IEEE double precision or 10^-17).
0292 //
0293 template <class T, class Policy>
0294 BOOST_MATH_GPU_ENABLED T igamma_temme_large(T a, T x, const Policy& pol, const boost::math::integral_constant<int, 53>&)
0295 {
0296    BOOST_MATH_STD_USING // ADL of std functions
0297    T sigma = (x - a) / a;
0298    T phi = -boost::math::log1pmx(sigma, pol);
0299    T y = a * phi;
0300    T z = sqrt(2 * phi);
0301    if(x < a)
0302       z = -z;
0303
0304    T workspace[10];
0305
0306    BOOST_MATH_STATIC const T C0[] = {
0307       static_cast<T>(-0.33333333333333333L),
0308       static_cast<T>(0.083333333333333333L),
0309       static_cast<T>(-0.014814814814814815L),
0310       static_cast<T>(0.0011574074074074074L),
0311       static_cast<T>(0.0003527336860670194L),
0312       static_cast<T>(-0.00017875514403292181L),
0313       static_cast<T>(0.39192631785224378e-4L),
0314       static_cast<T>(-0.21854485106799922e-5L),
0315       static_cast<T>(-0.185406221071516e-5L),
0316       static_cast<T>(0.8296711340953086e-6L),
0317       static_cast<T>(-0.17665952736826079e-6L),
0318       static_cast<T>(0.67078535434014986e-8L),
0319       static_cast<T>(0.10261809784240308e-7L),
0320       static_cast<T>(-0.43820360184533532e-8L),
0321       static_cast<T>(0.91476995822367902e-9L),
0322    };
0323    workspace[0] = tools::evaluate_polynomial(C0, z);
0324
0325    BOOST_MATH_STATIC const T C1[] = {
0326       static_cast<T>(-0.0018518518518518519L),
0327       static_cast<T>(-0.0034722222222222222L),
0328       static_cast<T>(0.0026455026455026455L),
0329       static_cast<T>(-0.00099022633744855967L),
0330       static_cast<T>(0.00020576131687242798L),
0331       static_cast<T>(-0.40187757201646091e-6L),
0332       static_cast<T>(-0.18098550334489978e-4L),
0333       static_cast<T>(0.76491609160811101e-5L),
0334       static_cast<T>(-0.16120900894563446e-5L),
0335       static_cast<T>(0.46471278028074343e-8L),
0336       static_cast<T>(0.1378633446915721e-6L),
0337       static_cast<T>(-0.5752545603517705e-7L),
0338       static_cast<T>(0.11951628599778147e-7L),
0339    };
0340    workspace[1] = tools::evaluate_polynomial(C1, z);
0341
0342    BOOST_MATH_STATIC const T C2[] = {
0343       static_cast<T>(0.0041335978835978836L),
0344       static_cast<T>(-0.0026813271604938272L),
0345       static_cast<T>(0.00077160493827160494L),
0346       static_cast<T>(0.20093878600823045e-5L),
0347       static_cast<T>(-0.00010736653226365161L),
0348       static_cast<T>(0.52923448829120125e-4L),
0349       static_cast<T>(-0.12760635188618728e-4L),
0350       static_cast<T>(0.34235787340961381e-7L),
0351       static_cast<T>(0.13721957309062933e-5L),
0352       static_cast<T>(-0.6298992138380055e-6L),
0353       static_cast<T>(0.14280614206064242e-6L),
0354    };
0355    workspace[2] = tools::evaluate_polynomial(C2, z);
0356
0357    BOOST_MATH_STATIC const T C3[] = {
0358       static_cast<T>(0.00064943415637860082L),
0359       static_cast<T>(0.00022947209362139918L),
0360       static_cast<T>(-0.00046918949439525571L),
0361       static_cast<T>(0.00026772063206283885L),
0362       static_cast<T>(-0.75618016718839764e-4L),
0363       static_cast<T>(-0.23965051138672967e-6L),
0364       static_cast<T>(0.11082654115347302e-4L),
0365       static_cast<T>(-0.56749528269915966e-5L),
0366       static_cast<T>(0.14230900732435884e-5L),
0367    };
0368    workspace[3] = tools::evaluate_polynomial(C3, z);
0369
0370    BOOST_MATH_STATIC const T C4[] = {
0371       static_cast<T>(-0.0008618882909167117L),
0372       static_cast<T>(0.00078403922172006663L),
0373       static_cast<T>(-0.00029907248030319018L),
0374       static_cast<T>(-0.14638452578843418e-5L),
0375       static_cast<T>(0.66414982154651222e-4L),
0376       static_cast<T>(-0.39683650471794347e-4L),
0377       static_cast<T>(0.11375726970678419e-4L),
0378    };
0379    workspace[4] = tools::evaluate_polynomial(C4, z);
0380
0381    BOOST_MATH_STATIC const T C5[] = {
0382       static_cast<T>(-0.00033679855336635815L),
0383       static_cast<T>(-0.69728137583658578e-4L),
0384       static_cast<T>(0.00027727532449593921L),
0385       static_cast<T>(-0.00019932570516188848L),
0386       static_cast<T>(0.67977804779372078e-4L),
0387       static_cast<T>(0.1419062920643967e-6L),
0388       static_cast<T>(-0.13594048189768693e-4L),
0389       static_cast<T>(0.80184702563342015e-5L),
0390       static_cast<T>(-0.22914811765080952e-5L),
0391    };
0392    workspace[5] = tools::evaluate_polynomial(C5, z);
0393
0394    BOOST_MATH_STATIC const T C6[] = {
0395       static_cast<T>(0.00053130793646399222L),
0396       static_cast<T>(-0.00059216643735369388L),
0397       static_cast<T>(0.00027087820967180448L),
0398       static_cast<T>(0.79023532326603279e-6L),
0399       static_cast<T>(-0.81539693675619688e-4L),
0400       static_cast<T>(0.56116827531062497e-4L),
0401       static_cast<T>(-0.18329116582843376e-4L),
0402    };
0403    workspace[6] = tools::evaluate_polynomial(C6, z);
0404
0405    BOOST_MATH_STATIC const T C7[] = {
0406       static_cast<T>(0.00034436760689237767L),
0407       static_cast<T>(0.51717909082605922e-4L),
0408       static_cast<T>(-0.00033493161081142236L),
0409       static_cast<T>(0.0002812695154763237L),
0410       static_cast<T>(-0.00010976582244684731L),
0411    };
0412    workspace[7] = tools::evaluate_polynomial(C7, z);
0413
0414    BOOST_MATH_STATIC const T C8[] = {
0415       static_cast<T>(-0.00065262391859530942L),
0416       static_cast<T>(0.00083949872067208728L),
0417       static_cast<T>(-0.00043829709854172101L),
0418    };
0419    workspace[8] = tools::evaluate_polynomial(C8, z);
0420    workspace[9] = static_cast<T>(-0.00059676129019274625L);
0421
0422    T result = tools::evaluate_polynomial<10, T, T>(workspace, 1/a);
0423    result *= exp(-y) / sqrt(2 * constants::pi<T>() * a);
0424    if(x < a)
0425       result = -result;
0426
0427    #ifdef BOOST_MATH_HAS_NVRTC
0428    if (boost::math::is_same_v<T, float>)
0429    {
0430       result += ::erfcf(::sqrtf(y)) / 2;
0431    }
0432    else
0433    {
0434       result += ::erfc(::sqrt(y)) / 2;
0435    }
0436    #else
0437    result += boost::math::erfc(sqrt(y), pol) / 2;
0438    #endif
0439
0440    return result;
0441 }
0442 //
0443 // This one is accurate for 24-bit mantissa's
0444 // (IEEE float precision, or 10^-8)
0445 //
0446 template <class T, class Policy>
0447 BOOST_MATH_GPU_ENABLED T igamma_temme_large(T a, T x, const Policy& pol, const boost::math::integral_constant<int, 24>&)
0448 {
0449    BOOST_MATH_STD_USING // ADL of std functions
0450    T sigma = (x - a) / a;
0451    T phi = -boost::math::log1pmx(sigma, pol);
0452    T y = a * phi;
0453    T z = sqrt(2 * phi);
0454    if(x < a)
0455       z = -z;
0456
0457    T workspace[3];
0458
0459    BOOST_MATH_STATIC const T C0[] = {
0460       static_cast<T>(-0.333333333L),
0461       static_cast<T>(0.0833333333L),
0462       static_cast<T>(-0.0148148148L),
0463       static_cast<T>(0.00115740741L),
0464       static_cast<T>(0.000352733686L),
0465       static_cast<T>(-0.000178755144L),
0466       static_cast<T>(0.391926318e-4L),
0467    };
0468    workspace[0] = tools::evaluate_polynomial(C0, z);
0469
0470    BOOST_MATH_STATIC const T C1[] = {
0471       static_cast<T>(-0.00185185185L),
0472       static_cast<T>(-0.00347222222L),
0473       static_cast<T>(0.00264550265L),
0474       static_cast<T>(-0.000990226337L),
0475       static_cast<T>(0.000205761317L),
0476    };
0477    workspace[1] = tools::evaluate_polynomial(C1, z);
0478
0479    BOOST_MATH_STATIC const T C2[] = {
0480       static_cast<T>(0.00413359788L),
0481       static_cast<T>(-0.00268132716L),
0482       static_cast<T>(0.000771604938L),
0483    };
0484    workspace[2] = tools::evaluate_polynomial(C2, z);
0485
0486    T result = tools::evaluate_polynomial(workspace, 1/a);
0487    result *= exp(-y) / sqrt(2 * constants::pi<T>() * a);
0488    if(x < a)
0489       result = -result;
0490
0491    #ifdef BOOST_MATH_HAS_NVRTC
0492    if (boost::math::is_same_v<T, float>)
0493    {
0494       result += ::erfcf(::sqrtf(y)) / 2;
0495    }
0496    else
0497    {
0498       result += ::erfc(::sqrt(y)) / 2;
0499    }
0500    #else
0501    result += boost::math::erfc(sqrt(y), pol) / 2;
0502    #endif
0503
0504    return result;
0505 }
0506 //
0507 // And finally, a version for 113-bit mantissa's
0508 // (128-bit long doubles, or 10^-34).
0509 // Note this one has been optimised for a > 200
0510 // It's use for a < 200 is not recommended, that would
0511 // require many more terms in the polynomials.
0512 //
0513 #ifndef BOOST_MATH_HAS_GPU_SUPPORT
0514
0515 template <class T, class Policy>
0516 BOOST_MATH_GPU_ENABLED T igamma_temme_large(T a, T x, const Policy& pol, const boost::math::integral_constant<int, 113>&)
0517 {
0518    BOOST_MATH_STD_USING // ADL of std functions
0519    T sigma = (x - a) / a;
0520    T phi = -boost::math::log1pmx(sigma, pol);
0521    T y = a * phi;
0522    T z = sqrt(2 * phi);
0523    if(x < a)
0524       z = -z;
0525
0526    T workspace[14];
0527
0528    BOOST_MATH_STATIC const T C0[] = {
0529       BOOST_MATH_BIG_CONSTANT(T, 113, -0.333333333333333333333333333333333333),
0530       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0833333333333333333333333333333333333),
0531       BOOST_MATH_BIG_CONSTANT(T, 113, -0.0148148148148148148148148148148148148),
0532       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00115740740740740740740740740740740741),
0533       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0003527336860670194003527336860670194),
0534       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000178755144032921810699588477366255144),
0535       BOOST_MATH_BIG_CONSTANT(T, 113, 0.391926317852243778169704095630021556e-4),
0536       BOOST_MATH_BIG_CONSTANT(T, 113, -0.218544851067999216147364295512443661e-5),
0537       BOOST_MATH_BIG_CONSTANT(T, 113, -0.185406221071515996070179883622956325e-5),
0538       BOOST_MATH_BIG_CONSTANT(T, 113, 0.829671134095308600501624213166443227e-6),
0539       BOOST_MATH_BIG_CONSTANT(T, 113, -0.17665952736826079304360054245742403e-6),
0540       BOOST_MATH_BIG_CONSTANT(T, 113, 0.670785354340149858036939710029613572e-8),
0541       BOOST_MATH_BIG_CONSTANT(T, 113, 0.102618097842403080425739573227252951e-7),
0542       BOOST_MATH_BIG_CONSTANT(T, 113, -0.438203601845335318655297462244719123e-8),
0543       BOOST_MATH_BIG_CONSTANT(T, 113, 0.914769958223679023418248817633113681e-9),
0544       BOOST_MATH_BIG_CONSTANT(T, 113, -0.255141939949462497668779537993887013e-10),
0545       BOOST_MATH_BIG_CONSTANT(T, 113, -0.583077213255042506746408945040035798e-10),
0546       BOOST_MATH_BIG_CONSTANT(T, 113, 0.243619480206674162436940696707789943e-10),
0547       BOOST_MATH_BIG_CONSTANT(T, 113, -0.502766928011417558909054985925744366e-11),
0548       BOOST_MATH_BIG_CONSTANT(T, 113, 0.110043920319561347708374174497293411e-12),
0549       BOOST_MATH_BIG_CONSTANT(T, 113, 0.337176326240098537882769884169200185e-12),
0550       BOOST_MATH_BIG_CONSTANT(T, 113, -0.13923887224181620659193661848957998e-12),
0551       BOOST_MATH_BIG_CONSTANT(T, 113, 0.285348938070474432039669099052828299e-13),
0552       BOOST_MATH_BIG_CONSTANT(T, 113, -0.513911183424257261899064580300494205e-15),
0553       BOOST_MATH_BIG_CONSTANT(T, 113, -0.197522882943494428353962401580710912e-14),
0554       BOOST_MATH_BIG_CONSTANT(T, 113, 0.809952115670456133407115668702575255e-15),
0555       BOOST_MATH_BIG_CONSTANT(T, 113, -0.165225312163981618191514820265351162e-15),
0556       BOOST_MATH_BIG_CONSTANT(T, 113, 0.253054300974788842327061090060267385e-17),
0557       BOOST_MATH_BIG_CONSTANT(T, 113, 0.116869397385595765888230876507793475e-16),
0558       BOOST_MATH_BIG_CONSTANT(T, 113, -0.477003704982048475822167804084816597e-17),
0559       BOOST_MATH_BIG_CONSTANT(T, 113, 0.969912605905623712420709685898585354e-18),
0560    };
0561    workspace[0] = tools::evaluate_polynomial(C0, z);
0562
0563    BOOST_MATH_STATIC const T C1[] = {
0564       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00185185185185185185185185185185185185),
0565       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00347222222222222222222222222222222222),
0566       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0026455026455026455026455026455026455),
0567       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000990226337448559670781893004115226337),
0568       BOOST_MATH_BIG_CONSTANT(T, 113, 0.000205761316872427983539094650205761317),
0569       BOOST_MATH_BIG_CONSTANT(T, 113, -0.401877572016460905349794238683127572e-6),
0570       BOOST_MATH_BIG_CONSTANT(T, 113, -0.180985503344899778370285914867533523e-4),
0571       BOOST_MATH_BIG_CONSTANT(T, 113, 0.76491609160811100846374214980916921e-5),
0572       BOOST_MATH_BIG_CONSTANT(T, 113, -0.16120900894563446003775221882217767e-5),
0573       BOOST_MATH_BIG_CONSTANT(T, 113, 0.464712780280743434226135033938722401e-8),
0574       BOOST_MATH_BIG_CONSTANT(T, 113, 0.137863344691572095931187533077488877e-6),
0575       BOOST_MATH_BIG_CONSTANT(T, 113, -0.575254560351770496402194531835048307e-7),
0576       BOOST_MATH_BIG_CONSTANT(T, 113, 0.119516285997781473243076536699698169e-7),
0577       BOOST_MATH_BIG_CONSTANT(T, 113, -0.175432417197476476237547551202312502e-10),
0578       BOOST_MATH_BIG_CONSTANT(T, 113, -0.100915437106004126274577504686681675e-8),
0579       BOOST_MATH_BIG_CONSTANT(T, 113, 0.416279299184258263623372347219858628e-9),
0580       BOOST_MATH_BIG_CONSTANT(T, 113, -0.856390702649298063807431562579670208e-10),
0581       BOOST_MATH_BIG_CONSTANT(T, 113, 0.606721510160475861512701762169919581e-13),
0582       BOOST_MATH_BIG_CONSTANT(T, 113, 0.716249896481148539007961017165545733e-11),
0583       BOOST_MATH_BIG_CONSTANT(T, 113, -0.293318664377143711740636683615595403e-11),
0584       BOOST_MATH_BIG_CONSTANT(T, 113, 0.599669636568368872330374527568788909e-12),
0585       BOOST_MATH_BIG_CONSTANT(T, 113, -0.216717865273233141017100472779701734e-15),
0586       BOOST_MATH_BIG_CONSTANT(T, 113, -0.497833997236926164052815522048108548e-13),
0587       BOOST_MATH_BIG_CONSTANT(T, 113, 0.202916288237134247736694804325894226e-13),
0588       BOOST_MATH_BIG_CONSTANT(T, 113, -0.413125571381061004935108332558187111e-14),
0589       BOOST_MATH_BIG_CONSTANT(T, 113, 0.828651623988309644380188591057589316e-18),
0590       BOOST_MATH_BIG_CONSTANT(T, 113, 0.341003088693333279336339355910600992e-15),
0591       BOOST_MATH_BIG_CONSTANT(T, 113, -0.138541953028939715357034547426313703e-15),
0592       BOOST_MATH_BIG_CONSTANT(T, 113, 0.281234665322887466568860332727259483e-16),
0593    };
0594    workspace[1] = tools::evaluate_polynomial(C1, z);
0595
0596    BOOST_MATH_STATIC const T C2[] = {
0597       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0041335978835978835978835978835978836),
0598       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00268132716049382716049382716049382716),
0599       BOOST_MATH_BIG_CONSTANT(T, 113, 0.000771604938271604938271604938271604938),
0600       BOOST_MATH_BIG_CONSTANT(T, 113, 0.200938786008230452674897119341563786e-5),
0601       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000107366532263651605215391223621676297),
0602       BOOST_MATH_BIG_CONSTANT(T, 113, 0.529234488291201254164217127180090143e-4),
0603       BOOST_MATH_BIG_CONSTANT(T, 113, -0.127606351886187277133779191392360117e-4),
0604       BOOST_MATH_BIG_CONSTANT(T, 113, 0.34235787340961380741902003904747389e-7),
0605       BOOST_MATH_BIG_CONSTANT(T, 113, 0.137219573090629332055943852926020279e-5),
0606       BOOST_MATH_BIG_CONSTANT(T, 113, -0.629899213838005502290672234278391876e-6),
0607       BOOST_MATH_BIG_CONSTANT(T, 113, 0.142806142060642417915846008822771748e-6),
0608       BOOST_MATH_BIG_CONSTANT(T, 113, -0.204770984219908660149195854409200226e-9),
0609       BOOST_MATH_BIG_CONSTANT(T, 113, -0.140925299108675210532930244154315272e-7),
0610       BOOST_MATH_BIG_CONSTANT(T, 113, 0.622897408492202203356394293530327112e-8),
0611       BOOST_MATH_BIG_CONSTANT(T, 113, -0.136704883966171134992724380284402402e-8),
0612       BOOST_MATH_BIG_CONSTANT(T, 113, 0.942835615901467819547711211663208075e-12),
0613       BOOST_MATH_BIG_CONSTANT(T, 113, 0.128722524000893180595479368872770442e-9),
0614       BOOST_MATH_BIG_CONSTANT(T, 113, -0.556459561343633211465414765894951439e-10),
0615       BOOST_MATH_BIG_CONSTANT(T, 113, 0.119759355463669810035898150310311343e-10),
0616       BOOST_MATH_BIG_CONSTANT(T, 113, -0.416897822518386350403836626692480096e-14),
0617       BOOST_MATH_BIG_CONSTANT(T, 113, -0.109406404278845944099299008640802908e-11),
0618       BOOST_MATH_BIG_CONSTANT(T, 113, 0.4662239946390135746326204922464679e-12),
0619       BOOST_MATH_BIG_CONSTANT(T, 113, -0.990510576390690597844122258212382301e-13),
0620       BOOST_MATH_BIG_CONSTANT(T, 113, 0.189318767683735145056885183170630169e-16),
0621       BOOST_MATH_BIG_CONSTANT(T, 113, 0.885922187259112726176031067028740667e-14),
0622       BOOST_MATH_BIG_CONSTANT(T, 113, -0.373782039804640545306560251777191937e-14),
0623       BOOST_MATH_BIG_CONSTANT(T, 113, 0.786883363903515525774088394065960751e-15),
0624    };
0625    workspace[2] = tools::evaluate_polynomial(C2, z);
0626
0627    BOOST_MATH_STATIC const T C3[] = {
0628       BOOST_MATH_BIG_CONSTANT(T, 113, 0.000649434156378600823045267489711934156),
0629       BOOST_MATH_BIG_CONSTANT(T, 113, 0.000229472093621399176954732510288065844),
0630       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000469189494395255712128140111679206329),
0631       BOOST_MATH_BIG_CONSTANT(T, 113, 0.000267720632062838852962309752433209223),
0632       BOOST_MATH_BIG_CONSTANT(T, 113, -0.756180167188397641072538191879755666e-4),
0633       BOOST_MATH_BIG_CONSTANT(T, 113, -0.239650511386729665193314027333231723e-6),
0634       BOOST_MATH_BIG_CONSTANT(T, 113, 0.110826541153473023614770299726861227e-4),
0635       BOOST_MATH_BIG_CONSTANT(T, 113, -0.567495282699159656749963105701560205e-5),
0636       BOOST_MATH_BIG_CONSTANT(T, 113, 0.14230900732435883914551894470580433e-5),
0637       BOOST_MATH_BIG_CONSTANT(T, 113, -0.278610802915281422405802158211174452e-10),
0638       BOOST_MATH_BIG_CONSTANT(T, 113, -0.16958404091930277289864168795820267e-6),
0639       BOOST_MATH_BIG_CONSTANT(T, 113, 0.809946490538808236335278504852724081e-7),
0640       BOOST_MATH_BIG_CONSTANT(T, 113, -0.191111684859736540606728140872727635e-7),
0641       BOOST_MATH_BIG_CONSTANT(T, 113, 0.239286204398081179686413514022282056e-11),
0642       BOOST_MATH_BIG_CONSTANT(T, 113, 0.206201318154887984369925818486654549e-8),
0643       BOOST_MATH_BIG_CONSTANT(T, 113, -0.946049666185513217375417988510192814e-9),
0644       BOOST_MATH_BIG_CONSTANT(T, 113, 0.215410497757749078380130268468744512e-9),
0645       BOOST_MATH_BIG_CONSTANT(T, 113, -0.138882333681390304603424682490735291e-13),
0646       BOOST_MATH_BIG_CONSTANT(T, 113, -0.218947616819639394064123400466489455e-10),
0647       BOOST_MATH_BIG_CONSTANT(T, 113, 0.979099895117168512568262802255883368e-11),
0648       BOOST_MATH_BIG_CONSTANT(T, 113, -0.217821918801809621153859472011393244e-11),
0649       BOOST_MATH_BIG_CONSTANT(T, 113, 0.62088195734079014258166361684972205e-16),
0650       BOOST_MATH_BIG_CONSTANT(T, 113, 0.212697836327973697696702537114614471e-12),
0651       BOOST_MATH_BIG_CONSTANT(T, 113, -0.934468879151743333127396765626749473e-13),
0652       BOOST_MATH_BIG_CONSTANT(T, 113, 0.204536712267828493249215913063207436e-13),
0653    };
0654    workspace[3] = tools::evaluate_polynomial(C3, z);
0655
0656    BOOST_MATH_STATIC const T C4[] = {
0657       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000861888290916711698604702719929057378),
0658       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00078403922172006662747403488144228885),
0659       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000299072480303190179733389609932819809),
0660       BOOST_MATH_BIG_CONSTANT(T, 113, -0.146384525788434181781232535690697556e-5),
0661       BOOST_MATH_BIG_CONSTANT(T, 113, 0.664149821546512218665853782451862013e-4),
0662       BOOST_MATH_BIG_CONSTANT(T, 113, -0.396836504717943466443123507595386882e-4),
0663       BOOST_MATH_BIG_CONSTANT(T, 113, 0.113757269706784190980552042885831759e-4),
0664       BOOST_MATH_BIG_CONSTANT(T, 113, 0.250749722623753280165221942390057007e-9),
0665       BOOST_MATH_BIG_CONSTANT(T, 113, -0.169541495365583060147164356781525752e-5),
0666       BOOST_MATH_BIG_CONSTANT(T, 113, 0.890750753220530968882898422505515924e-6),
0667       BOOST_MATH_BIG_CONSTANT(T, 113, -0.229293483400080487057216364891158518e-6),
0668       BOOST_MATH_BIG_CONSTANT(T, 113, 0.295679413754404904696572852500004588e-10),
0669       BOOST_MATH_BIG_CONSTANT(T, 113, 0.288658297427087836297341274604184504e-7),
0670       BOOST_MATH_BIG_CONSTANT(T, 113, -0.141897394378032193894774303903982717e-7),
0671       BOOST_MATH_BIG_CONSTANT(T, 113, 0.344635804994648970659527720474194356e-8),
0672       BOOST_MATH_BIG_CONSTANT(T, 113, -0.230245171745280671320192735850147087e-12),
0673       BOOST_MATH_BIG_CONSTANT(T, 113, -0.394092330280464052750697640085291799e-9),
0674       BOOST_MATH_BIG_CONSTANT(T, 113, 0.186023389685045019134258533045185639e-9),
0675       BOOST_MATH_BIG_CONSTANT(T, 113, -0.435632300505661804380678327446262424e-10),
0676       BOOST_MATH_BIG_CONSTANT(T, 113, 0.127860010162962312660550463349930726e-14),
0677       BOOST_MATH_BIG_CONSTANT(T, 113, 0.467927502665791946200382739991760062e-11),
0678       BOOST_MATH_BIG_CONSTANT(T, 113, -0.214924647061348285410535341910721086e-11),
0679       BOOST_MATH_BIG_CONSTANT(T, 113, 0.490881561480965216323649688463984082e-12),
0680    };
0681    workspace[4] = tools::evaluate_polynomial(C4, z);
0682
0683    BOOST_MATH_STATIC const T C5[] = {
0684       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000336798553366358150308767592718210002),
0685       BOOST_MATH_BIG_CONSTANT(T, 113, -0.697281375836585777429398828575783308e-4),
0686       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00027727532449593920787336425196507501),
0687       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000199325705161888477003360405280844238),
0688       BOOST_MATH_BIG_CONSTANT(T, 113, 0.679778047793720783881640176604435742e-4),
0689       BOOST_MATH_BIG_CONSTANT(T, 113, 0.141906292064396701483392727105575757e-6),
0690       BOOST_MATH_BIG_CONSTANT(T, 113, -0.135940481897686932784583938837504469e-4),
0691       BOOST_MATH_BIG_CONSTANT(T, 113, 0.80184702563342015397192571980419684e-5),
0692       BOOST_MATH_BIG_CONSTANT(T, 113, -0.229148117650809517038048790128781806e-5),
0693       BOOST_MATH_BIG_CONSTANT(T, 113, -0.325247355129845395166230137750005047e-9),
0694       BOOST_MATH_BIG_CONSTANT(T, 113, 0.346528464910852649559195496827579815e-6),
0695       BOOST_MATH_BIG_CONSTANT(T, 113, -0.184471871911713432765322367374920978e-6),
0696       BOOST_MATH_BIG_CONSTANT(T, 113, 0.482409670378941807563762631738989002e-7),
0697       BOOST_MATH_BIG_CONSTANT(T, 113, -0.179894667217435153025754291716644314e-13),
0698       BOOST_MATH_BIG_CONSTANT(T, 113, -0.630619450001352343517516981425944698e-8),
0699       BOOST_MATH_BIG_CONSTANT(T, 113, 0.316241762877456793773762181540969623e-8),
0700       BOOST_MATH_BIG_CONSTANT(T, 113, -0.784092425369742929000839303523267545e-9),
0701    };
0702    workspace[5] = tools::evaluate_polynomial(C5, z);
0703
0704    BOOST_MATH_STATIC const T C6[] = {
0705       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00053130793646399222316574854297762391),
0706       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000592166437353693882864836225604401187),
0707       BOOST_MATH_BIG_CONSTANT(T, 113, 0.000270878209671804482771279183488328692),
0708       BOOST_MATH_BIG_CONSTANT(T, 113, 0.790235323266032787212032944390816666e-6),
0709       BOOST_MATH_BIG_CONSTANT(T, 113, -0.815396936756196875092890088464682624e-4),
0710       BOOST_MATH_BIG_CONSTANT(T, 113, 0.561168275310624965003775619041471695e-4),
0711       BOOST_MATH_BIG_CONSTANT(T, 113, -0.183291165828433755673259749374098313e-4),
0712       BOOST_MATH_BIG_CONSTANT(T, 113, -0.307961345060330478256414192546677006e-8),
0713       BOOST_MATH_BIG_CONSTANT(T, 113, 0.346515536880360908673728529745376913e-5),
0714       BOOST_MATH_BIG_CONSTANT(T, 113, -0.202913273960586037269527254582695285e-5),
0715       BOOST_MATH_BIG_CONSTANT(T, 113, 0.578879286314900370889997586203187687e-6),
0716       BOOST_MATH_BIG_CONSTANT(T, 113, 0.233863067382665698933480579231637609e-12),
0717       BOOST_MATH_BIG_CONSTANT(T, 113, -0.88286007463304835250508524317926246e-7),
0718       BOOST_MATH_BIG_CONSTANT(T, 113, 0.474359588804081278032150770595852426e-7),
0719       BOOST_MATH_BIG_CONSTANT(T, 113, -0.125454150207103824457130611214783073e-7),
0720    };
0721    workspace[6] = tools::evaluate_polynomial(C6, z);
0722
0723    BOOST_MATH_STATIC const T C7[] = {
0724       BOOST_MATH_BIG_CONSTANT(T, 113, 0.000344367606892377671254279625108523655),
0725       BOOST_MATH_BIG_CONSTANT(T, 113, 0.517179090826059219337057843002058823e-4),
0726       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000334931610811422363116635090580012327),
0727       BOOST_MATH_BIG_CONSTANT(T, 113, 0.000281269515476323702273722110707777978),
0728       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000109765822446847310235396824500789005),
0729       BOOST_MATH_BIG_CONSTANT(T, 113, -0.127410090954844853794579954588107623e-6),
0730       BOOST_MATH_BIG_CONSTANT(T, 113, 0.277444515115636441570715073933712622e-4),
0731       BOOST_MATH_BIG_CONSTANT(T, 113, -0.182634888057113326614324442681892723e-4),
0732       BOOST_MATH_BIG_CONSTANT(T, 113, 0.578769494973505239894178121070843383e-5),
0733       BOOST_MATH_BIG_CONSTANT(T, 113, 0.493875893393627039981813418398565502e-9),
0734       BOOST_MATH_BIG_CONSTANT(T, 113, -0.105953670140260427338098566209633945e-5),
0735       BOOST_MATH_BIG_CONSTANT(T, 113, 0.616671437611040747858836254004890765e-6),
0736       BOOST_MATH_BIG_CONSTANT(T, 113, -0.175629733590604619378669693914265388e-6),
0737    };
0738    workspace[7] = tools::evaluate_polynomial(C7, z);
0739
0740    BOOST_MATH_STATIC const T C8[] = {
0741       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000652623918595309418922034919726622692),
0742       BOOST_MATH_BIG_CONSTANT(T, 113, 0.000839498720672087279993357516764983445),
0743       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000438297098541721005061087953050560377),
0744       BOOST_MATH_BIG_CONSTANT(T, 113, -0.696909145842055197136911097362072702e-6),
0745       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00016644846642067547837384572662326101),
0746       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000127835176797692185853344001461664247),
0747       BOOST_MATH_BIG_CONSTANT(T, 113, 0.462995326369130429061361032704489636e-4),
0748       BOOST_MATH_BIG_CONSTANT(T, 113, 0.455790986792270771162749294232219616e-8),
0749       BOOST_MATH_BIG_CONSTANT(T, 113, -0.105952711258051954718238500312872328e-4),
0750       BOOST_MATH_BIG_CONSTANT(T, 113, 0.678334290486516662273073740749269432e-5),
0751       BOOST_MATH_BIG_CONSTANT(T, 113, -0.210754766662588042469972680229376445e-5),
0752    };
0753    workspace[8] = tools::evaluate_polynomial(C8, z);
0754
0755    BOOST_MATH_STATIC const T C9[] = {
0756       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000596761290192746250124390067179459605),
0757       BOOST_MATH_BIG_CONSTANT(T, 113, -0.720489541602001055908571930225015052e-4),
0758       BOOST_MATH_BIG_CONSTANT(T, 113, 0.000678230883766732836161951166000673426),
0759       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000640147526026275845100045652582354779),
0760       BOOST_MATH_BIG_CONSTANT(T, 113, 0.000277501076343287044992374518205845463),
0761       BOOST_MATH_BIG_CONSTANT(T, 113, 0.181970083804651510461686554030325202e-6),
0762       BOOST_MATH_BIG_CONSTANT(T, 113, -0.847950711706850318239732559632810086e-4),
0763       BOOST_MATH_BIG_CONSTANT(T, 113, 0.610519208250153101764709122740859458e-4),
0764       BOOST_MATH_BIG_CONSTANT(T, 113, -0.210739201834048624082975255893773306e-4),
0765    };
0766    workspace[9] = tools::evaluate_polynomial(C9, z);
0767
0768    BOOST_MATH_STATIC const T C10[] = {
0769       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00133244544948006563712694993432717968),
0770       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00191443849856547752650089885832852254),
0771       BOOST_MATH_BIG_CONSTANT(T, 113, 0.0011089369134596637339607446329267522),
0772       BOOST_MATH_BIG_CONSTANT(T, 113, 0.993240412264229896742295262075817566e-6),
0773       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000508745012930931989848393025305956774),
0774       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00042735056665392884328432271160040444),
0775       BOOST_MATH_BIG_CONSTANT(T, 113, -0.000168588537679107988033552814662382059),
0776    };
0777    workspace[10] = tools::evaluate_polynomial(C10, z);
0778
0779    BOOST_MATH_STATIC const T C11[] = {
0780       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00157972766073083495908785631307733022),
0781       BOOST_MATH_BIG_CONSTANT(T, 113, 0.000162516262783915816898635123980270998),
0782       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00206334210355432762645284467690276817),
0783       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00213896861856890981541061922797693947),
0784       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00101085593912630031708085801712479376),
0785    };
0786    workspace[11] = tools::evaluate_polynomial(C11, z);
0787
0788    BOOST_MATH_STATIC const T C12[] = {
0789       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00407251211951401664727281097914544601),
0790       BOOST_MATH_BIG_CONSTANT(T, 113, 0.00640336283380806979482363809026579583),
0791       BOOST_MATH_BIG_CONSTANT(T, 113, -0.00404101610816766177473974858518094879),
0792    };
0793    workspace[12] = tools::evaluate_polynomial(C12, z);
0794    workspace[13] = -0.0059475779383993002845382844736066323L;
0795
0796    T result = tools::evaluate_polynomial(workspace, T(1/a));
0797    result *= exp(-y) / sqrt(2 * constants::pi<T>() * a);
0798    if(x < a)
0799       result = -result;
0800
0801    result += boost::math::erfc(sqrt(y), pol) / 2;
0802
0803    return result;
0804 }
0805
0806 #endif
0807
0808 }  // namespace detail
0809 }  // namespace math
0810 }  // namespace math
0811
0812
0813 #endif // BOOST_MATH_DETAIL_IGAMMA_LARGE
0814