Mercurial > hg > vamp-build-and-test

annotate DEPENDENCIES/generic/include/boost/math/special_functions/detail/igamma_large.hpp @ 133:4acb5d8d80b6 tip

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