1 #include "openturns/SpecFunc.hxx"
2 #include "openturns/incgam.hxx"
3
4 BEGIN_NAMESPACE_OPENTURNS
5
6 namespace GammaFunctions
7 {
8 #define explow -300
9 #define exphigh 300
incgam(const Scalar a,const Scalar x,Scalar & p,Scalar & q,SignedInteger & ierr)10 void incgam(const Scalar a,
11 const Scalar x,
12 Scalar & p,
13 Scalar & q,
14 SignedInteger & ierr)
15 {
16 // -------------------------------------------------------------
17 // Calculation of the incomplete gamma functions ratios P(a,x)
18 // and Q(a,x).
19 // -------------------------------------------------------------
20 // Inputs:
21 // a , argument of the functions
22 // x , argument of the functions
23 // Outputs:
24 // p, function P(a,x)
25 // q, function Q(a,x)
26 // ierr , error flag
27 // ierr = 0, computation succesful
28 // ierr = 1, overflow/underflow problems. The function values
29 // (P(a,x) and Q(a,x)) are set to zero.
30 // ----------------------------------------------------------------------
31 // Authors:
32 // Amparo Gil (U. Cantabria, Santander, Spain)
33 // e-mail: amparo.gil@unican.es
34 // Javier Segura (U. Cantabria, Santander, Spain)
35 // e-mail: javier.segura@unican.es
36 // Nico M. Temme (CWI, Amsterdam, The Netherlands)
37 // e-mail: nico.temme@cwi.nl
38 // -------------------------------------------------------------
39 // References: "Efficient and accurate algorithms for
40 // the computation and inversion of the incomplete gamma function ratios",
41 // A. Gil, J. Segura and N.M. Temme, submitted to SIAM J Sci Comput
42 // -------------------------------------------------------------------
43 Scalar lnx = -1.0;
44 ierr = 0;
45 if (a > alpha(x))
46 {
47 const Scalar dp = dompart(a, x, false);
48 if (dp < 0.0)
49 {
50 ierr = 1;
51 p = 0.0;
52 q = 0.0;
53 return;
54 }
55 if ((x < 0.3 * a) || (a < 12)) p = ptaylor(a, x, dp);
56 else p = pqasymp(a, x, dp, true);
57 q = 1.0 - p;
58 return;
59 } // a > alpha(x)
60 if (x < SpecFunc::MinScalar) lnx = SpecFunc::LogMinScalar;
61 else lnx = log(x);
62 if (a < -SpecFunc::MinScalar / lnx)
63 {
64 q = 0.0;
65 p = 1.0;
66 return;
67 }
68 if (x < 1.0)
69 {
70 const Scalar dp = dompart(a, x, true);
71 if (dp < 0.0)
72 {
73 ierr = 1;
74 q = 0.0;
75 p = 0.0;
76 return;
77 }
78 q = qtaylor(a, x, dp);
79 p = 1.0 - q;
80 return;
81 } // x < 1.0
82 const Scalar dp = dompart(a, x, false);
83 if (dp < 0.0)
84 {
85 ierr = 1;
86 p = 0.0;
87 q = 0.0;
88 return;
89 }
90 if ((x > 2.35 * a) || (a < 12)) q = qfraction(a, x, dp);
91 else q = pqasymp(a, x, dp, false);
92 p = 1.0 - q;
93 } // incgam
94
invincgam(const Scalar a,const Scalar p,const Scalar q,Scalar & xr,SignedInteger & ierr)95 void invincgam(const Scalar a,
96 const Scalar p,
97 const Scalar q,
98 Scalar & xr,
99 SignedInteger & ierr)
100 {
101 // -------------------------------------------------------------
102 // invincgam computes xr in the equations P(a, xr) = p and Q(a, xr) = q
103 // with a as a given positive parameter.
104 // In most cases, we invert the equation with min(p, q)
105 // -------------------------------------------------------------
106 // Inputs:
107 // a , argument of the functions
108 // p, function P(a, x)
109 // q, function Q(a, x)
110 // Outputs:
111 // xr , soluction of the equations P(a, xr) = p and Q(a, xr) = q
112 // with a as a given positive parameter.
113 // ierr , error flag
114 // ierr = 0, computation succesful
115 // ierr = -1, overflow problem in the computation of one of the
116 // gamma factors before starting the Newton iteration.
117 // The initial approximation to the root is given
118 // as output.
119 // ierr = -2, the number of iterations in the Newton method
120 // reached the upper limit N = 15. The last value
121 // obtained for the root is given as output.
122 // ------------------------------------------------------------------
123 ierr = 0;
124 Bool pcase = true;
125 Scalar porq = p;
126 Scalar s = -1.0;
127 if (p > 0.5)
128 {
129 pcase = false;
130 porq = q;
131 s = 1.0;
132 } // p > 0.5
133 const Scalar logr = (1.0 / a) * (log(p) + SpecFunc::LogGamma(a + 1.0));
134 Scalar x = 0.0;
135 Bool m0 = true;
136 Scalar eta = 0.0;
137 if (logr < log(0.2 * (1.0 + a)))
138 {
139 const Scalar r = exp(logr);
140 const Scalar ap1 = a + 1.0;
141 const Scalar ap12 = ap1 * ap1;
142 const Scalar ap13 = ap1 * ap12;
143 const Scalar ap14 = ap12 * ap12;
144 const Scalar ap2 = a + 2.0;
145 const Scalar ap22 = ap2 * ap2;
146 const Scalar ck1 = 1.0 / ap1;
147 const Scalar ck2 = 0.5 * (3.0 * a + 5.0) / (ap12 * ap2);
148 const Scalar ck3 = (31.0 + a * (33.0 + 8.0 * a)) / (3.0 * (ap13 * ap2 * (a + 3.0)));
149 const Scalar ck4 = (2888.0 + a * (5661.0 + a * (3971.0 + a * (1179.0 + 125.0 * a)))) / (24.0 * (ap14 * ap22 * (a + 3.0) * (a + 4.0)));
150 x = r * (1.0 + r * (ck1 + r * (ck2 + r * (ck3 + r * ck4))));
151 } // logr < log(0.2 * (1 + a))
152 else if ((q < std::min(0.02, exp(-1.5 * a) / SpecFunc::Gamma(a))) && (a < 10.0))
153 {
154 const Scalar b = 1.0 - a;
155 eta = sqrt(-2.0 / a * log(q * gamstar(a) * SpecFunc::SQRT2PI / sqrt(a)));
156 x = a * lambdaeta(eta);
157 const Scalar L = log(x);
158 const Scalar r = 1.0 / x;
159 if ((a > 0.12) || (x > 5.0))
160 {
161 const Scalar ck0 = L - 1.0;
162 const Scalar ck1 = 1.0 + 1.5 * b + L * (-b - 1.0 + 0.5 * L);
163 const Scalar ck2 = -2.0 + (-4.0 - 11.0 / 6.0 * b) * b + (2.0 + (4.0 + b) * b + (-1.5 * b - 1.0 + L / 3.0) * L) * L;
164 const Scalar ck3 = 6.0 + (13.5 + (10.0 + 25.0 / 12.0 * b) * b) * b + (-6.0 + (-14.0 + (-9.5 - b) * b) * b + (3.0 + (7.0 + 3.0 * b) * b + (-1.0 - 11.0 / 6.0 * b + 0.25 * L) * L) * L) * L;
165 x += -L + b * r * (ck0 + r * (ck1 + r * (ck2 + r * ck3)));
166 } // (a > 0.12) || (x > 5.0)
167 else x += -L + b * r * (L - 1.0);
168 }
169 else if (std::abs(porq - 0.5) < 1.0e-5) x = a - 1.0 / 3.0 + (8.0 / 405.0 + 184.0 / 25515.0 / a) / a;
170 else if (std::abs(a - 1.0) < 1.0e-4) x = (pcase ? -log1p(-p) : -log(q));
171 else if (a < 1.0) x = (pcase ? exp((1.0 / a) * (log(porq) + SpecFunc::LogGamma(a + 1.0))) : exp((1.0 / a) * (log1p(-porq) + SpecFunc::LogGamma(a + 1.0))));
172 else
173 {
174 // a >= 1.0
175 m0 = false;
176 const Scalar r = inverfc(2.0 * porq);
177 eta = s * r / sqrt(a * 0.5);
178 eta += (eps1(eta) + (eps2(eta) + eps3(eta) / a) / a) / a;
179 x = a * lambdaeta(eta);
180 } // a >= 1.0
181 Scalar t = 1.0;
182 UnsignedInteger n = 1;
183 // Implementation of the high order Newton-like method;
184 while ((t > SpecFunc::ScalarEpsilon) && (n < 15))
185 {
186 Scalar dx = 0.0;
187 Scalar r = 0.0;
188 if (m0)
189 {
190 const Scalar dlnr = (1.0 - a) * log(x) + x + SpecFunc::LogGamma(a);
191 if (dlnr > SpecFunc::LogMaxScalar)
192 {
193 n = 20;
194 ierr = -1;
195 continue;
196 } // dlnr > log(giant)
197 r = exp(dlnr);
198 }
199 else r = x * gamstar(a) / (sqrt(a) * SpecFunc::ISQRT2PI * exp(-0.5 * a * eta * eta));
200 SignedInteger ierrf;
201 Scalar px = -1.0;
202 Scalar qx = -1.0;
203 incgam(a, x, px, qx, ierrf);
204 r = (pcase ? -r * (px - p) : r * (qx - q));
205 if (a <= 0.05) dx = r;
206 else
207 {
208 // a > 0.05
209 const Scalar ck1 = (x - a + 1.0) / (2.0 * x);
210 if (a <= 0.1) dx = r * (1.0 + r * ck1);
211 else
212 {
213 // a > 0.1
214 const Scalar ck2 = (1.0 + (-3.0 + 2.0 * a) * a + (4.0 - 4.0 * a + 2.0 * x) * x) / (6.0 * x * x);
215 dx = r * (1.0 + r * (ck1 + r * ck2));
216 } // a > 0.1
217 } // a > 0.05
218 x += dx;
219 t = std::abs(dx / x);
220 ++n;
221 } // (t > 1.0e-15) && (n < 15)
222 if (n == 15) ierr = -2;
223 xr = x;
224 } // invincgam
225
exmin1(const Scalar x)226 Scalar exmin1(const Scalar x)
227 {
228 // computes (exp(x)-1)/x;
229 if (std::abs(x) < 3.65e-8) return 1.0 + 0.5 * x;
230 return expm1(x) / x;
231 } // exmin1
232
lnec(const Scalar x)233 Scalar lnec(const Scalar x)
234 {
235 // x > -1; lnec: = ln1: = ln(1+x)-x
236 if (std::abs(x) < 1.3e-5) return x * x * (-0.5 + x * (1.0 / 3.0 - 0.25 * x));
237 return log1p(x) - x;
238 } // lnec
239
alpha(const Scalar x)240 Scalar alpha(const Scalar x)
241 {
242 if (x > 0.25) return x + 0.25;
243 if (x >= SpecFunc::MinScalar) return -0.6931 / log(x);
244 return -0.6931 / SpecFunc::LogMinScalar;
245 } // alpha
246
dompart(const Scalar a,const Scalar x,const Bool qt)247 Scalar dompart(const Scalar a,
248 const Scalar x,
249 const Bool qt)
250 {
251 // dompart is approx. of x^a * exp(-x) / gamma(a+1) ;
252 Scalar r = -1.0;
253 const Scalar lnx = log(x);
254 if (a <= 1.0) r = -x + a * lnx;
255 else
256 {
257 if (x == a) r = 0.0;
258 else
259 {
260 const Scalar la = x / a;
261 r = a * (1.0 - la + log(la));
262 }
263 r += -0.5 * log(6.2832 * a);
264 }
265 const Scalar dp = r < explow ? 0.0 : exp(r);
266 if (qt) return dp;
267 if ((a < 3.0) || (x < 0.2)) return exp(a * lnx - x) / SpecFunc::Gamma(a + 1.0);
268 const Scalar mu = (x - a) / a;
269 const Scalar c = lnec(mu);
270 if ((a * c) > SpecFunc::LogMaxScalar) return -100.0;
271 return exp(a * c) / (sqrt(a) * SpecFunc::SQRT2PI * gamstar(a));
272 } // dompart
273
chepolsum(const Scalar x,const Point & a)274 Scalar chepolsum(const Scalar x,
275 const Point & a)
276 {
277 //{a[0]/2+a[1]T1(x)+...a[n]Tn(x); series of Chebychev polynomials}
278 const UnsignedInteger n = a.getDimension() - 1;
279 if (n == 0) return 0.5 * a[0];
280 if (n == 1) return 0.5 * a[0] + x * a[1];
281 const Scalar tx = x + x;
282 Scalar r = a[n];
283 Scalar h = a[n - 1] + r * tx;
284 for (UnsignedInteger k = n - 2; k >= 1; --k)
285 {
286 const Scalar s = r;
287 r = h;
288 h = a[k] + r * tx - s;
289 }
290 return 0.5 * a[0] - r + h * x;
291 } // chepolsum
292
auxgam(const Scalar x)293 Scalar auxgam(const Scalar x)
294 {
295 // function g in 1/gamma(x+1) = 1+x*(x-1)*g(x), -1 <= x <= 1
296 if (x < 0.0) return -(1.0 + (1.0 + x) * (1.0 + x) * auxgam(1.0 + x)) / (1.0 - x);
297 Point dr(18);
298 dr[0] = -1.013609258009865776949;
299 dr[1] = 0.784903531024782283535e-1;
300 dr[2] = 0.67588668743258315530e-2;
301 dr[3] = -0.12790434869623468120e-2;
302 dr[4] = 0.462939838642739585e-4;
303 dr[5] = 0.43381681744740352e-5;
304 dr[6] = -0.5326872422618006e-6;
305 dr[7] = 0.172233457410539e-7;
306 dr[8] = 0.8300542107118e-9;
307 dr[9] = -0.10553994239968e-9;
308 dr[10] = 0.39415842851e-11;
309 dr[11] = 0.362068537e-13;
310 dr[12] = -0.107440229e-13;
311 dr[13] = 0.5000413e-15;
312 dr[14] = -0.62452e-17;
313 dr[15] = -0.5185e-18;
314 dr[16] = 0.347e-19;
315 dr[17] = -0.9e-21;
316 return chepolsum(2.0 * x - 1.0, dr);
317 } // auxgam
318
gamstar(const Scalar x)319 Scalar gamstar(const Scalar x)
320 {
321 // gamstar(x) = exp(SpecFunc::GammaCorrection(x)), x > 0; or
322 // gamma(x)/(exp(-x+(x-0.5)*ln(x))/sqrt(2pi)
323 if (x >= 3.0) return exp(SpecFunc::GammaCorrection(x));
324 if (x > 0.0) return SpecFunc::Gamma(x) / (exp(-x + (x - 0.5) * log(x)) * SpecFunc::SQRT2PI);
325 return SpecFunc::MaxScalar;
326 } // gamstar
327
fractio(const Scalar x,const Point & r,const Point & s)328 Scalar fractio(const Scalar x,
329 const Point & r,
330 const Point & s)
331 {
332 const UnsignedInteger n = r.getSize() - 1;
333 Scalar a = r[n];
334 Scalar b = 1.0;
335 for (SignedInteger k = n - 1; k >= 0; --k)
336 {
337 a = a * x + r[k];
338 b = b * x + s[k];
339 }
340 return a / b;
341 } // fractio
342
pqasymp(const Scalar a,const Scalar x,const Scalar dp,const Bool p)343 Scalar pqasymp(const Scalar a,
344 const Scalar x,
345 const Scalar dp,
346 const Bool p)
347 {
348 if (dp == 0.0) return (p ? 0.0 : 1.0);
349 const SignedInteger s = p ? -1 : 1;
350 const Scalar mu = x / a - 1.0;
351 Scalar y = -lnec(mu);
352 Scalar eta = y < 0.0 ? 0.0 : sqrt(2.0 * y);
353 y *= a;
354 Scalar v = sqrt(std::abs(y));
355 if (mu < 0.0)
356 {
357 eta = -eta;
358 v = -v;
359 }
360 const Scalar u = 0.5 * SpecFunc::ErfC(s * v);
361 v = s * exp(-y) * saeta(a, eta) / (SpecFunc::SQRT2PI * sqrt(a));
362 return u + v;
363 } // FUNCTION pqasymp;
364
saeta(const Scalar a,const Scalar eta)365 Scalar saeta(const Scalar a,
366 const Scalar eta)
367 {
368 Point fm(27);
369 fm[0] = 1.0;
370 fm[1] = -1.0 / 3.0;
371 fm[2] = 1.0 / 12.0;
372 fm[3] = -2.0 / 135.0;
373 fm[4] = 1.0 / 864.0;
374 fm[5] = 1.0 / 2835.0;
375 fm[6] = -139.0 / 777600.0;
376 fm[7] = 1.0 / 25515.0;
377 fm[8] = -571.0 / 261273600.0;
378 fm[9] = -281.0 / 151559100.0;
379 fm[10] = 8.29671134095308601e-7;
380 fm[11] = -1.76659527368260793e-7;
381 fm[12] = 6.70785354340149857e-9;
382 fm[13] = 1.02618097842403080e-8;
383 fm[14] = -4.38203601845335319e-9;
384 fm[15] = 9.14769958223679023e-10;
385 fm[16] = -2.55141939949462497e-11;
386 fm[17] = -5.83077213255042507e-11;
387 fm[18] = 2.43619480206674162e-11;
388 fm[19] = -5.02766928011417559e-12;
389 fm[20] = 1.10043920319561347e-13;
390 fm[21] = 3.37176326240098538e-13;
391 fm[22] = -1.39238872241816207e-13;
392 fm[23] = 2.85348938070474432e-14;
393 fm[24] = -5.13911183424257258e-16;
394 fm[25] = -1.97522882943494428e-15;
395 fm[26] = 8.09952115670456133e-16;
396 Point bm(27);
397 bm[25] = fm[26];
398 bm[24] = fm[25];
399 for (UnsignedInteger m = 24; m > 0; --m) bm[m - 1] = fm[m] + (m + 1) * bm[m + 1] / a;
400 Scalar s = bm[0];
401 Scalar t = s;
402 Scalar y = eta;
403 for (UnsignedInteger m = 1; m < 25; ++m)
404 {
405 if (std::abs(t / s) <= SpecFunc::ScalarEpsilon) break;
406 t = bm[m] * y;
407 s += t;
408 y *= eta;
409 }
410 return s / (1.0 + bm[1] / a);
411 } // saeta
412
qfraction(const Scalar a,const Scalar x,const Scalar dp)413 Scalar qfraction(const Scalar a,
414 const Scalar x,
415 const Scalar dp)
416 {
417 if (dp == 0.0) return 0.0;
418 Scalar p = 0.0;
419 Scalar q = (x - 1.0 - a) * (x + 1.0 - a);
420 Scalar r = 4.0 * (x + 1.0 - a);
421 Scalar s = 1.0 - a;
422 Scalar ro = 0.0;
423 Scalar t = 1.0;
424 Scalar g = 1.0;
425 while (std::abs(t / g) >= SpecFunc::ScalarEpsilon)
426 {
427 p += s;
428 q += r;
429 r += 8.0;
430 s += 2.0;
431 const Scalar tau = p * (1.0 + ro);
432 ro = tau / (q - tau);
433 t *= ro;
434 g += t;
435 }
436 return (a / (x + 1.0 - a)) * g * dp;
437 } // qfraction
438
qtaylor(const Scalar a,const Scalar x,const Scalar dp)439 Scalar qtaylor(const Scalar a,
440 const Scalar x,
441 const Scalar dp)
442 {
443 if (dp == 0.0) return 0.0;
444 const Scalar lnx = log(x);
445 Scalar r = a * lnx;
446 Scalar q(r * exmin1(r)); // {q = x^a - 1}
447 Scalar s(a * (1.0 - a) * auxgam(a)); // {s = 1 - 1 / Gamma(1 + a)}
448 q *= 1.0 - s;
449 Scalar u(s - q); // {u = 1 - x^a/Gamma(1+a)}
450 Scalar p = a * x;
451 Scalar t = 1.0;
452 Scalar v = 1.0;
453 q = a + 1.0;
454 r = a + 3.0;
455 while (std::abs(t / v) > SpecFunc::ScalarEpsilon)
456 {
457 p += x;
458 q += r;
459 r += 2.0;
460 t *= -p / q;
461 v += t;
462 }
463 v *= a * (1.0 - s) * exp((a + 1.0) * lnx) / (a + 1.0);
464 return u + v;
465 } // qtaylor
466
ptaylor(const Scalar a,const Scalar x,const Scalar dp)467 Scalar ptaylor(const Scalar a,
468 const Scalar x,
469 const Scalar dp)
470 {
471 if (dp == 0.0) return 0.0;
472 Scalar p = 1.0;
473 Scalar c = 1.0;
474 Scalar r = a;
475 while (c > p * SpecFunc::ScalarEpsilon)
476 {
477 r += 1.0;
478 c *= x / r;
479 p += c;
480 }
481 return p * dp;
482 } // ptaylor
483
eps1(const Scalar eta)484 Scalar eps1(const Scalar eta)
485 {
486 if (std::abs(eta) < 1.0)
487 {
488 Point ak(5);
489 ak[0] = -3.333333333438e-1;
490 ak[1] = -2.070740359969e-1;
491 ak[2] = -5.041806657154e-2;
492 ak[3] = -4.923635739372e-3;
493 ak[4] = -4.293658292782e-5;
494 Point bk(5);
495 bk[0] = 1.000000000000e+0;
496 bk[1] = 7.045554412463e-1;
497 bk[2] = 2.118190062224e-1;
498 bk[3] = 3.048648397436e-2;
499 bk[4] = 1.605037988091e-3;
500 return ratfun(eta, ak, bk);
501 }
502 return log(eta / (lambdaeta(eta) - 1.0)) / eta;
503 } // eps1
504
eps2(const Scalar eta)505 Scalar eps2(const Scalar eta)
506 {
507 if (eta < -5.0)
508 {
509 const Scalar x = eta * eta;
510 const Scalar lnmeta = log(-eta);
511 return (12.0 - x - 6.0 * (lnmeta * lnmeta)) / (12.0 * x * eta);
512 }
513 if (eta < -2.0)
514 {
515 Point ak(5);
516 ak[0] = -1.72847633523e-2;
517 ak[1] = -1.59372646475e-2;
518 ak[2] = -4.64910887221e-3;
519 ak[3] = -6.06834887760e-4;
520 ak[4] = -6.14830384279e-6;
521 Point bk(5);
522 bk[0] = 1.00000000000e+0;
523 bk[1] = 7.64050615669e-1;
524 bk[2] = 2.97143406325e-1;
525 bk[3] = 5.79490176079e-2;
526 bk[4] = 5.74558524851e-3;
527 return ratfun(eta, ak, bk);
528 }
529 if (eta < 2.0)
530 {
531 Point ak(5);
532 ak[0] = -1.72839517431e-2;
533 ak[1] = -1.46362417966e-2;
534 ak[2] = -3.57406772616e-3;
535 ak[3] = -3.91032032692e-4;
536 ak[4] = 2.49634036069e-6;
537 Point bk(5);
538 bk[0] = 1.00000000000e+0;
539 bk[1] = 6.90560400696e-1;
540 bk[2] = 2.49962384741e-1;
541 bk[3] = 4.43843438769e-2;
542 bk[4] = 4.24073217211e-3;
543 return ratfun(eta, ak, bk);
544 }
545 if (eta < 1000.0)
546 {
547 Point ak(5);
548 ak[0] = 9.99944669480e-1;
549 ak[1] = 1.04649839762e+2;
550 ak[2] = 8.57204033806e+2;
551 ak[3] = 7.31901559577e+2;
552 ak[4] = 4.55174411671e+1;
553 Point bk(5);
554 bk[0] = 1.00000000000e+0;
555 bk[1] = 1.04526456943e+2;
556 bk[2] = 8.23313447808e+2;
557 bk[3] = 3.11993802124e+3;
558 bk[4] = 3.97003311219e+3;
559 return ratfun(1.0 / eta, ak, bk) / (-12.0 * eta);
560 }
561 return -1.0 / (12.0 * eta);
562 } // eps2
563
eps3(const Scalar eta)564 Scalar eps3(const Scalar eta)
565 {
566 if (eta < -8.0)
567 {
568 const Scalar x = eta * eta;
569 const Scalar y = log(-eta) / eta;
570 return (-30.0 + eta * y * (6.0 * x * y * y - 12.0 + x)) / (12.0 * eta * x * x);
571 }
572 if (eta < -4.0)
573 {
574 Point ak(5);
575 ak[0] = 4.95346498136e-2;
576 ak[1] = 2.99521337141e-2;
577 ak[2] = 6.88296911516e-3;
578 ak[3] = 5.12634846317e-4;
579 ak[4] = -2.01411722031e-5;
580 Point bk(5);
581 bk[0] = 1.00000000000e+0;
582 bk[1] = 7.59803615283e-1;
583 bk[2] = 2.61547111595e-1;
584 bk[3] = 4.64854522477e-2;
585 bk[4] = 4.03751193496e-3;
586 return ratfun(eta, ak, bk) / (eta * eta);
587 }
588 if (eta < -2.0)
589 {
590 Point ak(5);
591 ak[0] = 4.52313583942e-3;
592 ak[1] = 1.20744920113e-3;
593 ak[2] = -7.89724156582e-5;
594 ak[3] = -5.04476066942e-5;
595 ak[4] = -5.35770949796e-6;
596 Point bk(5);
597 bk[0] = 1.00000000000e+0;
598 bk[1] = 9.12203410349e-1;
599 bk[2] = 4.05368773071e-1;
600 bk[3] = 9.01638932349e-2;
601 bk[4] = 9.48935714996e-3;
602 return ratfun(eta, ak, bk);
603 }
604 if (eta < 2.0)
605 {
606 Point ak(5);
607 ak[0] = 4.39937562904e-3;
608 ak[1] = 4.87225670639e-4;
609 ak[2] = -1.28470657374e-4;
610 ak[3] = 5.29110969589e-6;
611 ak[4] = 1.57166771750e-7;
612 Point bk(5);
613 bk[0] = 1.00000000000e+0;
614 bk[1] = 7.94435257415e-1;
615 bk[2] = 3.33094721709e-1;
616 bk[3] = 7.03527806143e-2;
617 bk[4] = 8.06110846078e-3;
618 return ratfun(eta, ak, bk);
619 }
620 if (eta < 10.0)
621 {
622 Point ak(5);
623 ak[0] = -1.14811912320e-3;
624 ak[1] = -1.12850923276e-1;
625 ak[2] = 1.51623048511e+0;
626 ak[3] = -2.18472031183e-1;
627 ak[4] = 7.30002451555e-2;
628 Point bk(5);
629 bk[0] = 1.00000000000e+0;
630 bk[1] = 1.42482206905e+1;
631 bk[2] = 6.97360396285e+1;
632 bk[3] = 2.18938950816e+2;
633 bk[4] = 2.77067027185e+2;
634 return ratfun(1.0 / eta, ak, bk) / (eta * eta);
635 }
636 if (eta < 100.0)
637 {
638 Point ak(5);
639 ak[0] = -1.45727889667e-4;
640 ak[1] = -2.90806748131e-1;
641 ak[2] = -1.33085045450e+1;
642 ak[3] = 1.99722374056e+2;
643 ak[4] = -1.14311378756e+1;
644 Point bk(5);
645 bk[0] = 1.00000000000e+0;
646 bk[1] = 1.39612587808e+2;
647 bk[2] = 2.18901116348e+3;
648 bk[3] = 7.11524019009e+3;
649 bk[4] = 4.55746081453e+4;
650 return ratfun(1.0 / eta, ak, bk) / (eta * eta);
651 }
652 return -log(eta) / (12.0 * eta * eta * eta);
653 } // eps3
654
lambdaeta(const Scalar eta)655 Scalar lambdaeta(const Scalar eta)
656 {
657 // lambdaeta is the positive number satisfying;
658 // eta^2/2 = lambda-1-ln(lambda);
659 // with sign(lambda-1) = sign(eta);
660 if (eta == 0.0) return 1.0;
661 const Scalar z = 1.0 + 0.5 * eta * eta;
662 return exp(-z - SpecFunc::LambertW(-exp(-z), eta < 0.0));
663 } // lambdaeta
664
invq(const Scalar x)665 Scalar invq(const Scalar x)
666 {
667 // Abramowitx & Stegun 26.2.23;
668 Scalar t = sqrt(-2 * log(x));
669 return t - (2.515517 + t * (0.802853 + t * 0.010328)) / (1.0 + t * (1.432788 + t * (0.189269 + t * 0.001308)));
670 } // invq
671
inverfc(const Scalar x)672 Scalar inverfc(const Scalar x)
673 {
674 if (x > 1.0) return -inverfc(2.0 - x);
675 const Scalar y0 = 0.70710678 * invq(0.5 * x);
676 const Scalar f = SpecFunc::ErfC(y0) - x;
677 const Scalar y02 = y0 * y0;
678 const Scalar fp = -M_2_SQRTPI * exp(-y02);
679 const Scalar c1 = -1.0 / fp;
680 const Scalar c2 = y0;
681 const Scalar c3 = (4.0 * y02 + 1.0) / 3.0;
682 const Scalar c4 = y0 * (12.0 * y02 + 7.0) / 6.0;
683 const Scalar c5 = (8.0 * y02 + 7.0) * (12.0 * y02 + 1.0) / 30.0;
684 const Scalar r = f * c1;
685 const Scalar h = r * (1.0 + r * (c2 + r * (c3 + r * (c4 + r * c5))));
686 return y0 + h;
687 } // inverfc
688
ratfun(const Scalar x,const Point & ak,const Point & bk)689 Scalar ratfun(const Scalar x,
690 const Point & ak,
691 const Point & bk)
692 {
693 return (ak[0] + x * (ak[1] + x * (ak[2] + x * (ak[3] + x * ak[4])))) / (bk[0] + x * (bk[1] + x * (bk[2] + x * (bk[3] + x * bk[4]))));
694 } // ratfun
695
696 } // GammaFunctions
697
698 END_NAMESPACE_OPENTURNS
699