gammaUpperRegularized function
Returns the upper incomplete regularized gamma function Q(a,x) = 1/Gamma(a) * int(exp(-t)t^(a-1),t=0..x) for real a > 0, x > 0.
Implementation
double gammaUpperRegularized(double a, double x) {
const double epsilon = 0.000000000000001;
const double big = 4503599627370496.0;
const double bigInv = 2.22044604925031308085e-16;
if (x < 1.0 || x <= a) {
return 1.0 - gammaLowerRegularized(a, x);
}
double ax = a * stdmath.log(x) - x - gammaLn(a);
if (ax < -709.78271289338399) {
return a < x ? 0.0 : 1.0;
}
ax = stdmath.exp(ax);
double t;
double y = 1 - a;
double z = x + y + 1;
double c = 0.0;
double pkm2 = 1.0;
double qkm2 = x;
double pkm1 = x + 1;
double qkm1 = z * x;
double ans = pkm1 / qkm1;
do {
c = c + 1;
y = y + 1;
z = z + 2;
double yc = y * c;
double pk = pkm1 * z - pkm2 * yc;
double qk = qkm1 * z - qkm2 * yc;
if (qk != 0) {
double r = pk / qk;
t = ((ans - r) / r).abs();
ans = r;
} else {
t = 1.0;
}
pkm2 = pkm1;
pkm1 = pk;
qkm2 = qkm1;
qkm1 = qk;
if (pk.abs() > big) {
pkm2 = pkm2 * bigInv;
pkm1 = pkm1 * bigInv;
qkm2 = qkm2 * bigInv;
qkm1 = qkm1 * bigInv;
}
} while (t > epsilon);
return ans * ax;
}