function alngam ( xvalue, ifault ) !*****************************************************************************80 ! !! alngam() computes the logarithm of the gamma function. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 13 January 2008 ! ! Author: ! ! Original FORTRAN77 version by Allan Macleod. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Allan Macleod, ! Algorithm AS 245, ! A Robust and Reliable Algorithm for the Logarithm of the Gamma Function, ! Applied Statistics, ! Volume 38, Number 2, 1989, pages 397-402. ! ! Parameters: ! ! Input, real ( kind = rk ) XVALUE, the argument of the Gamma function. ! ! Output, integer IFAULT, error flag. ! 0, no error occurred. ! 1, XVALUE is less than or equal to 0. ! 2, XVALUE is too big. ! ! Output, real ( kind = rk ) ALNGAM, the logarithm of the gamma function of X. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) alngam real ( kind = rk ), parameter :: alr2pi = 0.918938533204673D+00 integer ifault real ( kind = rk ), dimension ( 9 ) :: r1 = (/ & -2.66685511495D+00, & -24.4387534237D+00, & -21.9698958928D+00, & 11.1667541262D+00, & 3.13060547623D+00, & 0.607771387771D+00, & 11.9400905721D+00, & 31.4690115749D+00, & 15.2346874070D+00 /) real ( kind = rk ), dimension ( 9 ) :: r2 = (/ & -78.3359299449D+00, & -142.046296688D+00, & 137.519416416D+00, & 78.6994924154D+00, & 4.16438922228D+00, & 47.0668766060D+00, & 313.399215894D+00, & 263.505074721D+00, & 43.3400022514D+00 /) real ( kind = rk ), dimension ( 9 ) :: r3 = (/ & -2.12159572323D+05, & 2.30661510616D+05, & 2.74647644705D+04, & -4.02621119975D+04, & -2.29660729780D+03, & -1.16328495004D+05, & -1.46025937511D+05, & -2.42357409629D+04, & -5.70691009324D+02 /) real ( kind = rk ), dimension ( 5 ) :: r4 = (/ & 0.279195317918525D+00, & 0.4917317610505968D+00, & 0.0692910599291889D+00, & 3.350343815022304D+00, & 6.012459259764103D+00 /) real ( kind = rk ) x real ( kind = rk ) x1 real ( kind = rk ) x2 real ( kind = rk ), parameter :: xlge = 5.10D+05 real ( kind = rk ), parameter :: xlgst = 1.0D+30 real ( kind = rk ) xvalue real ( kind = rk ) y x = xvalue alngam = 0.0D+00 ! ! Check the input. ! if ( xlgst <= x ) then ifault = 2 return end if if ( x <= 0.0D+00 ) then ifault = 1 return end if ifault = 0 ! ! Calculation for 0 < X < 0.5 and 0.5 <= X < 1.5 combined. ! if ( x < 1.5D+00 ) then if ( x < 0.5D+00 ) then alngam = - log ( x ) y = x + 1.0D+00 ! ! Test whether X < machine epsilon. ! if ( y == 1.0D+00 ) then return end if else alngam = 0.0D+00 y = x x = ( x - 0.5D+00 ) - 0.5D+00 end if alngam = alngam + x * (((( & r1(5) * y & + r1(4) ) * y & + r1(3) ) * y & + r1(2) ) * y & + r1(1) ) / (((( & y & + r1(9) ) * y & + r1(8) ) * y & + r1(7) ) * y & + r1(6) ) return end if ! ! Calculation for 1.5 <= X < 4.0. ! if ( x < 4.0D+00 ) then y = ( x - 1.0D+00 ) - 1.0D+00 alngam = y * (((( & r2(5) * x & + r2(4) ) * x & + r2(3) ) * x & + r2(2) ) * x & + r2(1) ) / (((( & x & + r2(9) ) * x & + r2(8) ) * x & + r2(7) ) * x & + r2(6) ) ! ! Calculation for 4.0 <= X < 12.0. ! else if ( x < 12.0D+00 ) then alngam = (((( & r3(5) * x & + r3(4) ) * x & + r3(3) ) * x & + r3(2) ) * x & + r3(1) ) / (((( & x & + r3(9) ) * x & + r3(8) ) * x & + r3(7) ) * x & + r3(6) ) ! ! Calculation for 12.0 <= X. ! else y = log ( x ) alngam = x * ( y - 1.0D+00 ) - 0.5D+00 * y + alr2pi if ( x <= xlge ) then x1 = 1.0D+00 / x x2 = x1 * x1 alngam = alngam + x1 * ( ( & r4(3) * & x2 + r4(2) ) * & x2 + r4(1) ) / ( ( & x2 + r4(5) ) * & x2 + r4(4) ) end if end if return end subroutine gamma_log_values ( n_data, x, fx ) !*****************************************************************************80 ! !! gamma_log_values() returns some values of the Log Gamma function. ! ! Discussion: ! ! In Mathematica, the function can be evaluated by: ! ! Log[Gamma[x]] ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 14 August 2004 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Milton Abramowitz, Irene Stegun, ! Handbook of Mathematical Functions, ! National Bureau of Standards, 1964, ! ISBN: 0-486-61272-4, ! LC: QA47.A34. ! ! Stephen Wolfram, ! The Mathematica Book, ! Fourth Edition, ! Cambridge University Press, 1999, ! ISBN: 0-521-64314-7, ! LC: QA76.95.W65. ! ! Parameters: ! ! Input/output, integer N_DATA. The user sets N_DATA to 0 ! before the first call. On each call, the routine increments N_DATA by 1, ! and returns the corresponding data; when there is no more data, the ! output value of N_DATA will be 0 again. ! ! Output, real ( kind = rk ) X, the argument of the function. ! ! Output, real ( kind = rk ) FX, the value of the function. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: n_max = 20 real ( kind = rk ) fx real ( kind = rk ), save, dimension ( n_max ) :: fx_vec = (/ & 0.1524063822430784D+01, & 0.7966778177017837D+00, & 0.3982338580692348D+00, & 0.1520596783998375D+00, & 0.0000000000000000D+00, & -0.4987244125983972D-01, & -0.8537409000331584D-01, & -0.1081748095078604D+00, & -0.1196129141723712D+00, & -0.1207822376352452D+00, & -0.1125917656967557D+00, & -0.9580769740706586D-01, & -0.7108387291437216D-01, & -0.3898427592308333D-01, & 0.00000000000000000D+00, & 0.69314718055994530D+00, & 0.17917594692280550D+01, & 0.12801827480081469D+02, & 0.39339884187199494D+02, & 0.71257038967168009D+02 /) integer n_data real ( kind = rk ) x real ( kind = rk ), save, dimension ( n_max ) :: x_vec = (/ & 0.20D+00, & 0.40D+00, & 0.60D+00, & 0.80D+00, & 1.00D+00, & 1.10D+00, & 1.20D+00, & 1.30D+00, & 1.40D+00, & 1.50D+00, & 1.60D+00, & 1.70D+00, & 1.80D+00, & 1.90D+00, & 2.00D+00, & 3.00D+00, & 4.00D+00, & 10.00D+00, & 20.00D+00, & 30.00D+00 /) if ( n_data < 0 ) then n_data = 0 end if n_data = n_data + 1 if ( n_max < n_data ) then n_data = 0 x = 0.0D+00 fx = 0.0D+00 else x = x_vec(n_data) fx = fx_vec(n_data) end if return end function lngamma ( z, ier ) !*****************************************************************************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 ! Stirling's 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: ! ! 13 January 2008 ! ! Author: ! ! Original FORTRAN77 version by Alan Miller. ! FORTRAN90 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. ! ! Parameters: ! ! Input, real ( kind = rk ) Z, the argument of the Gamma function. ! ! Output, integer IER, error flag. ! 0, no error occurred. ! 1, Z is less than or equal to 0. ! ! Output, real ( kind = rk ) LNGAMMA, the logarithm of the gamma function ! of Z. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ), dimension ( 9 ) :: a = (/ & 0.9999999999995183D+00, & 676.5203681218835D+00, & - 1259.139216722289D+00, & 771.3234287757674D+00, & - 176.6150291498386D+00, & 12.50734324009056D+00, & - 0.1385710331296526D+00, & 0.9934937113930748D-05, & 0.1659470187408462D-06 /) integer ier integer j real ( kind = rk ) lngamma real ( kind = rk ), parameter :: lnsqrt2pi = 0.9189385332046727D+00 real ( kind = rk ) tmp real ( kind = rk ) z if ( z <= 0.0D+00 ) then ier = 1 lngamma = 0.0D+00 return end if ier = 0 lngamma = 0.0D+00 tmp = z + 7.0D+00 do j = 9, 2, -1 lngamma = lngamma + a(j) / tmp tmp = tmp - 1.0D+00 end do lngamma = lngamma + a(1) lngamma = log ( lngamma ) + lnsqrt2pi - ( z + 6.5D+00 ) + & ( z - 0.5D+00 ) * log ( z + 6.5D+00 ) return end