d9gmit.f
1.65 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
double precision function d9gmit (a, x, algap1, sgngam, alx)
c july 1977 edition. w. fullerton, c3, los alamos scientific lab.
c
c compute tricomi-s incomplete gamma function for small x.
c
double precision a, x, algap1, sgngam, alx, ae, aeps, algs, alg2,
1 bot, eps, fk, s, sgng2, t, te, d1mach, dlngam, dexp, dlog
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 (x.le.0.d0) call seteru (24hd9gmit x should be gt 0, 24, 1, 2)
c
ma = a + 0.5d0
if (a.lt.0.d0) ma = a - 0.5d0
aeps = a - dble(float(ma))
c
ae = a
if (a.lt.(-0.5d0)) ae = aeps
c
t = 1.d0
te = ae
s = t
do 20 k=1,200
fk = k
te = -x*te/fk
t = te/(ae+fk)
s = s + t
if (dabs(t).lt.eps*dabs(s)) go to 30
20 continue
call seteru (54hd9gmit no convergence in 200 terms of taylor-s se
1ries, 54, 2, 2)
c
30 if (a.ge.(-0.5d0)) algs = -algap1 + dlog(s)
if (a.ge.(-0.5d0)) go to 60
c
algs = -dlngam(1.d0+aeps) + dlog(s)
s = 1.0d0
m = -ma - 1
if (m.eq.0) go to 50
t = 1.0d0
do 40 k=1,m
t = x*t/(aeps-dble(float(m+1-k)))
s = s + t
if (dabs(t).lt.eps*dabs(s)) go to 50
40 continue
c
50 d9gmit = 0.0d0
algs = -dble(float(ma))*dlog(x) + algs
if (s.eq.0.d0 .or. aeps.eq.0.d0) go to 60
c
sgng2 = sgngam * dsign (1.0d0, s)
alg2 = -x - algap1 + dlog(dabs(s))
c
if (alg2.gt.bot) d9gmit = sgng2 * dexp(alg2)
if (algs.gt.bot) d9gmit = d9gmit + dexp(algs)
return
c
60 d9gmit = dexp (algs)
return
c
end