! 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 *********************