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
137
|
DOUBLE PRECISION FUNCTION brcomp(a,b,x,y)
C-----------------------------------------------------------------------
C EVALUATION OF X**A*Y**B/BETA(A,B)
C-----------------------------------------------------------------------
C .. Scalar Arguments ..
DOUBLE PRECISION a,b,x,y
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,gam1,gamln1,rlog1
EXTERNAL algdiv,alnrel,bcorr,betaln,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
brcomp = 0.0D0
IF (x.EQ.0.0D0 .OR. y.EQ.0.0D0) RETURN
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)
brcomp = exp(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
brcomp = exp(z)
IF (brcomp.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
brcomp = brcomp* (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 brcomp = a0*exp(z)* (1.0D0+gam1(b0))/t
RETURN
C
C ALGORITHM FOR B0 .GE. 8
C
120 u = gamln1(a0) + algdiv(a0,b0)
brcomp = a0*exp(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 = exp(- (a*u+b*v))
brcomp = const*sqrt(b*x0)*z*exp(-bcorr(a,b))
RETURN
END
|
