fortran66のブログ

fortran について書きます。

統計力学などでは、よく階乗の対数を Stirling の公式で近似しますが、どのくらいの精度で近似になっているのか見てみます。

ここで Stirling の公式とは、n が大きいところで、
\ln(n\!)\simeq n(\ln(n) - 1)
というものです。n のオーダーに比例する項だけが問題になる時は、ここまでで十分となります。

計算では、高次の項を [1/n^7] のオーダーまでとった場合の差分も求めています。大体倍精度の場合は、この程度までとれば十分のはずです。

MODULE m_fact
IMPLICIT NONE
INTEGER, PARAMETER :: kd = SELECTED_REAL_KIND(8)
CONTAINS
REAL(kd) FUNCTION fact(n)
IMPLICIT NONE
INTEGER, INTENT(IN) :: n
INTEGER :: i
IF (n < 0) STOP
fact = 1.0_kd
DO i = 1, n
 fact = fact * i
END DO
RETURN
END FUNCTION fact
!---------------------------------------
REAL(kd) FUNCTION stirling(n)
IMPLICIT NONE
INTEGER, INTENT(IN) :: n
stirling = n * ( LOG(REAL(n, kd)) - 1.0_kd)
END FUNCTION stirling
!---------------------------------------
REAL(kd) FUNCTION stirling2(n)
IMPLICIT NONE
INTEGER, INTENT(IN) :: n
REAL(kd) :: pi, dn 
pi = 4.0_kd * ATAN(1.0_kd)
dn = REAL(n, kd)
stirling2 = 0.5_kd * LOG(2.0_kd * pi * n) &
          + 1.0_kd / (  12.0_kd * dn)     &
          - 1.0_kd / ( 360.0_kd * dn**3)  &
          + 1.0_kd / (1260.0_kd * dn**5)  &
          - 1.0_kd / (1680.0_kd * dn**7)
END FUNCTION stirling2
!---------------------------------------
END MODULE m_fact
!=======================================
PROGRAM factorial
USE m_fact
IMPLICIT NONE
INTEGER :: i
DO i = 1, 170
 PRINT '(i4, 2f20.14, es25.14)', i, LOG(fact(i)), stirling(i), &
                      LOG(fact(i)) - stirling(i) - stirling2(i)
END DO
STOP
END PROGRAM factorial


読み方は、スターリングの公式だが・・・w
持ってない読んでない。

こっちは読んだ。江青女史は半纏足www 
でも女子高生の方がいいよね!