-
Notifications
You must be signed in to change notification settings - Fork 0
/
gamma.f
152 lines (146 loc) · 5.65 KB
/
gamma.f
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
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
c-----------------------------------------------------------------------
c file gamma.f
c a few useful utilities.
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c code organization.
c-----------------------------------------------------------------------
c 1. gamma.
c 2. rgr.
c 3. rgam.
c-----------------------------------------------------------------------
c subprogram 0. gamma.
c module declarations.
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c delcarations.
c-----------------------------------------------------------------------
MODULE gamma_mod
USE local_mod
IMPLICIT NONE
CONTAINS
c-----------------------------------------------------------------------
c subprogram 1. gamma.
c evaluates gamma function.
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c declarations.
c-----------------------------------------------------------------------
FUNCTION gamma(x)
REAL(r8) :: gamma,x
c-----------------------------------------------------------------------
c computations.
c-----------------------------------------------------------------------
gamma=1d0/rgr(x)
c-----------------------------------------------------------------------
c terminate.
c-----------------------------------------------------------------------
RETURN
END FUNCTION gamma
c-----------------------------------------------------------------------
c subprogram 2. rgr.
c evaluates reciprocal of gamma function.
c CACM Algorithm 80, William Holsten, UCSD.
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c declarations.
c-----------------------------------------------------------------------
FUNCTION rgr(x1)
REAL(r8) :: rgr,x1
REAL(r8) :: x,y
c-----------------------------------------------------------------------
c compute special value for x=0.
c-----------------------------------------------------------------------
x=x1
IF(x == 0.)THEN
rgr=0.
RETURN
ENDIF
c-----------------------------------------------------------------------
c compute special values for x=1.
c-----------------------------------------------------------------------
IF(x == 1.)THEN
rgr=1.
RETURN
ENDIF
c-----------------------------------------------------------------------
c for x greater than or equal to 1, reduce to the interval 0 to 1.
c-----------------------------------------------------------------------
IF(x >= 1.)THEN
y=1.
DO while(x > 1.)
x=x-1.
y=y*x
ENDDO
c-----------------------------------------------------------------------
c complete the case of x greater than or equal to 1.
c-----------------------------------------------------------------------
IF(x == 1.)THEN
rgr=1./y
ELSE
rgr=rgam(x)/y
ENDIF
RETURN
ENDIF
c-----------------------------------------------------------------------
c treat the case of x = -1.
c-----------------------------------------------------------------------
IF(x == -1.)THEN
rgr=0.
RETURN
ENDIF
c-----------------------------------------------------------------------
c treat the case of x greater than -1.
c-----------------------------------------------------------------------
IF(x > -1.)THEN
rgr=rgam(x)
RETURN
ENDIF
c-----------------------------------------------------------------------
c treat all other cases.
c-----------------------------------------------------------------------
y=x
DO while(x < -2.)
x=x+1
y=y*x
ENDDO
rgr=rgam(x)*y
c-----------------------------------------------------------------------
c terminate.
c-----------------------------------------------------------------------
RETURN
END FUNCTION rgr
c-----------------------------------------------------------------------
c subprogram 3. rgam.
c polynomial approximation.
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c declarations.
c-----------------------------------------------------------------------
FUNCTION rgam(x)
REAL(r8) :: rgam,x
INTEGER :: i
REAL(r8) :: b(0:13),z
c-----------------------------------------------------------------------
c data statements.
c-----------------------------------------------------------------------
data b/
$ 1.000000000000d0,-0.422784335092d0,-0.233093736365d0,
$ 0.191091101162d0,-0.024552490887d0,-0.017645242118d0,
$ 0.008023278113d0,-0.000804341335d0,-0.000360851496d0,
$ 0.000145624324d0,-0.000017527917d0,-0.000002625721d0,
$ 0.000001328554d0,-0.000000181220d0/
c-----------------------------------------------------------------------
c computations.
c-----------------------------------------------------------------------
z=b(13)
DO i=12,0,-1
z=z*x+b(i)
ENDDO
c-----------------------------------------------------------------------
c terminate.
c-----------------------------------------------------------------------
rgam=z*x*(x+1)
RETURN
END FUNCTION rgam
END MODULE gamma_mod