function binom (n, m)
  c june 1977 edition.   w. fullerton, c3, los alamos scientific lab.
  c
        external  alnrel, 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