% Gamma(z) computes the Gamma function using the integral representation
%            /oo  z-1  -t
% Gamma(z) = |   t    e   dt
%            /0
% The integral is evaluated using the midpoint method. It is terminated
% when the integrand drops below a given precision.
% A minimum cutoff Tmin is set in case the integrand is 0 at t=0
% If z is an integer, the factorial function is used.
% If z is negative, the recursion formula Gamma(z+1) = zGamma(z) is used
function G = Gamma(z)
Tmin = 10;    % minimum cutoff for integral
dt   = 1e-6;  % grid spacing for integral

if (z > 0)
  if (floor(z) == z)  % this means that z is an integer
    error('Gamma function undefined for non-positive integers')
  else
    G = Gamma(z)/z;  % use the recursion formula
end

G = 0;
t = 0;     % start at midpoint of first interval
if ( (t < Tmin) || (fn > prec) )
    fn = t^(z-1)*exp^-t;           % evaluate integrand
    G  = G + fn*dt;
end