!     Last change:  DOS  29 Jul 2000   12:08 pm
!                                       *** copyright 2000 ***
! *** filename pdgama.f95               *** John F. Monahan **
!                                       **********************
      program ppdgama
!  test program for incomplete gamma function
      implicit none
      real(kind=8) x,a,p,pp,r,q,qq
      real(kind=8) pdgama,qdgama
      integer i,j
!
  20  format(' Test of incomplete gamma function '/3x,'a',8x,'x',8x,  &
     &       'pdgama',4x,'sum',9x,'qdgama',8x,'sum')
  21  format(2f8.4,4f12.8)
  22  format(2f14.6,2e20.8)
!      
      open( unit=6, file='pdgama.out' )
!
      write(6,20)
!     
!    first test on integer values of a & compare to analytic result
!
      do j = 1,20
      x = real(j,8)/5.d0                !KIND
      pp = 1.
      qq = 0.
      r = exp(-x)
      do i = 1,5
      a = real(i,8)                     !KIND
      pp = pp - r
      qq = qq + r
      r = r*x/a
      p = pdgama(a,x)
      q = qdgama(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='pdgama.tst' )
      do i = 1,30
      a = real(i*i,8)/10.d0             !KIND
      do j = 1,40
      x = a*real(j,8)/16.d0             !KIND
      p = pdgama(a,x)
      q = qdgama(a,x)
      write(8,22) a,x,p,q
      end do  ! loop on j
      end do  ! loop on i
      stop
      end program ppdgama
      real(kind=8) function pdgama(a,x)
!  computes incomplete gamma function
!
!  integral from 0 to x of  (x**(a-1) ) * exp(-x) dx / gamma(a)
!
!  if x <= 0 then pdgama = 0
!  if a <= 0 then pdgama = -1
!  if convergence problems then pdgama = -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(kind=8), intent(in) :: a,x
      real(kind=8), parameter :: eps = 1.d-12
      real(kind=8) t,r,st,ak,fk
      real(kind=8) dlgama
      integer k
!    first flag some out of range values
      pdgama = -1.d0
      if( a .le. 0. ) return
      pdgama = 0.d0
      if( x .le. 0.d0 ) return
!
!    triage -- three methods depending on a and x
!
      if( a .gt. x +.25d0 ) then
!
!   monotone series for little g
!
      t = 1.d0/a
      st = t
      do k = 1,100             ! iteration max is 100
      t = t*x/(a+real(k,8))             !KIND
      st = st + t
      if( abs(t) .le. st*eps ) exit
      end do  !  loop on k
      if( k .le. 100 ) then
      pdgama = st*exp(a*log(x)-x-dlgama(a))
                      else
!
!    failure to converge -- error
!
      pdgama = -2.d0
                      end if  ! if ( k .le. 100 )
      return
                          end if  ! if( a .gt. x+.25 )
      if( x .gt. 1.5d0 ) then
!
!    Gautschi's Region III
!    rewrite continued fraction as power series
!
      r = 0.d0
      t = 1.d0
      st = 1.d0
      do k = 1,100             ! iteration max is 100
      fk = real(k,8)                    !KIND
      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. 100 ) then
      pdgama = 1.d0 - st*exp(a*log(x)-x-dlgama(a))/(x+1.-a)
                      else
!
!    failure to converge -- error
!
      pdgama = -2.d0
                      end if  ! if ( k .le. 100 )
      return
                          end if  ! if( x .gt. 1.5 )
!
!    Gautschi's Region I
!    integration by parts and alternating series expansion (4.10)
!
      st = 1.d0
      r = 1.d0
      do k = 1,100             ! iteration max is 100
!                *** notice k+1 here
      fk = real(k+1,8)                  !KIND
      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. 100 ) then
      pdgama = (1.+1./a-st*x)*exp( a*log(x) - dlgama(a) ) / (a+1.)
                      else
!
!    failure to converge -- error
!
      pdgama = -2.d0
                      end if  ! if ( k .le. 100 )
      return
      end function pdgama
      real(kind=8) function qdgama(a,x)
!  computes incomplete gamma function
!
!  integral from x to infinity of  (x**(a-1) ) * exp(-x) dx / gamma(a)
!
!  if x <= 0 then qdgama = 1
!  if a <= 0 then qdgama = -1
!  if convergence problems then qdgama = -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(kind=8), intent(in) :: a,x
      real(kind=8), parameter :: eps = 1.d-12
      real(kind=8) t,r,st,ak,fk
      real(kind=8) dlgama
      integer k
!    first flag some out of range values
      qdgama = -1.d0
      if( a .le. 0. ) return
      qdgama = 1.d0
      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,100             ! iteration max is 100
      t = t*x/(a+real(k,8))             !KIND
      st = st + t
      if( abs(t) .le. st*eps ) exit
      end do  !  loop on k
      if( k .le. 100 ) then
      qdgama = 1.d0 - st*exp(a*log(x)-x-dlgama(a))
                      else
!
!    failure to converge -- error
!
      qdgama = -2.d0
                      end if  ! if ( k .le. 100 )
      return
                          end if  ! if( a .gt. x+.25 )
      if( x .gt. 1.5d0 ) then
!
!    Gautschi's Region III
!    rewrite continued fraction as power series
!
      r = 0.d0
      t = 1.d0
      st = 1.d0
      do k = 1,100             ! iteration max is 100
      fk = real(k,8)                    !KIND
      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. 100 ) then
      qdgama = st*exp(a*log(x)-x-dlgama(a))/(x+1.-a)
                      else
!
!    failure to converge -- error
!
      qdgama = -2.d0
                      end if  ! if ( k .le. 100 )
      return
                          end if  ! if( x .gt. 1.5 )
!
!    Gautschi's Region I
!    integration by parts and alternating series expansion (4.10)
!
      st = 1.d0
      r = 1.d0
      do k = 1,100             ! iteration max is 100
!                *** notice k+1 here
      fk = real(k+1,8)                  !KIND
      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. 100 ) then
      qdgama = 1. - (1.+1./a-st*x)*exp( a*log(x) - dlgama(a) ) / (a+1.)
                      else
!
!    failure to converge -- error
!
      qdgama = -2.d0
                      end if  ! if ( k .le. 100 )
      return
      end function qdgama
      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 pdgama.f95        *********************