*DECK DFAC
      DOUBLE PRECISION FUNCTION DFAC (N)
C***BEGIN PROLOGUE  DFAC
C***PURPOSE  Compute the factorial function.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C1
C***TYPE      DOUBLE PRECISION (FAC-S, DFAC-D)
C***KEYWORDS  FACTORIAL, FNLIB, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C DFAC(N) calculates the double precision factorial for integer
C argument N.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  D9LGMC, DGAMLM, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C***END PROLOGUE  DFAC
      DOUBLE PRECISION FACN(31), SQ2PIL, X, XMAX, XMIN,  D9LGMC
      SAVE FACN, SQ2PIL, NMAX
      DATA FACN  (  1) / +.1000000000 0000000000 0000000000 000 D+1    /
      DATA FACN  (  2) / +.1000000000 0000000000 0000000000 000 D+1    /
      DATA FACN  (  3) / +.2000000000 0000000000 0000000000 000 D+1    /
      DATA FACN  (  4) / +.6000000000 0000000000 0000000000 000 D+1    /
      DATA FACN  (  5) / +.2400000000 0000000000 0000000000 000 D+2    /
      DATA FACN  (  6) / +.1200000000 0000000000 0000000000 000 D+3    /
      DATA FACN  (  7) / +.7200000000 0000000000 0000000000 000 D+3    /
      DATA FACN  (  8) / +.5040000000 0000000000 0000000000 000 D+4    /
      DATA FACN  (  9) / +.4032000000 0000000000 0000000000 000 D+5    /
      DATA FACN  ( 10) / +.3628800000 0000000000 0000000000 000 D+6    /
      DATA FACN  ( 11) / +.3628800000 0000000000 0000000000 000 D+7    /
      DATA FACN  ( 12) / +.3991680000 0000000000 0000000000 000 D+8    /
      DATA FACN  ( 13) / +.4790016000 0000000000 0000000000 000 D+9    /
      DATA FACN  ( 14) / +.6227020800 0000000000 0000000000 000 D+10   /
      DATA FACN  ( 15) / +.8717829120 0000000000 0000000000 000 D+11   /
      DATA FACN  ( 16) / +.1307674368 0000000000 0000000000 000 D+13   /
      DATA FACN  ( 17) / +.2092278988 8000000000 0000000000 000 D+14   /
      DATA FACN  ( 18) / +.3556874280 9600000000 0000000000 000 D+15   /
      DATA FACN  ( 19) / +.6402373705 7280000000 0000000000 000 D+16   /
      DATA FACN  ( 20) / +.1216451004 0883200000 0000000000 000 D+18   /
      DATA FACN  ( 21) / +.2432902008 1766400000 0000000000 000 D+19   /
      DATA FACN  ( 22) / +.5109094217 1709440000 0000000000 000 D+20   /
      DATA FACN  ( 23) / +.1124000727 7776076800 0000000000 000 D+22   /
      DATA FACN  ( 24) / +.2585201673 8884976640 0000000000 000 D+23   /
      DATA FACN  ( 25) / +.6204484017 3323943936 0000000000 000 D+24   /
      DATA FACN  ( 26) / +.1551121004 3330985984 0000000000 000 D+26   /
      DATA FACN  ( 27) / +.4032914611 2660563558 4000000000 000 D+27   /
      DATA FACN  ( 28) / +.1088886945 0418352160 7680000000 000 D+29   /
      DATA FACN  ( 29) / +.3048883446 1171386050 1504000000 000 D+30   /
      DATA FACN  ( 30) / +.8841761993 7397019545 4361600000 000 D+31   /
      DATA FACN  ( 31) / +.2652528598 1219105863 6308480000 000 D+33   /
      DATA SQ2PIL / 0.9189385332 0467274178 0329736405 62 D0 /
      DATA NMAX / 0 /
C***FIRST EXECUTABLE STATEMENT  DFAC
      IF (NMAX.NE.0) GO TO 10
      CALL DGAMLM (XMIN, XMAX)
      NMAX = XMAX - 1.D0
C
 10   IF (N .LT. 0) CALL XERMSG ('SLATEC', 'DFAC',
     +   'FACTORIAL OF NEGATIVE INTEGER UNDEFINED', 1, 2)
C
      IF (N.LE.30) DFAC = FACN(N+1)
      IF (N.LE.30) RETURN
C
      IF (N .GT. NMAX) CALL XERMSG ('SLATEC', 'DFAC',
     +   'N SO BIG FACTORIAL(N) OVERFLOWS', 2, 2)
C
      X = N + 1
      DFAC = EXP ((X-0.5D0)*LOG(X) - X + SQ2PIL + D9LGMC(X) )
C
      RETURN
      END