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
|
DOUBLE PRECISION FUNCTION bpser(a,b,x,eps)
C-----------------------------------------------------------------------
C POWER SERIES EXPANSION FOR EVALUATING IX(A,B) WHEN B .LE. 1
C OR B*X .LE. 0.7. EPS IS THE TOLERANCE USED.
C-----------------------------------------------------------------------
C .. Scalar Arguments ..
DOUBLE PRECISION a,b,eps,x
C ..
C .. Local Scalars ..
DOUBLE PRECISION a0,apb,b0,c,n,sum,t,tol,u,w,z
INTEGER i,m
C ..
C .. External Functions ..
DOUBLE PRECISION algdiv,betaln,gam1,gamln1
EXTERNAL algdiv,betaln,gam1,gamln1
C ..
C .. Intrinsic Functions ..
INTRINSIC abs,dble,dlog,dmax1,dmin1,exp
C ..
C .. Executable Statements ..
C
bpser = 0.0D0
IF (x.EQ.0.0D0) RETURN
C-----------------------------------------------------------------------
C COMPUTE THE FACTOR X**A/(A*BETA(A,B))
C-----------------------------------------------------------------------
a0 = dmin1(a,b)
IF (a0.LT.1.0D0) GO TO 10
z = a*dlog(x) - betaln(a,b)
bpser = exp(z)/a
GO TO 100
10 b0 = dmax1(a,b)
IF (b0.GE.8.0D0) GO TO 90
IF (b0.GT.1.0D0) GO TO 40
C
C PROCEDURE FOR A0 .LT. 1 AND B0 .LE. 1
C
bpser = x**a
IF (bpser.EQ.0.0D0) RETURN
C
apb = a + b
IF (apb.GT.1.0D0) GO TO 20
z = 1.0D0 + gam1(apb)
GO TO 30
20 u = dble(a) + dble(b) - 1.D0
z = (1.0D0+gam1(u))/apb
C
30 c = (1.0D0+gam1(a))* (1.0D0+gam1(b))/z
bpser = bpser*c* (b/apb)
GO TO 100
C
C PROCEDURE FOR A0 .LT. 1 AND 1 .LT. B0 .LT. 8
C
40 u = gamln1(a0)
m = b0 - 1.0D0
IF (m.LT.1) GO TO 60
c = 1.0D0
DO 50 i = 1,m
b0 = b0 - 1.0D0
c = c* (b0/ (a0+b0))
50 CONTINUE
u = dlog(c) + u
C
60 z = a*dlog(x) - u
b0 = b0 - 1.0D0
apb = a0 + b0
IF (apb.GT.1.0D0) GO TO 70
t = 1.0D0 + gam1(apb)
GO TO 80
70 u = dble(a0) + dble(b0) - 1.D0
t = (1.0D0+gam1(u))/apb
80 bpser = exp(z)* (a0/a)* (1.0D0+gam1(b0))/t
GO TO 100
C
C PROCEDURE FOR A0 .LT. 1 AND B0 .GE. 8
C
90 u = gamln1(a0) + algdiv(a0,b0)
z = a*dlog(x) - u
bpser = (a0/a)*exp(z)
100 IF (bpser.EQ.0.0D0 .OR. a.LE.0.1D0*eps) RETURN
C-----------------------------------------------------------------------
C COMPUTE THE SERIES
C-----------------------------------------------------------------------
sum = 0.0D0
n = 0.0D0
c = 1.0D0
tol = eps/a
110 n = n + 1.0D0
c = c* (0.5D0+ (0.5D0-b/n))*x
w = c/ (a+n)
sum = sum + w
IF (abs(w).GT.tol) GO TO 110
bpser = bpser* (1.0D0+a*sum)
RETURN
END
|
