double precision function dfac (n)
c june 1977 edition.   w. fullerton, c3, los alamos scientific lab.
      double precision facn(31), sq2pil, x, xmax, xmin,  d9lgmc,
     1  dexp, dlog
      external d9lgmc, dexp, dlog
c
      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   /
c
      data sq2pil / 0.9189385332 0467274178 0329736405 62 d0 /
      data nmax / 0 /
c
      if (nmax.ne.0) go to 10
      call d9gaml (xmin, xmax)
      nmax = xmax - 1.d0
c
 10   if (n.lt.0) call seteru (
     1  47hdfac    factorial of negative integer undefined, 47, 1, 2)
c
      if (n.le.30) dfac = facn(n+1)
      if (n.le.30) return
c
      if (n.gt.nmax) call seteru (
     1  39hdfac    n so big factorial(n) overflows, 39, 2, 2)
c
      x = n + 1
      dfac = dexp ((x-0.5d0)*dlog(x) - x + sq2pil + d9lgmc(x) )
c
      return
      end