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
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
|
DOUBLE PRECISION FUNCTION brcmp1(mu,a,b,x,y)
C-----------------------------------------------------------------------
C EVALUATION OF EXP(MU) * (X**A*Y**B/BETA(A,B))
C-----------------------------------------------------------------------
C .. Scalar Arguments ..
DOUBLE PRECISION a,b,x,y
INTEGER mu
C ..
C .. Local Scalars ..
DOUBLE PRECISION a0,apb,b0,c,const,e,h,lambda,lnx,lny,t,u,v,x0,y0,
+ z
INTEGER i,n
C ..
C .. External Functions ..
DOUBLE PRECISION algdiv,alnrel,bcorr,betaln,esum,gam1,gamln1,rlog1
EXTERNAL algdiv,alnrel,bcorr,betaln,esum,gam1,gamln1,rlog1
C ..
C .. Intrinsic Functions ..
INTRINSIC abs,dble,dlog,dmax1,dmin1,exp,sqrt
C ..
C .. Data statements ..
C-----------------
C CONST = 1/SQRT(2*PI)
C-----------------
DATA const/.398942280401433D0/
C ..
C .. Executable Statements ..
C
a0 = dmin1(a,b)
IF (a0.GE.8.0D0) GO TO 130
C
IF (x.GT.0.375D0) GO TO 10
lnx = dlog(x)
lny = alnrel(-x)
GO TO 30
10 IF (y.GT.0.375D0) GO TO 20
lnx = alnrel(-y)
lny = dlog(y)
GO TO 30
20 lnx = dlog(x)
lny = dlog(y)
C
30 z = a*lnx + b*lny
IF (a0.LT.1.0D0) GO TO 40
z = z - betaln(a,b)
brcmp1 = esum(mu,z)
RETURN
C-----------------------------------------------------------------------
C PROCEDURE FOR A .LT. 1 OR B .LT. 1
C-----------------------------------------------------------------------
40 b0 = dmax1(a,b)
IF (b0.GE.8.0D0) GO TO 120
IF (b0.GT.1.0D0) GO TO 70
C
C ALGORITHM FOR B0 .LE. 1
C
brcmp1 = esum(mu,z)
IF (brcmp1.EQ.0.0D0) RETURN
C
apb = a + b
IF (apb.GT.1.0D0) GO TO 50
z = 1.0D0 + gam1(apb)
GO TO 60
50 u = dble(a) + dble(b) - 1.D0
z = (1.0D0+gam1(u))/apb
C
60 c = (1.0D0+gam1(a))* (1.0D0+gam1(b))/z
brcmp1 = brcmp1* (a0*c)/ (1.0D0+a0/b0)
RETURN
C
C ALGORITHM FOR 1 .LT. B0 .LT. 8
C
70 u = gamln1(a0)
n = b0 - 1.0D0
IF (n.LT.1) GO TO 90
c = 1.0D0
DO 80 i = 1,n
b0 = b0 - 1.0D0
c = c* (b0/ (a0+b0))
80 CONTINUE
u = dlog(c) + u
C
90 z = z - u
b0 = b0 - 1.0D0
apb = a0 + b0
IF (apb.GT.1.0D0) GO TO 100
t = 1.0D0 + gam1(apb)
GO TO 110
100 u = dble(a0) + dble(b0) - 1.D0
t = (1.0D0+gam1(u))/apb
110 brcmp1 = a0*esum(mu,z)* (1.0D0+gam1(b0))/t
RETURN
C
C ALGORITHM FOR B0 .GE. 8
C
120 u = gamln1(a0) + algdiv(a0,b0)
brcmp1 = a0*esum(mu,z-u)
RETURN
C-----------------------------------------------------------------------
C PROCEDURE FOR A .GE. 8 AND B .GE. 8
C-----------------------------------------------------------------------
130 IF (a.GT.b) GO TO 140
h = a/b
x0 = h/ (1.0D0+h)
y0 = 1.0D0/ (1.0D0+h)
lambda = a - (a+b)*x
GO TO 150
140 h = b/a
x0 = 1.0D0/ (1.0D0+h)
y0 = h/ (1.0D0+h)
lambda = (a+b)*y - b
C
150 e = -lambda/a
IF (abs(e).GT.0.6D0) GO TO 160
u = rlog1(e)
GO TO 170
160 u = e - dlog(x/x0)
C
170 e = lambda/b
IF (abs(e).GT.0.6D0) GO TO 180
v = rlog1(e)
GO TO 190
180 v = e - dlog(y/y0)
C
190 z = esum(mu,- (a*u+b*v))
brcmp1 = const*sqrt(b*x0)*z*exp(-bcorr(a,b))
RETURN
END
|
