I tried the accuracy of three Stirling's approximations in python

I've been worried about the accuracy of Stirling's approximation for a long time, but I tried it. It's a simple verification, but the results are amazing. The accuracy is so different depending on the calculation method.

The simplest approximation is

$ \log n! = n \log n - n +O(\log n)$ It will be.

For the gamma function $ \ Gamma (z) $, $ \ gamma (n) = (n-1)! $ holds for the positive integer $ n $. The following equation holds as a computer approximation of the gamma function.

2\ln \gamma(z) \sim \ln 2\pi -\ln z + z \left ( 2\ln z+\ln \left(z \sinh \frac{1}{z}+\frac{1}{810z^6} \right)-2\right ) It has a precision of 8 decimal places. Also, as a simpler approximation

\ln \gamma(z) \sim 0.5\left(\ln 2\pi -\ln z\right) + z \left ( \ln \left(z+ \frac{1}{12z-\frac{1}{10z}}-1\right )\right) there is.

Try these precisions in Python.

%matplotlib inline
from scipy.special import factorial
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import statsmodels.api as sm
import itertools
for i in range(1,10,1):
    f=factorial(i)
    ff=1
    for ii in range(1,i+1,1):
        ff=ff*ii
    fff=np.exp(i*np.log(i)-i)
    j=i+1
    ffff=0.5*(np.log(2*np.pi)-np.log(j))+j*(np.log(j+1/(12*j-1/10/j))-1)
    fffff=0.5*(np.log(2*np.pi)-np.log(j)+j*(2*np.log(j)+np.log(j*np.sinh(1/j)+1/810/j**6)-2))
    print("{0:10d}{1:10.0f}{2: 10.0f}{3:15.8f}{4:18.8f}{5:18.8f}".format(i,f,ff,fff,np.exp(ffff),np.exp(fffff)))

         1         1         1     0.36787944        0.99998309        1.00000267
         2         2         2     0.54134113        1.99999509        2.00000037
         3         6         6     1.34425085        5.99999637        6.00000016
         4        24        24     4.68880356       23.99999516       24.00000014
         5       120       120    21.05608437      119.99999017      120.00000020
         6       720       720   115.64866155      719.99997253      720.00000040
         7      5040      5040   750.97400956     5039.99990097     5040.00000112
         8     40320     40320  5628.12896825    40319.99955910    40320.00000396
         9    362880    362880 47811.48664666   362879.99765201   362880.00001709

The accuracy of these results was amazing. Let's increase the number of digits of $ n $.

for i in range(100,1000,100):
    f=factorial(i)
    ff=1
    for ii in range(1,i+1,1):
        ff=ff*ii
    fff=np.exp(i*np.log(i)-i)
    j=i+1
    ffff=0.5*(np.log(2*np.pi)-np.log(j))+j*(np.log(j+1/(12*j-1/10/j))-1)
    fffff=0.5*(np.log(2*np.pi)-np.log(j)+j*(2*np.log(j)+np.log(j*np.sinh(1/j)+1/810/j**6)-2))
    print(i,f,np.exp(ffff),np.exp(fffff))

100 9.332621544394415e+157 9.332621544394225e+157 9.332621544394755e+157
200 inf inf inf
300 inf inf inf
400 inf inf inf
500 inf inf inf
600 inf inf inf
700 inf inf inf
800 inf inf inf
900 inf inf inf

ff and fff are omitted because they cannot be output properly. Up to this point, the calculation limits of the calculation approximation formula are shown. But with ln Factorial, you can calculate farther.

for i in range(1000,10000,1000):
    f=factorial(i)
    ff=1
    for ii in range(1,i+1,1):
        ff=ff*ii
    fff=np.exp(i*np.log(i)-i)
    j=i+1
    ffff=0.5*(np.log(2*np.pi)-np.log(j))+j*(np.log(j+1/(12*j-1/10/j))-1)
    fffff=0.5*(np.log(2*np.pi)-np.log(j)+j*(2*np.log(j)+np.log(j*np.sinh(1/j)+1/810/j**6)-2))
    print("{0:10d}{1: 18.8f}{2:18.8f}".format(i,ffff,fffff))

      1000     5912.12817849     5912.12817849
      2000    13206.52435051    13206.52435051
      3000    21024.02485305    21024.02485305
      4000    29181.26454459    29181.26454459
      5000    37591.14350888    37591.14350888
      6000    46202.35719906    46202.35719906
      7000    54981.00377941    54981.00377941
      8000    63902.98711266    63902.98711266
      9000    72950.29014459    72950.29014459

Calculations that take the logarithm of the factorial can be calculated up to a fairly large $ n $. This can greatly extend the computational limit, such as when using multiplicity with very small probabilities.

reference: Wiki: Stirling's approximation