Copyright (c) 2001 by
Paul Godfrey
Intersil Corp.
Palm Bay, FL 32940
(321) 7295359
pgodfrey@intersil.com
A note on the computation of the convergent
Lanczos complex Gamma approximation.

Abstract
The convergent approximation of the Gamma function as developed by Lanczos
is examined. It is found that the coefficients in the series expansion
can be written as the product of 4 matrices. This allows greater accuracy
in computing the function. The coefficients of various length expansions
were computed and the best approximation was found. An 11 term expansion
was found that provides an accuracy of 13 significant digits everywhere in
the complex plane.
I. Lanczos (1) has derived a convergent approximation for the Gamma
function valid in the complex half plane Re(z)>0. That, coupled with the
reflection formula
pi*z
(z)! = 
z!*sin(pi*z)
allows the Gamma function to be determined anywhere in the complex plane.
In (1) Lanczos gives equations 32 through 35 to determine the coefficients
in the expansion. This is a convergent approximation instead of the more
widely known asymptotic approximation of Stirling (2).
Error bounds are rigously derived in (3) by Spouge and Lanczos gives error
bounds for g=1 through 5 as well. The greatest accuracy is given as
e < 2*10^10
for g=5 using a 7 term expansion.
The equations that Lanczos (1) gives;
(32) z! = sqrt(2*pi)*(z+g+1/2)^(z+1/2)*e^(z+g+1/2)*Ag(z)
1 z z*(z1)
(33) Ag(z) = p0* + p1* + p2* + ...
2 z + 1 (z+1)*(z+2)
(34) pk = sum_from a=0 to k( C(2k,2a)*Fg(a) )
(35) Fg(a) = sqrt(2)/pi * (a1/2)!*(a+g+1/2)^(a+1/2)*e^(a+g+1/2)
are subject to roundoff error if computed as written by Lanczos.
For accuracy better than the best error (2*e10) that Lanczos gives,
more thought must be given as to how to compute equations 32 through 35.
The widespread availability of programs such as Matlab, Maple, and
Mathematica allowing double precision (or better) accuracy as well as
vector operations provides motivation for a highly accurate, vectorized
computation. Henceforth, "a" will be a vector of length k+2 and g will
be a small, integer scalar. Worked examples will use g=5, the largest
used in (1) and in (4).
Equation (35) can be written as
(2*a1)!! e^(a+g+1/2)
(36) Fg(a) = sqrt(2/pi) * 
2^a (a+g+1/2)^(a+1/2)
Fg(a) will be a vector of size (k+2,1)
It is worth noting that Fg(a) is the only floating point calculation
required. For g=5 we have:
T
Fg(a) ~= [83.24 16.01 7.02 4.10 2.78 2.06 1.63] = f
Equation (34) is recognized as an inner product multiply operation or
row*vector multiply. The coefficients C(2k,2a) are the coefficients of the
Chebyshev polynomials as found in the even columns of Table 22.3 of (2).
For g=5 this becomes:
[0.5 0 0 0 0 0 0]
[1 2 0 0 0 0 0]
[ 1 8 8 0 0 0 0]
C(2k,2a) = [1 18 48 32 0 0 0] = C
[ 1 32 160 256 128 0 0]
[1 50 400 1120 1280 512 0]
[ 1 72 840 3584 6912 6144 2048]
Thus, pk = C*f
Equation (33) must first be expanded into it's constituent partial
fractions, that is:
z(z1) 2 6
 = 1 +  + 
(z+1)*(z+2) z + 1 z + 2
etc.
The partial fraction expansion Lanczos accomplishes with the aid of the
Binomial coefficients as found in Table 24.1 of (2). The resulting
expansion actually uses scaled Binomial coefficients. Accordingly, we
compute (33) as a diagonal scaling matrix multiplied by a Binomial matrix
as:
Ag(z) = Z*D*B*pk
where
1 1 1
Z = [1    ...]
z + 1 z + 2 z + 3
and
D = diag([1 1 6 30 140 630 2772])
the diagonal entries of D are, essentially, Sloane sequence A002457 (7),
which gives:
(2n+2)!
D(n,n) =  n=2,...
2*n!(n+1)!
and finally,
[1 1 1 1 1 1 1]
[0 1 2 3 4 5 6]
[0 0 1 4 10 20 35]
B = [0 0 0 1 6 21 56]
[0 0 0 0 1 8 36]
[0 0 0 0 0 1 10]
[0 0 0 0 0 0 1]
where the rows of B come from the odd columns of Table 24.1 of (2) with
alternating signs.
Ag(z) is thus found to be:
Ag(z) = Z*D*B*C*f
which has dimensions (for g=5) of (1,7)*(7,7)*(7,7)*(7,7)*(7,1)
or, in terms of coefficients
c = D*B*C*f
which need be precomputed only once.
We compute D*B*C first and find that all of it's elements are integers
with the exception of the (1,1) element which is 1/2.
[ 0.5 42 560 2688 5760 5632 2048 ]
[ 21 882 7840 28224 48384 39424 12288 ]
[ 420 23520 235200 903168 1612800 1351680 430080 ]
D*B*C = [ 2520 158760 1693440 6773760 12441600 10644480 3440640 ]
[ 6300 423360 4704000 19353600 36288000 31539200 10321920 ]
[ 6930 485100 5544000 23284800 44352000 39029760 12902400 ]
[ 2772 199584 2328480 9934848 19160064 17031168 5677056 ]
The factoring of the Lanczos equations into a product of 4 precomputed
integer matrices allows the elimination of most, but not all, of the
roundoff errors that would have occurred in computing the coefficients if
we had not done so. Even so, the final computation of the coefficients
c = (D*B*C)*f
needs to be done in quadruple precision to maintain double precision
accuracy in c. For g=5 the above will generate the coefficients found
in (1) and (4).
II. The coefficients for g=5 through 13 were computed and the resultant
Gamma functions compared. The values of z ranged from z=0.5 through z=171
and were compared with the values of Maple's Gamma function. An overall
error performance was computed as the RSS (square root of the sum of the
relative errors squared). The best performance was exhibited by g=9 with
a relative error less than 10^13. This is 1000 times more accurate than
best results shown by Lanczos. The 11 coefficients are:
c = [ 1.000000000000000174663 ]
[ 5716.400188274341379136 ]
[14815.30426768413909044 ]
[ 14291.49277657478554025 ]
[ 6348.160217641458813289 ]
[ 1301.608286058321874105 ]
[ 108.1767053514369634679 ]
[ 2.605696505611755827729 ]
[ 0.7423452510201416151527e2]
[ 0.5384136432509564062961e7]
[ 0.4023533141268236372067e8]
Further accuracy improvements would require an implementation using
better than the standard double precision computation.
It is interesting to observe that
24*sum(c) ~= 12*g^2 + 23 = Integer
This provides a quick check on c.
This is, perhaps, more useful than the pk check suggested by Lanczos.
An unproven conjecture is that this holds for all g = positive integer.
For a practical Gamma function, one could always fall back on the
definition of the factorial function to provide exact integer results
whenever z is a positive integer. If z is an integer plus 1/2 then exact
results can be computed in terms of sqrt(pi) as well.
Also, the Digamma or Psi function can also benefit from these coefficients
to provide a very accurate implementation suitable for complex
arguments. (2)
Implementations in Matlab of the Lanczos complex Gamma function and
of the complex Psi function can be found at:
http://www.mathworks.com/support/ftp/mathsv5.shtml
Finally, it is interesting to note the the less accurate coefficients
derived by Spouge in (3) can also be derived as
(z+n)*z!
Ag(z) =  z=n, n=1,2,...
sqrt(2*pi)*(z+g+1/2)^(z+1/2)*e^(z+g+1/2)
which turns into an interesting, but tedious, exercise in L'Hopital's rule.
Conclusion
The convergent expansion of the Gamma function due to Lanczos is examined.
It is found that the computation of the coefficients in the series
expansion can be computed accurately by factoring the coefficients into a
product of 4 matrices. An 11 term expansion was found that provides a
relative accuracy of 10^{13} everywhere in the complex plane.
(1) Lanczos, C. J. "A Precision Approximation of the Gamma Function"
J. SIAM Numer. Anal. Ser. B, Vol. 1 1964. pp. 8696
(2) Abramowitz and Stegun. "Handbook of Mathematical Functions"
9th edition. Dover Publications. 1972. pp. 256, 795, 828
(3) Spouge, J.L. "Computation of the gamma, digamma, and trigamma
functions", SIAM J. Numer. Anal. 31 (1994), no. 3, 931944.
(4) Press, et al. "Numerical Recipes" Cambridge University. Section 6.1
(5) Luke, Y. "The Special Functions and their Approximations" Vol. 1
Academic Press, 1969. pp. 29
(6) Zhang, S. "Computation of Special Functions" John Wiley & Sons. 1996
Ch. 3
(7) Online Encyclopedia of Integers Sequences. Sequence A002457
oeis.org.
(8) Online Inverse Symbolic Calculator
http://www.cecm.sfu.ca/projects/ISC/
(9) Ng, E.W. "A Comparison of Computational Methods and Algorithms for
the Complex Gamma Function" ACM Trans. Math. Software 1 (1975), 5670.
