binom.f
1.36 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
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