1CS REAL FUNCTION ALGAMA(X) 2 DOUBLE PRECISION FUNCTION DLGAMA(X) 3C---------------------------------------------------------------------- 4C 5C This routine calculates the LOG(GAMMA) function for a positive real 6C argument X. Computation is based on an algorithm outlined in 7C references 1 and 2. The program uses rational functions that 8C theoretically approximate LOG(GAMMA) to at least 18 significant 9C decimal digits. The approximation for X > 12 is from reference 10C 3, while approximations for X < 12.0 are similar to those in 11C reference 1, but are unpublished. The accuracy achieved depends 12C on the arithmetic system, the compiler, the intrinsic functions, 13C and proper selection of the machine-dependent constants. 14C 15C 16C********************************************************************* 17C********************************************************************* 18C 19C Explanation of machine-dependent constants 20C 21C beta - radix for the floating-point representation 22C maxexp - the smallest positive power of beta that overflows 23C XBIG - largest argument for which LN(GAMMA(X)) is representable 24C in the machine, i.e., the solution to the equation 25C LN(GAMMA(XBIG)) = beta**maxexp 26C XINF - largest machine representable floating-point number; 27C approximately beta**maxexp. 28C EPS - The smallest positive floating-point number such that 29C 1.0+EPS .GT. 1.0 30C FRTBIG - Rough estimate of the fourth root of XBIG 31C 32C 33C Approximate values for some important machines are: 34C 35C beta maxexp XBIG 36C 37C CRAY-1 (S.P.) 2 8191 9.62E+2461 38C Cyber 180/855 39C under NOS (S.P.) 2 1070 1.72E+319 40C IEEE (IBM/XT, 41C SUN, etc.) (S.P.) 2 128 4.08E+36 42C IEEE (IBM/XT, 43C SUN, etc.) (D.P.) 2 1024 2.55D+305 44C IBM 3033 (D.P.) 16 63 4.29D+73 45C VAX D-Format (D.P.) 2 127 2.05D+36 46C VAX G-Format (D.P.) 2 1023 1.28D+305 47C 48C 49C XINF EPS FRTBIG 50C 51C CRAY-1 (S.P.) 5.45E+2465 7.11E-15 3.13E+615 52C Cyber 180/855 53C under NOS (S.P.) 1.26E+322 3.55E-15 6.44E+79 54C IEEE (IBM/XT, 55C SUN, etc.) (S.P.) 3.40E+38 1.19E-7 1.42E+9 56C IEEE (IBM/XT, 57C SUN, etc.) (D.P.) 1.79D+308 2.22D-16 2.25D+76 58C IBM 3033 (D.P.) 7.23D+75 2.22D-16 2.56D+18 59C VAX D-Format (D.P.) 1.70D+38 1.39D-17 1.20D+9 60C VAX G-Format (D.P.) 8.98D+307 1.11D-16 1.89D+76 61C 62C************************************************************** 63C************************************************************** 64C 65C Error returns 66C 67C The program returns the value XINF for X .LE. 0.0 or when 68C overflow would occur. The computation is believed to 69C be free of underflow and overflow. 70C 71C 72C Intrinsic functions required are: 73C 74C LOG 75C 76C 77C References: 78C 79C 1) W. J. Cody and K. E. Hillstrom, 'Chebyshev Approximations for 80C the Natural Logarithm of the Gamma Function,' Math. Comp. 21, 81C 1967, pp. 198-203. 82C 83C 2) K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May, 84C 1969. 85C 86C 3) Hart, Et. Al., Computer Approximations, Wiley and sons, New 87C York, 1968. 88C 89C 90C Authors: W. J. Cody and L. Stoltz 91C Argonne National Laboratory 92C 93C Latest modification: June 16, 1988 94C 95C---------------------------------------------------------------------- 96* 97* Some modifs from Bruno (25 Feb 2005): 98* - return Nan (in place of XINF) if x <= 0 99* - return Inf (in place of XINF) if x > XBIG 100* 101* In fact an error indicator should be returned in these cases to prevent 102* the user... 103* 104 INTEGER I 105CS REAL 106 DOUBLE PRECISION 107 1 C,CORR,D1,D2,D4,EPS,FRTBIG,FOUR,HALF,ONE,PNT68,P1,P2,P4, 108 2 Q1,Q2,Q4,RES,SQRTPI,THRHAL,TWELVE,TWO,X,XBIG,XDEN,XINF, 109 3 XM1,XM2,XM4,XNUM,Y,YSQ,ZERO 110 DIMENSION C(7),P1(8),P2(8),P4(8),Q1(8),Q2(8),Q4(8) 111C---------------------------------------------------------------------- 112C Mathematical constants 113C---------------------------------------------------------------------- 114CS DATA ONE,HALF,TWELVE,ZERO/1.0E0,0.5E0,12.0E0,0.0E0/, 115CS 1 FOUR,THRHAL,TWO,PNT68/4.0E0,1.5E0,2.0E0,0.6796875E0/, 116CS 2 SQRTPI/0.9189385332046727417803297E0/ 117 DATA ONE,HALF,TWELVE,ZERO/1.0D0,0.5D0,12.0D0,0.0D0/, 118 1 FOUR,THRHAL,TWO,PNT68/4.0D0,1.5D0,2.0D0,0.6796875D0/, 119 2 SQRTPI/0.9189385332046727417803297D0/ 120C---------------------------------------------------------------------- 121C Machine dependent parameters 122C---------------------------------------------------------------------- 123CS DATA XBIG,XINF,EPS,FRTBIG/4.08E36,3.401E38,1.19E-7,1.42E9/ 124 DATA XBIG,XINF,EPS,FRTBIG/2.55D305,1.79D308,2.22D-16,2.25D76/ 125C---------------------------------------------------------------------- 126C Numerator and denominator coefficients for rational minimax 127C approximation over (0.5,1.5). 128C---------------------------------------------------------------------- 129CS DATA D1/-5.772156649015328605195174E-1/ 130CS DATA P1/4.945235359296727046734888E0,2.018112620856775083915565E2, 131CS 1 2.290838373831346393026739E3,1.131967205903380828685045E4, 132CS 2 2.855724635671635335736389E4,3.848496228443793359990269E4, 133CS 3 2.637748787624195437963534E4,7.225813979700288197698961E3/ 134CS DATA Q1/6.748212550303777196073036E1,1.113332393857199323513008E3, 135CS 1 7.738757056935398733233834E3,2.763987074403340708898585E4, 136CS 2 5.499310206226157329794414E4,6.161122180066002127833352E4, 137CS 3 3.635127591501940507276287E4,8.785536302431013170870835E3/ 138 DATA D1/-5.772156649015328605195174D-1/ 139 DATA P1/4.945235359296727046734888D0,2.018112620856775083915565D2, 140 1 2.290838373831346393026739D3,1.131967205903380828685045D4, 141 2 2.855724635671635335736389D4,3.848496228443793359990269D4, 142 3 2.637748787624195437963534D4,7.225813979700288197698961D3/ 143 DATA Q1/6.748212550303777196073036D1,1.113332393857199323513008D3, 144 1 7.738757056935398733233834D3,2.763987074403340708898585D4, 145 2 5.499310206226157329794414D4,6.161122180066002127833352D4, 146 3 3.635127591501940507276287D4,8.785536302431013170870835D3/ 147C---------------------------------------------------------------------- 148C Numerator and denominator coefficients for rational minimax 149C Approximation over (1.5,4.0). 150C---------------------------------------------------------------------- 151CS DATA D2/4.227843350984671393993777E-1/ 152CS DATA P2/4.974607845568932035012064E0,5.424138599891070494101986E2, 153CS 1 1.550693864978364947665077E4,1.847932904445632425417223E5, 154CS 2 1.088204769468828767498470E6,3.338152967987029735917223E6, 155CS 3 5.106661678927352456275255E6,3.074109054850539556250927E6/ 156CS DATA Q2/1.830328399370592604055942E2,7.765049321445005871323047E3, 157CS 1 1.331903827966074194402448E5,1.136705821321969608938755E6, 158CS 2 5.267964117437946917577538E6,1.346701454311101692290052E7, 159CS 3 1.782736530353274213975932E7,9.533095591844353613395747E6/ 160 DATA D2/4.227843350984671393993777D-1/ 161 DATA P2/4.974607845568932035012064D0,5.424138599891070494101986D2, 162 1 1.550693864978364947665077D4,1.847932904445632425417223D5, 163 2 1.088204769468828767498470D6,3.338152967987029735917223D6, 164 3 5.106661678927352456275255D6,3.074109054850539556250927D6/ 165 DATA Q2/1.830328399370592604055942D2,7.765049321445005871323047D3, 166 1 1.331903827966074194402448D5,1.136705821321969608938755D6, 167 2 5.267964117437946917577538D6,1.346701454311101692290052D7, 168 3 1.782736530353274213975932D7,9.533095591844353613395747D6/ 169C---------------------------------------------------------------------- 170C Numerator and denominator coefficients for rational minimax 171C Approximation over (4.0,12.0). 172C---------------------------------------------------------------------- 173CS DATA D4/1.791759469228055000094023E0/ 174CS DATA P4/1.474502166059939948905062E4,2.426813369486704502836312E6, 175CS 1 1.214755574045093227939592E8,2.663432449630976949898078E9, 176CS 2 2.940378956634553899906876E10,1.702665737765398868392998E11, 177CS 3 4.926125793377430887588120E11,5.606251856223951465078242E11/ 178CS DATA Q4/2.690530175870899333379843E3,6.393885654300092398984238E5, 179CS 2 4.135599930241388052042842E7,1.120872109616147941376570E9, 180CS 3 1.488613728678813811542398E10,1.016803586272438228077304E11, 181CS 4 3.417476345507377132798597E11,4.463158187419713286462081E11/ 182 DATA D4/1.791759469228055000094023D0/ 183 DATA P4/1.474502166059939948905062D4,2.426813369486704502836312D6, 184 1 1.214755574045093227939592D8,2.663432449630976949898078D9, 185 2 2.940378956634553899906876D10,1.702665737765398868392998D11, 186 3 4.926125793377430887588120D11,5.606251856223951465078242D11/ 187 DATA Q4/2.690530175870899333379843D3,6.393885654300092398984238D5, 188 2 4.135599930241388052042842D7,1.120872109616147941376570D9, 189 3 1.488613728678813811542398D10,1.016803586272438228077304D11, 190 4 3.417476345507377132798597D11,4.463158187419713286462081D11/ 191C---------------------------------------------------------------------- 192C Coefficients for minimax approximation over (12, INF). 193C---------------------------------------------------------------------- 194CS DATA C/-1.910444077728E-03,8.4171387781295E-04, 195CS 1 -5.952379913043012E-04,7.93650793500350248E-04, 196CS 2 -2.777777777777681622553E-03,8.333333333333333331554247E-02, 197CS 3 5.7083835261E-03/ 198 DATA C/-1.910444077728D-03,8.4171387781295D-04, 199 1 -5.952379913043012D-04,7.93650793500350248D-04, 200 2 -2.777777777777681622553D-03,8.333333333333333331554247D-02, 201 3 5.7083835261D-03/ 202C---------------------------------------------------------------------- 203 Y = X 204 IF ((Y .GT. ZERO) .AND. (Y .LE. XBIG)) THEN 205 IF (Y .LE. EPS) THEN 206 RES = -LOG(Y) 207 ELSE IF (Y .LE. THRHAL) THEN 208C---------------------------------------------------------------------- 209C EPS .LT. X .LE. 1.5 210C---------------------------------------------------------------------- 211 IF (Y .LT. PNT68) THEN 212 CORR = -LOG(Y) 213 XM1 = Y 214 ELSE 215 CORR = ZERO 216 XM1 = (Y - HALF) - HALF 217 END IF 218 IF ((Y .LE. HALF) .OR. (Y .GE. PNT68)) THEN 219 XDEN = ONE 220 XNUM = ZERO 221 DO 140 I = 1, 8 222 XNUM = XNUM*XM1 + P1(I) 223 XDEN = XDEN*XM1 + Q1(I) 224 140 CONTINUE 225 RES = CORR + (XM1 * (D1 + XM1*(XNUM/XDEN))) 226 ELSE 227 XM2 = (Y - HALF) - HALF 228 XDEN = ONE 229 XNUM = ZERO 230 DO 220 I = 1, 8 231 XNUM = XNUM*XM2 + P2(I) 232 XDEN = XDEN*XM2 + Q2(I) 233 220 CONTINUE 234 RES = CORR + XM2 * (D2 + XM2*(XNUM/XDEN)) 235 END IF 236 ELSE IF (Y .LE. FOUR) THEN 237C---------------------------------------------------------------------- 238C 1.5 .LT. X .LE. 4.0 239C---------------------------------------------------------------------- 240 XM2 = Y - TWO 241 XDEN = ONE 242 XNUM = ZERO 243 DO 240 I = 1, 8 244 XNUM = XNUM*XM2 + P2(I) 245 XDEN = XDEN*XM2 + Q2(I) 246 240 CONTINUE 247 RES = XM2 * (D2 + XM2*(XNUM/XDEN)) 248 ELSE IF (Y .LE. TWELVE) THEN 249C---------------------------------------------------------------------- 250C 4.0 .LT. X .LE. 12.0 251C---------------------------------------------------------------------- 252 XM4 = Y - FOUR 253 XDEN = -ONE 254 XNUM = ZERO 255 DO 340 I = 1, 8 256 XNUM = XNUM*XM4 + P4(I) 257 XDEN = XDEN*XM4 + Q4(I) 258 340 CONTINUE 259 RES = D4 + XM4*(XNUM/XDEN) 260 ELSE 261C---------------------------------------------------------------------- 262C Evaluate for argument .GE. 12.0, 263C---------------------------------------------------------------------- 264 RES = ZERO 265 IF (Y .LE. FRTBIG) THEN 266 RES = C(7) 267 YSQ = Y * Y 268 DO 450 I = 1, 6 269 RES = RES / YSQ + C(I) 270 450 CONTINUE 271 END IF 272 RES = RES/Y 273 CORR = LOG(Y) 274 RES = RES + SQRTPI - HALF*CORR 275 RES = RES + Y*(CORR-ONE) 276 END IF 277 ELSE 278C---------------------------------------------------------------------- 279C Return for bad arguments 280C---------------------------------------------------------------------- 281* modif from Bruno (see comment at the beginning) 282* RES = XINF 283 if (X .le. 0.d0) then 284 CALL returnananfortran(RES) 285 else 286C this means that X > XBIG and so that log(gamma) overflows 287C bad trick to get Inf 288 RES = 2*XINF 289 endif 290* end modif 291 END IF 292C---------------------------------------------------------------------- 293C Final adjustments and return 294C---------------------------------------------------------------------- 295CS ALGAMA = RES 296 DLGAMA = RES 297 RETURN 298C ---------- Last line of DLGAMA ---------- 299 END 300 301 302