REAL FUNCTION GAMDWC ( X )
C      
C APPROXIMATION OF THE GAMMA FUNCTION USING DAVID W. CANTRELL'S
C CONVERGENT EXPANSION POSTED IN NEWSGROUP SCI.MATH 2001-11-05
C "A new convergent expansion for the gamma function"
C
C Change history:
C 07.11.01 ALL CONSTANTS WITH 20 DECIMAL DIGITS TO ENABLE
C          CONVERSION TO DOUBLE PRECISION
C          exponentiation now performed in double precision
C 06.11.01 INITIAL VERSION
C Author: Hugo Pfoertner
C Contact Info: http://www.pfoertner.org/
C
      REAL X
C
C LOCAL VARIABLES
      REAL SQPP, C1, C2, C3, C4, CF, HSG, Z, PP5, XT
      DOUBLEPRECISION DEINV
C
C CONSTANTS
      PARAMETER ( DEINV = 0.36787944117144232160D+0 )
C SQRT(2*PI)
      PARAMETER ( SQPP = 2.5066282746310005024E+0 )
C
C COEFFICIENTS FOR CONTINUED FRACTION
      PARAMETER ( C1 = 0.71251855517070757051E+0,
     &            C2 = 5.4503922554948506451E+0,
     &            C3 = 0.33106364977586039145E+0,
     &            C4 = 2.9914210678640645176E+0 )
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
C STATEMENT FUNCTIONS
C
C CONTINUED FRACTION
      CF(Z) = 1.0 /
     &          ( C1 * Z + 1.0 /
     &                     ( C2 * Z + 1.0 /
     &                                ( C3 * Z + 1.0 / ( C4 * Z ) ) ) )
C
C HALF-SHIFT GAMMA EXPANSION,
C EXPONENTIATION PERFORMED IN DOUBLE PRECISION
      HSG (Z) = SQPP *
     &        SNGL(( DBLE(Z) * DEINV )**Z ) /
     &                   ( 1.0 + 1.0 / ( 24.0*Z - 0.5 + CF(Z) ) )
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7..
C
C Check for range
      IF ( X .GT. 34.5 ) STOP 'OVERFLOW IN GAMMA'
C
C POCHHAMMER POLYNOMIAL
      PP5 = X * (X+1.0) * (X+2.0) * (X+3.0) * (X+4.0)
C
C RECURRENCE FORMULA
      XT = X + 4.5
      GAMDWC = HSG ( XT ) / PP5
C
      RETURN
C END OF FUNCTION GAMDWC
      END