gcf.f90

FORTRANvisualfortran常用数值算法集及源码

SUBROUTINE gcf(gammcf,a,x,gln)
INTEGER ITMAX
REAL a,gammcf,gln,x,EPS,FPMIN
PARAMETER (ITMAX=100,EPS=3.e-7,FPMIN=1.e-30)
!USES gammln
INTEGER i
REAL an,b,c,d,del,h,gammln
gln=gammln(a)
b=x+1.-a
c=1./FPMIN
d=1./b
h=d
do i=1,ITMAX
  an=-i*(i-a)
  b=b+2.
  d=an*d+b
  if(abs(d)<FPMIN) d=FPMIN
  c=b+an/c
  if(abs(c)<FPMIN) c=FPMIN
  d=1./d
  del=d*c
  h=h*del
  if(abs(del-1.)<EPS) then
    gammcf=exp(-x+a*log(x)-gln)*h
    return
  endif
end do
pause 'a too large, ITMAX too small in gcf'
END SUBROUTINE gcf