function [f] = gammax(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
