From octave-sources-request at bevo dot che dot wisc dot edu Thu Mar 21 01:52:30 2002
Subject: Re: Better Gamma implementation
From: Paul Kienzle <pkienzle at jazz dot ncnr dot nist dot gov>
To: Timothy Reluga <treluga at amath dot washington dot edu>,    octave-sources@bevo.che.wisc.edu
Cc: pgodfrey at intersil dot com
Date: Thu, 21 Mar 2002 02:52:22 -0500

See also:

	http://www.octave.org/mailing-lists/octave-sources/1996/3

You can add it to octave-forge (http://octave.sf.net) if John does
not wish to include it in the main octave distribution.

Paul Kienzle
pkienzle at users dot sf dot net

On Thu, Mar 21, 2002 at 12:35:07AM -0600, Timothy Reluga wrote:
> 
> 	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
>