統計力学などでは、よく階乗の対数を Stirling の公式で近似しますが、どのくらいの精度で近似になっているのか見てみます。
ここで Stirling の公式とは、n が大きいところで、
というものです。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
でも女子高生の方がいいよね!