d9gmic.f
1.88 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
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
double precision function d9gmic (a, x, alx)
c july 1977 edition. w. fullerton, c3, los alamos scientific lab.
c
c compute the complementary incomplete gamma function for a near
c a negative integer and for small x.
c
double precision a, x, alx, alng, bot, eps, euler, fk, fkp1, fm,
1 s, sgng, t, te, d1mach, dlngam, dexp, dlog
external d1mach, dlngam
data euler / 0.5772156649 0153286060 6512090082 40 d0 /
data eps, bot / 2*0.d0 /
c
if (eps.ne.0.d0) go to 10
eps = 0.5d0*d1mach(3)
bot = dlog (d1mach(1))
c
10 if (a.gt.0.d0) call seteru (
1 41hd9gmic a must be near a negative integer, 41, 2, 2)
if (x.le.0.d0) call seteru (25hd9gmic x must be gt zero, 25, 3,2)
c
m = -(a - 0.5d0)
fm = m
c
te = 1.0d0
t = 1.0d0
s = t
do 20 k=1,200
fkp1 = k + 1
te = -x*te/(fm+fkp1)
t = te/fkp1
s = s + t
if (dabs(t).lt.eps*s) go to 30
20 continue
call seteru (57hd9gmic no convergence in 200 terms of continued f
1raction, 57, 4, 2)
c
30 d9gmic = -alx - euler + x*s/(fm+1.0d0)
if (m.eq.0) return
c
if (m.eq.1) d9gmic = -d9gmic - 1.d0 + 1.d0/x
if (m.eq.1) return
c
te = fm
t = 1.d0
s = t
mm1 = m - 1
do 40 k=1,mm1
fk = k
te = -x*te/fk
t = te/(fm-fk)
s = s + t
if (dabs(t).lt.eps*dabs(s)) go to 50
40 continue
c
50 do 60 k=1,m
d9gmic = d9gmic + 1.0d0/dble(float(k))
60 continue
c
sgng = 1.0d0
if (mod(m,2).eq.1) sgng = -1.0d0
alng = dlog(d9gmic) - dlngam(fm+1.d0)
c
d9gmic = 0.d0
if (alng.gt.bot) d9gmic = sgng * dexp(alng)
if (s.ne.0.d0) d9gmic = d9gmic +
1 dsign (dexp(-fm*alx+dlog(dabs(s)/fm)), s)
c
if (d9gmic.eq.0.d0 .and. s.eq.0.d0) call seteru (
1 25hd9gmic result underflows, 25, 1, 0)
return
c
end