dbinom.f
1.58 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
50
51
52
53
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