CS REAL FUNCTION ALGAMA(X) CD DOUBLE PRECISION FUNCTION DLGAMA(X) C---------------------------------------------------------------------- C C This routine calculates the LOG(GAMMA) function for a positive real C argument X. Computation is based on an algorithm outlined in C references 1 and 2. The program uses rational functions that C theoretically approximate LOG(GAMMA) to at least 18 significant C decimal digits. The approximation for X > 12 is from reference C 3, while approximations for X < 12.0 are similar to those in C reference 1, but are unpublished. The accuracy achieved depends C on the arithmetic system, the compiler, the intrinsic functions, C and proper selection of the machine-dependent constants. C C C********************************************************************* C********************************************************************* C C Explanation of machine-dependent constants C C beta - radix for the floating-point representation C maxexp - the smallest positive power of beta that overflows C XBIG - largest argument for which LN(GAMMA(X)) is representable C in the machine, i.e., the solution to the equation C LN(GAMMA(XBIG)) = beta**maxexp C XINF - largest machine representable floating-point number; C approximately beta**maxexp. C EPS - The smallest positive floating-point number such that C 1.0+EPS .GT. 1.0 C FRTBIG - Rough estimate of the fourth root of XBIG C C C Approximate values for some important machines are: C C beta maxexp XBIG C C CRAY-1 (S.P.) 2 8191 9.62E+2461 C Cyber 180/855 C under NOS (S.P.) 2 1070 1.72E+319 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 2 128 4.08E+36 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 2 1024 2.55D+305 C IBM 3033 (D.P.) 16 63 4.29D+73 C VAX D-Format (D.P.) 2 127 2.05D+36 C VAX G-Format (D.P.) 2 1023 1.28D+305 C C C XINF EPS FRTBIG C C CRAY-1 (S.P.) 5.45E+2465 7.11E-15 3.13E+615 C Cyber 180/855 C under NOS (S.P.) 1.26E+322 3.55E-15 6.44E+79 C IEEE (IBM/XT, C SUN, etc.) (S.P.) 3.40E+38 1.19E-7 1.42E+9 C IEEE (IBM/XT, C SUN, etc.) (D.P.) 1.79D+308 2.22D-16 2.25D+76 C IBM 3033 (D.P.) 7.23D+75 2.22D-16 2.56D+18 C VAX D-Format (D.P.) 1.70D+38 1.39D-17 1.20D+9 C VAX G-Format (D.P.) 8.98D+307 1.11D-16 1.89D+76 C C************************************************************** C************************************************************** C C Error returns C C The program returns the value XINF for X .LE. 0.0 or when C overflow would occur. The computation is believed to C be free of underflow and overflow. C C C Intrinsic functions required are: C C LOG C C C References: C C 1) W. J. Cody and K. E. Hillstrom, 'Chebyshev Approximations for C the Natural Logarithm of the Gamma Function,' Math. Comp. 21, C 1967, pp. 198-203. C C 2) K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May, C 1969. C C 3) Hart, Et. Al., Computer Approximations, Wiley and sons, New C York, 1968. C C C Authors: W. J. Cody and L. Stoltz C Argonne National Laboratory C C Latest modification: June 16, 1988 C C---------------------------------------------------------------------- INTEGER I CS REAL CD DOUBLE PRECISION 1 C,CORR,D1,D2,D4,EPS,FRTBIG,FOUR,HALF,ONE,PNT68,P1,P2,P4, 2 Q1,Q2,Q4,RES,SQRTPI,THRHAL,TWELVE,TWO,X,XBIG,XDEN,XINF, 3 XM1,XM2,XM4,XNUM,Y,YSQ,ZERO DIMENSION C(7),P1(8),P2(8),P4(8),Q1(8),Q2(8),Q4(8) C---------------------------------------------------------------------- C Mathematical constants C---------------------------------------------------------------------- CS DATA ONE,HALF,TWELVE,ZERO/1.0E0,0.5E0,12.0E0,0.0E0/, CS 1 FOUR,THRHAL,TWO,PNT68/4.0E0,1.5E0,2.0E0,0.6796875E0/, CS 2 SQRTPI/0.9189385332046727417803297E0/ CD DATA ONE,HALF,TWELVE,ZERO/1.0D0,0.5D0,12.0D0,0.0D0/, CD 1 FOUR,THRHAL,TWO,PNT68/4.0D0,1.5D0,2.0D0,0.6796875D0/, CD 2 SQRTPI/0.9189385332046727417803297D0/ C---------------------------------------------------------------------- C Machine dependent parameters C---------------------------------------------------------------------- CS DATA XBIG,XINF,EPS,FRTBIG/4.08E36,3.401E38,1.19E-7,1.42E9/ CD DATA XBIG,XINF,EPS,FRTBIG/2.55D305,1.79D308,2.22D-16,2.25D76/ C---------------------------------------------------------------------- C Numerator and denominator coefficients for rational minimax C approximation over (0.5,1.5). C---------------------------------------------------------------------- CS DATA D1/-5.772156649015328605195174E-1/ CS DATA P1/4.945235359296727046734888E0,2.018112620856775083915565E2, CS 1 2.290838373831346393026739E3,1.131967205903380828685045E4, CS 2 2.855724635671635335736389E4,3.848496228443793359990269E4, CS 3 2.637748787624195437963534E4,7.225813979700288197698961E3/ CS DATA Q1/6.748212550303777196073036E1,1.113332393857199323513008E3, CS 1 7.738757056935398733233834E3,2.763987074403340708898585E4, CS 2 5.499310206226157329794414E4,6.161122180066002127833352E4, CS 3 3.635127591501940507276287E4,8.785536302431013170870835E3/ CD DATA D1/-5.772156649015328605195174D-1/ CD DATA P1/4.945235359296727046734888D0,2.018112620856775083915565D2, CD 1 2.290838373831346393026739D3,1.131967205903380828685045D4, CD 2 2.855724635671635335736389D4,3.848496228443793359990269D4, CD 3 2.637748787624195437963534D4,7.225813979700288197698961D3/ CD DATA Q1/6.748212550303777196073036D1,1.113332393857199323513008D3, CD 1 7.738757056935398733233834D3,2.763987074403340708898585D4, CD 2 5.499310206226157329794414D4,6.161122180066002127833352D4, CD 3 3.635127591501940507276287D4,8.785536302431013170870835D3/ C---------------------------------------------------------------------- C Numerator and denominator coefficients for rational minimax C Approximation over (1.5,4.0). C---------------------------------------------------------------------- CS DATA D2/4.227843350984671393993777E-1/ CS DATA P2/4.974607845568932035012064E0,5.424138599891070494101986E2, CS 1 1.550693864978364947665077E4,1.847932904445632425417223E5, CS 2 1.088204769468828767498470E6,3.338152967987029735917223E6, CS 3 5.106661678927352456275255E6,3.074109054850539556250927E6/ CS DATA Q2/1.830328399370592604055942E2,7.765049321445005871323047E3, CS 1 1.331903827966074194402448E5,1.136705821321969608938755E6, CS 2 5.267964117437946917577538E6,1.346701454311101692290052E7, CS 3 1.782736530353274213975932E7,9.533095591844353613395747E6/ CD DATA D2/4.227843350984671393993777D-1/ CD DATA P2/4.974607845568932035012064D0,5.424138599891070494101986D2, CD 1 1.550693864978364947665077D4,1.847932904445632425417223D5, CD 2 1.088204769468828767498470D6,3.338152967987029735917223D6, CD 3 5.106661678927352456275255D6,3.074109054850539556250927D6/ CD DATA Q2/1.830328399370592604055942D2,7.765049321445005871323047D3, CD 1 1.331903827966074194402448D5,1.136705821321969608938755D6, CD 2 5.267964117437946917577538D6,1.346701454311101692290052D7, CD 3 1.782736530353274213975932D7,9.533095591844353613395747D6/ C---------------------------------------------------------------------- C Numerator and denominator coefficients for rational minimax C Approximation over (4.0,12.0). C---------------------------------------------------------------------- CS DATA D4/1.791759469228055000094023E0/ CS DATA P4/1.474502166059939948905062E4,2.426813369486704502836312E6, CS 1 1.214755574045093227939592E8,2.663432449630976949898078E9, CS 2 2.940378956634553899906876E10,1.702665737765398868392998E11, CS 3 4.926125793377430887588120E11,5.606251856223951465078242E11/ CS DATA Q4/2.690530175870899333379843E3,6.393885654300092398984238E5, CS 2 4.135599930241388052042842E7,1.120872109616147941376570E9, CS 3 1.488613728678813811542398E10,1.016803586272438228077304E11, CS 4 3.417476345507377132798597E11,4.463158187419713286462081E11/ CD DATA D4/1.791759469228055000094023D0/ CD DATA P4/1.474502166059939948905062D4,2.426813369486704502836312D6, CD 1 1.214755574045093227939592D8,2.663432449630976949898078D9, CD 2 2.940378956634553899906876D10,1.702665737765398868392998D11, CD 3 4.926125793377430887588120D11,5.606251856223951465078242D11/ CD DATA Q4/2.690530175870899333379843D3,6.393885654300092398984238D5, CD 2 4.135599930241388052042842D7,1.120872109616147941376570D9, CD 3 1.488613728678813811542398D10,1.016803586272438228077304D11, CD 4 3.417476345507377132798597D11,4.463158187419713286462081D11/ C---------------------------------------------------------------------- C Coefficients for minimax approximation over (12, INF). C---------------------------------------------------------------------- CS DATA C/-1.910444077728E-03,8.4171387781295E-04, CS 1 -5.952379913043012E-04,7.93650793500350248E-04, CS 2 -2.777777777777681622553E-03,8.333333333333333331554247E-02, CS 3 5.7083835261E-03/ CD DATA C/-1.910444077728D-03,8.4171387781295D-04, CD 1 -5.952379913043012D-04,7.93650793500350248D-04, CD 2 -2.777777777777681622553D-03,8.333333333333333331554247D-02, CD 3 5.7083835261D-03/ C---------------------------------------------------------------------- Y = X IF ((Y .GT. ZERO) .AND. (Y .LE. XBIG)) THEN IF (Y .LE. EPS) THEN RES = -LOG(Y) ELSE IF (Y .LE. THRHAL) THEN C---------------------------------------------------------------------- C EPS .LT. X .LE. 1.5 C---------------------------------------------------------------------- IF (Y .LT. PNT68) THEN CORR = -LOG(Y) XM1 = Y ELSE CORR = ZERO XM1 = (Y - HALF) - HALF END IF IF ((Y .LE. HALF) .OR. (Y .GE. PNT68)) THEN XDEN = ONE XNUM = ZERO DO 140 I = 1, 8 XNUM = XNUM*XM1 + P1(I) XDEN = XDEN*XM1 + Q1(I) 140 CONTINUE RES = CORR + (XM1 * (D1 + XM1*(XNUM/XDEN))) ELSE XM2 = (Y - HALF) - HALF XDEN = ONE XNUM = ZERO DO 220 I = 1, 8 XNUM = XNUM*XM2 + P2(I) XDEN = XDEN*XM2 + Q2(I) 220 CONTINUE RES = CORR + XM2 * (D2 + XM2*(XNUM/XDEN)) END IF ELSE IF (Y .LE. FOUR) THEN C---------------------------------------------------------------------- C 1.5 .LT. X .LE. 4.0 C---------------------------------------------------------------------- XM2 = Y - TWO XDEN = ONE XNUM = ZERO DO 240 I = 1, 8 XNUM = XNUM*XM2 + P2(I) XDEN = XDEN*XM2 + Q2(I) 240 CONTINUE RES = XM2 * (D2 + XM2*(XNUM/XDEN)) ELSE IF (Y .LE. TWELVE) THEN C---------------------------------------------------------------------- C 4.0 .LT. X .LE. 12.0 C---------------------------------------------------------------------- XM4 = Y - FOUR XDEN = -ONE XNUM = ZERO DO 340 I = 1, 8 XNUM = XNUM*XM4 + P4(I) XDEN = XDEN*XM4 + Q4(I) 340 CONTINUE RES = D4 + XM4*(XNUM/XDEN) ELSE C---------------------------------------------------------------------- C Evaluate for argument .GE. 12.0, C---------------------------------------------------------------------- RES = ZERO IF (Y .LE. FRTBIG) THEN RES = C(7) YSQ = Y * Y DO 450 I = 1, 6 RES = RES / YSQ + C(I) 450 CONTINUE END IF RES = RES/Y CORR = LOG(Y) RES = RES + SQRTPI - HALF*CORR RES = RES + Y*(CORR-ONE) END IF ELSE C---------------------------------------------------------------------- C Return for bad arguments C---------------------------------------------------------------------- RES = XINF END IF C---------------------------------------------------------------------- C Final adjustments and return C---------------------------------------------------------------------- CS ALGAMA = RES CD DLGAMA = RES RETURN C ---------- Last line of DLGAMA ---------- END