function [ value, ier ] = lngamma ( z )

%*****************************************************************************80
%
%% lngamma() computes Log(Gamma(X)) using a Lanczos approximation.
%
%  Discussion:
%
%    This algorithm is not part of the Applied Statistics algorithms.
%    It is slower but gives 14 or more significant decimal digits
%    accuracy, except around X = 1 and X = 2.   The Lanczos series from
%    which this algorithm is derived is interesting in that it is a
%    convergent series approximation for the gamma function, whereas
%    the familiar series due to De Moivre (and usually wrongly called
%    the Stirling approximation) is only an asymptotic approximation, as
%    is the true and preferable approximation due to Stirling.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 January 2008
%
%  Author:
%
%    Original FORTRAN77 version by Alan Miller.
%    MATLAB  version by John Burkardt.
%
%  Reference:
%
%    Cornelius Lanczos,
%    A precision approximation of the gamma function,
%    SIAM Journal on Numerical Analysis, B,
%    Volume 1, 1964, pages 86-96.
%
%  Input:
%
%    real Z, the argument of the Gamma function.
%
%  Output:
%
%    real VALUE, the logarithm of the gamma function of Z.
%
%    integer IER, error flag.
%    0, no error occurred.
%    1, Z is less than or equal to 0.
%
  a = [ ...
         0.9999999999995183, ...
       676.5203681218835, ...
    - 1259.139216722289, ...
       771.3234287757674, ...
     - 176.6150291498386, ...
        12.50734324009056, ...
       - 0.1385710331296526, ...
         0.9934937113930748E-05, ...
         0.1659470187408462E-06 ];

  lnsqrt2pi = 0.9189385332046727;

  if ( z <= 0.0 )
    ier = 1;
    value = 0.0;
    return
  end

  ier = 0;

  value = 0.0;
  tmp = z + 7.0;

  for j = 9 : -1 : 2
    value = value + a(j) / tmp;
    tmp = tmp - 1.0;
  end

  value = value + a(1);
  value = log ( value ) + lnsqrt2pi - ( z + 6.5 ) + ...
    ( z - 0.5 ) * log ( z + 6.5 );

  return
end

