From octave-sources-request at bevo dot che dot wisc dot edu Thu Mar 21 00:36:31 2002
Subject: Better Gamma implementation
From: Timothy Reluga <treluga at amath dot washington dot edu>
To: <octave-sources at bevo dot che dot wisc dot edu>
cc: <pgodfrey at intersil dot com>
Date: Thu, 21 Mar 2002 00:35:07 -0600


	Having recently had need for accurate evaluation of the gamma
function on the complex unit circle, I was disappointed to discover
Octave's current implementation doesn't handle complex arguements.
	I'm not sure how the current implementation works.  It sounds like
it may be legacy fortran code.  As such, there was no easy resolution to
my problem there.
	In the archives, I've found some old implementations based on
Stirling's asymptotic expansion.  Unfortunately, this is a divergent
series, and unable to provide arbitrary accuracy.
	Looking for alternatives, I've found Viktor Toth's web page
<http://www.rskey.org/gamma.htm> where he discusses using a Lanczos
Approximation in the complex domain.  Even better, Paul Godfrey has
explored the accuracy of the Lanczos approximation and already implemented
it in an octave-compatible fashion. See
<http://winnie.fit.edu/~gabdo/paulbio.html>.  I've emailed Mr. Godfrey,
and he has agreed to allow his implementation to be included with Octave,
in the maintainers are so inclined. If there is anything I can do to
facilitate it's incorporation into the next development release, let me
know.

						Tim Reluga

-------------------begin Gamma.m---------------

function [f] = Gamma(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 at intersil dot 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

%Lanczos approximation for the complex plane
%calculated using vpa digits(256)
%the best set of coeffs was selected from
%a solution space of g=0 to 32 with 1 to 32 terms
%these coeffs really give superb performance
%of 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