Gerard,

I read your info about the Gamma function
and hope that the enclosed may be of some help.

The first document is a Matlab implementation
of the complex Gamma function good to 13 digits
everywhere in the complex plane.

The second document explains how to EASILY compute
the exact Lanczos coefficients with minimum roundoff error.

Regards,
    Paul [Godfrey]

function [f] = gamma(z)
% GAMMA  Gamma function valid in the entire complex plane.
%        Gives exact results for integer arguments.
%        Relative accuracy is better than 10^-13.
%        This routine uses an excellent Lanczos series
%        approximation for the complex Gamma function.
%
%usage: [f] = gamma(z)
%        z may be complex and of any size.
%        Also  n! = prod(1:n) = gamma(n+1)
%        
%tested under version 5.3.1
%
%References: C. Lanczos, SIAM JNA  1, 1964. pp. 86-96
%            Y. Luke, "The Special ... approximations", 1969 pp. 29-31
%            Y. Luke, "Algorithms ... functions", 1977
%            J. Spouge,  SIAM JNA 31, 1994. pp. 931
%            W. Press,  "Numerical Recipes"
%            S. Chang, "Computation of special functions", 1996
%
%see also:   GAMMA GAMMALN GAMMAINC PSI
%see also:   mhelp GAMMA

%Paul Godfrey
%pgodfrey@intersil.com
%11-15-00

twopi=pi+pi;
[row, col] = size(z);
z=z(:);
zz=z;

f = 0.*z; % reserve space in advance

p=find(real(z)<0);
if ~isempty(p)
   z(p)=-z(p);
end

%Lanczos approximation for the complex plane
%calculated with vpa digits(128)
c = [   1.000000000000000174663;
     5716.400188274341379136;
   -14815.30426768413909044;
    14291.49277657478554025;
    -6348.160217641458813289;
     1301.608286058321874105;
     -108.1767053514369634679;
        2.605696505611755827729;
       -0.7423452510201416151527e-2;
        0.5384136432509564062961e-7;
       -0.4023533141268236372067e-8];

%Matlab's own Gamma is based upon one by W.J. Cody from Oct. 12, 1989

% c is calculated from c=D*B*C*F
% where D is [1  -Sloane's sequence A002457],
% fact(2n+2)/(2fact(n)fact(n+1)), diagonal
% and B is rows from the odd  cols of A & Stegun table 24.1 (Binomial),
% unit Upper triangular
% and C is cols from the even cols of A & Stegun table 22.3 (Cheby),
% C(1)=1/2, Lower triangular
% and F is sqrt(2)/pi*gamma(z+0.5).*exp(z+g+0.5).*(z+g+0.5).^-(z+0.5)
% gamma(z+0.5) can be computed using factorials,2^z, and sqrt(pi).
% A & Stegun formulas 6.1.8 and 6.1.12)
%
% for z=0:(g+1) where g=9 for this example. Num. Recipes used g=5
% gamma functions were made for g=5 to g=13.
% g=9 gave the best error performance
% for n=1:171. This accuracy is comparable to Matlab's
% real only Gamma function.

%for example, for g=5 we have dimension (g+2) as
%
%D=diag([1 -1 -6 -30 -140 -630 -2772])
%
%B=[1  1  1  1  1   1   1;
%   0  1 -2  3 -4   5  -6;
%   0  0  1 -4 10 -20  35;
%   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]
%
%C=[0.5  0    0     0     0     0    0;
%   -1   2    0     0     0     0    0;
%    1  -8    8     0     0     0    0;
%   -1  18  -48    32     0     0    0;
%    1 -32  160  -256   128     0    0;
%   -1  50 -400  1120 -1280   512    0;
%    1 -72  840 -3584  6912 -6144 2048]
%
%
%F=[ 83.2488738328781364;
%    16.0123164052516814;
%     7.0235516528280364;
%     4.1065601620725384;
%     2.7864466488340058;
%     2.0690586753686773;
%     1.6309293967598249]
%
%so  c = (D*(B*C))*F
%
%this will generate the Num. Recp. values
%(if you used vpa for calculating F)
%notice that (D*B*C) always contains only integers
%(except for the 1/2 value)

%Note: 24*sum(c) ~= Integer = 12*g*g+23 if you calculated c correctly
%for this example g=5 so 24*sum(c) = 322.99... ~= 12*5*5+23 = 323

%Spouge's approximate formula for the c's can be calculated from applying
%L'Hopitals rule to the partial fraction expansion of the sum portion of
%Lanczos's formula. These values are easyer to calculate than D*B*C*F
%but are not as accurate. (An approx. of an approx.)


g=9;
%Num Recipes used g=9
%for a lower order approximation
t=z+g;
s=0;
for k=g+2:-1:2
    s=s+c(k)./t;
    t=t-1;
end
s=s+c(1);
ss=(z+g-0.5);
s=log(s.*sqrt(twopi)) + (z-0.5).*log(ss)-ss;

LogofGamma = s;
f = exp(LogofGamma);

if ~isempty(p)
   f(p)=-pi./(zz(p).*f(p).*sin(pi*zz(p)));
end

p=find(round(zz)==zz & imag(zz)==0 & real(zz)<=0);
if ~isempty(p)
   f(p)=Inf;
end

%make results exact for integer arguments
p=find(round(zz)==zz & imag(zz)==0 & real(zz)>0 );
if ~isempty(p)
   zmax=max(zz(p));
   pp=1;
   for zint=1:zmax
       p=find(zz==zint);
       if ~isempty(p)
          f(p)=pp;
       end
       pp=pp*zint;
   end
end

f=reshape(f,row,col);

return

%A demo of this routine is:
clc
clear all
close all
x=-4:1/16:4.5;
y=-4:1/16:4;
[X,Y]=meshgrid(x,y);
z=X+i*Y;
f=gamma(z);
p=find(abs(f)>10);
f(p)=10;

mesh(x,y,abs(f),phase(f));
view([-40 30]);
rotate3d;
return

%to see the relative accuracy for real n
clc
clear all
close all
warning off
format long g

n=[0.5:0.5:171.5];
n=n(:);

which gamma
lg=gamma(n);%Lanczos complex gamma
wd=pwd;
%get ready to use the other Gamma function
cd ../
which gamma
mg=gamma(n);%Matlab's gamma
cd(wd)

%let's consider Maple's Gamma to be the true and accurate standard
G=mfun('gamma',n);

elg=abs((G-lg)./G);
meanelg=mean(elg);
stdelg =svd( elg);
eelg=log10(elg);

emg=abs((G-mg)./G);
meanemg=mean(emg);
stdemg =svd( emg);
eemg=log10(emg);

figure(1)
plot(n,eelg,'bh')
hold on;
plot(n,eemg,'r*')
grid on

xlabel('n')
ylabel('Relative accuracy compared to Maple Gamma')
title('Matlab real Gamma in red, Lanczos complex Gamma in blue')

MeanErrorAndStd=[meanemg stdemg;
                 meanelg stdelg]

return

Copyright (c) 2001 by

Paul Godfrey
Intersil Corp.
Palm Bay, FL 32940
(321) 729-5359
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*(z-1)
(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 * (a-1/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*e-10) 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*a-1)!! 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(z-1)           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 pre-computed 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.7423452510201416151527e-2]
    [     0.5384136432509564062961e-7]
    [    -0.4023533141268236372067e-8]


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. 86-96
(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, 931-944.
(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), 56-70.