boost/math/special_functions/detail/igamma_large.hpp

   1 //  Copyright John Maddock 2006.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5 //
   6 // This file implements the asymptotic expansions of the incomplete
   7 // gamma functions P(a, x) and Q(a, x), used when a is large and
   8 // x ~ a.
   9 //
  10 // The primary reference is:
  11 //
  12 // "The Asymptotic Expansion of the Incomplete Gamma Functions"
  13 // N. M. Temme.
  14 // Siam J. Math Anal. Vol 10 No 4, July 1979, p757.
  15 //
  16 // A different way of evaluating these expansions,
  17 // plus a lot of very useful background information is in:
  18 //
  19 // "A Set of Algorithms For the Incomplete Gamma Functions."
  20 // N. M. Temme.
  21 // Probability in the Engineering and Informational Sciences,
  22 // 8, 1994, 291.
  23 //
  24 // An alternative implementation is in:
  25 //
  26 // "Computation of the Incomplete Gamma Function Ratios and their Inverse."
  27 // A. R. Didonato and A. H. Morris.
  28 // ACM TOMS, Vol 12, No 4, Dec 1986, p377.
  29 //
  30 // There are various versions of the same code below, each accurate
  31 // to a different precision.  To understand the code, refer to Didonato
  32 // and Morris, from Eq 17 and 18 onwards.
  33 //
  34 // The coefficients used here are not taken from Didonato and Morris:
  35 // the domain over which these expansions are used is slightly different
  36 // to theirs, and their constants are not quite accurate enough for
  37 // 128-bit long double's.  Instead the coefficients were calculated
  38 // using the methods described by Temme p762 from Eq 3.8 onwards.
  39 // The values obtained agree with those obtained by Didonato and Morris
  40 // (at least to the first 30 digits that they provide).
  41 // At double precision the degrees of polynomial required for full
  42 // machine precision are close to those recomended to Didonato and Morris,
  43 // but of course many more terms are needed for larger types.
  44 //
  45 #ifndef BOOST_MATH_DETAIL_IGAMMA_LARGE
  46 #define BOOST_MATH_DETAIL_IGAMMA_LARGE
  47
  48 #ifdef _MSC_VER
  49 #pragma once
  50 #endif
  51
  52 namespace boost{ namespace math{ namespace detail{
  53
  54 // This version will never be called (at runtime), it's a stub used
  55 // when T is unsuitable to be passed to these routines:
  56 //
  57 template <class T, class Policy>
  58 inline T igamma_temme_large(T, T, const Policy& /* pol */, mpl::int_<0> const *)
  59 {
  60    // stub function, should never actually be called
  61    BOOST_ASSERT(0);
  62    return 0;
  63 }
  64 //
  65 // This version is accurate for up to 64-bit mantissa's,
  66 // (80-bit long double, or 10^-20).
  67 //
  68 template <class T, class Policy>
  69 T igamma_temme_large(T a, T x, const Policy& pol, mpl::int_<64> const *)
  70 {
  71    BOOST_MATH_STD_USING // ADL of std functions
  72    T sigma = (x - a) / a;
  73    T phi = -boost::math::log1pmx(sigma, pol);
  74    T y = a * phi;
  75    T z = sqrt(2 * phi);
  76    if(x < a)
  77       z = -z;
  78
  79    T workspace[13];
  80
  81    static const T C0[] = {
  82       -0.333333333333333333333L,
  83       0.0833333333333333333333L,
  84       -0.0148148148148148148148L,
  85       0.00115740740740740740741L,
  86       0.000352733686067019400353L,
  87       -0.0001787551440329218107L,
  88       0.39192631785224377817e-4L,
  89       -0.218544851067999216147e-5L,
  90       -0.18540622107151599607e-5L,
  91       0.829671134095308600502e-6L,
  92       -0.176659527368260793044e-6L,
  93       0.670785354340149858037e-8L,
  94       0.102618097842403080426e-7L,
  95       -0.438203601845335318655e-8L,
  96       0.914769958223679023418e-9L,
  97       -0.255141939949462497669e-10L,
  98       -0.583077213255042506746e-10L,
  99       0.243619480206674162437e-10L,
 100       -0.502766928011417558909e-11L,
 101    };
 102    workspace[0] = tools::evaluate_polynomial(C0, z);
 103
 104    static const T C1[] = {
 105       -0.00185185185185185185185L,
 106       -0.00347222222222222222222L,
 107       0.00264550264550264550265L,
 108       -0.000990226337448559670782L,
 109       0.000205761316872427983539L,
 110       -0.40187757201646090535e-6L,
 111       -0.18098550334489977837e-4L,
 112       0.764916091608111008464e-5L,
 113       -0.161209008945634460038e-5L,
 114       0.464712780280743434226e-8L,
 115       0.137863344691572095931e-6L,
 116       -0.575254560351770496402e-7L,
 117       0.119516285997781473243e-7L,
 118       -0.175432417197476476238e-10L,
 119       -0.100915437106004126275e-8L,
 120       0.416279299184258263623e-9L,
 121       -0.856390702649298063807e-10L,
 122    };
 123    workspace[1] = tools::evaluate_polynomial(C1, z);
 124
 125    static const T C2[] = {
 126       0.00413359788359788359788L,
 127       -0.00268132716049382716049L,
 128       0.000771604938271604938272L,
 129       0.200938786008230452675e-5L,
 130       -0.000107366532263651605215L,
 131       0.529234488291201254164e-4L,
 132       -0.127606351886187277134e-4L,
 133       0.342357873409613807419e-7L,
 134       0.137219573090629332056e-5L,
 135       -0.629899213838005502291e-6L,
 136       0.142806142060642417916e-6L,
 137       -0.204770984219908660149e-9L,
 138       -0.140925299108675210533e-7L,
 139       0.622897408492202203356e-8L,
 140       -0.136704883966171134993e-8L,
 141    };
 142    workspace[2] = tools::evaluate_polynomial(C2, z);
 143
 144    static const T C3[] = {
 145       0.000649434156378600823045L,
 146       0.000229472093621399176955L,
 147       -0.000469189494395255712128L,
 148       0.000267720632062838852962L,
 149       -0.756180167188397641073e-4L,
 150       -0.239650511386729665193e-6L,
 151       0.110826541153473023615e-4L,
 152       -0.56749528269915965675e-5L,
 153       0.142309007324358839146e-5L,
 154       -0.278610802915281422406e-10L,
 155       -0.169584040919302772899e-6L,
 156       0.809946490538808236335e-7L,
 157       -0.191111684859736540607e-7L,
 158    };
 159    workspace[3] = tools::evaluate_polynomial(C3, z);
 160
 161    static const T C4[] = {
 162       -0.000861888290916711698605L,
 163       0.000784039221720066627474L,
 164       -0.000299072480303190179733L,
 165       -0.146384525788434181781e-5L,
 166       0.664149821546512218666e-4L,
 167       -0.396836504717943466443e-4L,
 168       0.113757269706784190981e-4L,
 169       0.250749722623753280165e-9L,
 170       -0.169541495365583060147e-5L,
 171       0.890750753220530968883e-6L,
 172       -0.229293483400080487057e-6L,
 173    };
 174    workspace[4] = tools::evaluate_polynomial(C4, z);
 175
 176    static const T C5[] = {
 177       -0.000336798553366358150309L,
 178       -0.697281375836585777429e-4L,
 179       0.000277275324495939207873L,
 180       -0.000199325705161888477003L,
 181       0.679778047793720783882e-4L,
 182       0.141906292064396701483e-6L,
 183       -0.135940481897686932785e-4L,
 184       0.801847025633420153972e-5L,
 185       -0.229148117650809517038e-5L,
 186    };
 187    workspace[5] = tools::evaluate_polynomial(C5, z);
 188
 189    static const T C6[] = {
 190       0.000531307936463992223166L,
 191       -0.000592166437353693882865L,
 192       0.000270878209671804482771L,
 193       0.790235323266032787212e-6L,
 194       -0.815396936756196875093e-4L,
 195       0.561168275310624965004e-4L,
 196       -0.183291165828433755673e-4L,
 197       -0.307961345060330478256e-8L,
 198       0.346515536880360908674e-5L,
 199       -0.20291327396058603727e-5L,
 200       0.57887928631490037089e-6L,
 201    };
 202    workspace[6] = tools::evaluate_polynomial(C6, z);
 203
 204    static const T C7[] = {
 205       0.000344367606892377671254L,
 206       0.517179090826059219337e-4L,
 207       -0.000334931610811422363117L,
 208       0.000281269515476323702274L,
 209       -0.000109765822446847310235L,
 210       -0.127410090954844853795e-6L,
 211       0.277444515115636441571e-4L,
 212       -0.182634888057113326614e-4L,
 213       0.578769494973505239894e-5L,
 214    };
 215    workspace[7] = tools::evaluate_polynomial(C7, z);
 216
 217    static const T C8[] = {
 218       -0.000652623918595309418922L,
 219       0.000839498720672087279993L,
 220       -0.000438297098541721005061L,
 221       -0.696909145842055197137e-6L,
 222       0.000166448466420675478374L,
 223       -0.000127835176797692185853L,
 224       0.462995326369130429061e-4L,
 225    };
 226    workspace[8] = tools::evaluate_polynomial(C8, z);
 227
 228    static const T C9[] = {
 229       -0.000596761290192746250124L,
 230       -0.720489541602001055909e-4L,
 231       0.000678230883766732836162L,
 232       -0.0006401475260262758451L,
 233       0.000277501076343287044992L,
 234    };
 235    workspace[9] = tools::evaluate_polynomial(C9, z);
 236
 237    static const T C10[] = {
 238       0.00133244544948006563713L,
 239       -0.0019144384985654775265L,
 240       0.00110893691345966373396L,
 241    };
 242    workspace[10] = tools::evaluate_polynomial(C10, z);
 243
 244    static const T C11[] = {
 245       0.00157972766073083495909L,
 246       0.000162516262783915816899L,
 247       -0.00206334210355432762645L,
 248       0.00213896861856890981541L,
 249       -0.00101085593912630031708L,
 250    };
 251    workspace[11] = tools::evaluate_polynomial(C11, z);
 252
 253    static const T C12[] = {
 254       -0.00407251211951401664727L,
 255       0.00640336283380806979482L,
 256       -0.00404101610816766177474L,
 257    };
 258    workspace[12] = tools::evaluate_polynomial(C12, z);
 259
 260    T result = tools::evaluate_polynomial(workspace, 1/a);
 261    result *= exp(-y) / sqrt(2 * constants::pi<T>() * a);
 262    if(x < a)
 263       result = -result;
 264
 265    result += boost::math::erfc(sqrt(y), pol) / 2;
 266
 267    return result;
 268 }
 269 //
 270 // This one is accurate for 53-bit mantissa's
 271 // (IEEE double precision or 10^-17).
 272 //
 273 template <class T, class Policy>
 274 T igamma_temme_large(T a, T x, const Policy& pol, mpl::int_<53> const *)
 275 {
 276    BOOST_MATH_STD_USING // ADL of std functions
 277    T sigma = (x - a) / a;
 278    T phi = -boost::math::log1pmx(sigma, pol);
 279    T y = a * phi;
 280    T z = sqrt(2 * phi);
 281    if(x < a)
 282       z = -z;
 283
 284    T workspace[10];
 285
 286    static const T C0[] = {
 287       static_cast<T>(-0.33333333333333333L),
 288       static_cast<T>(0.083333333333333333L),
 289       static_cast<T>(-0.014814814814814815L),
 290       static_cast<T>(0.0011574074074074074L),
 291       static_cast<T>(0.0003527336860670194L),
 292       static_cast<T>(-0.00017875514403292181L),
 293       static_cast<T>(0.39192631785224378e-4L),
 294       static_cast<T>(-0.21854485106799922e-5L),
 295       static_cast<T>(-0.185406221071516e-5L),
 296       static_cast<T>(0.8296711340953086e-6L),
 297       static_cast<T>(-0.17665952736826079e-6L),
 298       static_cast<T>(0.67078535434014986e-8L),
 299       static_cast<T>(0.10261809784240308e-7L),
 300       static_cast<T>(-0.43820360184533532e-8L),
 301       static_cast<T>(0.91476995822367902e-9L),
 302    };
 303    workspace[0] = tools::evaluate_polynomial(C0, z);
 304
 305    static const T C1[] = {
 306       static_cast<T>(-0.0018518518518518519L),
 307       static_cast<T>(-0.0034722222222222222L),
 308       static_cast<T>(0.0026455026455026455L),
 309       static_cast<T>(-0.00099022633744855967L),
 310       static_cast<T>(0.00020576131687242798L),
 311       static_cast<T>(-0.40187757201646091e-6L),
 312       static_cast<T>(-0.18098550334489978e-4L),
 313       static_cast<T>(0.76491609160811101e-5L),
 314       static_cast<T>(-0.16120900894563446e-5L),
 315       static_cast<T>(0.46471278028074343e-8L),
 316       static_cast<T>(0.1378633446915721e-6L),
 317       static_cast<T>(-0.5752545603517705e-7L),
 318       static_cast<T>(0.11951628599778147e-7L),
 319    };
 320    workspace[1] = tools::evaluate_polynomial(C1, z);
 321
 322    static const T C2[] = {
 323       static_cast<T>(0.0041335978835978836L),
 324       static_cast<T>(-0.0026813271604938272L),
 325       static_cast<T>(0.00077160493827160494L),
 326       static_cast<T>(0.20093878600823045e-5L),
 327       static_cast<T>(-0.00010736653226365161L),
 328       static_cast<T>(0.52923448829120125e-4L),
 329       static_cast<T>(-0.12760635188618728e-4L),
 330       static_cast<T>(0.34235787340961381e-7L),
 331       static_cast<T>(0.13721957309062933e-5L),
 332       static_cast<T>(-0.6298992138380055e-6L),
 333       static_cast<T>(0.14280614206064242e-6L),
 334    };
 335    workspace[2] = tools::evaluate_polynomial(C2, z);
 336
 337    static const T C3[] = {
 338       static_cast<T>(0.00064943415637860082L),
 339       static_cast<T>(0.00022947209362139918L),
 340       static_cast<T>(-0.00046918949439525571L),
 341       static_cast<T>(0.00026772063206283885L),
 342       static_cast<T>(-0.75618016718839764e-4L),
 343       static_cast<T>(-0.23965051138672967e-6L),
 344       static_cast<T>(0.11082654115347302e-4L),
 345       static_cast<T>(-0.56749528269915966e-5L),
 346       static_cast<T>(0.14230900732435884e-5L),
 347    };
 348    workspace[3] = tools::evaluate_polynomial(C3, z);
 349
 350    static const T C4[] = {
 351       static_cast<T>(-0.0008618882909167117L),
 352       static_cast<T>(0.00078403922172006663L),
 353       static_cast<T>(-0.00029907248030319018L),
 354       static_cast<T>(-0.14638452578843418e-5L),
 355       static_cast<T>(0.66414982154651222e-4L),
 356       static_cast<T>(-0.39683650471794347e-4L),
 357       static_cast<T>(0.11375726970678419e-4L),
 358    };
 359    workspace[4] = tools::evaluate_polynomial(C4, z);
 360
 361    static const T C5[] = {
 362       static_cast<T>(-0.00033679855336635815L),
 363       static_cast<T>(-0.69728137583658578e-4L),
 364       static_cast<T>(0.00027727532449593921L),
 365       static_cast<T>(-0.00019932570516188848L),
 366       static_cast<T>(0.67977804779372078e-4L),
 367       static_cast<T>(0.1419062920643967e-6L),
 368       static_cast<T>(-0.13594048189768693e-4L),
 369       static_cast<T>(0.80184702563342015e-5L),
 370       static_cast<T>(-0.22914811765080952e-5L),
 371    };
 372    workspace[5] = tools::evaluate_polynomial(C5, z);
 373
 374    static const T C6[] = {
 375       static_cast<T>(0.00053130793646399222L),
 376       static_cast<T>(-0.00059216643735369388L),
 377       static_cast<T>(0.00027087820967180448L),
 378       static_cast<T>(0.79023532326603279e-6L),
 379       static_cast<T>(-0.81539693675619688e-4L),
 380       static_cast<T>(0.56116827531062497e-4L),
 381       static_cast<T>(-0.18329116582843376e-4L),
 382    };
 383    workspace[6] = tools::evaluate_polynomial(C6, z);
 384
 385    static const T C7[] = {
 386       static_cast<T>(0.00034436760689237767L),
 387       static_cast<T>(0.51717909082605922e-4L),
 388       static_cast<T>(-0.00033493161081142236L),
 389       static_cast<T>(0.0002812695154763237L),
 390       static_cast<T>(-0.00010976582244684731L),
 391    };
 392    workspace[7] = tools::evaluate_polynomial(C7, z);
 393
 394    static const T C8[] = {
 395       static_cast<T>(-0.00065262391859530942L),
 396       static_cast<T>(0.00083949872067208728L),
 397       static_cast<T>(-0.00043829709854172101L),
 398    };
 399    workspace[8] = tools::evaluate_polynomial(C8, z);
 400    workspace[9] = static_cast<T>(-0.00059676129019274625L);
 401
 402    T result = tools::evaluate_polynomial(workspace, 1/a);
 403    result *= exp(-y) / sqrt(2 * constants::pi<T>() * a);
 404    if(x < a)
 405       result = -result;
 406
 407    result += boost::math::erfc(sqrt(y), pol) / 2;
 408
 409    return result;
 410 }
 411 //
 412 // This one is accurate for 24-bit mantissa's
 413 // (IEEE float precision, or 10^-8)
 414 //
 415 template <class T, class Policy>
 416 T igamma_temme_large(T a, T x, const Policy& pol, mpl::int_<24> const *)
 417 {
 418    BOOST_MATH_STD_USING // ADL of std functions
 419    T sigma = (x - a) / a;
 420    T phi = -boost::math::log1pmx(sigma, pol);
 421    T y = a * phi;
 422    T z = sqrt(2 * phi);
 423    if(x < a)
 424       z = -z;
 425
 426    T workspace[3];
 427
 428    static const T C0[] = {
 429       static_cast<T>(-0.333333333L),
 430       static_cast<T>(0.0833333333L),
 431       static_cast<T>(-0.0148148148L),
 432       static_cast<T>(0.00115740741L),
 433       static_cast<T>(0.000352733686L),
 434       static_cast<T>(-0.000178755144L),
 435       static_cast<T>(0.391926318e-4L),
 436    };
 437    workspace[0] = tools::evaluate_polynomial(C0, z);
 438
 439    static const T C1[] = {
 440       static_cast<T>(-0.00185185185L),
 441       static_cast<T>(-0.00347222222L),
 442       static_cast<T>(0.00264550265L),
 443       static_cast<T>(-0.000990226337L),
 444       static_cast<T>(0.000205761317L),
 445    };
 446    workspace[1] = tools::evaluate_polynomial(C1, z);
 447
 448    static const T C2[] = {
 449       static_cast<T>(0.00413359788L),
 450       static_cast<T>(-0.00268132716L),
 451       static_cast<T>(0.000771604938L),
 452    };
 453    workspace[2] = tools::evaluate_polynomial(C2, z);
 454
 455    T result = tools::evaluate_polynomial(workspace, 1/a);
 456    result *= exp(-y) / sqrt(2 * constants::pi<T>() * a);
 457    if(x < a)
 458       result = -result;
 459
 460    result += boost::math::erfc(sqrt(y), pol) / 2;
 461
 462    return result;
 463 }
 464 //
 465 // And finally, a version for 113-bit mantissa's
 466 // (128-bit long doubles, or 10^-34).
 467 // Note this one has been optimised for a > 200
 468 // It's use for a < 200 is not recomended, that would
 469 // require many more terms in the polynomials.
 470 //
 471 template <class T, class Policy>
 472 T igamma_temme_large(T a, T x, const Policy& pol, mpl::int_<113> const *)
 473 {
 474    BOOST_MATH_STD_USING // ADL of std functions
 475    T sigma = (x - a) / a;
 476    T phi = -boost::math::log1pmx(sigma, pol);
 477    T y = a * phi;
 478    T z = sqrt(2 * phi);
 479    if(x < a)
 480       z = -z;
 481
 482    T workspace[14];
 483
 484    static const T C0[] = {
 485       -0.333333333333333333333333333333333333L,
 486       0.0833333333333333333333333333333333333L,
 487       -0.0148148148148148148148148148148148148L,
 488       0.00115740740740740740740740740740740741L,
 489       0.0003527336860670194003527336860670194L,
 490       -0.000178755144032921810699588477366255144L,
 491       0.391926317852243778169704095630021556e-4L,
 492       -0.218544851067999216147364295512443661e-5L,
 493       -0.185406221071515996070179883622956325e-5L,
 494       0.829671134095308600501624213166443227e-6L,
 495       -0.17665952736826079304360054245742403e-6L,
 496       0.670785354340149858036939710029613572e-8L,
 497       0.102618097842403080425739573227252951e-7L,
 498       -0.438203601845335318655297462244719123e-8L,
 499       0.914769958223679023418248817633113681e-9L,
 500       -0.255141939949462497668779537993887013e-10L,
 501       -0.583077213255042506746408945040035798e-10L,
 502       0.243619480206674162436940696707789943e-10L,
 503       -0.502766928011417558909054985925744366e-11L,
 504       0.110043920319561347708374174497293411e-12L,
 505       0.337176326240098537882769884169200185e-12L,
 506       -0.13923887224181620659193661848957998e-12L,
 507       0.285348938070474432039669099052828299e-13L,
 508       -0.513911183424257261899064580300494205e-15L,
 509       -0.197522882943494428353962401580710912e-14L,
 510       0.809952115670456133407115668702575255e-15L,
 511       -0.165225312163981618191514820265351162e-15L,
 512       0.253054300974788842327061090060267385e-17L,
 513       0.116869397385595765888230876507793475e-16L,
 514       -0.477003704982048475822167804084816597e-17L,
 515       0.969912605905623712420709685898585354e-18L,
 516    };
 517    workspace[0] = tools::evaluate_polynomial(C0, z);
 518
 519    static const T C1[] = {
 520       -0.00185185185185185185185185185185185185L,
 521       -0.00347222222222222222222222222222222222L,
 522       0.0026455026455026455026455026455026455L,
 523       -0.000990226337448559670781893004115226337L,
 524       0.000205761316872427983539094650205761317L,
 525       -0.401877572016460905349794238683127572e-6L,
 526       -0.180985503344899778370285914867533523e-4L,
 527       0.76491609160811100846374214980916921e-5L,
 528       -0.16120900894563446003775221882217767e-5L,
 529       0.464712780280743434226135033938722401e-8L,
 530       0.137863344691572095931187533077488877e-6L,
 531       -0.575254560351770496402194531835048307e-7L,
 532       0.119516285997781473243076536699698169e-7L,
 533       -0.175432417197476476237547551202312502e-10L,
 534       -0.100915437106004126274577504686681675e-8L,
 535       0.416279299184258263623372347219858628e-9L,
 536       -0.856390702649298063807431562579670208e-10L,
 537       0.606721510160475861512701762169919581e-13L,
 538       0.716249896481148539007961017165545733e-11L,
 539       -0.293318664377143711740636683615595403e-11L,
 540       0.599669636568368872330374527568788909e-12L,
 541       -0.216717865273233141017100472779701734e-15L,
 542       -0.497833997236926164052815522048108548e-13L,
 543       0.202916288237134247736694804325894226e-13L,
 544       -0.413125571381061004935108332558187111e-14L,
 545       0.828651623988309644380188591057589316e-18L,
 546       0.341003088693333279336339355910600992e-15L,
 547       -0.138541953028939715357034547426313703e-15L,
 548       0.281234665322887466568860332727259483e-16L,
 549    };
 550    workspace[1] = tools::evaluate_polynomial(C1, z);
 551
 552    static const T C2[] = {
 553       0.0041335978835978835978835978835978836L,
 554       -0.00268132716049382716049382716049382716L,
 555       0.000771604938271604938271604938271604938L,
 556       0.200938786008230452674897119341563786e-5L,
 557       -0.000107366532263651605215391223621676297L,
 558       0.529234488291201254164217127180090143e-4L,
 559       -0.127606351886187277133779191392360117e-4L,
 560       0.34235787340961380741902003904747389e-7L,
 561       0.137219573090629332055943852926020279e-5L,
 562       -0.629899213838005502290672234278391876e-6L,
 563       0.142806142060642417915846008822771748e-6L,
 564       -0.204770984219908660149195854409200226e-9L,
 565       -0.140925299108675210532930244154315272e-7L,
 566       0.622897408492202203356394293530327112e-8L,
 567       -0.136704883966171134992724380284402402e-8L,
 568       0.942835615901467819547711211663208075e-12L,
 569       0.128722524000893180595479368872770442e-9L,
 570       -0.556459561343633211465414765894951439e-10L,
 571       0.119759355463669810035898150310311343e-10L,
 572       -0.416897822518386350403836626692480096e-14L,
 573       -0.109406404278845944099299008640802908e-11L,
 574       0.4662239946390135746326204922464679e-12L,
 575       -0.990510576390690597844122258212382301e-13L,
 576       0.189318767683735145056885183170630169e-16L,
 577       0.885922187259112726176031067028740667e-14L,
 578       -0.373782039804640545306560251777191937e-14L,
 579       0.786883363903515525774088394065960751e-15L,
 580    };
 581    workspace[2] = tools::evaluate_polynomial(C2, z);
 582
 583    static const T C3[] = {
 584       0.000649434156378600823045267489711934156L,
 585       0.000229472093621399176954732510288065844L,
 586       -0.000469189494395255712128140111679206329L,
 587       0.000267720632062838852962309752433209223L,
 588       -0.756180167188397641072538191879755666e-4L,
 589       -0.239650511386729665193314027333231723e-6L,
 590       0.110826541153473023614770299726861227e-4L,
 591       -0.567495282699159656749963105701560205e-5L,
 592       0.14230900732435883914551894470580433e-5L,
 593       -0.278610802915281422405802158211174452e-10L,
 594       -0.16958404091930277289864168795820267e-6L,
 595       0.809946490538808236335278504852724081e-7L,
 596       -0.191111684859736540606728140872727635e-7L,
 597       0.239286204398081179686413514022282056e-11L,
 598       0.206201318154887984369925818486654549e-8L,
 599       -0.946049666185513217375417988510192814e-9L,
 600       0.215410497757749078380130268468744512e-9L,
 601       -0.138882333681390304603424682490735291e-13L,
 602       -0.218947616819639394064123400466489455e-10L,
 603       0.979099895117168512568262802255883368e-11L,
 604       -0.217821918801809621153859472011393244e-11L,
 605       0.62088195734079014258166361684972205e-16L,
 606       0.212697836327973697696702537114614471e-12L,
 607       -0.934468879151743333127396765626749473e-13L,
 608       0.204536712267828493249215913063207436e-13L,
 609    };
 610    workspace[3] = tools::evaluate_polynomial(C3, z);
 611
 612    static const T C4[] = {
 613       -0.000861888290916711698604702719929057378L,
 614       0.00078403922172006662747403488144228885L,
 615       -0.000299072480303190179733389609932819809L,
 616       -0.146384525788434181781232535690697556e-5L,
 617       0.664149821546512218665853782451862013e-4L,
 618       -0.396836504717943466443123507595386882e-4L,
 619       0.113757269706784190980552042885831759e-4L,
 620       0.250749722623753280165221942390057007e-9L,
 621       -0.169541495365583060147164356781525752e-5L,
 622       0.890750753220530968882898422505515924e-6L,
 623       -0.229293483400080487057216364891158518e-6L,
 624       0.295679413754404904696572852500004588e-10L,
 625       0.288658297427087836297341274604184504e-7L,
 626       -0.141897394378032193894774303903982717e-7L,
 627       0.344635804994648970659527720474194356e-8L,
 628       -0.230245171745280671320192735850147087e-12L,
 629       -0.394092330280464052750697640085291799e-9L,
 630       0.186023389685045019134258533045185639e-9L,
 631       -0.435632300505661804380678327446262424e-10L,
 632       0.127860010162962312660550463349930726e-14L,
 633       0.467927502665791946200382739991760062e-11L,
 634       -0.214924647061348285410535341910721086e-11L,
 635       0.490881561480965216323649688463984082e-12L,
 636    };
 637    workspace[4] = tools::evaluate_polynomial(C4, z);
 638
 639    static const T C5[] = {
 640       -0.000336798553366358150308767592718210002L,
 641       -0.697281375836585777429398828575783308e-4L,
 642       0.00027727532449593920787336425196507501L,
 643       -0.000199325705161888477003360405280844238L,
 644       0.679778047793720783881640176604435742e-4L,
 645       0.141906292064396701483392727105575757e-6L,
 646       -0.135940481897686932784583938837504469e-4L,
 647       0.80184702563342015397192571980419684e-5L,
 648       -0.229148117650809517038048790128781806e-5L,
 649       -0.325247355129845395166230137750005047e-9L,
 650       0.346528464910852649559195496827579815e-6L,
 651       -0.184471871911713432765322367374920978e-6L,
 652       0.482409670378941807563762631738989002e-7L,
 653       -0.179894667217435153025754291716644314e-13L,
 654       -0.630619450001352343517516981425944698e-8L,
 655       0.316241762877456793773762181540969623e-8L,
 656       -0.784092425369742929000839303523267545e-9L,
 657    };
 658    workspace[5] = tools::evaluate_polynomial(C5, z);
 659
 660    static const T C6[] = {
 661       0.00053130793646399222316574854297762391L,
 662       -0.000592166437353693882864836225604401187L,
 663       0.000270878209671804482771279183488328692L,
 664       0.790235323266032787212032944390816666e-6L,
 665       -0.815396936756196875092890088464682624e-4L,
 666       0.561168275310624965003775619041471695e-4L,
 667       -0.183291165828433755673259749374098313e-4L,
 668       -0.307961345060330478256414192546677006e-8L,
 669       0.346515536880360908673728529745376913e-5L,
 670       -0.202913273960586037269527254582695285e-5L,
 671       0.578879286314900370889997586203187687e-6L,
 672       0.233863067382665698933480579231637609e-12L,
 673       -0.88286007463304835250508524317926246e-7L,
 674       0.474359588804081278032150770595852426e-7L,
 675       -0.125454150207103824457130611214783073e-7L,
 676    };
 677    workspace[6] = tools::evaluate_polynomial(C6, z);
 678
 679    static const T C7[] = {
 680       0.000344367606892377671254279625108523655L,
 681       0.517179090826059219337057843002058823e-4L,
 682       -0.000334931610811422363116635090580012327L,
 683       0.000281269515476323702273722110707777978L,
 684       -0.000109765822446847310235396824500789005L,
 685       -0.127410090954844853794579954588107623e-6L,
 686       0.277444515115636441570715073933712622e-4L,
 687       -0.182634888057113326614324442681892723e-4L,
 688       0.578769494973505239894178121070843383e-5L,
 689       0.493875893393627039981813418398565502e-9L,
 690       -0.105953670140260427338098566209633945e-5L,
 691       0.616671437611040747858836254004890765e-6L,
 692       -0.175629733590604619378669693914265388e-6L,
 693    };
 694    workspace[7] = tools::evaluate_polynomial(C7, z);
 695
 696    static const T C8[] = {
 697       -0.000652623918595309418922034919726622692L,
 698       0.000839498720672087279993357516764983445L,
 699       -0.000438297098541721005061087953050560377L,
 700       -0.696909145842055197136911097362072702e-6L,
 701       0.00016644846642067547837384572662326101L,
 702       -0.000127835176797692185853344001461664247L,
 703       0.462995326369130429061361032704489636e-4L,
 704       0.455790986792270771162749294232219616e-8L,
 705       -0.105952711258051954718238500312872328e-4L,
 706       0.678334290486516662273073740749269432e-5L,
 707       -0.210754766662588042469972680229376445e-5L,
 708    };
 709    workspace[8] = tools::evaluate_polynomial(C8, z);
 710
 711    static const T C9[] = {
 712       -0.000596761290192746250124390067179459605L,
 713       -0.720489541602001055908571930225015052e-4L,
 714       0.000678230883766732836161951166000673426L,
 715       -0.000640147526026275845100045652582354779L,
 716       0.000277501076343287044992374518205845463L,
 717       0.181970083804651510461686554030325202e-6L,
 718       -0.847950711706850318239732559632810086e-4L,
 719       0.610519208250153101764709122740859458e-4L,
 720       -0.210739201834048624082975255893773306e-4L,
 721    };
 722    workspace[9] = tools::evaluate_polynomial(C9, z);
 723
 724    static const T C10[] = {
 725       0.00133244544948006563712694993432717968L,
 726       -0.00191443849856547752650089885832852254L,
 727       0.0011089369134596637339607446329267522L,
 728       0.993240412264229896742295262075817566e-6L,
 729       -0.000508745012930931989848393025305956774L,
 730       0.00042735056665392884328432271160040444L,
 731       -0.000168588537679107988033552814662382059L,
 732    };
 733    workspace[10] = tools::evaluate_polynomial(C10, z);
 734
 735    static const T C11[] = {
 736       0.00157972766073083495908785631307733022L,
 737       0.000162516262783915816898635123980270998L,
 738       -0.00206334210355432762645284467690276817L,
 739       0.00213896861856890981541061922797693947L,
 740       -0.00101085593912630031708085801712479376L,
 741    };
 742    workspace[11] = tools::evaluate_polynomial(C11, z);
 743
 744    static const T C12[] = {
 745       -0.00407251211951401664727281097914544601L,
 746       0.00640336283380806979482363809026579583L,
 747       -0.00404101610816766177473974858518094879L,
 748    };
 749    workspace[12] = tools::evaluate_polynomial(C12, z);
 750    workspace[13] = -0.0059475779383993002845382844736066323L;
 751
 752    T result = tools::evaluate_polynomial(workspace, 1/a);
 753    result *= exp(-y) / sqrt(2 * constants::pi<T>() * a);
 754    if(x < a)
 755       result = -result;
 756
 757    result += boost::math::erfc(sqrt(y), pol) / 2;
 758
 759    return result;
 760 }
 761
 762
 763 }  // namespace detail
 764 }  // namespace math
 765 }  // namespace math
 766
 767
 768 #endif // BOOST_MATH_DETAIL_IGAMMA_LARGE
 769