伽玛函数（Gamma函数）gamma.m

function [f] = gamma1(z) % GAMMA Gamma function valid in the entire complex plane. % Accuracy is 15 significant digits along the real axis % and 13 significant digits elsewhere. % This routine uses a superb Lanczos series % approximation for the complex Gamma function. % % z may be complex and of any size. % Also n! = prod(1:n) = gamma(n+1) % %usage: [f] = gamma(z) % %tested on versions 6.0 and 5.3.1 under Sun Solaris 5.5.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-944 % W. Press, "Numerical Recipes" % S. Chang, "Computation of special functions", 1996 % W. J. Cody "An Overview of Software Development for Special % Functions", 1975 % %see also: GAMMA GAMMALN GAMMAINC PSI %see also: mhelp GAMMA % %Paul Godfrey %pgodfrey@intersil.com %http://winnie.fit.edu/~gabdo/gamma.txt %Sept 11, 2001 siz = 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 % 15 sig. digits for 0<=real(z)<=171 % coeffs should sum to about g*g/2+23/24 g=607/128; % best results when 4<=g<=5 c = [ 0.99999999999999709182; 57.156235665862923517; -59.597960355475491248; 14.136097974741747174; -0.49191381609762019978; .33994649984811888699e-4; .46523628927048575665e-4; -.98374475304879564677e-4; .15808870322491248884e-3; -.21026444172410488319e-3; .21743961811521264320e-3; -.16431810653676389022e-3; .84418223983852743293e-4; -.26190838401581408670e-4; .36899182659531622704e-5]; %Num Recipes used g=5 with 7 terms %for a less effective approximation z=z-1; zh =z+0.5; zgh=zh+g; %trick for avoiding FP overflow above z=141 zp=zgh.^(zh*0.5); ss=0.0; for pp=size(c,1)-1:-1:1 ss=ss+c(pp+1)./(z+pp); end %sqrt(2Pi) sq2pi= 2.5066282746310005024157652848110; f=(sq2pi*(c(1)+ss)).*((zp.*exp(-zgh)).*zp); f(z==0 | z==1) = 1.0; %adjust for negative real parts if ~isempty(p) f(p)=-pi./(zz(p).*f(p).*sin(pi*zz(p))); end %adjust for negative poles p=find(round(zz)==zz & imag(zz)==0 & real(zz)<=0); if ~isempty(p) f(p)=Inf; end f=reshape(f,siz); return