#! /usr/bin/env python
#
def alogam ( x ):

#*****************************************************************************80
#
## ALOGAM computes the logarithm of the Gamma function.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    23 December 2015
#
#  Author:
#
#    Original FORTRAN77 version by Malcolm Pike, David Hill.
#    Python version by John Burkardt.
#
#  Reference:
#
#    Malcolm Pike, David Hill,
#    Algorithm 291:
#    Logarithm of Gamma Function,
#    Communications of the ACM,
#    Volume 9, Number 9, September 1966, page 684.
#
#  Parameters:
#
#    Input, real X, the argument of the Gamma function.
#    X should be greater than 0.
#
#    Output, real ALOGAM, the logarithm of the Gamma
#    function of X.
#
  import numpy as np

  if ( x <= 0.0 ):
    value = 0.0
    return value

  y = x

  if ( x < 7.0 ):

    f = 1.0
    z = y

    while ( z < 7.0 ):
      f = f * z
      z = z + 1.0

    y = z
    f = - np.log ( f )

  else:

    f = 0.0

  z = 1.0 / y / y

  value = f + ( y - 0.5 ) * np.log ( y ) - y \
    + 0.918938533204673 + \
    ((( \
    - 0.000595238095238   * z \
    + 0.000793650793651 ) * z \
    - 0.002777777777778 ) * z \
    + 0.083333333333333 ) / y

  return value

def alogam_test ( ):

#*****************************************************************************80
#
## ALOGAM_TEST tests ALOGAM.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    23 December 2015
#
#  Author:
#
#    John Burkardt
#
  import platform
  from gamma_log_values import gamma_log_values

  print ( '' )
  print ( 'ALOGAM_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  ALOGAM evaluates the logarithm of the gamma function.' )
  print ( '' )
  print ( '      X            GAMMA_LOG(X)            GAMMA_LOG(X)' )
  print ( '      X            tabulated               computed' )
  print ( '' )

  n_data = 0

  while ( True ):

    n_data, x, fx1 = gamma_log_values ( n_data )

    if ( n_data == 0 ):
      break

    fx2 = alogam ( x )

    print ( '  %12f  %24.16g  %24.16g' % ( x, fx1, fx2 ) )
#
#  Terminate.
#
  print ( '' )
  print ( 'ALOGAM_TEST:' )
  print ( '  Normal end of execution.' )

  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  alogam_test ( )
  timestamp ( )