gamic.f
2.89 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
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
function gamic (a, x)
c july 1977 edition. w. fullerton, c3, los alamos scientific lab.
c
c evaluate the complementary incomplete gamma function
c
c gamic = integral from t = x to infinity of exp(-t) * t**(a-1.) .
c
c gamic is evaluated for arbitrary real values of a and for non-negative
c values of x (even though gamic is defined for x .lt. 0.0), except that
c for x = 0 and a .le. 0.0, gamic is undefined.
c
c a slight deterioration of 2 or 3 digits accuracy will occur when
c gamic is very large or very small in absolute value, because log-
c arithmic variables are used. also, if the parameter a is very close
c to a negative integer (but not a negative integer), there is a loss
c of accuracy, which is reported if the result is less than half
c machine precision.
c
c ref. -- w. gautschi, an evaluation procedure for incomplete gamma
c functions, acm trans. math. software.
c
external alngam, r1mach, r9gmic, r9gmit, r9lgic,
1 r9lgit
data eps, sqeps, alneps, bot / 4*0.0 /
c
if (eps.ne.0.) go to 10
eps = 0.5*r1mach(3)
sqeps = sqrt(r1mach(4))
alneps = -alog(r1mach(3))
bot = alog(r1mach(1))
c
10 if (x.lt.0.0) call seteru (21hgamic x is negative, 21, 2, 2)
c
if (x.gt.0.0) go to 20
if (a.le.0.0) call seteru (
1 46hgamic x = 0 and a le 0 so gamic is undefined, 46, 3, 2)
c
gamic = exp (alngam(a+1.0) - alog(a))
return
c
20 alx = alog(x)
sga = 1.0
if (a.ne.0.0) sga = sign (1.0, a)
ma = a + 0.5*sga
aeps = a - float(ma)
c
izero = 0
if (x.ge.1.0) go to 60
c
if (a.gt.0.5 .or. abs(aeps).gt.0.001) go to 50
fm = -ma
e = 2.0
if (fm.gt.1.0) e = 2.0*(fm+2.0)/(fm*fm-1.0)
e = e - alx*x**(-0.001)
if (e*abs(aeps).gt.eps) go to 50
c
gamic = r9gmic (a, x, alx)
return
c
50 call algams (a+1.0, algap1, sgngam)
gstar = r9gmit (a, x, algap1, sgngam, alx)
if (gstar.eq.0.0) izero = 1
if (gstar.ne.0.0) alngs = alog (abs(gstar))
if (gstar.ne.0.0) sgngs = sign (1.0, gstar)
go to 70
c
60 if (a.lt.x) gamic = exp (r9lgic(a, x, alx))
if (a.lt.x) return
c
sgngam = 1.0
algap1 = alngam (a+1.0)
sgngs = 1.0
alngs = r9lgit (a, x, algap1)
c
c evaluation of gamic(a,x) in terms of tricomi-s incomplete gamma fn.
c
70 h = 1.0
if (izero.eq.1) go to 80
c
t = a*alx + alngs
if (t.gt.alneps) go to 90
if (t.gt.(-alneps)) h = 1.0 - sgngs*exp(t)
c
if (abs(h).lt.sqeps) call erroff
if (abs(h).lt.sqeps) call seteru (
1 32hgamic result lt half precision, 32, 1, 1)
c
80 sgng = sign (1.0, h) * sga * sgngam
t = alog(abs(h)) + algap1 - alog(abs(a))
if (t.lt.bot) call erroff
gamic = sgng * exp(t)
return
c
90 sgng = -sgngs * sga * sgngam
t = t + algap1 - alog(abs(a))
if (t.lt.bot) call erroff
gamic = sgng * exp(t)
return
c
end