function GM = gm(z)
% GM 	Complete gamma function.
%
%	Y = GM(X) Returns the complete gamma function of X 
%    	NOTE: For integers: gm(n) = (n-1)! with gm(0) = inf and gm(1) = 1
%
%    	GM (with no input arguments) invokes the following example:
%	% Find the gamma function for a few values 
%	% >>yy = [0 1 2 3;0.5 1.5 2.5 3.5]
%    	  >>gy = gm(yy)


% ADSP Toolbox: Version 2.0 
% For use with "Analog and Digital Signal Processing", 2nd Ed.
% Published by PWS Publishing Co.
%
% Ashok Ambardar, EE Dept. MTU, Houghton, MI 49931, USA
% http://www.ee.mtu/faculty/akambard.html
% e-mail: akambard@mtu.edu
% Copyright (c) 1998


if nargin==0,help gm,disp('Strike a key to see results of the example')
pause,yy=[0 1 2 3;0.5 1.5 2.5 3.5],gy=gm(yy),return,end

if exist('gamma')==2,
i=find(z <= 0 & rem(z,1)== 0);
z(i)=ones(size(i));
GM=eval('gamma(z)');
GM(i)=inf*z(i);
return,end

tol=2e-16;
x=(1-z).*(z<0)+z.*(z>=0);gmb=21;
A=exp(x*log(gmb)-gmb);[m,n]=size(x);
GM=ones(m,n);for j=1:m*n
if z(j)~=1 & z(j)~=2
if (z(j) < 0 & rem(z(j),1)== 0),GM(j) = -inf;
elseif abs(z(j))>172 | z(j)==0,GM(j)=inf;
elseif (z(j)<24) & (z(j)>0) & (rem(z(j),1)==0),GM(j)=prod(2:z(j)-1);
else
GM(j)=A(j)*(ghgf(1,1,[1 1 1 1],[1+x(j) 1 1 1],gmb,tol)/x(j) +...
ghgf(2,0,[1-x(j) 1 1 1],[1 1 1 1],-1/gmb,tol)/gmb);
if z(j)<0,GM(j)=pi/(GM(j)*sin(pi*z(j)));end
end
end
end