home | index | units | counting | geometry | algebra | trigonometry & functions | calculus analysis | sets & logic | number theory | recreational | misc | nomenclature & history | physics

# Implementing the Gamma Function

Copyright © 2000-2001 Paul Godfrey : The text below is material copyrighted by Paul Godfrey, privately communicated to us and reproduced here by permission.

## Lanczos Implementation of the Gamma Function

```Gerard,

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]```
2001-06-07
 ```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. ```Copyright © 2000-2001, Paul Godfrey