gammapInv function
Implementation
double gammapInv(num p, num a) {
const epsilon = 1.0e-8;
final a1 = a - 1.0;
final gln = gammaLn(a);
var x = 0.0;
var afac = 0.0;
var lna1 = 0.0;
if (p >= 1.0) {
return max(100, a + 100 * sqrt(a));
} else if (p <= 0) {
return 0.0;
} else if (a > 1.0) {
lna1 = log(a1);
afac = exp(a1 * (lna1 - 1.0) - gln);
final pp = p < 0.5 ? p : 1.0 - p;
final t = sqrt(-2 * log(pp));
x = (2.30753 + t * 0.27061) / (1.0 + t * (0.99229 + t * 0.04481)) - t;
if (p < 0.5) {
x = -x;
}
x = max(1.0e-3,
a * pow(1.0 - 1.0 / (9.0 * a) - x / (3.0 * sqrt(a)), 3).toDouble());
} else {
final t = 1.0 - a * (0.253 + a * 0.12);
if (p < t) {
x = pow(p / t, 1.0 / a).toDouble();
} else {
x = 1.0 - log(1.0 - (p - t) / (1.0 - t));
}
}
for (var j = 0; j < 12; j++) {
if (x <= 0.0) {
return 0.0;
}
final err = lowRegGamma(a, x) - p;
var t = a > 1.0
? afac * exp(-(x - a1) + a1 * (log(x) - lna1))
: exp(-x + a1 * log(x) - gln);
final u = err / t;
x -= (t = u / (1.0 - 0.5 * min(1.0, u * ((a - 1.0) / x - 1.0))));
if (x <= 0.0) {
x = 0.5 * (x + t);
}
if (t.abs() < epsilon * x) {
break;
}
}
return x;
}