      FUNCTION FAC (N)
      DIMENSION FACN(26)
      SAVE FACN, SQ2PIL, NMAX
      DATA FACN( 1) / 1.0E0 /
      DATA FACN( 2) / 1.0E0 /
      DATA FACN( 3) / 2.0E0 /
      DATA FACN( 4) / 6.0E0 /
      DATA FACN( 5) / 24.0E0 /
      DATA FACN( 6) / 120.0E0 /
      DATA FACN( 7) / 720.0E0 /
      DATA FACN( 8) / 5040.0E0 /
      DATA FACN( 9) / 40320.0E0 /
      DATA FACN(10) / 362880.0E0 /
      DATA FACN(11) / 3628800.0E0 /
      DATA FACN(12) / 39916800.0E0 /
      DATA FACN(13) / 479001600.0E0 /
      DATA FACN(14) / 6227020800.0E0 /
      DATA FACN(15) / 87178291200.0E0 /
      DATA FACN(16) / 1307674368000.0E0 /
      DATA FACN(17) / 20922789888000.0E0 /
      DATA FACN(18) / 355687428096000.0E0 /
      DATA FACN(19) / 6402373705728000.0E0 /
      DATA FACN(20) /  .12164510040883200E18 /
      DATA FACN(21) /  .24329020081766400E19 /
      DATA FACN(22) /  .51090942171709440E20 /
      DATA FACN(23) /  .11240007277776077E22 /
      DATA FACN(24) /  .25852016738884977E23 /
      DATA FACN(25) /  .62044840173323944E24 /
      DATA FACN(26) /  .15511210043330986E26 /
      DATA SQ2PIL / 0.9189385332 0467274E0/
      DATA NMAX / 0 /
C***FIRST EXECUTABLE STATEMENT  FAC
      IF (NMAX.NE.0) GO TO 10
      CALL GAMLIM (XMIN, XMAX)
      NMAX = XMAX - 1.
C
 10   IF (N .LT. 0) CALL XERMSG ('SLATEC', 'FAC',
     +   'FACTORIAL OF NEGATIVE INTEGER UNDEFINED', 1, 2)
C
      IF (N.LE.25) FAC = FACN(N+1)
      IF (N.LE.25) RETURN
C
      IF (N .GT. NMAX) CALL XERMSG ('SLATEC', 'FAC',
     +   'N SO BIG FACTORIAL(N) OVERFLOWS', 2, 2)
C
      X = N + 1
      FAC = EXP ( (X-0.5)*LOG(X) - X + SQ2PIL + R9LGMC(X) )
C
      RETURN
      END
