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