Algorithm to calculate lower incomplete gamma function

Question

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.

Answer 1

It seems that the algorithm is based on the approximation of the following continued fraction:

(The image from the Boost C++ Library documentation)

To approximate such a fraction one can use the following recursion:

The result of approximation is value of the upper incomplete gamma function. Lower incomplete gamma function can be derived using the following equation:

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