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  d1mach, d9lgmc,  dlnrel
        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