factrl.f90

来自「FORTRANvisualfortran常用数值算法集及源码」· F90 代码 · 共 23 行

FUNCTION factrl(n)
INTEGER n
REAL factrl
!USES gammln
INTEGER j,ntop
REAL a(33),gammln
SAVE ntop,a
DATA ntop,a(1)/0,1./
if (n<0) then
  pause 'negative factorial in factrl'
else if (n<=ntop) then
  factrl=a(n+1)
else if (n<=32) then
  do j=ntop+1,n
    a(j+1)=j*a(j)
  end do
  ntop=n
  factrl=a(n+1)
else
  factrl=exp(gammln(n+1.))
endif
END FUNCTION factrl