PROGRAM stirling
C------------------------------------------------------------------------
C STIRLING2.FOR
C Program to test the accuracy of Stirling's approximation.
C The formula is:
C n! = [sqrt(2*pi*n)][(n/e)**n]
C Form given by Stephenson (p.241) is
C ln(n!) ~ 0.5*ln(2*pi)-n+(n+0.5)*ln(n) 
C   
C The output is n; approx. ln(n!); ln(n!); error ratio; ln(gamma(n-1))
C------------------------------------------------------------------------
      REAL*8       DKLN10
      PARAMETER(DKLN10 = 2.30258509299D0)              ! A & S p.274
      REAL*8       DK2PI
      PARAMETER(DK2PI = 6.28318530718D0)               ! A & S p.3
      INTEGER*4    nn, jj, jline
      REAL*8       lnnfac, lnstirl, xx, rnn
      INTEGER*4    ICHAN1, OCHAN1, NLINEMAX
      PARAMETER(ICHAN1 = 11)
      PARAMETER(OCHAN1 = 12)
      PARAMETER(NLINEMAX = 100)
      CHARACTER*12 strfile1, strfile2
      REAL*8       asgamm(1:NLINEMAX), ferr, f5

      strfile1 = 'A.S.p274.txt'         ! Abramovich & Stegun, p.274.
      strfile2 = 'STIRLING.txt'
      OPEN(ICHAN1, FILE=strfile1, ERR=80, STATUS='OLD')
      DO 10 jline = 1, NLINEMAX
          ! using free-format READ seems to chop the data at each
          ! item of whitespace, so tabs are okay...
          READ(ICHAN1, *, END=12) nn, asgamm(jline)
10    CONTINUE
12    jj = 0
      CLOSE(ICHAN1)
      OPEN(OCHAN1, FILE=strfile2, ERR=80, STATUS='UNKNOWN')
      lnnfac = 0.D0
      xx = 0.D0
      ferr = 0.D0                       ! skip ferr, f5 for n = 1
      f5 = 0.D0
      WRITE(6, *) '  n       Stirling             n!'
     &    //'    frac. error        A. & S.'
      WRITE(OCHAN1, *) '  n       Stirling             n!'
     &    //'    frac. error        A. & S.'
      DO 20 nn = 1, 50
          rnn = DFLOAT(nn)
          xx = DLOG(rnn)
          lnnfac = lnnfac+xx
          lnstirl = DLOG(DK2PI)/2.D0-rnn+(rnn+0.5D0)*xx
          IF( nn .GT. 1 ) THEN
              ferr = lnstirl/lnnfac-1.D0
              f5  = DKLN10*asgamm(nn+1)
          ENDIF
          WRITE(6, 90)      nn, lnstirl, lnnfac, ferr, f5
          WRITE(OCHAN1, 90) nn, lnstirl, lnnfac, ferr, f5
          ! prettify the output...
          jj = jj+1
          IF( jj .EQ. 10 ) THEN
              jj = 0
              WRITE(6, *)
          ENDIF
20    CONTINUE
      CLOSE(OCHAN1)
80    STOP
90    FORMAT(I4, 4(F15.9))
      END