スターリングの公式等の近似精度をn!とlog(n!)の両方で検証してみた

はじめに

数値計算において、非常に大きな$n!$や$\log n!$などの階乗の計算は時間がかかる場合があります。そのため、階乗の計算にはスターリングの公式などを利用して近似的に求めることがあります。スターリングの公式を高次の項まで展開したスターリング級数もあり、より精度を向上させることができます。
本記事では、スターリング級数の精度と計算速度について、数値計算を通して検証していきます。また、Wikipedia情報ですが、スターリングの公式以外の近似式があり、それらとの比較も行っていきます。

スターリングの公式

スターリングの公式は、階乗$n!$を近似するための公式です。一般的な形式は

n! \sim \sqrt{2πn}\left(\frac{n}{e}\right )^n

となります。ここで、$π$は円周率、$e$はネイビア数です。

ただし、この公式自体だけでは精度が十分でないことがあり、追加の項を含めたスターリング級数は、

n!\sim {\sqrt {2\pi n}}\left({\frac {n}{e}}\right)^{n}\left(1+{\frac {1}{12n}}+{\frac {1}{288n^{2}}}-{\frac {139}{51840n^{3}}}-{\frac {571}{2488320n^{4}}+\frac{163879}{209018880n^{5}}+\frac{5246819}{75246796800n^{6}}}-\frac{4483131259}{902961561600n^{7}}+\cdots \right)

となります。(この級数は、分子: https://oeis.org/A001163 、分母: https://oeis.org/A001164 に載っておりますので興味がある方はご覧ください。)

また、対数の場合のスターリングの公式は、

\ln n!\sim n\ln n-n

となり、スターリング級数は、

\ln n!\sim n\ln n-n+{\frac {1}{2}}\ln 2\pi n+{\frac {1}{12n}}-{\frac {1}{360n^{3}}}+{\frac {1}{1260n^{5}}}-{\frac {1}{1680n^{7}}}+{\frac {1}{1188n^{9}}}-{\frac {691}{360360n^{11}}}+\cdots

と表されます。
この対数のスターリング級数は、以下から求められますので、ご参考にしてください。

\ln \Gamma (z)={\frac {\log 2\pi }{2}}-z+\left(z-{\frac {1}{2}}\right)\log z+\sum _{k=1}^{m}{\frac {B_{2k}}{(2k)(2k-1)z^{2k-1}}}+O\left(z^{-2m}\right)

ここで、$\Gamma (z)$は、ガンマ関数、$B$はベルヌーイ数です。また、このガンマ関数の式をテイラー展開することで、上式のn!のスターリング級数の各項を計算することができます。詳しくは

をご参照ください。

その他の近似式

スターリングの公式以外に、比較的高速で精度の伴った数値計算が行えるような式として、Wikipediaに載っていることもあり、以下2つ含め本記事の後半で比較します。

Robert H. (2002)の近似式

Robert H. (2002)によって考案された近似式で、

n! \approx {\sqrt {\frac {2\pi }{n+1}}}\left({\frac {n+1}{e}}{\sqrt {(n+1)\sinh {\frac {1}{n+1}}+{\frac {1}{810(n+1)^{6}}}}}\right)^{n+1}

で与えられ、対数で

\ln n!\approx \frac{1}{2} \left(\ln(2\pi )- \ln (n+1)+(n+1)\left(2\ln (n+1)+\ln \left((n+1)\sinh {\frac {1}{n+1}}+{\frac {1}{810(n+1)^{6}}}\right)-2\right) \right)

と表されます。

詳しくは、

を参照ください。

Gergő N. (2007)の近似式

Gergő N. (2007)によって考案された近似式で、

n!\approx {\sqrt {\frac {2\pi }{n+1}}}\left({\frac {1}{e}}\left(n+1+{\frac {1}{12(n+1)-{\frac {1}{10(n+1)}}}}\right)\right)^{n+1}

で与えられ、対数で

\ln n!\approx {\tfrac {1}{2}}\left(\ln(2\pi )-\ln (n+1)\right)+(n+1)\left(\ln \left((n+1)+{\frac {1}{12(n+1)-{\frac {1}{10(n+1)}}}}\right)-1\right)

と表されます。

詳しくは、

を参照ください。

数値実験

Google Colabで作成した本記事のコードは、こちらにあります。

数値実験では、次数を変えた場合の精度をグラフで比較します。速度を計測します。

各種インポート

import numpy as np
import math
import matplotlib.pyplot as plt
import time

スターリング級数 n!

まず、n!のスターリング級数の関数を書きます。

n!のスターリング級数

def calc_stirling_series(n):
    stirling_series = np.zeros((8,))
    stirling_series[0] = math.sqrt(2 * math.pi * n) * (n / math.e) ** n
    stirling_series[1] = stirling_series[0] + stirling_series[0] * 1 / (12 * n)
    stirling_series[2] = stirling_series[1] + stirling_series[0] * 1 / (288 * n ** 2)
    stirling_series[3] = stirling_series[2] - stirling_series[0] * 139 / (51840 * n ** 3)
    stirling_series[4] = stirling_series[3] - stirling_series[0] * 571 / (2488320 * n ** 4)
    stirling_series[5] = stirling_series[4] + stirling_series[0] * 163879 / (209018880 * n ** 5)
    stirling_series[6] = stirling_series[5] + stirling_series[0] * 5246819 / (75246796800 * n ** 6)
    stirling_series[7] = stirling_series[6] + stirling_series[0] * 4483131259 / (902961561600 * n ** 7)
    return stirling_series

n=2~100まで1刻みでグラフに表示します。ここで、n = 0は、スターリング級数では分母に0を代入できないので使っていません。また、n = 1は、分母をn!で割って差分を表示するときに、ln1!=0で割ることになってしまうので省きました。

n!をグラフに表示

n_range = list(range(2, 100))
factorials = np.zeros((len(n_range), 1))
stirling_seriess = np.zeros((len(n_range), 8))

for i, n in enumerate(n_range):
    factorials[i] = math.factorial(n)
    stirling_seriess[i] = calc_stirling_series(n)


fig, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(16, 6))
stirling_legend = ['term 1', 'term 2', 'term 3', 'term 4', 'term 5', 'term 6', 'term 7', 'term 8']

ax0.set_title('stirling_series vs. factorial')
ax0.plot(n_range, factorials)
ax0.plot(n_range, stirling_seriess)
ax0.set_yscale('log')
ax0.set_xlabel('n')
ax0.set_ylabel('value')
ax0.legend(['factorial']+stirling_legend)

ax1.set_title('abs(stirling_series - factorial)')
ax1.plot(n_range, abs(factorials-stirling_seriess))
ax1.set_yscale('log')
ax1.set_xlabel('n')
ax1.set_ylabel('value')
ax1.legend(stirling_legend)

ax2.set_title('abs(stirling_series - factorial) / factorial')
ax2.plot(n_range, abs(factorials-stirling_seriess)/factorials)
ax2.set_yscale('log')
ax2.set_xlabel('n')
ax2.set_ylabel('value')
ax2.legend(stirling_legend)

plt.show()

図の説明
左: n!とスターリング級数で近似した場合の結果
中央: n!とスターリング級数の近似の差分の絶対値 (差分)
右: n!とスターリング級数の近似の差分の絶対値をn!で割ったもの (差分の増減率)
判例: スターリング級数の項を変えた時の比較

図を見ると、左図では全ての級数に対してn!に重なっていることが確認できますが、中央図のように差分を見ると、やはり級数の項が増えるほどn!に近づく傾向が見られます。
n!は、指数関数的に増加するため、中央図だけではどの程度の割合の誤差を含んでいるのかわかりにくいので、n!で差分を割った右図を見てみと、n!が小さいところでは、誤差が占める割合が大きいが、nが大きくなると誤差の占める割合が小さくなることがわかります。
また、項が6, 7, 8ではほとんど同じ分布になっていて縦に振動しているように見えますが、これはfloat64型の精度が15桁程度であるため、丸め込まれて同じような分布になってしまっていると考えられます。

スターリング級数 ln(n!)

同様に対数でも比較します。

ln(n!)のグラフに書くためのコード

ln(n!)のグラフ表示

n_range = list(range(2, 100))
log_factorials = np.zeros((len(n_range), 1))
log_stirling_seriess = np.zeros((len(n_range), 8))

for i, n in enumerate(n_range):
    log_factorials[i] = math.log(math.factorial(n))
    log_stirling_seriess[i] = calc_log_stirling_series(n)


fig, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(16, 6))
stirling_legend = ['term 1', 'term 2', 'term 3', 'term 4', 'term 5', 'term 6', 'term 7','term 8']

ax0.set_title('stirling_series vs. factorial')
ax0.plot(n_range, log_factorials)
ax0.plot(n_range, log_stirling_seriess)
ax0.set_yscale('log')
ax0.set_xlabel('n')
ax0.set_ylabel('value')
ax0.legend(['factorial']+stirling_legend)

ax1.set_title('abs(stirling_series - factorial)')
ax1.plot(n_range, abs(log_factorials-log_stirling_seriess))
ax1.set_yscale('log')
ax1.set_xlabel('n')
ax1.set_ylabel('value')
ax1.legend(stirling_legend)

ax2.set_title('abs(stirling_series - factorial) / factorial')
ax2.plot(n_range, abs(log_factorials-log_stirling_seriess)/log_factorials)
ax2.set_yscale('log')
ax2.set_xlabel('n')
ax2.set_ylabel('value')
ax2.legend(stirling_legend)

plt.show()

グラフの説明はn!の場合と同じため割愛します。n!で全て話してしまって、特に話すことはないです。

その他の近似式との比較

結論から申し上げると、精度としてはRobertさんがスターリング級数の4項目と5項目の間くらいで、Gergőさんが4項目程度の結果になりました。

その他の近似式 n!

n!の近似式のグラフに書くためのコード

n!の近似式

def calc_factorial_approx_by_Robert(n):
    # Suggested by Robert H. (2002)
    # citing: rskey.org/CMS/index.php/the-library/11
    factorial_approx = math.sqrt(2 * math.pi / (n+1)) * (((n+1) / math.e) * math.sqrt((n+1) * math.sinh(1 / (n+1)) + (1 / (810 * (n-1) ** 6)))) ** (n+1)
    return factorial_approx

def calc_factorial_approx_by_Gergo(n):
    # Suggested by Gergő N. (2007)
    # Reference: Nemes, Gergő (2010), "New asymptotic expansion for the Gamma function", Archiv der Mathematik, 95 (2): 161–169, doi:10.1007/s00013-010-0146-9
    factorial_approx = math.sqrt(2 * math.pi / (n+1)) * ((n+1 + (1 / (12 * (n+1) - 1 / (10 * (n+1))))) / math.e) ** (n+1)
    return factorial_approx


n_range = list(range(2, 100))
factorials = np.zeros((len(n_range), 1))
stirling_seriess = np.zeros((len(n_range), 8))
factorial_approx_by_Robert = np.zeros((len(n_range), 1))
factorial_approx_by_Gergo = np.zeros((len(n_range), 1))

for i, n in enumerate(n_range):
    factorials[i] = math.factorial(n)
    stirling_seriess[i] = calc_stirling_series(n)
    factorial_approx_by_Robert[i] = calc_factorial_approx_by_Robert(n)
    factorial_approx_by_Robert[i] = calc_factorial_approx_by_Robert(n)
    factorial_approx_by_Gergo[i] = calc_factorial_approx_by_Gergo(n)


fig, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(16, 6))
stirling_terms = [4, 5]

ax0.set_title('approx_formula vs. factorial')
ax0.plot(n_range, factorials, label='factorial')
ax0.plot(n_range, stirling_seriess[:, stirling_terms[0]], label=f'{stirling_terms[0]=}')
ax0.plot(n_range, stirling_seriess[:, stirling_terms[1]], label=f'{stirling_terms[1]=}')
ax0.plot(n_range, factorial_approx_by_Robert, label='Robert')
ax0.plot(n_range, factorial_approx_by_Gergo, label='Gergo')
ax0.set_yscale('log')
ax0.set_xlabel('n')
ax0.set_ylabel('value')
ax0.legend()

ax1.set_title('abs(approx_formula - factorial)')
ax1.plot(n_range, abs(factorials-stirling_seriess)[:, stirling_terms[0]], label=f'{stirling_terms[0]=}')
ax1.plot(n_range, abs(factorials-stirling_seriess)[:, stirling_terms[1]], label=f'{stirling_terms[1]=}')
ax1.plot(n_range, abs(factorials-factorial_approx_by_Robert), label='Robert')
ax1.plot(n_range, abs(factorials-factorial_approx_by_Gergo), label='Gergo')
ax1.set_yscale('log')
ax1.set_xlabel('n')
ax1.set_ylabel('value')
ax1.legend()

ax2.set_title('abs(approx_formula - factorial) / factorial')
ax2.plot(n_range, abs((factorials-stirling_seriess)/factorials)[:, stirling_terms[0]], label=f'{stirling_terms[0]=}')
ax2.plot(n_range, abs((factorials-stirling_seriess)/factorials)[:, stirling_terms[1]], label=f'{stirling_terms[1]=}')
ax2.plot(n_range, abs(factorials-factorial_approx_by_Robert)/factorials, label='Robert')
ax2.plot(n_range, abs(factorials-factorial_approx_by_Gergo)/factorials, label='Gergo')
ax2.set_yscale('log')
ax2.set_xlabel('n')
ax2.set_ylabel('value')
ax2.legend()

plt.show()

その他の近似式 ln(n!)

ln(n!)の近似式のグラフに書くためのコード

ln(n!)の近似式

def calc_log_factorial_approx_by_Robert(n):
    # Suggested by Robert H. (2002)
    # Reference: rskey.org/CMS/index.php/the-library/11
    log_factorial_approx = 0.5 * (math.log(2 * math.pi) - math.log(n+1)) + 0.5 * (n+1) * (2 * math.log(n+1) + math.log((n+1) * math.sinh(1 / (n+1)) + (1 / (810 * (n+1) ** 6))) - 2)
    return log_factorial_approx

def calc_log_factorial_approx_by_Gergo(n):
    # Suggested by Gergő N. (2007)
    # Reference: Nemes, Gergő (2010), "New asymptotic expansion for the Gamma function", Archiv der Mathematik, 95 (2): 161–169, doi:10.1007/s00013-010-0146-9
    log_factorial_approx = 0.5 * (math.log(2 * math.pi) - math.log(n+1)) + (n+1) * (math.log((n+1) + (1 / (12 * (n+1) - 1 / (10 * (n+1))))) - 1)
    return log_factorial_approx


n_range = list(range(2, 100))
log_factorials = np.zeros((len(n_range), 1))
stirling_seriess = np.zeros((len(n_range), 8))
log_factorial_approx_by_Robert = np.zeros((len(n_range), 1))
log_factorial_approx_by_Gergo = np.zeros((len(n_range), 1))

for i, n in enumerate(n_range):
    log_factorials[i] = math.log(math.factorial(n))
    log_stirling_seriess[i] = calc_log_stirling_series(n)
    log_factorial_approx_by_Robert[i] = calc_log_factorial_approx_by_Robert(n)
    log_factorial_approx_by_Robert[i] = calc_log_factorial_approx_by_Robert(n)
    log_factorial_approx_by_Gergo[i] = calc_log_factorial_approx_by_Gergo(n)


fig, (ax0, ax1, ax2) = plt.subplots(1, 3, figsize=(16, 6))
stirling_terms = [4, 5]

ax0.set_title('approx_formula vs. factorial')
ax0.plot(n_range, log_factorials, label='factorial')
ax0.plot(n_range, log_stirling_seriess[:, stirling_terms[0]], label=f'{stirling_terms[0]=}')
ax0.plot(n_range, log_stirling_seriess[:, stirling_terms[1]], label=f'{stirling_terms[1]=}')
ax0.plot(n_range, log_factorial_approx_by_Robert, label='Robert')
ax0.plot(n_range, log_factorial_approx_by_Gergo, label='Gergo')
ax0.set_yscale('log')
ax0.set_xlabel('n')
ax0.set_ylabel('value')
ax0.legend()

ax1.set_title('abs(approx_formula - factorial)')
ax1.plot(n_range, abs(log_factorials-log_stirling_seriess)[:, stirling_terms[0]], label=f'{stirling_terms[0]=}')
ax1.plot(n_range, abs(log_factorials-log_stirling_seriess)[:, stirling_terms[1]], label=f'{stirling_terms[1]=}')
ax1.plot(n_range, abs(log_factorials-log_factorial_approx_by_Robert), label='Robert')
ax1.plot(n_range, abs(log_factorials-log_factorial_approx_by_Gergo), label='Gergo')
ax1.set_yscale('log')
ax1.set_xlabel('n')
ax1.set_ylabel('value')
ax1.legend()

ax2.set_title('abs(approx_formula - factorial) / factorial')
ax2.plot(n_range, abs((log_factorials-log_stirling_seriess)/log_factorials)[:, stirling_terms[0]], label=f'{stirling_terms[0]=}')
ax2.plot(n_range, abs((log_factorials-log_stirling_seriess)/log_factorials)[:, stirling_terms[1]], label=f'{stirling_terms[1]=}')
ax2.plot(n_range, abs(log_factorials-log_factorial_approx_by_Robert)/log_factorials, label='Robert')
ax2.plot(n_range, abs(log_factorials-log_factorial_approx_by_Gergo)/log_factorials, label='Gergo')
ax2.set_yscale('log')
ax2.set_xlabel('n')
ax2.set_ylabel('value')
ax2.legend()

plt.show()

処理速度

処理速度も見てみましょう。
一応、先にコード書いたスターリング級数が逐次代入していたので、平等に比較するために一つ一つを分けて関数に書いておきます。一応、先に書いたRobertさんとGergőさんの近似式も改めて書いておきます。

ここで、processing_func_speed(func, n, iteration)は、処理速度を測るための関数で、調べたい関数func、その階乗n、試行回数iterationを入れて計測するための関数です。

処理速度に使う関数

# 処理時間計測用の関数
def processing_func_speed(func, n, iteration):
    start_time = time.time()
    for _ in range(iteration):
        func(n)
    end_time = time.time()
    elapsed_avg_time = (end_time - start_time) / iteration
    return elapsed_avg_time

# n!用の関数
def calc_factorial(n):
    factorial = math.factorial(n)
    return factorial

def calc_stirling_series_term_1(n):
    stirling = math.sqrt(2 * math.pi * n) * (n / math.e) ** n
    return stirling

def calc_stirling_series_term_2(n):
    stirling = math.sqrt(2 * math.pi * n) * (n / math.e) ** n * (1 + 1 / (12 * n))
    return stirling

def calc_stirling_series_term_3(n):
    stirling = math.sqrt(2 * math.pi * n) * (n / math.e) ** n * (1 + 1 / (12 * n)\
            + 1 / (288 * n ** 2))
    return stirling

def calc_stirling_series_term_4(n):
    stirling = math.sqrt(2 * math.pi * n) * (n / math.e) ** n * (1 + 1 / (12 * n)\
            + 1 / (288 * n ** 2) - 139 / (51840 * n ** 3))
    return stirling

def calc_stirling_series_term_5(n):
    stirling = math.sqrt(2 * math.pi * n) * (n / math.e) ** n * (1 + 1 / (12 * n)\
            + 1 / (288 * n ** 2) - 139 / (51840 * n ** 3) - 571 / (2488320 * n ** 4))
    return stirling
def calc_stirling_series_term_6(n):
    stirling = math.sqrt(2 * math.pi * n) * (n / math.e) ** n * (1 + 1 / (12 * n)\
            + 1 / (288 * n ** 2) - 139 / (51840 * n ** 3) - 571 / (2488320 * n ** 4)\
            + 163879 / (209018880 * n ** 5))
    return stirling

def calc_stirling_series_term_7(n):
    stirling = math.sqrt(2 * math.pi * n) * (n / math.e) ** n * (1 + 1 / (12 * n)\
            + 1 / (288 * n ** 2) - 139 / (51840 * n ** 3) - 571 / (2488320 * n ** 4)\
            + 163879 / (209018880 * n ** 5) + 5246819 / (75246796800 * n ** 6))
    return stirling

def calc_factorial_approx_by_Robert(n):
    # Suggested by Robert H. (2002)
    # citing: rskey.org/CMS/index.php/the-library/11
    factorial_approx = math.sqrt(2 * math.pi / (n+1)) * (((n+1) / math.e) * math.sqrt((n+1) * math.sinh(1 / (n+1)) + (1 / (810 * (n-1) ** 6)))) ** (n+1)
    return factorial_approx

def calc_factorial_approx_by_Gergo(n):
    # Suggested by Gergő N. (2007)
    # Reference: Nemes, Gergő (2010), "New asymptotic expansion for the Gamma function", Archiv der Mathematik, 95 (2): 161–169, doi:10.1007/s00013-010-0146-9
    factorial_approx = math.sqrt(2 * math.pi / (n+1)) * ((n+1 + (1 / (12 * (n+1) - 1 / (10 * (n+1))))) / math.e) ** (n+1)
    return factorial_approx

# log(n!)用の関数
def calc_log_factorial(n):
    log_factorial = math.log(math.factorial(n))
    return log_factorial

def calc_log_stirling_series_term_1(n):
    log_stirling = n * math.log(n) - n
    return log_stirling

def calc_log_stirling_series_term_2(n):
    log_stirling = n * math.log(n) - n + 0.5 * math.log(2 * math.pi * n)
    return log_stirling

def calc_log_stirling_series_term_3(n):
    log_stirling = n * math.log(n) - n + 0.5 * math.log(2 * math.pi * n)\
            + 1 / (12 * n)
    return log_stirling

def calc_log_stirling_series_term_4(n):
    log_stirling = n * math.log(n) - n + 0.5 * math.log(2 * math.pi * n)\
            + 1 / (12 * n) - 1 / (360 * n**3)
    return log_stirling

def calc_log_stirling_series_term_5(n):
    log_stirling = n * math.log(n) - n + 0.5 * math.log(2 * math.pi * n)\
            + 1 / (12 * n) - 1 / (360 * n**3) + 1 / (1260 * n**5)
    return log_stirling

def calc_log_stirling_series_term_6(n):
    log_stirling = n * math.log(n) - n + 0.5 * math.log(2 * math.pi * n)\
            + 1 / (12 * n) - 1 / (360 * n**3) + 1 / (1260 * n**5)\
            - 1 / (1680 * n**7)
    return log_stirling

def calc_log_stirling_series_term_7(n):
    log_stirling = n * math.log(n) - n + 0.5 * math.log(2 * math.pi * n)\
            + 1 / (12 * n) - 1 / (360 * n**3) + 1 / (1260 * n**5)\
            - 1 / (1680 * n**7) + 1 / (1188 * n**9)
    return log_stirling

def calc_log_factorial_approx_by_Robert(n):
    # Suggested by Robert H. (2002)
    # Reference: rskey.org/CMS/index.php/the-library/11
    log_factorial_approx = 0.5 * (math.log(2 * math.pi) - math.log(n+1)) + 0.5 * (n+1) * (2 * math.log(n+1) + math.log((n+1) * math.sinh(1 / (n+1)) + (1 / (810 * (n+1) ** 6))) - 2)
    return log_factorial_approx

def calc_log_factorial_approx_by_Gergo(n):
    # Suggested by Gergő N. (2007)
    # Reference: Nemes, Gergő (2010), "New asymptotic expansion for the Gamma function", Archiv der Mathematik, 95 (2): 161–169, doi:10.1007/s00013-010-0146-9
    log_factorial_approx = 0.5 * (math.log(2 * math.pi) - math.log(n+1)) + (n+1) * (math.log((n+1) + (1 / (12 * (n+1) - 1 / (10 * (n+1))))) - 1)
    return log_factorial_approx

調べる階乗は、n_range = [2, 4, 8, 16, 32, 64, 128]、試行回数はiteration = 3000000で、計測してみます。

n!の近似式の処理速度

n!の処理速度を計測するためのコード

n!の近似式の処理速度比較

funcs = [calc_factorial, calc_stirling_series_term_1, 
         calc_stirling_series_term_2, calc_stirling_series_term_3,
         calc_stirling_series_term_4, calc_stirling_series_term_5,
         calc_stirling_series_term_6, calc_stirling_series_term_7,
         calc_factorial_approx_by_Robert, calc_factorial_approx_by_Gergo]
function_legend = [func.__name__ for func in funcs]

n_range = [2**n for n in range(1, 8)]
times = np.zeros((len(n_range), len(funcs)))
iteration = 300000

for i, n in enumerate(n_range):
    for j, func in enumerate(funcs):
        times[i, j] = processing_func_speed(func, n, iteration)

fig, ax = plt.subplots()
ax.plot(n_range, times, '.-')
ax.legend(function_legend, loc='upper left', bbox_to_anchor=(1, 1))
ax.set_ylim(0,)
ax.set_xlabel('n')
ax.set_ylabel('time')
plt.show()

純粋な階乗を計算式factorialは、nに比例して計算時間が大きくなっています。
Robertさんの近似式の精度は、スタリーング級数のtermが4, 5の間でしたが、n=16, 32, 128では、それよりも処理時間が小さくなっていますね。
また、Gergőさんの近似式も、スタリーング級数のtermが4あたりでしたが、こちらもn=16, 32, 128では、それよりも処理時間が小さくなっていますね。

n=64でのRobertさんとGergőさんの処理速度の2倍程度の増加の原因が気になりますが、概ね計算速度が向上しているのがわかります。

ln(n!)の近似式の処理速度

ln(n!)の処理速度を計測するためのコード

ln(n!)の近似式の処理速度比較

funcs = [calc_log_factorial, calc_log_stirling_series_term_1, 
         calc_log_stirling_series_term_2, calc_log_stirling_series_term_3,
         calc_log_stirling_series_term_4, calc_log_stirling_series_term_5,
         calc_log_stirling_series_term_6, calc_log_stirling_series_term_7,
         calc_log_factorial_approx_by_Robert, calc_log_factorial_approx_by_Gergo]
function_legend = [func.__name__ for func in funcs]

n_range = [2**n for n in range(1, 8)]
times = np.zeros((len(n_range), len(funcs)))
iteration = 300000

for i, n in enumerate(n_range):
    for j, func in enumerate(funcs):
        times[i, j] = processing_func_speed(func, n, iteration)

fig, ax = plt.subplots()
ax.plot(n_range, times, '.-')
ax.legend(function_legend, loc='upper left', bbox_to_anchor=(1, 1))
ax.set_ylim(0,)
ax.set_xlabel('n')
ax.set_ylabel('time')
plt.show()

Robertさんの近似式の精度は、スタリーング級数のtermが4, 5あたりと比較すると、それよりも多く時間がかかってしまっています。
また、Gergőさんの近似式も、スタリーング級数のtermが4あたりと比較すると、n=16, 32, 128では、ほとんど処理速度が変わらず、n=64に至っては、処理速度の低下が見られます。

まとめ

スターリング級数とその他の近似式について紹介しました。この記事が目的に合わせた適切な式を選択の一助になれば幸いです。

参考資料

スターリング級数
分子:

分母: