!     Last change:  DOS  29 Jul 2000   12:20 pm
!                                       *** copyright 2000 ***
! *** filename slgamk.f95               *** John F. Monahan **
!                                       **********************
      program pslgamk
!  test of slgamk -- log( k ! )
      implicit none
      real(kind=8) dlgk,fact,dlgfac,slgfac
      real(kind=8) slgamk
      integer k
!               
  20  format('  Test of SLGAMK'/3x,'k',8x,'k!',8x,'log(k!)',8x,'SLGAMK')
  21  format(i4,f14.2,4x,2f14.8)
  22  format(/'   Second Test'                              &
     &   /3x,'k',6x,'log(k)',4x,'SLGAMK(k)-SLGAMK(k-1)')
  23  format(i4,4x,2f14.8)
!                   
      open( unit=6, file='slgamk.out' )
!                   for small values, give the comparison
      write(6,20)
      fact = 1.d0
      do k = 1,10
      fact = fact*real(k,8)             !KIND
      dlgfac = log(fact)
      slgfac = slgamk(k)
      write(6,21) k,fact,dlgfac,slgfac
      end do
!                   for larger values, compare difference
      write(6,22)
      do k = 8,80,4
      dlgk = log(real(k,8))             !KIND
      dlgfac = slgamk(k) - slgamk(k-1)
      write(6,23) k,dlgk,dlgfac
      end do
      stop
      end program pslgamk

      REAL(KIND=8) FUNCTION SLGAMK(K)
!  PRODUCES NATURAL LOGARITHM OF K!                                     
!  FOR SMALL VALUES OF K -- USES TABLE
!  FOR LARGER VALUES -- USES STIRLING'S APPROXIMATION
!
!  J F MONAHAN (SEPT, 1990) DEPT OF STATISTICS, NCSU, RALEIGH, NC USA   
!      RECODED OCTOBER 1999, APRIL 2000 FOR FORTRAN 95
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: K
      REAL(KIND=8), DIMENSION(36) :: TAB
      REAL(KIND=8) F,DLF
      REAL(KIND=8), PARAMETER :: CHL2PI = .918938533205D0 ! log(2pi)/2
!
      DATA TAB/   0.D0, 0.6931471806D0, 1.7917594692D0, 3.1780538303D0&
     &, 4.7874917428D0, 6.5792512120D0, 8.5251613611D0,10.6046029027D0&
     &,12.8018274801D0,15.1044125731D0,17.5023078459D0,19.9872144957D0&
     &,22.5521638531D0,25.1912211827D0,27.8992713838D0,30.6718601061D0&
     &,33.5050734501D0,36.3954452080D0,39.3398841872D0,42.3356164608D0&
     &,45.3801388985D0,48.4711813518D0,51.6066755678D0,54.7847293981D0&
     &,58.0036052230D0,61.2617017610D0,64.5575386270D0,67.8897431372D0&
     &,71.2570389672D0,74.6582363488D0,78.0922235533D0,81.5579594561D0&
     &,85.0544670176D0,88.5808275422D0,92.1361756037D0,95.7196945421D0/
!                                                                       
      SLGAMK = 0.D0
      IF( (K.EQ.0) .OR. (K.EQ.1) ) RETURN                               
      IF( K .LE. 36 ) THEN
      SLGAMK = TAB(K)                                                   
      RETURN                                                            
!                      USE STIRLING'S FORMULA                           
                      ELSE
      F = DBLE(K)
      DLF = LOG(F)
      SLGAMK = ( (F+.5D0)*DLF - F )                         &
     &          + ( CHL2PI + ( 1.D0 - 1.D0/(30.D0*F*F) )/(12.D0*F) )
                      END IF ! ( K .LE. 36 )
!
      RETURN
      END FUNCTION SLGAMK
! *** end of filename slgamk.f95        *********************