Stirling’s formula and its asymptotic series

Shows that sqrt(2 pi n) (n/e)^n approximates n! with relative error about 1/(12n), and that Stirling’s full series for log Gamma(x) is asymptotic: at small x, adding terms first helps and then hurts.

import math

import matplotlib.pyplot as plt
import numpy as np

from mathematicskit.special_functions import log_gamma_function, stirling_factorial, stirling_log_gamma

The leading term: n! ~ sqrt(2 pi n) (n/e)^n

for n in (1, 5, 10, 50, 100):
    ratio = math.factorial(n) / stirling_factorial(float(n))
    print(f"n = {n:3d}: n!/Stirling = {ratio:.6f}, 1 + 1/(12n) = {1 + 1 / (12 * n):.6f}")

n =   1: n!/Stirling = 1.084438, 1 + 1/(12n) = 1.083333
n =   5: n!/Stirling = 1.016784, 1 + 1/(12n) = 1.016667
n =  10: n!/Stirling = 1.008365, 1 + 1/(12n) = 1.008333
n =  50: n!/Stirling = 1.001668, 1 + 1/(12n) = 1.001667
n = 100: n!/Stirling = 1.000834, 1 + 1/(12n) = 1.000833

Error of the truncated series versus number of terms

terms = np.arange(0, 16)
fig, ax = plt.subplots()
for x in (1.0, 2.0, 5.0):
    errors = [abs(stirling_log_gamma(x, terms=int(k)) - log_gamma_function(x)) for k in terms]
    ax.semilogy(terms, errors, "o-", label=f"x = {x:g}")
    print(f"x = {x:g}: best with {int(np.argmin(errors))} terms, error {min(errors):.2e}")
ax.set_xlabel("number of Bernoulli correction terms")
ax.set_ylabel(r"$|\,\mathrm{series} - \ln\Gamma(x)\,|$")
ax.set_title("Stirling's series is asymptotic, not convergent")
ax.legend()

x = 1: best with 3 terms, error 2.88e-04
x = 2: best with 7 terms, error 3.90e-07
x = 5: best with 15 terms, error 8.88e-16

<matplotlib.legend.Legend object at 0x11c4ff4d0>