REAL FUNCTION fGAMMA(X) !D DOUBLE PRECISION FUNCTION fGAMMA(X) !---------------------------------------------------------------------- ! ! This routine calculates the GAMMA function for a real argument X. ! Computation is based on an algorithm outlined in reference 1. ! The program uses rational functions that approximate the GAMMA ! function to at least 20 significant decimal digits. Coefficients ! for the approximation over the interval (1,2) are unpublished. ! Those for the approximation for X .GE. 12 are from reference 2. ! The accuracy achieved depends on the arithmetic system, the ! compiler, the intrinsic functions, and proper selection of the ! machine-dependent constants. ! ! !******************************************************************* !******************************************************************* ! ! Explanation of machine-dependent constants ! ! beta - radix for the floating-point representation ! maxexp - the smallest positive power of beta that overflows ! XBIG - the largest argument for which GAMMA(X) is representable ! in the machine, i.e., the solution to the equation ! GAMMA(XBIG) = beta**maxexp ! XINF - the largest machine representable floating-point number; ! approximately beta**maxexp ! EPS - the smallest positive floating-point number such that ! 1.0+EPS .GT. 1.0 ! XMININ - the smallest positive floating-point number such that ! 1/XMININ is machine representable ! ! Approximate values for some important machines are: ! ! beta maxexp XBIG ! ! CRAY-1 (S.P.) 2 8191 966.961 ! Cyber 180/855 ! under NOS (S.P.) 2 1070 177.803 ! IEEE (IBM/XT, ! SUN, etc.) (S.P.) 2 128 35.040 ! IEEE (IBM/XT, ! SUN, etc.) (D.P.) 2 1024 171.624 ! IBM 3033 (D.P.) 16 63 57.574 ! VAX D-Format (D.P.) 2 127 34.844 ! VAX G-Format (D.P.) 2 1023 171.489 ! ! XINF EPS XMININ ! ! CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466 ! Cyber 180/855 ! under NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294 ! IEEE (IBM/XT, ! SUN, etc.) (S.P.) 3.40E+38 1.19E-7 1.18E-38 ! IEEE (IBM/XT, ! SUN, etc.) (D.P.) 1.79D+308 2.22D-16 2.23D-308 ! IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76 ! VAX D-Format (D.P.) 1.70D+38 1.39D-17 5.88D-39 ! VAX G-Format (D.P.) 8.98D+307 1.11D-16 1.12D-308 ! !******************************************************************* !******************************************************************* ! ! Error returns ! ! The program returns the value XINF for singularities or ! when overflow would occur. The computation is believed ! to be free of underflow and overflow. ! ! ! Intrinsic functions required are: ! ! INT, DBLE, EXP, LOG, REAL, SIN ! ! ! References: "An Overview of Software Development for Special ! Functions", W. J. Cody, Lecture Notes in Mathematics, ! 506, Numerical Analysis Dundee, 1975, G. A. Watson ! (ed.), Springer Verlag, Berlin, 1976. ! ! Computer Approximations, Hart, Et. Al., Wiley and ! sons, New York, 1968. ! ! Latest modification: October 12, 1989 ! ! Authors: W. J. Cody and L. Stoltz ! Applied Mathematics Division ! Argonne National Laboratory ! Argonne, IL 60439 ! !---------------------------------------------------------------------- implicit none INTEGER I,N LOGICAL PARITY REAL C,CONV,EPS,FACT,HALF,ONE,P,PI,Q,RES,SQRTPI,SUM,TWELVE, & TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO DIMENSION C(7),P(8),Q(8) !---------------------------------------------------------------------- ! Mathematical constants !---------------------------------------------------------------------- DATA ONE,HALF,TWELVE,TWO,ZERO/1.0E0,0.5E0,12.0E0,2.0E0,0.0E0/, & SQRTPI/0.9189385332046727417803297E0/, & PI/3.1415926535897932384626434E0/ !D DATA ONE,HALF,TWELVE,TWO,ZERO/1.0D0,0.5D0,12.0D0,2.0D0,0.0D0/, !D 1 SQRTPI/0.9189385332046727417803297D0/, !D 2 PI/3.1415926535897932384626434D0/ !---------------------------------------------------------------------- ! Machine dependent parameters !---------------------------------------------------------------------- DATA XBIG,XMININ,EPS/35.040E0,1.18E-38,1.19E-7/, & XINF/3.4E38/ !D DATA XBIG,XMININ,EPS/171.624D0,2.23D-308,2.22D-16/, !D 1 XINF/1.79D308/ !---------------------------------------------------------------------- ! Numerator and denominator coefficients for rational minimax ! approximation over (1,2). !---------------------------------------------------------------------- DATA P/-1.71618513886549492533811E+0,2.47656508055759199108314E+1,& -3.79804256470945635097577E+2,6.29331155312818442661052E+2,& 8.66966202790413211295064E+2,-3.14512729688483675254357E+4,& -3.61444134186911729807069E+4,6.64561438202405440627855E+4/ DATA Q/-3.08402300119738975254353E+1,3.15350626979604161529144E+2,& -1.01515636749021914166146E+3,-3.10777167157231109440444E+3,& 2.25381184209801510330112E+4,4.75584627752788110767815E+3,& -1.34659959864969306392456E+5,-1.15132259675553483497211E+5/ !D DATA P/-1.71618513886549492533811D+0,2.47656508055759199108314D+1, !D 1 -3.79804256470945635097577D+2,6.29331155312818442661052D+2, !D 2 8.66966202790413211295064D+2,-3.14512729688483675254357D+4, !D 3 -3.61444134186911729807069D+4,6.64561438202405440627855D+4/ !D DATA Q/-3.08402300119738975254353D+1,3.15350626979604161529144D+2, !D 1 -1.01515636749021914166146D+3,-3.10777167157231109440444D+3, !D 2 2.25381184209801510330112D+4,4.75584627752788110767815D+3, !D 3 -1.34659959864969306392456D+5,-1.15132259675553483497211D+5/ !---------------------------------------------------------------------- ! Coefficients for minimax approximation over (12, INF). !---------------------------------------------------------------------- DATA C/-1.910444077728E-03,8.4171387781295E-04, & -5.952379913043012E-04,7.93650793500350248E-04, & -2.777777777777681622553E-03,8.333333333333333331554247E-02, & 5.7083835261E-03/ !D DATA C/-1.910444077728D-03,8.4171387781295D-04, !D 1 -5.952379913043012D-04,7.93650793500350248D-04, !D 2 -2.777777777777681622553D-03,8.333333333333333331554247D-02, !D 3 5.7083835261D-03/ !---------------------------------------------------------------------- ! Statement functions for conversion between integer and float !---------------------------------------------------------------------- CONV(I) = REAL(I) !D CONV(I) = DBLE(I) PARITY = .FALSE. FACT = ONE N = 0 Y = X IF (Y .LE. ZERO) THEN !---------------------------------------------------------------------- ! Argument is negative !---------------------------------------------------------------------- Y = -X Y1 = AINT(Y) RES = Y - Y1 IF (RES .NE. ZERO) THEN IF (Y1 .NE. AINT(Y1*HALF)*TWO) PARITY = .TRUE. FACT = -PI / SIN(PI*RES) Y = Y + ONE ELSE RES = XINF GO TO 900 END IF END IF !---------------------------------------------------------------------- ! Argument is positive !---------------------------------------------------------------------- IF (Y .LT. EPS) THEN !---------------------------------------------------------------------- ! Argument .LT. EPS !---------------------------------------------------------------------- IF (Y .GE. XMININ) THEN RES = ONE / Y ELSE RES = XINF GO TO 900 END IF ELSE IF (Y .LT. TWELVE) THEN Y1 = Y IF (Y .LT. ONE) THEN !---------------------------------------------------------------------- ! 0.0 .LT. argument .LT. 1.0 !---------------------------------------------------------------------- Z = Y Y = Y + ONE ELSE !---------------------------------------------------------------------- ! 1.0 .LT. argument .LT. 12.0, reduce argument if necessary !---------------------------------------------------------------------- N = INT(Y) - 1 Y = Y - CONV(N) Z = Y - ONE END IF !---------------------------------------------------------------------- ! Evaluate approximation for 1.0 .LT. argument .LT. 2.0 !---------------------------------------------------------------------- XNUM = ZERO XDEN = ONE DO 260 I = 1, 8 XNUM = (XNUM + P(I)) * Z XDEN = XDEN * Z + Q(I) 260 CONTINUE RES = XNUM / XDEN + ONE IF (Y1 .LT. Y) THEN !---------------------------------------------------------------------- ! Adjust result for case 0.0 .LT. argument .LT. 1.0 !---------------------------------------------------------------------- RES = RES / Y1 ELSE IF (Y1 .GT. Y) THEN !---------------------------------------------------------------------- ! Adjust result for case 2.0 .LT. argument .LT. 12.0 !---------------------------------------------------------------------- DO 290 I = 1, N RES = RES * Y Y = Y + ONE 290 CONTINUE END IF ELSE !---------------------------------------------------------------------- ! Evaluate for argument .GE. 12.0, !---------------------------------------------------------------------- IF (Y .LE. XBIG) THEN YSQ = Y * Y SUM = C(7) DO 350 I = 1, 6 SUM = SUM / YSQ + C(I) 350 CONTINUE SUM = SUM/Y - Y + SQRTPI SUM = SUM + (Y-HALF)*LOG(Y) RES = EXP(SUM) ELSE RES = XINF GO TO 900 END IF END IF !---------------------------------------------------------------------- ! Final adjustments and return !---------------------------------------------------------------------- IF (PARITY) RES = -RES IF (FACT .NE. ONE) RES = FACT / RES 900 fGAMMA = RES RETURN ! ---------- Last line of fGAMMA ---------- END