! Last change: DOS 29 Jul 2000 11:52 am ! *** copyright 2000 *** ! *** filename dgamma.f95 *** John F. Monahan ** ! ********************** PROGRAM PDGAMMA ! DEMONSTRATION/TEST OF DGAMMA AND DLGAMA ! DGAMMA -- COMPUTES GAMMA FUNCTION ! DLGAMA -- COMPUTES NATURAL LOG OF GAMMA FUNCTION ! DOUBLE PRECISION(KIND=8 OR 2 -- DEPENDING ON COMPILER) VERSIONS ! ! COMPUTE TABLES TO MATCH ABRAMOWITZ AND STEGUN TABLES 6.1, 6.3, 6.4 ! IMPLICIT NONE REAL(KIND=8) X,Y,Z,AL10 REAL(KIND=8) DGAMMA,DLGAMA INTEGER I ! 21 FORMAT(F6.2,2F20.15) 22 FORMAT(F6.2,2E25.15) 25 FORMAT('TABLE 6.1 -- GAMMA FROM 1 TO 2') 26 FORMAT(' X DGAMMA(X) DLGAMA(X) ') 27 FORMAT(//'TABLE 6.3 -- GAMMA FOR INTEGER VALUES') 28 FORMAT(//'TABLE 6.4 -- LOG10(GAMMA(X))') 29 FORMAT(' X DLGAMA(X) LOG10(GAMMA(X))') ! OPEN( UNIT=6, FILE='dgamma.out' ) ! TABLE 6.1 HAS GAMMA FOR X FROM 1 TO 2 WRITE(6,25) WRITE(6,26) DO I = 1,21 X = REAL( 19 + I ,8)/20.D0 !KIND Y = DGAMMA(X) Z = DLGAMA(X) WRITE(6,21) X,Y,Z END DO ! LOOP ON I ! TABLE 6.3 HAS GAMMA FOR INTEGERS 1 TO 50 WRITE(6,27) WRITE(6,26) DO I = 1,50 X = REAL(I,8) !KIND Y = DGAMMA( X ) Z = DLGAMA( X ) WRITE(6,22) X,Y,Z END DO ! LOOP ON I ! TABLE 6.4 HAS LOG10(GAMMA(X)) FOR INTEGER 1 TO 101 WRITE(6,28) WRITE(6,29) AL10 = LOG(10.D0) DO I = 1,101,4 X = REAL(I,8) !KIND Y = DLGAMA(X) Z = Y / AL10 WRITE(6,21) X,Y,Z END DO ! LOOP ON I STOP END PROGRAM PDGAMMA REAL(KIND=8) FUNCTION DGAMMA(X) ! COMPUTES GAMMA FUNCTION FOR POSITIVE ARGUMENTS ! X .LE. 0. RETURNS -1.0 ! X .GE. 35.0 RETURNS GAMMA(35) SINCE BIGGER OVERFLOWS ! USE DLGAMA FOR LARGE ARGUMENTS IMPLICIT NONE REAL(KIND=8), INTENT(IN) :: X REAL(KIND=8), DIMENSION(7) :: P ! COEFS FOR RATIONAL FUNCTION REAL(KIND=8), DIMENSION(6) :: Q REAL(KIND=8), DIMENSION(3) :: PT ! COEFS FOR STIRLING CORRECTION REAL(KIND=8), DIMENSION(2) :: QT REAL(KIND=8) XN,A,V REAL(KIND=8), PARAMETER :: LSQR2PI = .918938533204672742D0 DATA P/ .378601050348257245d4, .207745979389418732d4, & & .893581804523749814d3, .222112396168011795d3, & & .489543462279099381d2, .61260674503360843d1, & & .778079585613300576d0 / DATA Q/ .378601050348257187d4, .476793860503687915d3, & & -.867230987531102994d3, .835500586679197696d2, & & .507884753288954097d2, -.134004147857813483d2 / DATA PT/ .288119283935546015d0, .498030766924499634d0, & & .691561607375687d-1 / DATA QT/ .345743140722674507d1, .609161691641660296d1 / ! USES METHODS 5230,5401 FROM 'COMPUTER APPROXIMATIONS' ! ED. BY JF HART, SIAM SERIES, WILEY(1968) ! 27 MAR 74 JFM W/RWS,RM ! ** MODIFICATIONS JUNE 1988, 1993 ** ! RECODED OCTOBER 1999 FOR FORTRAN 95 IF( X .LE. 0.D0 ) THEN DGAMMA = -1.D0 ! NOT DESIGNED FOR NEGATIVE ARGUMENTS RETURN ENDIF ! FOR LARGE X (>8) USE STIRLING ROUTE IF( X .GT. 8.D0 ) THEN ! LARGE X -- USE STIRLING WITH CORRECTION XN = X * X ! CORRECTION FOR STIRLING V = (PT(1)+(PT(2)+PT(3)/XN)/XN) / (QT(1)+(QT(2)+1.D0/XN)/XN) V = (X-0.5D0)*LOG(X) - X + V/X + LSQR2PI DGAMMA = EXP( V ) ELSE ! FOR SMALL X, USE REPRODUCTION FORMULA ! TO GET NUMBER BETWEEN 2 AND 3 V = 1.D0 XN = X DO WHILE( XN .LT. 2.D0 ) V = V/XN XN = XN + 1.D0 END DO ! WHILE( XN .LT. 2.D0 ) DO WHILE( XN .GT. 3.D0 ) XN = XN - 1.D0 V = V * XN END DO ! WHILE( XN .GT. 3.D0 ) A = XN-2.D0 ! RATIONAL FUNCTION APPROXIMATION DGAMMA=(P(1)+A*(P(2)+A*(P(3)+A*(P(4)+A*(P(5)+A*(P(6)+A*P(7)))))))& & * V / (Q(1)+A*(Q(2)+A*(Q(3)+A*(Q(4)+A*(Q(5)+A*(Q(6)+A)))))) ENDIF RETURN END FUNCTION DGAMMA REAL(KIND=8) FUNCTION DLGAMA(X) ! COMPUTES NATURAL LOG OF GAMMA FUNCTION FOR POSITIVE ARGUMENTS ! X .LE. 0. RETURNS -1.0 IMPLICIT NONE REAL(KIND=8), INTENT(IN) :: X REAL(KIND=8), DIMENSION(7) :: P ! COEFS FOR RATIONAL FUNCTION REAL(KIND=8), DIMENSION(6) :: Q REAL(KIND=8), DIMENSION(3) :: PT ! COEFS FOR STIRLING CORRECTION REAL(KIND=8), DIMENSION(2) :: QT REAL(KIND=8) XN,A,V REAL(KIND=8), PARAMETER :: LSQR2PI = .918938533204672742D0 DATA P/ .378601050348257245d4, .207745979389418732d4, & & .893581804523749814d3, .222112396168011795d3, & & .489543462279099381d2, .61260674503360843d1, & & .778079585613300576d0 / DATA Q/ .378601050348257187d4, .476793860503687915d3, & & -.867230987531102994d3, .835500586679197696d2, & & .507884753288954097d2, -.134004147857813483d2 / DATA PT/ .288119283935546015d0, .498030766924499634d0, & & .691561607375687d-1 / DATA QT/ .345743140722674507d1, .609161691641660296d1 / ! USES METHODS 5230,5401 FROM 'COMPUTER APPROXIMATIONS' ! ED. BY JF HART, SIAM SERIES, WILEY(1968) ! 27 MAR 74 JFM W/RWS,RM ! ** MODIFICATIONS JUNE 1988, 1993 ** ! RECODED OCTOBER 1999 FOR FORTRAN 95 IF( X .LE. 0.D0 ) THEN DLGAMA = -1.D0 ! NOT DESIGNED FOR NEGATIVE ARGUMENTS RETURN ENDIF ! FOR LARGE X (>8) USE STIRLING ROUTE IF( X .GT. 8.D0 ) THEN ! LARGE X -- USE STIRLING WITH CORRECTION XN = X * X ! CORRECTION FOR STIRLING V = (PT(1)+(PT(2)+PT(3)/XN)/XN) / (QT(1)+(QT(2)+1.D0/XN)/XN) DLGAMA = (X-0.5D0)*LOG(X) - X + V/X + LSQR2PI ELSE ! FOR SMALL X, USE REPRODUCTION FORMULA ! TO GET NUMBER BETWEEN 2 AND 3 V = 1.D0 XN = X DO WHILE( XN .LT. 2.D0 ) V = V/XN XN = XN + 1.D0 END DO ! WHILE( XN .LT. 2.D0 ) DO WHILE( XN .GT. 3.D0 ) XN = XN - 1.D0 V = V * XN END DO ! WHILE( XN .GT. 3.D0 ) A = XN-2.D0 ! RATIONAL FUNCTION APPROXIMATION V =(P(1)+A*(P(2)+A*(P(3)+A*(P(4)+A*(P(5)+A*(P(6)+A*P(7))))))) & & * V / (Q(1)+A*(Q(2)+A*(Q(3)+A*(Q(4)+A*(Q(5)+A*(Q(6)+A)))))) DLGAMA = LOG(V) ENDIF RETURN END FUNCTION DLGAMA ! *** end of filename dgamma.f95 *********************