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
|
DOUBLE PRECISION FUNCTION bfrac(a,b,x,y,lambda,eps)
C-----------------------------------------------------------------------
C CONTINUED FRACTION EXPANSION FOR IX(A,B) WHEN A,B .GT. 1.
C IT IS ASSUMED THAT LAMBDA = (A + B)*Y - B.
C-----------------------------------------------------------------------
C .. Scalar Arguments ..
DOUBLE PRECISION a,b,eps,lambda,x,y
C ..
C .. Local Scalars ..
DOUBLE PRECISION alpha,an,anp1,beta,bn,bnp1,c,c0,c1,e,n,p,r,r0,s,
+ t,w,yp1
C ..
C .. External Functions ..
DOUBLE PRECISION brcomp
EXTERNAL brcomp
C ..
C .. Intrinsic Functions ..
INTRINSIC abs
C ..
C .. Executable Statements ..
C--------------------
bfrac = brcomp(a,b,x,y)
IF (bfrac.EQ.0.0D0) RETURN
C
c = 1.0D0 + lambda
c0 = b/a
c1 = 1.0D0 + 1.0D0/a
yp1 = y + 1.0D0
C
n = 0.0D0
p = 1.0D0
s = a + 1.0D0
an = 0.0D0
bn = 1.0D0
anp1 = 1.0D0
bnp1 = c/c1
r = c1/c
C
C CONTINUED FRACTION CALCULATION
C
10 n = n + 1.0D0
t = n/a
w = n* (b-n)*x
e = a/s
alpha = (p* (p+c0)*e*e)* (w*x)
e = (1.0D0+t)/ (c1+t+t)
beta = n + w/s + e* (c+n*yp1)
p = 1.0D0 + t
s = s + 2.0D0
C
C UPDATE AN, BN, ANP1, AND BNP1
C
t = alpha*an + beta*anp1
an = anp1
anp1 = t
t = alpha*bn + beta*bnp1
bn = bnp1
bnp1 = t
C
r0 = r
r = anp1/bnp1
IF (abs(r-r0).LE.eps*r) GO TO 20
C
C RESCALE AN, BN, ANP1, AND BNP1
C
an = an/bnp1
bn = bn/bnp1
anp1 = r
bnp1 = 1.0D0
GO TO 10
C
C TERMINATION
C
20 bfrac = bfrac*r
RETURN
END
|
