1
CS REAL FUNCTION ALGAMA(X)
2
DOUBLE PRECISION FUNCTION DLGAMA(X)
3
C----------------------------------------------------------------------
5
C This routine calculates the LOG(GAMMA) function for a positive real
6
C argument X. Computation is based on an algorithm outlined in
7
C references 1 and 2. The program uses rational functions that
8
C theoretically approximate LOG(GAMMA) to at least 18 significant
9
C decimal digits. The approximation for X > 12 is from reference
10
C 3, while approximations for X < 12.0 are similar to those in
11
C reference 1, but are unpublished. The accuracy achieved depends
12
C on the arithmetic system, the compiler, the intrinsic functions,
13
C and proper selection of the machine-dependent constants.
16
C*********************************************************************
17
C*********************************************************************
19
C Explanation of machine-dependent constants
21
C beta - radix for the floating-point representation
22
C maxexp - the smallest positive power of beta that overflows
23
C XBIG - largest argument for which LN(GAMMA(X)) is representable
24
C in the machine, i.e., the solution to the equation
25
C LN(GAMMA(XBIG)) = beta**maxexp
26
C XINF - largest machine representable floating-point number;
27
C approximately beta**maxexp.
28
C EPS - The smallest positive floating-point number such that
30
C FRTBIG - Rough estimate of the fourth root of XBIG
33
C Approximate values for some important machines are:
37
C CRAY-1 (S.P.) 2 8191 9.62E+2461
39
C under NOS (S.P.) 2 1070 1.72E+319
41
C SUN, etc.) (S.P.) 2 128 4.08E+36
43
C SUN, etc.) (D.P.) 2 1024 2.55D+305
44
C IBM 3033 (D.P.) 16 63 4.29D+73
45
C VAX D-Format (D.P.) 2 127 2.05D+36
46
C VAX G-Format (D.P.) 2 1023 1.28D+305
51
C CRAY-1 (S.P.) 5.45E+2465 7.11E-15 3.13E+615
53
C under NOS (S.P.) 1.26E+322 3.55E-15 6.44E+79
55
C SUN, etc.) (S.P.) 3.40E+38 1.19E-7 1.42E+9
57
C SUN, etc.) (D.P.) 1.79D+308 2.22D-16 2.25D+76
58
C IBM 3033 (D.P.) 7.23D+75 2.22D-16 2.56D+18
59
C VAX D-Format (D.P.) 1.70D+38 1.39D-17 1.20D+9
60
C VAX G-Format (D.P.) 8.98D+307 1.11D-16 1.89D+76
62
C**************************************************************
63
C**************************************************************
67
C The program returns the value XINF for X .LE. 0.0 or when
68
C overflow would occur. The computation is believed to
69
C be free of underflow and overflow.
72
C Intrinsic functions required are:
79
C 1) W. J. Cody and K. E. Hillstrom, 'Chebyshev Approximations for
80
C the Natural Logarithm of the Gamma Function,' Math. Comp. 21,
83
C 2) K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May,
86
C 3) Hart, Et. Al., Computer Approximations, Wiley and sons, New
90
C Authors: W. J. Cody and L. Stoltz
91
C Argonne National Laboratory
93
C Latest modification: June 16, 1988
95
C----------------------------------------------------------------------
99
1 C,CORR,D1,D2,D4,EPS,FRTBIG,FOUR,HALF,ONE,PNT68,P1,P2,P4,
100
2 Q1,Q2,Q4,RES,SQRTPI,THRHAL,TWELVE,TWO,X,XBIG,XDEN,XINF,
101
3 XM1,XM2,XM4,XNUM,Y,YSQ,ZERO
102
DIMENSION C(7),P1(8),P2(8),P4(8),Q1(8),Q2(8),Q4(8)
103
C----------------------------------------------------------------------
104
C Mathematical constants
105
C----------------------------------------------------------------------
106
CS DATA ONE,HALF,TWELVE,ZERO/1.0E0,0.5E0,12.0E0,0.0E0/,
107
CS 1 FOUR,THRHAL,TWO,PNT68/4.0E0,1.5E0,2.0E0,0.6796875E0/,
108
CS 2 SQRTPI/0.9189385332046727417803297E0/
109
DATA ONE,HALF,TWELVE,ZERO/1.0D0,0.5D0,12.0D0,0.0D0/,
110
1 FOUR,THRHAL,TWO,PNT68/4.0D0,1.5D0,2.0D0,0.6796875D0/,
111
2 SQRTPI/0.9189385332046727417803297D0/
112
C----------------------------------------------------------------------
113
C Machine dependent parameters
114
C----------------------------------------------------------------------
115
CS DATA XBIG,XINF,EPS,FRTBIG/4.08E36,3.401E38,1.19E-7,1.42E9/
116
DATA XBIG,XINF,EPS,FRTBIG/2.55D305,1.79D308,2.22D-16,2.25D76/
117
C----------------------------------------------------------------------
118
C Numerator and denominator coefficients for rational minimax
119
C approximation over (0.5,1.5).
120
C----------------------------------------------------------------------
121
CS DATA D1/-5.772156649015328605195174E-1/
122
CS DATA P1/4.945235359296727046734888E0,2.018112620856775083915565E2,
123
CS 1 2.290838373831346393026739E3,1.131967205903380828685045E4,
124
CS 2 2.855724635671635335736389E4,3.848496228443793359990269E4,
125
CS 3 2.637748787624195437963534E4,7.225813979700288197698961E3/
126
CS DATA Q1/6.748212550303777196073036E1,1.113332393857199323513008E3,
127
CS 1 7.738757056935398733233834E3,2.763987074403340708898585E4,
128
CS 2 5.499310206226157329794414E4,6.161122180066002127833352E4,
129
CS 3 3.635127591501940507276287E4,8.785536302431013170870835E3/
130
DATA D1/-5.772156649015328605195174D-1/
131
DATA P1/4.945235359296727046734888D0,2.018112620856775083915565D2,
132
1 2.290838373831346393026739D3,1.131967205903380828685045D4,
133
2 2.855724635671635335736389D4,3.848496228443793359990269D4,
134
3 2.637748787624195437963534D4,7.225813979700288197698961D3/
135
DATA Q1/6.748212550303777196073036D1,1.113332393857199323513008D3,
136
1 7.738757056935398733233834D3,2.763987074403340708898585D4,
137
2 5.499310206226157329794414D4,6.161122180066002127833352D4,
138
3 3.635127591501940507276287D4,8.785536302431013170870835D3/
139
C----------------------------------------------------------------------
140
C Numerator and denominator coefficients for rational minimax
141
C Approximation over (1.5,4.0).
142
C----------------------------------------------------------------------
143
CS DATA D2/4.227843350984671393993777E-1/
144
CS DATA P2/4.974607845568932035012064E0,5.424138599891070494101986E2,
145
CS 1 1.550693864978364947665077E4,1.847932904445632425417223E5,
146
CS 2 1.088204769468828767498470E6,3.338152967987029735917223E6,
147
CS 3 5.106661678927352456275255E6,3.074109054850539556250927E6/
148
CS DATA Q2/1.830328399370592604055942E2,7.765049321445005871323047E3,
149
CS 1 1.331903827966074194402448E5,1.136705821321969608938755E6,
150
CS 2 5.267964117437946917577538E6,1.346701454311101692290052E7,
151
CS 3 1.782736530353274213975932E7,9.533095591844353613395747E6/
152
DATA D2/4.227843350984671393993777D-1/
153
DATA P2/4.974607845568932035012064D0,5.424138599891070494101986D2,
154
1 1.550693864978364947665077D4,1.847932904445632425417223D5,
155
2 1.088204769468828767498470D6,3.338152967987029735917223D6,
156
3 5.106661678927352456275255D6,3.074109054850539556250927D6/
157
DATA Q2/1.830328399370592604055942D2,7.765049321445005871323047D3,
158
1 1.331903827966074194402448D5,1.136705821321969608938755D6,
159
2 5.267964117437946917577538D6,1.346701454311101692290052D7,
160
3 1.782736530353274213975932D7,9.533095591844353613395747D6/
161
C----------------------------------------------------------------------
162
C Numerator and denominator coefficients for rational minimax
163
C Approximation over (4.0,12.0).
164
C----------------------------------------------------------------------
165
CS DATA D4/1.791759469228055000094023E0/
166
CS DATA P4/1.474502166059939948905062E4,2.426813369486704502836312E6,
167
CS 1 1.214755574045093227939592E8,2.663432449630976949898078E9,
168
CS 2 2.940378956634553899906876E10,1.702665737765398868392998E11,
169
CS 3 4.926125793377430887588120E11,5.606251856223951465078242E11/
170
CS DATA Q4/2.690530175870899333379843E3,6.393885654300092398984238E5,
171
CS 2 4.135599930241388052042842E7,1.120872109616147941376570E9,
172
CS 3 1.488613728678813811542398E10,1.016803586272438228077304E11,
173
CS 4 3.417476345507377132798597E11,4.463158187419713286462081E11/
174
DATA D4/1.791759469228055000094023D0/
175
DATA P4/1.474502166059939948905062D4,2.426813369486704502836312D6,
176
1 1.214755574045093227939592D8,2.663432449630976949898078D9,
177
2 2.940378956634553899906876D10,1.702665737765398868392998D11,
178
3 4.926125793377430887588120D11,5.606251856223951465078242D11/
179
DATA Q4/2.690530175870899333379843D3,6.393885654300092398984238D5,
180
2 4.135599930241388052042842D7,1.120872109616147941376570D9,
181
3 1.488613728678813811542398D10,1.016803586272438228077304D11,
182
4 3.417476345507377132798597D11,4.463158187419713286462081D11/
183
C----------------------------------------------------------------------
184
C Coefficients for minimax approximation over (12, INF).
185
C----------------------------------------------------------------------
186
CS DATA C/-1.910444077728E-03,8.4171387781295E-04,
187
CS 1 -5.952379913043012E-04,7.93650793500350248E-04,
188
CS 2 -2.777777777777681622553E-03,8.333333333333333331554247E-02,
189
CS 3 5.7083835261E-03/
190
DATA C/-1.910444077728D-03,8.4171387781295D-04,
191
1 -5.952379913043012D-04,7.93650793500350248D-04,
192
2 -2.777777777777681622553D-03,8.333333333333333331554247D-02,
194
C----------------------------------------------------------------------
196
IF ((Y .GT. ZERO) .AND. (Y .LE. XBIG)) THEN
199
ELSE IF (Y .LE. THRHAL) THEN
200
C----------------------------------------------------------------------
201
C EPS .LT. X .LE. 1.5
202
C----------------------------------------------------------------------
203
IF (Y .LT. PNT68) THEN
208
XM1 = (Y - HALF) - HALF
210
IF ((Y .LE. HALF) .OR. (Y .GE. PNT68)) THEN
214
XNUM = XNUM*XM1 + P1(I)
215
XDEN = XDEN*XM1 + Q1(I)
217
RES = CORR + (XM1 * (D1 + XM1*(XNUM/XDEN)))
219
XM2 = (Y - HALF) - HALF
223
XNUM = XNUM*XM2 + P2(I)
224
XDEN = XDEN*XM2 + Q2(I)
226
RES = CORR + XM2 * (D2 + XM2*(XNUM/XDEN))
228
ELSE IF (Y .LE. FOUR) THEN
229
C----------------------------------------------------------------------
230
C 1.5 .LT. X .LE. 4.0
231
C----------------------------------------------------------------------
236
XNUM = XNUM*XM2 + P2(I)
237
XDEN = XDEN*XM2 + Q2(I)
239
RES = XM2 * (D2 + XM2*(XNUM/XDEN))
240
ELSE IF (Y .LE. TWELVE) THEN
241
C----------------------------------------------------------------------
242
C 4.0 .LT. X .LE. 12.0
243
C----------------------------------------------------------------------
248
XNUM = XNUM*XM4 + P4(I)
249
XDEN = XDEN*XM4 + Q4(I)
251
RES = D4 + XM4*(XNUM/XDEN)
253
C----------------------------------------------------------------------
254
C Evaluate for argument .GE. 12.0,
255
C----------------------------------------------------------------------
257
IF (Y .LE. FRTBIG) THEN
261
RES = RES / YSQ + C(I)
266
RES = RES + SQRTPI - HALF*CORR
267
RES = RES + Y*(CORR-ONE)
270
C----------------------------------------------------------------------
271
C Return for bad arguments
272
C----------------------------------------------------------------------
275
C----------------------------------------------------------------------
276
C Final adjustments and return
277
C----------------------------------------------------------------------
281
C ---------- Last line of DLGAMA ----------
added
removed
routines/calelm/dlgama.f
