gamma.f
3.73 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
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
function gamma (x)
c jan 1984 edition. w. fullerton, c3, los alamos scientific lab.
dimension gcs(23)
external csevl, inits, r1mach, r9lgmc
c
data gcs ( 1) / .0085711955 90989331e0/
data gcs ( 2) / .0044153813 24841007e0/
data gcs ( 3) / .0568504368 1599363e0/
data gcs ( 4) /-.0042198353 96418561e0/
data gcs ( 5) / .0013268081 81212460e0/
data gcs ( 6) /-.0001893024 529798880e0/
data gcs ( 7) / .0000360692 532744124e0/
data gcs ( 8) /-.0000060567 619044608e0/
data gcs ( 9) / .0000010558 295463022e0/
data gcs (10) /-.0000001811 967365542e0/
data gcs (11) / .0000000311 772496471e0/
data gcs (12) /-.0000000053 542196390e0/
data gcs (13) / .0000000009 193275519e0/
data gcs (14) /-.0000000001 577941280e0/
data gcs (15) / .0000000000 270798062e0/
data gcs (16) /-.0000000000 046468186e0/
data gcs (17) / .0000000000 007973350e0/
data gcs (18) /-.0000000000 001368078e0/
data gcs (19) / .0000000000 000234731e0/
data gcs (20) /-.0000000000 000040274e0/
data gcs (21) / .0000000000 000006910e0/
data gcs (22) /-.0000000000 000001185e0/
data gcs (23) / .0000000000 000000203e0/
c
data pi /3.14159 26535 89793 24e0/
c sq2pil is alog (sqrt (2.*pi) )
data sq2pil /0.91893 85332 04672 74e0/
data ngcs, xmin, xmax, xsml, dxrel /0, 4*0.0 /
c
if (ngcs.ne.0) go to 10
c
c ---------------------------------------------------------------------
c initialize. find legal bounds for x, and determine the number of
c terms in the series required to attain an accuracy ten times better
c than machine precision.
c
ngcs = inits (gcs, 23, 0.1*r1mach(3))
c
call r9gaml (xmin, xmax)
xsml = exp (amax1 (alog (r1mach(1)), -alog(r1mach(2))) + 0.01)
dxrel = sqrt (r1mach(4))
c
c ---------------------------------------------------------------------
c finish initialization. start evaluating gamma(x).
c
10 y = abs(x)
if (y.gt.10.0) go to 50
c
c compute gamma(x) for abs(x) .le. 10.0. reduce interval and
c find gamma(1+y) for 0. .le. y .lt. 1. first of all.
c
n = x
if (x.lt.0.) n = n - 1
y = x - float(n)
n = n - 1
gamma = 0.9375 + csevl(2.*y-1., gcs, ngcs)
if (n.eq.0) return
c
if (n.gt.0) go to 30
c
c compute gamma(x) for x .lt. 1.
c
n = -n
if (x.eq.0.) call seteru (14hgamma x is 0, 14, 4, 2)
if (x.lt.0. .and. x+float(n-2).eq.0.) call seteru (
1 31hgamma x is a negative integer, 31, 4, 2)
if (x.lt.(-0.5) .and. abs((x-aint(x-0.5))/x).lt.dxrel) call
1 seteru (68hgamma answer lt half precision because x too near n
2egative integer, 68, 1, 1)
if (y.lt.xsml) call seteru (
1 54hgamma x is so close to 0.0 that the result overflows,
2 54, 5, 2)
c
do 20 i=1,n
gamma = gamma / (x+float(i-1))
20 continue
return
c
c gamma(x) for x .ge. 2.
c
30 do 40 i=1,n
gamma = (y+float(i))*gamma
40 continue
return
c
c compute gamma(x) for abs(x) .gt. 10.0. recall y = abs(x).
c
50 if (x.gt.xmax) call seteru (32hgamma x so big gamma overflows,
1 32, 3, 2)
c
gamma = 0.
if (x.lt.xmin) call seteru (35hgamma x so small gamma underflows
1 , 35, 2, 0)
if (x.lt.xmin) return
c
gamma = exp((y-0.5)*alog(y) - y + sq2pil + r9lgmc(y) )
if (x.gt.0.) return
c
if (abs((x-aint(x-0.5))/x).lt.dxrel) call seteru (
1 61hgamma answer lt half precision, x too near negative integer
2 , 61, 1, 1)
c
sinpiy = sin (pi*y)
if (sinpiy.eq.0.) call seteru (
1 31hgamma x is a negative integer, 31, 4, 2)
c
gamma = -pi / (y*sinpiy*gamma)
c
return
end