!     Last change:  DOS  29 Jul 2000   12:11 pm
!                                       *** copyright 2000 ***
! *** filename pgama.f95                *** John F. Monahan **
!                                       **********************
      program ppgama
!  test program for incomplete gamma function
      implicit none
      real x,a,p,pp,r,q,qq
      real pgama,qgama
      integer i,j
!
  20  format(' Test of incomplete gamma function '/3x,'a',8x,'x',8x,  &
     &       'pgama',4x,'sum',9x,'qgama',8x,'sum')
  21  format(2f8.4,4f12.8)
  22  format(2f14.6,2e20.8)
!      
      open( unit=6, file='pgama.out' )
!
      write(6,20)
!     
!    first test on integer values of a & compare to analytic result
!
      do j = 1,20
      x = real(j)/5.
      pp = 1.
      qq = 0.
      r = exp(-x)
      do i = 1,5
      a = real(i)
      pp = pp - r
      qq = qq + r
      r = r*x/real(i)
      p = pgama(a,x)
      q = qgama(a,x)
      write(6,21) a,x,p,pp,q,qq
      end do  ! loop on i
      end do  ! loop on j
!
!    second test -- get a lot of values -- put in file for comparison
!
      open( unit=8, file='pgama.tst' )
      do i = 1,30
      a = real(i*i)/10.
      do j = 1,40
      x = a*real(j)/16.
      p = pgama(a,x)
      q = qgama(a,x)
      write(8,22) a,x,p,q
      end do  ! loop on j
      end do  ! loop on i
      stop
      end program ppgama
      real function pgama(a,x)
!  computes incomplete gamma function
!
!  integral from 0 to x of  (x**(a-1) ) * exp(-x) dx / gamma(a)
!
!  if x <= 0 then pgama = 0
!  if a <= 0 then pgama = -1
!  if convergence problems then pgama = -2
!
!  follows methods for regions I, III and tail as described in
!  Walter Gautschi (1979) "A Computational Procedure for Incomplete
!   Gamma Functions," ACM TOMS, Volume 5, pp. 466-481.
!
! J F Monahan (June 1997) Dept of Stat, NC State U, Raleigh, NC USA
!   recoded February, April 2000 for Fortran 95
      implicit none
      real, intent(in) :: a,x
      real, parameter :: eps = 1.e-6
      real t,r,st,ak,fk
      real algama
      integer k
!    first flag some out of range values
      pgama = -1.
      if( a .le. 0. ) return
      pgama = 0.
      if( x .le. 0. ) return
!
!    triage -- three methods depending on a and x
!
      if( a .gt. x +.25 ) then
!
!   monotone series for little g
!
      t = 1./a
      st = t
      do k = 1,50             ! iteration max is 50
      t = t*x/(a+real(k))
      st = st + t
      if( abs(t) .le. st*eps ) exit
      end do  !  loop on k
      if( k .le. 50 ) then
      pgama = st*exp(a*log(x)-x-algama(a))
                      else
!
!    failure to converge -- error
!
      pgama = -2.
                      end if  ! if ( k .le. 50 )
      return
                          end if  ! if( a .gt. x+.25 )
      if( x .gt. 1.5 ) then
!
!    Gautschi's Region III
!    rewrite continued fraction as power series
!
      r = 0.
      t = 1.
      st = 1.
      do k = 1,50             ! iteration max is 50
      fk = real(k)
      ak = fk*(a-fk)/( (x+fk+fk-1.-a) * (x+fk+fk+1.-a) )
      r = -ak*(1.+r)/( 1.+ak*(1.+r) )
      t = t*r
      st = st + t
      if( abs(t) .le. st*eps ) exit
      end do  ! loop on k
      if( k .le. 50 ) then
      pgama = 1. - st*exp(a*log(x)-x-algama(a))/(x+1.-a)
                      else
!
!    failure to converge -- error
!
      pgama = -2.
                      end if  ! if ( k .le. 50 )
      return
                          end if  ! if( x .gt. 1.5 )
!
!    Gautschi's Region I
!    integration by parts and alternating series expansion (4.10)
!
      st = 1.
      r = 1.
      do k = 1,50             ! iteration max is 50
!                *** notice k+1 here
      fk = real(k+1)
      r = -x*r/fk
      t = (a+1.) * r / (a+fk)
      st = st + t
      if( abs(t) .le. st*eps ) exit
      end do  ! loop on k
      if( k .le. 50 ) then
      pgama = (1.+1./a-st*x)*exp( a*log(x) - algama(a) ) / (a+1.)
                      else
!
!    failure to converge -- error
!
      pgama = -2.
                      end if  ! if ( k .le. 50 )
      return
      end function pgama
      real function qgama(a,x)
!  computes incomplete gamma function
!
!  integral from x to infinity of  (x**(a-1) ) * exp(-x) dx / gamma(a)
!
!  if x <= 0 then qgama = 1
!  if a <= 0 then qgama = -1
!  if convergence problems then qgama = -2
!
!  follows methods for regions I, III and tail as described in
!  Walter Gautschi (1979) "A Computational Procedure for Incomplete
!   Gamma Functions," ACM TOMS, Volume 5, pp. 466-481.
!
! J F Monahan (June 1997) Dept of Stat, NC State U, Raleigh, NC USA
!   recoded February, April 2000 for Fortran 95
      implicit none
      real, intent(in) :: a,x
      real, parameter :: eps = 1.e-6
      real t,r,st,ak,fk
      real algama
      integer k
!    first flag some out of range values
      qgama = -1.
      if( a .le. 0. ) return
      qgama = 1.
      if( x .le. 0. ) return
!
!    triage -- three methods depending on a and x
!
      if( a .gt. x +.25 ) then
!
!   monotone series for little g
!
      t = 1./a
      st = t
      do k = 1,50             ! iteration max is 50
      t = t*x/(a+real(k))
      st = st + t
      if( abs(t) .le. st*eps ) exit
      end do  !  loop on k
      if( k .le. 50 ) then
      qgama = 1. - st*exp(a*log(x)-x-algama(a))
                      else
!
!    failure to converge -- error
!
      qgama = -2.
                      end if  ! if ( k .le. 50 )
      return
                          end if  ! if( a .gt. x+.25 )
      if( x .gt. 1.5 ) then
!
!    Gautschi's Region III
!    rewrite continued fraction as power series
!
      r = 0.
      t = 1.
      st = 1.
      do k = 1,50             ! iteration max is 50
      fk = real(k)
      ak = fk*(a-fk)/( (x+fk+fk-1.-a) * (x+fk+fk+1.-a) )
      r = -ak*(1.+r)/( 1.+ak*(1.+r) )
      t = t*r
      st = st + t
      if( abs(t) .le. st*eps ) exit
      end do  ! loop on k
      if( k .le. 50 ) then
      qgama = st*exp(a*log(x)-x-algama(a))/(x+1.-a)
                      else
!
!    failure to converge -- error
!
      qgama = -2.
                      end if  ! if ( k .le. 50 )
      return
                          end if  ! if( x .gt. 1.5 )
!
!    Gautschi's Region I
!    integration by parts and alternating series expansion (4.10)
!
      st = 1.
      r = 1.
      do k = 1,50             ! iteration max is 50
!                *** notice k+1 here
      fk = real(k+1)
      r = -x*r/fk
      t = (a+1.) * r / (a+fk)
      st = st + t
      if( abs(t) .le. st*eps ) exit
      end do  ! loop on k
      if( k .le. 50 ) then
      qgama = 1. - (1.+1./a-st*x)*exp( a*log(x) - algama(a) ) / (a+1.)
                      else
!
!    failure to converge -- error
!
      qgama = -2.
                      end if  ! if ( k .le. 50 )
      return
      end function qgama
      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 pgama.f95         *********************