double precision function dbinom (n, m)
c june 1977 edition.   w. fullerton, c3, los alamos scientific lab.
      double precision corr, fintmx, sq2pil, xk, xn, xnk, dint, d9lgmc,
     1  d1mach, dexp, dlnrel, dlog
      external alog, d1mach, d9lgmc, dexp, dint, dlnrel, dlog
      real bilnmx
c
      data sq2pil / 0.9189385332 0467274178 0329736405 62 d0 /
      data bilnmx, fintmx / 0.0, 0.0d0 /
c
      if (bilnmx.ne.0.0) go to 10
      bilnmx = dlog(d1mach(2)) - 0.0001d0
      fintmx = 0.9d0/d1mach(3)
c
 10   if (n.lt.0 .or. m.lt.0) call seteru (
     1  22hdbinom  n or m lt zero, 22, 1, 2)
      if (n.lt.m) call seteru (14hdbinom  n lt m, 14, 2, 2)
c
      k = min0 (m, n-m)
      if (k.gt.20) go to 30
      if (float(k)*alog(amax0(n,1)).gt.bilnmx) go to 30
c
      dbinom = 1.0d0
      if (k.eq.0) return
      do 20 i=1,k
        xn = n - i + 1
        xk = i
        dbinom = dbinom * (xn/xk)
 20   continue
c
      if (dbinom.lt.fintmx) dbinom = dint (dbinom+0.5d0)
      return
c
c if k.lt.9, approx is not valid and answer is close to the overflow lim
 30   if (k.lt.9) call seteru(51hdbinom  result overflows because n and/
     1or m too big, 51, 3, 2)
c
      xn = n + 1
      xk = k + 1
      xnk = n - k + 1
c
      corr = d9lgmc(xn) - d9lgmc(xk) - d9lgmc(xnk)
      dbinom = xk*dlog(xnk/xk) - xn*dlnrel(-(xk-1.0d0)/xn)
     1  -0.5d0*dlog(xn*xnk/xk) + 1.0d0 - sq2pil + corr
c
      if (dbinom.gt.dble(bilnmx)) call seteru (
     1  51hdbinom  result overflows because n and/or m too big, 51, 3,2)
c
      dbinom = dexp (dbinom)
      if (dbinom.lt.fintmx) dbinom = dint (dbinom+0.5d0)
c
      return
      end