function binom (n, m)
c june 1977 edition.   w. fullerton, c3, los alamos scientific lab.
c
      external aint, alnrel, alog, exp, r1mach, r9lgmc
      data sq2pil / 0.9189385332 0467274e0 /
      data bilnmx, fintmx / 0.0, 0.0 /
c
      if (bilnmx.ne.0.0) go to 10
      bilnmx = alog (r1mach(2))
      fintmx = 0.9/r1mach(3)
c
 10   if (n.lt.0 .or. m.lt.0) call seteru (
     1  22hbinom   n or m lt zero, 22, 1, 2)
      if (n.lt.m) call seteru (14hbinom   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
      binom = 1.
      if (k.eq.0) return
c
      do 20 i=1,k
        binom = binom * float(n-i+1)/float(i)
 20   continue
c
      if (binom.lt.fintmx) binom = aint (binom+0.5)
      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(51hbinom   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 = r9lgmc(xn) - r9lgmc(xk) - r9lgmc(xnk)
      binom = xk*alog(xnk/xk) - xn*alnrel(-(xk-1.)/xn)
     1  - 0.5*alog(xn*xnk/xk) + 1.0 - sq2pil + corr
c
      if (binom.gt.bilnmx) call seteru (
     1  51hbinom   result overflows because n and/or m too big, 51, 3,2)
c
      binom = exp (binom)
      if (binom.lt.fintmx) binom = aint (binom+0.5)
c
      return
      end