Efficient and accurate algorithms for the computation and inversion of the incomplete gamma function ratios

Algorithms for the numerical evaluation of the incomplete gamma function ratios $P(a,x)=\gamma(a,x)/\Gamma(a)$ and $Q(a,x)=\Gamma(a,x)/\Gamma(a)$ are described for positive values of $a$ and $x$. Also, inversion methods are given for solving the equations $P(a,x)=p$, $Q(a,x)=q$, with $0<p,q<1$. Both the direct computation and the inversion of the incomplete gamma function ratios are used in many problems in statistics and applied probability. The analytical approach from earlier literature is summarized and new initial estimates are derived for starting the inversion algorithms. The performance of the associated software to our algorithms (the Fortran 90 module {\bf IncgamFI}) is analyzed and compared with earlier published algorithms.


Introduction.
As it is well known, the chi-squared distribution (or its equivalent, the incomplete gamma integral) plays a key role in many applied probability problems. The incomplete gamma functions are defined by γ(a, x) = where we assume that a and x are positive. The ratios P (a, x), Q(a, x) are the standard chi-squared probability functions P (χ 2 |ν) and Q(χ 2 |ν) with parameters a = ν/2 and x = χ 2 /2.
Not only the direct computation but also the inversion of cumulative distribution functions is an important topic in statistics, probability theory, communication theory and econometrics, in particular for computing percentage points of the gamma and beta distributions, which have the chi-square, F , and Student's t-distributions as special forms. In the tails of these distributions the numerical inversion is not very easy, and for these standard distributions asymptotic formulas are available.
As an example of application in Telecom Engineering, consider a communication system with one transmit and r receive antennas operating over a flat Rayleigh fading MIMO channel. For a given communication rate Φ, the outage probability can be expressed as [6] P (1×r) out = 1 Γ(r) γ r, where S and N are the signal and noise power at the detection moment, respectively, and Γ and γ are the gamma and incomplete gamma functions. On the other hand, for a communication system with t transmit and one receive antennas one has that the outage probability is given by [6]: Then, in order to express the communication rate Φ in terms of the desired outage probability P out , one has to invert the incomplete gamma function.
In this paper, we present numerical algorithms for the computation and inversion of the incomplete gamma function ratios. A Fortran 90 version of the algorithms is made available at our website 1 .
In the numerical algorithms for the computation of the ratios, both P (a, x) and Q(a, x) are computed. First the primary function (the smaller of the two) is computed, and next the other one by using P (a, x) + Q(a, x) = 1. (1.5) In particular for large values of a, x we have a transition at a ∼ x, with P (a, x) 1 2 when a x, Q(a, x) 1 2 when a x.
(1. 6) In the next section we use more refined relations for small values of the parameters. The algorithms are partly based on [3], where also negative values of a are considered for the the pair {Γ(a, z), x −a P (a, x)}, the second element being analytic at x = 0. For applications in probability theory and mathematical statistics we prefer working with the ratios. Also, when a ∼ x (both large) the second element becomes approximately 1 2 x −a , and underflow may occur, in which case P (a, x) cannot be computed. The element Γ(a, x) may also become too large to handle.
For the large parameter case, in particular when a ∼ x, we do not use Gautschi's approach (continued fractions), but use the method given in [8]; see also [4, §8.3] and [10, §5.2]. This method is based on uniform asymptotic expansions of the incomplete gamma functions, see [7]. In [2] these uniform expansions are also used. These authors use several expansions of coefficients for the critical region a ∼ x, and we will explain that a more efficient expansion can be used.
For the numerical inversion we consider the equations for a given value of a and we compute x. For several cases we give new initial estimates for starting a numerical inversion process. The approach for large a is based on asymptotic methods developed in [9]; see also [4, §10.3.1] and [10, §6].
The resulting algorithms and software significantly improve both the accuracy and ranges of computation of the algorithm presented in [2].
2. Methods of computation. We describe the methods and the domains in the (x, a) quarter plane where they are used. First we define a function for separating the (a, x) quarter plane in two parts for assigning the primary function, that is, the function P (a, x) or Q(a, x) that has to be computed first.
As a minor modification of the function introduced in [3] we use (2.1) Then, the primary function is The function α(x) has the same asymptotic behavior as the one given in [3] for small x. The changeover point x = 1 2 is slightly better. At x = a = 1 4 (changeover point in [3]) we have P (a, x) ≃ 0.74, while at the point x = a = 1 2 we have P (a, x) ≃ 0.68.
In many representations we have to deal with the function see for example representation (2.11). For large values of a straightforward computation of Γ(a + 1) will cause overflow. When x is also large, say x ∼ a, D(a, x) may be computable, as can be seen when using Stirling's formula for the gamma function. We write D(a, x) in the form 2π/a a a e −a , a > 0, (2.5) and The quantity η arises in other cases in this paper, and when we take the square root in (2.6) we assume that for λ > 0: sign(η) = sign(λ− 1). Then, we have η ∼ λ− 1 if λ ∼ 1. The function Γ * (a) has the asymptotic expansion (Stirling series) This function is included in our software package. For testing the algorithms it is useful to use the recurrence relations P (a + 1, x) = P (a, x) − D(a, x), Q(a + 1, x) = Q(a, x) + D(a, x). (2.8) For large values of a and x we can use a scaled version by writing p(a, x) = P (a, x)/D(a, x), q(a, x) = Q(a, x)/D(a, x), x > 0, (2.9) and these functions satisfy the recursion x a + 1 p(a + 1, x) = p(a, x) − 1, x a + 1 q(a + 1, x) = q(a, x) + 1. The domains of computation are established following a compromise between efficiency and accuracy: when two methods provide the same accuracy in a certain parameter region, the selection of one method or another will depend on the efficiency of each of these methods. The recurrence relations (2.8) and (2.10) will provide numerical checks for testing the accuracy of the resulting algorithm in all regions of the (x, a)-plane.  Figure 2.1.
The expansion is where we use the Pochhammer symbol (a + 1) n = Γ(a + n + 1)/Γ(a + 1). The series converges for all a and x, and the rate of convergence improves as a/x → ∞. The terms of the series are decreasing, because we apply this expansion when a > x.
To have an idea about the number of terms needed for obtaining a certain error ε after truncating the series we write x n (a + 1) n = S n0 (a, x) + R n0 (a, x), and we compute the smallest n = n 0 that satisfies (2.14) In Table 2.1 we show the smallest number n satisfying (2.14) for ε = 10 −15 and several values of a and x. For large a we see that the number of terms needed (n) becomes constant. This can be understood by observing that the left-hand side of (2.14) becomes roughly λ n , λ = x/a, if a ≫ n. Because for large values of a we use the method of this section for P (a, x) only if λ ≤ 3/10, see Figure 2.1, it follows that not more than 30 terms are needed for a ≤ 10000.
The remainder R n0 in (2.13) can be written in the form which is again an incomplete gamma function. The function u(t) = t − (a + n 0 − 1) ln t is monotonic in (0, x], and we can integrate with respect to u, giving The expansions of this section will be used for . This domain is indicated by QT in Figure 2.1. For details and discussion we refer to [3, §4.1].
For this case we use the expansion Straightforward use of the relation Q(a, x) = 1 − P (a, x) should be avoided when a is small, and we write For the first term we have available an algorithm to compute the function g(a) in the representation The second term can be computed by using an expansion of 1 − x a = 1 − exp(a ln(x)).
in which we can use the series in (2.20) by skipping the term with n = 0.
. . . , The front factor term (x+1−a) is not causing problems, because we use the continued fraction for x ≥ a. When a = 1, 2, 3, . . . the fraction is terminating. Several algorithms are available for the numerical evaluation of continued fractions. See, for example, [4, §6.6]. Gautschi used a conversion into an infinite series of the form In Table 2 The continued fraction is an excellent alternative for the asymptotic expansion (2.29) 2.5. P (a, x), Q(a, x): uniform asymptotic expansion. The domain of computation of this section is indicated by UA in Figure 2.1. For more details on the used method we refer to [8]; see also [4, §8.3] and [10, §5.2]. We summarize the main steps for constructing an algorithm for this case.
We use the representations [7] Q(a, x) = 1 2 erfc(η a/2) + R a (η), the complementary error function. The quantity η is defined in (2.6), with again Although analytical expressions for the coefficients C n (η) are available, these representations are difficult to evaluate numerically for small values of η, that is near the transition x ∼ a. In [2] power series expansions of the coefficients C n (η) for n = 0, 1, 2, . . . , 9 for small values of η are used.
The first values of d n are To describe the algorithm, we choose a positive integer N , put β N +2 = β N +1 = 0, and compute the sequence from the recurrence relation (2.35). This recursion is stable in the backward direction. Because as an approximation for S a (η). We use the approximation in (2.40) for computing the incomplete gamma functions in IEEE double precision for a ≥ 12 and |η| ≤ 1. We need the storage of 25 coefficients d n , and in the series in (2.40). For a = 12 we need 25 terms; as a increases the convergence in the algorithm improves and we need a fewer number of terms.
The value η = −1 corresponds to λ = 0.30 . . . , and the value η = 1 to λ = 2.35 . . . . In Figure 2.1 we show the area indicated by UA in the (x, a) quarter-plane where we can apply the algorithm to obtain IEEE double precision.
3. Inversion methods. Several approaches are available in the (statistical) literature for computing the inverse of cumulative distribution functions, where often a first approximation of x is constructed, based on asymptotic estimates, but this first approximation may not be reliable. Higher approximations can be obtained by numerical inversion techniques, which require evaluation of the incomplete gamma functions. This may be rather time consuming, especially when a is large.
In [2] the inversion is considered also; we are using different methods based on analytic inversion of power series and asymptotic expansions. In this way it is clear how the first steps in the inversion method are taken. We use Newton methods when reliable starting values are available.
We solve the equations for x, with a as a given positive parameter. We consider several cases, which are schematically indicated in  p, a). We assume that the user provides both p and q, which is important when min(p, q) is small.

Small values of p.
When p is small we use the series in (2.20), and write the inversion problem as where we assume that r is small. Inverting this relation, we obtain x = r + ∞ n=2 c k r k , and the first few coefficients are c 2 = 1 a + 1 , c 3 = 3a + 5 2(a + 1) 2 (a + 2) , c 4 = 8a 2 + 33a + 31 3(a + 1) 3 (a + 2)(a + 3) , c 5 = 125a 4 + 1179a 3 + 3971a 2 + 5661a + 2888 24(1 + a) 4 (a + 2) 2 (a + 3)(a + 4) . (3. 3) It appears that c k = O((a + 1) −k+1 for large values of a, and from numerical experiments by checking several values of p, q, a which were relevant for these cases, we conclude that if r < 0.2(1 + a), that is, p < (0.2(1 + a)) a /Γ(1 + a) we can obtain 4 digits accuracy in x with the coefficients shown in (3.3). This is enough for starting a Newton method for obtaining higher accuracy. This method also works when a is small, because in that case r ∼ exp( 1 a ln p) becomes small for all fixed p ∈ (0, 1).

Small values of q.
When q is small we use the asymptotic expansion in (2.29). A first approximation x 0 of x is obtained from the equation Higher approximations of x are obtained in the form (3.5) These coefficients, as well as those given in (3.3), are obtained by symbolic computation. Maple codes for computing these or more coefficients can be obtained at our website 2 .
This method works with rather small values of q (large values of x 0 ). When we assume that 0 < a < 10 and q < e − 1 2 a /Γ(a + 1) we obtain about 4 digits accuracy in x, which is enough for starting a Newton method for obtaining higher accuracy. For larger values of a and q, the method of §3.4 can be used.

Small values of a.
We consider the inversion of P (a, x) = p for a ∈ (0, 1). For these values of a we cannot derive expansions for obtaining a reliable starting value, unless p is small in which case we can use the result of §3.1. Also, the method of the next section cannot be used in this case.

(3.8)
Then the solution x of P (a, x) = p with a ∈ (0, 1) satisfies x l < x < x u . The same results hold for the inversion of Q(a, x) = q for a ∈ (0, 1) when p is replaced with 1 − q. These bounds of x can be used for starting values for the Newton method. Of the two possibilities, the best option is x l because the Newton method necessarily converges from this starting value. The reason is that P (a, x) is an increasing function with negative second derivative; elementary graphical arguments show that if an starting value x 0 < x is chosen, then the Newton iteration x n+1 = x n − P (a, x n )/P ′ (a, x n ) produces a monotonically increasing sequence which is bounded by x, and therefore a converging sequence. The same is true for Q(a, x) because it is decreasing and with positive second derivative.
For this case, the initial approximation may be inaccurate, however convergence is certain and not very expensive, as numerical experiments show.

Large values of a.
We perform the inversion of (3.1) with respect to the parameter η by using the representations (2.30). Afterwards we have to compute λ and x from the relation for η in (2.6) and λ = x/a. We concentrate on the second equation in (3.1). For details we refer to [9]; see also [4, §10.3.1] and [10, §6].
We rewrite the inversion problem in the form 1 2 erfc η a/2 + R a (η) = q, q ∈ [0, 1], (3.9) which is equivalent to the second equation in (3.1), and we denote the solution of the above equation by η(q, a).
To start the procedure we consider R a (η) in (3.9) as a perturbation, and we define the number η 0 = η 0 (q, a) as the real number that satisfies the equation For large values of a the value η defined by (3.9) can be approximated by η 0 : η(q, a) = η 0 (q, a) + ε(η 0 , a), (3.11) and it is possible to expand as a → ∞. The coefficients ε j (η 0 , a) can be written explicitly as functions of η 0 , the first coefficient being where λ 0 follows from inverting (2.6) with η replaced by η 0 . From numerical tests it follows that we can obtain 3 or 4 significant digits when using (3.12) with 4 terms for a ≥ 1 and q ∈ (0, 1). This is enough to start a stable Newton method. Remark 1. We start the inversion of P (a, x) = p with the equation and with the first approximation η 0 (p, a) = −η 0 (q, a) (cf. (3.10)), with p + q = 1. The results for P (a, x) = p then follow from the results of this section with η 0 (q, a) replaced by −η 0 (p, a), throughout.

Complementary error function.
The inversion of this function is the first step in the large a asymptotic inversion method for the incomplete gamma function ratios. We summarize results from [11, §7.17]; for more details, see [4, §10.2].
More coefficients in the expansion can be found in [5], but they can easily be obtained by using computer algebra and formal manipulation of power series.
For small values of x we have an asymptotic expansion. Let t, α, and β be defined by t = 2 πx 2 , α = 1 ln t , β = ln(ln t). (3.16) Then we have the expansion (3.17) The first coefficients x k are given by (3.18) In §4.1 we use a higher-order Newton process to obtain better approximations. For other methods of the inversion of the error functions, we refer to [5], where coefficients of the Maclaurin expansion of y = inverf x, the inverse of x = erf y, are given, with Chebyshev coefficients for an expansion on the y-interval [−0.8, 0.8]. For small values of x (not smaller than 10 −300 ) high-precision coefficients of Chebyshev expansions are given for the numerical evaluation of y = inverfc x. For rational Chebyshev (near-minimax) approximations for y = inverfc x, we refer to [1], where xvalues are considered in the x-interval [10 −10000 , 1], with relative errors ranging down to 10 −23 . An asymptotic formula for the region x → 0 is also given.

4.
High order Newton-like methods. We describe the known method (see [4, §10.6] and [10, §6.8]) for constructing Newton-like methods of high order, with details for the inversion of the complementary error function and the incomplete gamma functions. For certain values of the parameters we will use these methods in the inversion algorithms.
The method is based on inverting the power series in h: where ζ is a zero of f and ζ 0 is an approximation. Also, f k = f (k) (ζ 0 ). We expand, assuming f (ζ 0 ) is small, and find, when f 1 = 0, When we neglect in (4.2) the coefficients c k with k ≥ 2, we obtain Newton's Rule, When f (z) satisfies a simple ordinary differential equation, the higher derivatives can be replaced by combinations of lower derivatives.

Incomplete gamma function.
When starting with an initial value x 0 > 0 for the inversion of f (x) = P (a, x) − p, the coefficients c k can be derived from the integral representation of P (a, x). This gives For much smaller values of a, this fourth order method has some convergence problems. However, as commented before, the plain Newton method converges with certainty. For values of a smaller than 0.05 the plain Newton is used in the algorithm.

Associated Software and Testing.
A Fortran 90 module (IncgamFI) implementing the algorithms can be obtained at our website 4 . The module includes the public routines incgam, for the computation of P (a, x) and Q(a, x), and invincgam, for the computation of x in the equations P (a, x) = p and Q(a, x) = q with a as a given positive parameter; p and q are also inputs of this routine. As a test of the routine incgam, we compute the maximum relative errors for the relations given in equation (2.8) using 10 6 and 10 7 random points in the following two regions of the (x, a)-plane, respectively. We obtain: 1. These errors constitute the accuracy claim of our algorithm incgam. Complementarily, we have performed tests in larger ranges by using the scaled expressions given in (2.9) and considering the relation (2.10) in the regions where the function D(a, x) (2.3) can be explicitly factored out. The use of scaled functions allows to test these methods for large values of x and a. We consider first the cases of the Taylor expansion and the continued fraction method. Using 10 7 and 10 8 random points in the following two regions of the (x, a)-plane (excluding the points where asymptotic expansions are used, see Figure 2.1), we obtain as maximum relative errors: 1. 2πa Erfc(η a/2) + S a (η) , where Erfc x = e x 2 erfc x, D(a, x) is defined in (2.3), and S a (η) denotes the series as in (2.32). In this way we can test the algorithm for the scaled function q(a, x) = Q(a, x)/D(a, x) by using the relation in (2.10). Similarly for p(a, x) = Q(a, x)/D(a, x). Using 10 7 random points in the region (0, 10 4 ] × (0, 10 4 ] of the (x, a)-plane, the maximum relative error obtained when computing the scaled functions using uniform asymptotic expansions is 4 10 −14 .
The results obtained with the scaled expressions confirm the stability and accuracy of the methods used for computing P (a, x) and Q(a, x).
For the inversion algorithm, testing is made by checking that the composition of the functions with their inverse is the identity: from the values of a and x = x in we compute p = P (a, x in ) and q = Q(a, x in ) and from the values of p, q and a we obtain x out from the inversion algorithm; we compare the values x in and x out and compute the relative errors. As a first test, we generate 10 7 random points in the region (x, a) ∈ (0, 100] × (0, 100]. In the test we exclude the points where problems related to the underflow limit in double precision arithmetic when computing the incomplete gamma function ratios appear (see Figure 5.1). The maximum relative error is 1.42 10 −11 . It is important to note that this value is obtained for a point (x, a) near the region of underflow problems: the value of the function P (a, x) at that point was ∼ 10 −297 . Some loss of accuracy is expected in these cases.
The accuracy of the initial estimates discussed in sections (3.1), (3.2) and (3.4) is illustrated in Figure 5.2, where relative distances between the initial estimate (x ini ) and the true value (x) are plotted. 10 4 random points have been considered in the plane (x, a) ∈ (0, 100) × (1, 100). As can be seen, the poorest estimate in the test is located at a relative distance less than 5 10 −3 to the real value.
The number of iterations used in the inversion algorithm is also tested. Figure  5.3 shows the number of iterations used in the region (x, a) ∈ (0, 100] × (0, 100] for computing the x values in the equations p = P (a, x) and q = Q(a, x) within an accuracy of 10 −12 − 10 −14 . As can be seen, 2 or 3 iterations are enough in most of the points of (x, a)-plane.
As a final comment, our module IncgamFI clearly improves the algorithm provided in [2], having its associated Fortran 77 routines GRATIO and GAMINV (single precision routines) both a more limited range of validity and accuracy than  our algorithms. As an illustration, Figure 5.4 shows the points where the inversion routine GAMINV fails when performing the same test as used for our inversion algorithm. A point is plotted when the relative accuracy was greater than 0.1. The routine GRATIO was used for the direct computation of the functions.