gamma1.f90
References to this file elsewhere.
1 ! REV PROCESSED 213 LINES OF CODE. PROGRAM DONE.
2 REAL FUNCTION GAMMA1(X)
3 !D DOUBLE PRECISION FUNCTION DGAMMA(X)
4 !----------------------------------------------------------------------
5 !
6 ! THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X.
7 ! COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN REFERENCE 1.
8 ! THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA
9 ! FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS. COEFFICIENTS
10 ! FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED.
11 ! THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM REFERENCE 2.
12 ! THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC SYSTEM, THE
13 ! COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER SELECTION OF THE
14 ! MACHINE-DEPENDENT CONSTANTS.
15 !
16 !
17 !*******************************************************************
18 !*******************************************************************
19 !
20 ! EXPLANATION OF MACHINE-DEPENDENT CONSTANTS
21 !
22 ! BETA - RADIX FOR THE FLOATING-POINT REPRESENTATION
23 ! MAXEXP - THE SMALLEST POSITIVE POWER OF BETA THAT OVERFLOWS
24 ! XBIG - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE
25 ! IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION
26 ! GAMMA(XBIG) = BETA**MAXEXP
27 ! XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER;
28 ! APPROXIMATELY BETA**MAXEXP
29 ! EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
30 ! 1.0+EPS .GT. 1.0
31 ! XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT
32 ! 1/XMININ IS MACHINE REPRESENTABLE
33 !
34 ! APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE:
35 !
36 ! BETA MAXEXP XBIG
37 !
38 ! CRAY-1 (S.P.) 2 8191 966.961
39 ! CYBER 180/855
40 ! UNDER NOS (S.P.) 2 1070 177.803
41 ! IEEE (IBM/XT,
42 ! SUN, ETC.) (S.P.) 2 128 35.040
43 ! IEEE (IBM/XT,
44 ! SUN, ETC.) (D.P.) 2 1024 171.624
45 ! IBM 3033 (D.P.) 16 63 57.574
46 ! VAX D-FORMAT (D.P.) 2 127 34.844
47 ! VAX G-FORMAT (D.P.) 2 1023 171.489
48 !
49 ! XINF EPS XMININ
50 !
51 ! CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466
52 ! CYBER 180/855
53 ! UNDER NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294
54 ! IEEE (IBM/XT,
55 ! SUN, ETC.) (S.P.) 3.40E+38 1.19E-7 1.18E-38
56 ! IEEE (IBM/XT,
57 ! SUN, ETC.) (D.P.) 1.79D+308 2.22D-16 2.23D-308
58 ! IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76
59 ! VAX D-FORMAT (D.P.) 1.70D+38 1.39D-17 5.88D-39
60 ! VAX G-FORMAT (D.P.) 8.98D+307 1.11D-16 1.12D-308
61 !
62 !*******************************************************************
63 !*******************************************************************
64 !
65 ! ERROR RETURNS
66 !
67 ! THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR
68 ! WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED
69 ! TO BE FREE OF UNDERFLOW AND OVERFLOW.
70 !
71 !
72 ! INTRINSIC FUNCTIONS REQUIRED ARE:
73 !
74 ! INT, DBLE, EXP, LOG, REAL, SIN
75 !
76 !
77 ! REFERENCES: AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL
78 ! FUNCTIONS W. J. CODY, LECTURE NOTES IN MATHEMATICS,
79 ! 506, NUMERICAL ANALYSIS DUNDEE, 1975, G. A. WATSON
80 ! (ED.), SPRINGER VERLAG, BERLIN, 1976.
81 !
82 ! COMPUTER APPROXIMATIONS, HART, ET. AL., WILEY AND
83 ! SONS, NEW YORK, 1968.
84 !
85 ! LATEST MODIFICATION: OCTOBER 12, 1989
86 !
87 ! AUTHORS: W. J. CODY AND L. STOLTZ
88 ! APPLIED MATHEMATICS DIVISION
89 ! ARGONNE NATIONAL LABORATORY
90 ! ARGONNE, IL 60439
91 !
92 !----------------------------------------------------------------------
93 INTEGER I,N
94 LOGICAL PARITY
95 !D DOUBLE PRECISION
96 REAL :: &
97 C,CONV,EPS,FACT,HALF,ONE,P,PI,Q,RES,SQRTPI,SUM,TWELVE, &
98 TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO
99 DIMENSION C(7),P(8),Q(8)
100 !----------------------------------------------------------------------
101 ! MATHEMATICAL CONSTANTS
102 !----------------------------------------------------------------------
103 DATA ONE,HALF,TWELVE,TWO,ZERO/1.0E0,0.5E0,12.0E0,2.0E0,0.0E0/, &
104 SQRTPI/0.9189385332046727417803297E0/, &
105 PI/3.1415926535897932384626434E0/
106 !D DATA ONE,HALF,TWELVE,TWO,ZERO/1.0D0,0.5D0,12.0D0,2.0D0,0.0D0/,
107 !D 1 SQRTPI/0.9189385332046727417803297D0/,
108 !D 2 PI/3.1415926535897932384626434D0/
109 !----------------------------------------------------------------------
110 ! MACHINE DEPENDENT PARAMETERS
111 !----------------------------------------------------------------------
112 DATA XBIG,XMININ,EPS/35.040E0,1.18E-38,1.19E-7/, &
113 XINF/3.4E38/
114 !D DATA XBIG,XMININ,EPS/171.624D0,2.23D-308,2.22D-16/,
115 !D 1 XINF/1.79D308/
116 !----------------------------------------------------------------------
117 ! NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX
118 ! APPROXIMATION OVER (1,2).
119 !----------------------------------------------------------------------
120 DATA P/-1.71618513886549492533811E+0,2.47656508055759199108314E+1, &
121 -3.79804256470945635097577E+2,6.29331155312818442661052E+2, &
122 8.66966202790413211295064E+2,-3.14512729688483675254357E+4, &
123 -3.61444134186911729807069E+4,6.64561438202405440627855E+4/
124 DATA Q/-3.08402300119738975254353E+1,3.15350626979604161529144E+2, &
125 -1.01515636749021914166146E+3,-3.10777167157231109440444E+3, &
126 2.25381184209801510330112E+4,4.75584627752788110767815E+3, &
127 -1.34659959864969306392456E+5,-1.15132259675553483497211E+5/
128 !D DATA P/-1.71618513886549492533811D+0,2.47656508055759199108314D+1,
129 !D 1 -3.79804256470945635097577D+2,6.29331155312818442661052D+2,
130 !D 2 8.66966202790413211295064D+2,-3.14512729688483675254357D+4,
131 !D 3 -3.61444134186911729807069D+4,6.64561438202405440627855D+4/
132 !D DATA Q/-3.08402300119738975254353D+1,3.15350626979604161529144D+2,
133 !D 1 -1.01515636749021914166146D+3,-3.10777167157231109440444D+3,
134 !D 2 2.25381184209801510330112D+4,4.75584627752788110767815D+3,
135 !D 3 -1.34659959864969306392456D+5,-1.15132259675553483497211D+5/
136 !----------------------------------------------------------------------
137 ! COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF).
138 !----------------------------------------------------------------------
139 DATA C/-1.910444077728E-03,8.4171387781295E-04, &
140 -5.952379913043012E-04,7.93650793500350248E-04, &
141 -2.777777777777681622553E-03,8.333333333333333331554247E-02, &
142 5.7083835261E-03/
143 !D DATA C/-1.910444077728D-03,8.4171387781295D-04,
144 !D 1 -5.952379913043012D-04,7.93650793500350248D-04,
145 !D 2 -2.777777777777681622553D-03,8.333333333333333331554247D-02,
146 !D 3 5.7083835261D-03/
147 !----------------------------------------------------------------------
148 ! STATEMENT FUNCTIONS FOR CONVERSION BETWEEN INTEGER AND FLOAT
149 !----------------------------------------------------------------------
150 CONV(I) = REAL(I)
151 !D CONV(I) = DBLE(I)
152 PARITY=.FALSE.
153 FACT=ONE
154 N=0
155 Y=X
156 IF(Y.LE.ZERO)THEN
157 !----------------------------------------------------------------------
158 ! ARGUMENT IS NEGATIVE
159 !----------------------------------------------------------------------
160 Y=-X
161 Y1=AINT(Y)
162 RES=Y-Y1
163 IF(RES.NE.ZERO)THEN
164 IF(Y1.NE.AINT(Y1*HALF)*TWO)PARITY=.TRUE.
165 FACT=-PI/SIN(PI*RES)
166 Y=Y+ONE
167 ELSE
168 RES=XINF
169 GOTO 900
170 ENDIF
171 ENDIF
172 !----------------------------------------------------------------------
173 ! ARGUMENT IS POSITIVE
174 !----------------------------------------------------------------------
175 IF(Y.LT.EPS)THEN
176 !----------------------------------------------------------------------
177 ! ARGUMENT .LT. EPS
178 !----------------------------------------------------------------------
179 IF(Y.GE.XMININ)THEN
180 RES=ONE/Y
181 ELSE
182 RES=XINF
183 GOTO 900
184 ENDIF
185 ELSEIF(Y.LT.TWELVE)THEN
186 Y1=Y
187 IF(Y.LT.ONE)THEN
188 !----------------------------------------------------------------------
189 ! 0.0 .LT. ARGUMENT .LT. 1.0
190 !----------------------------------------------------------------------
191 Z=Y
192 Y=Y+ONE
193 ELSE
194 !----------------------------------------------------------------------
195 ! 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY
196 !----------------------------------------------------------------------
197 N=INT(Y)-1
198 Y=Y-CONV(N)
199 Z=Y-ONE
200 ENDIF
201 !----------------------------------------------------------------------
202 ! EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0
203 !----------------------------------------------------------------------
204 XNUM=ZERO
205 XDEN=ONE
206 DO 260 I=1,8
207 XNUM=(XNUM+P(I))*Z
208 XDEN=XDEN*Z+Q(I)
209 260 CONTINUE
210 RES=XNUM/XDEN+ONE
211 IF(Y1.LT.Y)THEN
212 !----------------------------------------------------------------------
213 ! ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0
214 !----------------------------------------------------------------------
215 RES=RES/Y1
216 ELSEIF(Y1.GT.Y)THEN
217 !----------------------------------------------------------------------
218 ! ADJUST RESULT FOR CASE 2.0 .LT. ARGUMENT .LT. 12.0
219 !----------------------------------------------------------------------
220 DO 290 I=1,N
221 RES=RES*Y
222 Y=Y+ONE
223 290 CONTINUE
224 ENDIF
225 ELSE
226 !----------------------------------------------------------------------
227 ! EVALUATE FOR ARGUMENT .GE. 12.0,
228 !----------------------------------------------------------------------
229 IF(Y.LE.XBIG)THEN
230 YSQ=Y*Y
231 SUM=C(7)
232 DO 350 I=1,6
233 SUM=SUM/YSQ+C(I)
234 350 CONTINUE
235 SUM=SUM/Y-Y+SQRTPI
236 SUM=SUM+(Y-HALF)*LOG(Y)
237 RES=EXP(SUM)
238 ELSE
239 RES=XINF
240 GOTO 900
241 ENDIF
242 ENDIF
243 !----------------------------------------------------------------------
244 ! FINAL ADJUSTMENTS AND RETURN
245 !----------------------------------------------------------------------
246 IF(PARITY)RES=-RES
247 IF(FACT.NE.ONE)RES=FACT/RES
248 900 GAMMA1=RES
249 !D900 DGAMMA = RES
250 ! ---------- LAST LINE OF GAMMA ----------
251 END FUNCTION GAMMA1
252