# Algorithm to calculate lower incomplete gamma function

I'm trying to understand the implementation of the algoritm here. Please see GammaLowerRegularized(a, x) function.

I understand 1st part of the function for x <= 1 || x <= a. The authors use Kummer's confluent hypergeometric function M(1,a,x), because it converges under the specified conditions.

But don't understand the rest of the algorithm. On each iteration the following ratio is calculated:

(p2*z - p3*yc) / (q2*z - q3*yc)

Here is the code fragment:

        public static double GammaLowerRegularized(double a, double x)
{
const double epsilon = 0.000000000000001;
const double big = 4503599627370496.0;
const double bigInv = 2.22044604925031308085e-16;

double ax = (a*Math.Log(x)) - x - GammaLn(a);

// ...Skipped...

int c = 0;
double y = 1 - a;
double z = x + y + 1;

double p3 = 1;
double q3 = x;
double p2 = x + 1;
double q2 = z*x;
double ans = p2/q2;

double error;

do
{
c++;
y += 1;
z += 2;
double yc = y*c;

double p = (p2*z) - (p3*yc);
double q = (q2*z) - (q3*yc);

if (q != 0)
{
double nextans = p/q;
error = Math.Abs((ans - nextans)/nextans);
ans = nextans;
}
else
{
// zero div, skip
error = 1;
}

// shift
p3 = p2;
p2 = p;
q3 = q2;
q2 = q;

// normalize fraction when the numerator becomes large
if (Math.Abs(p) > big)
{
p3 *= bigInv;
p2 *= bigInv;
q3 *= bigInv;
q2 *= bigInv;
}
}
while (error > epsilon);

return 1d - (Math.Exp(ax)*ans);
}


Could you please explain 1) what series is calculated in this function and 2) what approximation approach is used? Or maybe you can suggest another good implementation?

The difference between lower incomplete gamma function and regularized lower incomplete gamma function doesn't matter for me, because the algorithms are almost identical.

LowerGamma(a, x) = Gamma(a) - UpperGamma(a, x)