C 階乗の近似計算
C メインプログラム
      PROGRAM KAIJOU
      DOUBLE PRECISION STRLNG,FCTRL,STRLG2,STRLG3,STRLGL
      WRITE(*,*) 'Give N : '
      READ(*,*)N
      WRITE(*,200)'Value of direct multiplication : ',FCTRL(N)
C      WRITE(*,200)'Value of the Stirling formula  : ',STRLNG(N+1)
C      WRITE(*,200)'Value of the 2nd asymptotic    : ',STRLG2(N+1)
      WRITE(*,200)'Value of the 3rd asymptotic    : ',STRLG3(N+1)
C      WRITE(*,300)'Striling through log           : ',
C     +    STRLGL(N+1,IEX),' E',IEX
  200 FORMAT(1H ,A33,D25.18)
  300 FORMAT(1H ,A33,F21.18,A2,I4)
      END
C
C 階乗を定義で計算
C
      DOUBLE PRECISION FUNCTION FCTRL(N)
      FCTRL=1.0D0
      DO 100 I=2,N
          FCTRL=FCTRL*I
  100 CONTINUE
      RETURN
      END
C
C Stirling の公式を使う
C
      DOUBLE PRECISION FUNCTION STRLNG(N)
      DOUBLE PRECISION PI,E
      PI=3.14159265358979323846
      E =2.71828182845904523536
      STRLNG=DSQRT(PI*2/N)*(DBLE(N)/E)**N
C ここを次のように書くと，N^N の部分が一足先にオーバーフローしてしまい，
C 直接 N! を計算する方法に負けてしまう．
C     STRLNG=DSQRT(PI*2/N)*DEXP(-DBLE(N))*DBLE(N)**N
      RETURN
      END
C
C Stirling の公式に漸近第２項を追加したもの
C
      DOUBLE PRECISION FUNCTION STRLG2(N)
      DOUBLE PRECISION PI,E
      PARAMETER (PI=3.14159265358979323846D0,E=2.71828182845904523536D0)
      STRLG2=DSQRT(PI*2/N)*(DBLE(N)/E)**(N)*DEXP(1.0D0/12/N)
      RETURN
      END
C
C Stirling の公式に漸近第３項まで追加したもの
C
      DOUBLE PRECISION FUNCTION STRLG3(N)
      REAL*8 PI,E
      PARAMETER (PI=3.14159265358979323846D0,E=2.71828182845904523536D0)
      PRINT *,PI
      PRINT *,PI*2/N
      PRINT *,DSQRT(PI*2/N)
      PRINT *,(DBLE(N)/E)**N
      PRINT *,DEXP(1.0D0/12/N-1.0D0/360/N/N/N)
      STRLG3=DSQRT(PI*2/N)*(DBLE(N)/E)**N
     +       *DEXP(1.0D0/12/N-1.0D0/360/N/N/N)
      RETURN
      END
C
C Stirling の公式の log を用いてオーバーフローしないようにしたもの
C
      DOUBLE PRECISION FUNCTION STRLGL(N,IEX)
      DOUBLE PRECISION PI,E,DL10E,DL
      PI=3.14159265358979323846
      E =2.71828182845904523536
      DL10E=0.43429448190325182765
      DL=DLOG(PI*2/N)/2+DLOG(DBLE(N)/E)*N
     +       +(DBLE(1)/12/N-DBLE(1)/360/N/N/N)
      DL=DL*DL10E
      IEX=IDINT(DL)+1
      STRLGL=DBLE(10)**(DL-IEX)
      RETURN
      END