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