スターリングの定理による階乗計算の正確さ

2024 1/05

2024.01.05 2024.01.05

$\displaystyle n!=\varGamma(n+1)\sim\sqrt{2\pi n}\Bigl(\frac{n}{e}\Bigr)^{n}(n\ge1)$

とその前に、スターリングの定理を計算するときにテイラー展開をしましたが、ここをもっと求める字数を増やすとさらに精度が増すことが考えられます。そこで、次の通り計算できます。

$n!\simeq\sqrt{2\pi n}\left(\dfrac{n}{e}\right)^n$

$n!\simeq\sqrt{2\pi n}\left(\dfrac{n}{e}\right)^n\left(1+\dfrac{1}{12n}\right)$

$n!\simeq\sqrt{2\pi n}\left(\dfrac{n}{e}\right)^n\left(1+\dfrac{1}{12n}+\dfrac{1}{288n^2}\right)$

$n!\simeq\sqrt{2\pi n}\left(\dfrac{n}{e}\right)^n\left(1+\dfrac{1}{12n}+\dfrac{1}{288n^2}-\dfrac{139}{51840n^3}\right)$

そこで、これらの式でどこまで階乗の値に迫れるか、調べてみます。

import math
import numpy as np
n = np.arange(1, 21)
stirling = np.sqrt(2*np.pi*n)*(n/np.e)**n
stirling2 = stirling*(1+1/(12*n))
stirling3 = stirling*(1+1/(12*n) + 1/(288*n**2))
stirling4 = stirling*(1+1/(12*n) + 1/(288*n**2)-139/(51840*n**3))
stirling5 = stirling*(1+1/(12*n) + 1/(288*n**2)-139/(51840*n**3) - 571/(2488320*n**4))
print('--+----------------------+----------------------+----------------------+----------------------+--------------------+')
print('n |   2pi*(n)*(n/e)**n   +-------- +1/12n-------+------ +1/288n^2------+-- -139/51840n^3------+----- factorial ----+')
print('--+----------------------+----------------------+----------------------+----------------------+--------------------+')

for i in n:
    print(f'{i:>2d}|{int(round(stirling[i-1],0)):>22d}|{int(round(stirling2[i-1],0)):>22d}|{int(round(stirling3[i-1],0)):>22d}|{int(round(stirling4[i-1],0)):>22d}|{math.factorial(i):>20d}|')

2列目から5列目はスターリングの定理によって計算したもので、右に行くほど精度が上がると考えられます。一番右側の列はfactorial関数を使って計算したものです。ちょっと見、かなり精度が高そうです。そこで、factorial関数で求めた正しい値に対するパーセンテージを計算します。

n=np.arange(1, 31)
stirling=np.sqrt(2*np.pi*n)*(n/np.e)**n
stirling2=stirling*(1+1/(12*n))
stirling3=stirling*(1+1/(12*n)+1/(288*n**2))
stirling4=stirling*(1+1/(12*n)+1/(288*n**2)-139/(51840*n**3))
stirling5=stirling*(1+1/(12*n)+1/(288*n**2)-139/(51840*n**3)-571/(2488320*n**4))
print('--+--------------+--------------+--------------+--------------+-----------------------------------+')
print('n |2pi*n*(n/e)**n+- +1/12n------+- +1/288n^2---+ -139/51840n^3+----       factorial          -----+')
print('--+--------------+--------------+--------------+--------------+-----------------------------------+')

for i in n:
    fact=math.factorial(i)
    print(f'{i:>2d}|{stirling[i-1]/fact:>14.8%}|{stirling2[i-1]/fact:>14.8%}|{stirling3[i-1]/fact:>14.8%}|{stirling4[i-1]/fact:>14.8%}|{factorial(i):>35.0f}|')