!     Last change:  DOS  29 Jul 2000   11:57 am
!                                       *** copyright 2000 ***
! *** filename gamma.f95                *** John F. Monahan **
!                                       **********************
      PROGRAM PGAMMA
!  DEMONSTRATION/TEST OF GAMMA AND ALGAMA
!   GAMMA -- COMPUTES GAMMA FUNCTION
!   ALGAMA -- COMPUTES NATURAL LOG OF GAMMA FUNCTION
!
!  COMPUTE TABLES TO MATCH ABRAMOWITZ AND STEGUN TABLES 6.1, 6.3, 6.4
!
      IMPLICIT NONE
      REAL X,Y,Z,AL10
      REAL GAMMA,ALGAMA
      INTEGER I
!
  21  FORMAT(F6.2,2F15.7)                                               
  22  FORMAT(F6.2,2E20.8)
  25  FORMAT('TABLE 6.1 -- GAMMA FROM 1 TO 2')
  26  FORMAT('    X          GAMMA(X)        ALGAMA(X) ')
  27  FORMAT(//'TABLE 6.3 -- GAMMA FOR INTEGER VALUES')
  28  FORMAT(//'TABLE 6.4 -- LOG10(GAMMA(X))')
  29  FORMAT('    X        ALGAMA(X)      LOG10(GAMMA(X))')
!
      OPEN( UNIT=6, FILE='gamma.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 )/20.
      Y = GAMMA(X)
      Z = ALGAMA(X) 
      WRITE(6,21) X,Y,Z                                                    
      END DO  ! LOOP ON I
!                TABLE 6.3 HAS GAMMA FOR INTEGERS 1 TO 35 (OVERFLOW)
      WRITE(6,27)
      WRITE(6,26)
      DO I = 1,35
      X = REAL(I) 
      Y = GAMMA( X )
      Z = ALGAMA( 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.)
      DO I = 1,101,4
      X = REAL(I)                                               
      Y = ALGAMA(X)
      Z = Y / AL10
      WRITE(6,21) X,Y,Z                                                   
      END DO  ! LOOP ON I
      STOP
      END PROGRAM PGAMMA
      REAL FUNCTION GAMMA(X)
!  COMPUTES GAMMA FUNCTION FOR POSITIVE ARGUMENTS
!   X .LE. 0.  RETURNS -1.0
!   X .GE. 35.0 RETURNS GAMMA(35) SINCE BIGGER OVERFLOWS
!  USE ALGAMA FOR LARGE ARGUMENTS
      IMPLICIT NONE
      REAL, INTENT(IN) :: X
      REAL, DIMENSION(5) :: P,Q         ! COEFS FOR RATIONAL FUNCTION
      REAL, DIMENSION(3) :: R           ! COEFS FOR STIRLING CORRECTION
      REAL XN,A,V
      REAL, PARAMETER :: LSQR2PI = .918938533
      DATA P/ -.510046056E1, -.368309363E2, -.664664448E2,  &
     &        -.195692802E3, -.185072903E1 /
      DATA Q/ -.114321842E2,  .368848969E2,  .16269403E2,   &
     &        -.195692802E3, 1.0 /
      DATA R/ -.277765545E-2, .833333330E-1, .77783067E-3 /      
!   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. ) THEN
      GAMMA = -1.0            ! NOT DESIGNED FOR NEGATIVE ARGUMENTS
      RETURN
                        ENDIF
!                          FOR LARGE X (>8) USE STIRLING ROUTE
      IF( X .GT. 8. ) THEN
!                           LARGE X -- USE STIRLING WITH CORRECTION
      XN = X * X
!                           CORRECTION FOR STIRLING
      V = R(2) + ( R(3)/XN + R(1) )/XN                                  
      V = (X-0.5)*LOG(X) - X + V/X + LSQR2PI
      GAMMA = EXP( V )
                      ELSE
!                          FOR SMALL X, USE REPRODUCTION FORMULA
!                           TO GET NUMBER BETWEEN 2 AND 3
      V = 1.0                                                             
      XN = X                                                              
      DO WHILE( XN .LT. 2. )
      V = V/XN
      XN = XN + 1.
      END DO  ! WHILE( XN .LT. 2. )
      DO WHILE( XN .GT. 3. )
      XN = XN - 1.
      V = V * XN
      END DO  ! WHILE( XN .GT. 3. )
      A = XN-2.0
!                           RATIONAL FUNCTION APPROXIMATION
      GAMMA = V*(P(4)+A*(P(3)+A*(P(2)+A*(P(1)+A*P(5))))) /  &
     &   (Q(4)+A*(Q(3)+A*(Q(2)+A*(Q(1)+A*Q(5)))))
                      ENDIF
      RETURN
      END FUNCTION GAMMA
      REAL FUNCTION ALGAMA(X)
!  COMPUTES NATURAL LOG OF GAMMA FUNCTION FOR POSITIVE ARGUMENTS
!   X .LE. 0.  RETURNS -1.0
      IMPLICIT NONE
      REAL, INTENT(IN) :: X
      REAL, DIMENSION(5) :: P,Q         ! COEFS FOR RATIONAL FUNCTION
      REAL, DIMENSION(3) :: R           ! COEFS FOR STIRLING CORRECTION
      REAL XN,A,V
      REAL, PARAMETER :: LSQR2PI = .918938533
      DATA P/ -.510046056E1, -.368309363E2, -.664664448E2,  &
     &        -.195692802E3, -.185072903E1 /
      DATA Q/ -.114321842E2,  .368848969E2,  .16269403E2,   &
     &        -.195692802E3, 1.0 /
      DATA R/ -.277765545E-2, .833333330E-1, .77783067E-3 /      
!   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. ) THEN
      ALGAMA = -1.0            ! NOT DESIGNED FOR NEGATIVE ARGUMENTS
      RETURN
                        ENDIF
!                          FOR LARGE X (>8) USE STIRLING ROUTE
      IF( X .GT. 8. ) THEN
!                           LARGE X -- USE STIRLING WITH CORRECTION
      XN = X * X
!                           CORRECTION FOR STIRLING
      V = R(2) + ( R(3)/XN + R(1) )/XN                                  
      ALGAMA = (X-0.5)*LOG(X) - X + V/X + LSQR2PI
                      ELSE
!                          FOR SMALL X, USE REPRODUCTION FORMULA
!                           TO GET NUMBER BETWEEN 2 AND 3
      V = 1.0                                                             
      XN = X                                                              
      DO WHILE( XN .LT. 2. )
      V = V/XN
      XN = XN + 1.
      END DO  ! WHILE( XN .LT. 2. )
      DO WHILE( XN .GT. 3. )
      XN = XN - 1.
      V = V * XN
      END DO  ! WHILE( XN .GT. 3. )
      A = XN-2.0
!                           RATIONAL FUNCTION APPROXIMATION
      V = V*(P(4)+A*(P(3)+A*(P(2)+A*(P(1)+A*P(5))))) /  &
     &   (Q(4)+A*(Q(3)+A*(Q(2)+A*(Q(1)+A*Q(5)))))
      ALGAMA = LOG(V)
                      ENDIF
      RETURN
      END FUNCTION ALGAMA
! *** end of filename gamma.f95         *********************