This site is supported by donations to The OEIS Foundation.

# Approximations to the Factorial Function.

KEYWORDS: Factorial function, Gamma function, Stirling's formula, Stieltjes' formula, continued fraction approximation.

 Formula Numerator Denominator Stirling A001163 A001164 De Moivre A182935 A144618 Stieltjes A005146 A005147 Lanczos A090674 A090675 Nemes A181855 A181856 NemesG A182912 A182913 Wehmeier A182916 A182917 Gosper A182919 A182920 / New / A182914 A182915 / New / A277000 A277001 / New / A277002 A277003

## Introduction

The factorial of n is defined as

${\displaystyle n!=\int _{0}^{\infty }\ x^{n}\ e^{-x}}$.

The most famous approximation to this function is the De Moivre-Stirling formula

${\displaystyle S(n)={\sqrt {2\pi }}\ n^{n+1/2}\ e^{-n}}$.

To accelerate the asymptotic convergence of this formula several correction functions σ(n) were proposed.

${\displaystyle S^{*}(n)={\sqrt {2\pi }}\ n^{n+1/2}\ e^{-n}\ \sigma (n)}$.

For example the most simple one is σ(n) = (1 + 1/(12n)).

The formula we want to look at here starts from a different base structure than the De Moivre-Stirling formula. The base structure is more uniform; it does not mix n and n + 1/2.

${\displaystyle T(n)={\sqrt {2\pi }}\ (n+1/2)^{n+1/2}\ e^{-n-1/2}}$.

Again this formula has to be adjusted to become useful. Here we propose a correction function τ(n) which keeps the base structure of a half-shifted argument.

${\displaystyle T^{*}(n)={\sqrt {2\pi }}\ (n+1/2)^{n+1/2}\ e^{-n-1/2}\ \tau (n+1/2)^{n+1/2}}$

The details of this correction function will be discussed below. The formula was, to the best of my knowledge, introduced for the first time by the present author on Jan 9, 2007 on his homepage [1]. It turns out that the efficiency of the new formula is comparable to Stieltjes' famous continued fraction approximation to the factorial function.

## Stieltjes' approximation

In 1894 T. J. Stieltjes proposed an interesting approach to the factorial function.

${\displaystyle n!=\exp \left(\log(2\pi )/2-n+\left(n+{\frac {1}{2}}\right)\log(n)+p(n)\right).}$

For p(n) Stieltjes developed a continued fraction

${\displaystyle p(n)={\cfrac {a_{0}}{n+{\cfrac {a_{1}}{n+{\cfrac {a_{2}}{n+{\cfrac {a_{3}}{n+\ddots }}}}}}}}}$

The first few coefficients an are in A005146 and A005147.

 n 0 1 2 3 4 an 1 / 12 1 / 30 53 / 210 195 / 371 22999 / 22737

First we have to compute all the constants involved. This is far from being trivial. Stieltjes himself groaned about it: "Le calcul est trés pénible ... la loi de ces nombres paraît extrêmement compliqué". Using the algorithms of Rutishauser and Akiyama-Tanigawa we can compute nowadays the coefficients much more easily than Stieltjes could. More on this can be found here. A short Maple program is here.

## The new approximation

We will also consider a continued fraction, this time for τ(n)

${\displaystyle \tau (n)=n\times {\cfrac {c_{0}}{n+{\cfrac {c_{1}}{n+{\cfrac {c_{2}}{n+{\cfrac {c_{3}}{n+\ddots }}}}}}}}}$

The first few coefficients cn are in A182914 and A182915. (Note in passing that 45360 is the smallest integer with exactly 100 divisors (A163816).)

 n 0 1 2 3 4 cn 1 1 / 24 3 / 80 18029 / 45360 6272051 / 14869008

Again now the question is how to compute these coefficients. Of course this amounts here to give a precise definition of our approximation.

## The coefficients of the new formula

c := array(0..20, [1/24]): w := 'w':

approx := proc(deg, n) local N, A, k, p;
N := n + 1/2; A := N;
for k from deg by -1 to 0 do A := N + c[k] / A od;
p := N^2/A; ln(2*Pi)/2 + N*(ln(p)-1); exp(%) end:

coeffCF := proc(deg) global w, c;
local R, k; R := $1..0; if w = 'w' then R :=$1..deg; w := deg
elif deg > w then R := $w+1..deg; w := deg fi; for k in R do c[k] := a; convert(asympt(approx(k,m)-m!,m,2*k+4),polynom): c[k] := solve(%=0,a); print(coeff: ,k,c[k]); od end:  Assuming the coefficients are already computed use the function 'approx', otherwise call the function 'coeffCF' first. The parameter 'deg' is the degree of the formula requested. For example: coeffCF(5); approx(3,10): evalf(%,20); 3628800.0000007818727 approx(5,10): evalf(%,20); 3628800.0000000002108  ## The Wehmeier formulas The technique used in the last section clearly can also be applied to other basic structures for the approximation of the factorial. I learned this approach from Stefan Wehmeier in the Usenet group de.sci.mathematik in March 2006 dsm1. Wehmeier used it to derive a generalization of Gosper's approximation to n!. dsm2. He proceeded as follows: W := array(0..20): w := 'w': wehmeier := proc(deg, n) local A, k; A := n; for k from deg by -1 to 0 do A := A+W[k]/n^k od; sqrt(2*Pi*A)*n^n*exp(-n) end: coeffW := proc(deg) global w, W; local R, k; R :=$0..0;
if w = 'w' then R := $0..deg; w := deg elif deg > w then R :=$w..deg; w := deg fi;

for k in R do
W[k] := a;
convert(asympt(wehmeier(k,m)-m!,m,k+2),polynom):
W[k] := solve(%=0,a); print(coeff: ,k,W[k]);
od end:


With the same calling conventions as above we get

coeffW(5);
wehmeier(3,10): evalf(%,20);  3628799.9727503853301
wehmeier(5,10): evalf(%,20);  3628800.0002087858324

The first few coefficients Wn are in A182916 and A182917.

 n 0 1 2 3 4 Wn 1 / 6 1 / 72 -31 / 6480 -139 / 155520 9871 / 6531840

## The Nemes-G formulas

This month (March 2011) Gergö Nemes published a new approximation formula, "More Accurate Approximations for the Gamma Function", Thai Journal of Mathematics, Vol. 9,(1), p. 21-28. Nemes starts, like Wehmeier, from  Gosper's approximation. But then he continues with the following basic structure: A(n+1/4)sqrt(2Pi(n+1/6))(n/e)^n. Again this expansion can be easily derived with the same technique used above; only minor adjustments are necessary.

H := array(0..20): w := 'w':

nemesG := proc(deg, n) local A,N,k; A:=1; N := n+1/4;
for k from deg by -1 to 2 do A := A + H[k] / N^k od;
sqrt(2*Pi*(n+1/6))*n^n*exp(-n)*A end:

coeffH := proc(deg) global w, H; local R, k;
R := $2..0; if w = 'w' then R :=$2..deg; w := deg
elif deg > w then R := $w..deg; w := deg fi; for k in R do H[k] := a; convert(asympt(nemesG(k,m)-m!,m,k+1),polynom): H[k] := solve(%=0,a); print(coeff: ,k,H[k]); od end:  With the same calling conventions as above we get coeffH(5); nemesG(3,10): evalf(%,20); 3628800.4061837742630 nemesG(5,10): evalf(%,20); 3628799.9981087549165  The first few coefficients Hn are in A182912 and A182913.  n 0 1 2 3 4 Hn 1 0 1 / 144 -1 / 12960 -257 / 207360 ## The Gosper formulas There is even a simpler way to generalize Gosper's approximation. The basic formula is A(n)sqrt(2Pi(n+1/6))(n/e)^n. The same Ansatz as in the Nemes formula but without a shift of the argument. G := array(0..20): w := 'w': gosper := proc(deg, n) local A, k; A := 1; for k from deg by -1 to 2 do A := A + G[k]/n^k od; A*sqrt(2*Pi*(n+1/6))*(n/exp(1))^n end: coeffG := proc(deg) global w, G; local R, k; R :=$2..0;
if w = 'w' then R := $2..deg; w := deg elif deg > w then R :=$w..deg; w := deg fi;

for k in R do
G[k] := a;
convert(asympt(gosper(k,m)-m!,m,k+1),polynom):
G[k] := solve(%=0,a); print(coeff: ,k,G[k]);
od end:


With the same calling conventions as above we get

coeffG(8);
gosper(3,10): evalf(%,20);  3628799.9289952224556
gosper(5,10): evalf(%,20);  3628800.0001794192645 

The first few coefficients Gn (the coefficients are not (yet) in the OEIS database).

 n 0 1 2 3 4 Gn 1 0 1 / 144 -23 / 6480 5 / 41472

## The Stirling formulas

Stirling's approximation is the most famous approximation to the factorial function — for historic reasons. Today there is no reason to use these formulas for computational purposes. Formulas which perform much better are known.

h := proc(k) option remember; local j; if(k=0,1,
/(1+1/(k+1))) end:

coeffStirling := proc(n) option remember;
h(2*n)*2^n*pochhammer(1/2,n) end:

 n 0 1 2 3 4 Sn 1 1 / 12 1 / 288 -139 / 51840 -571 / 2488320

## A small benchmark

First we rewrite the formulas in a more compact way.

new := proc(n) local C, N, A; N := n + 1/2;
C := [1/24,3/80,18029/45360,6272051/14869008];
A := N^2/(N+C[1]/(N+C[2]/(N+C[3]/(N+C[4]/N))));
sqrt(2*Pi)*exp(N*(ln(A)-1)) end:

wehmeier := proc(n) local W, N, A; N := n + 1/6;
W := [1/72,-31/6480,-139/155520,9871/6531840];
A := N+W[1]/n+W[2]/n^2+W[3]/n^3+W[4]/n^4;
sqrt(2*Pi*A)*n^n*exp(-n) end:

nemes := proc(n) local N, H, A; N := n + 1/4;
H := [1/144,-1/12960,-257/207360,-53/2612736];
A := 1+H[1]/N^2+H[2]/N^3+H[3]/N^4+H[4]/N^5;
A*sqrt(2*Pi*(n+1/6))*n^n*exp(-n) end:

gosper := proc(n) local A, G;
G := [1/144,-23/6480,5/41472,4939/6531840];
A := 1+G[1]/n^2+G[2]/n^3+G[3]/n^4+G[4]/n^5;
A*sqrt(2*Pi*(n+1/6))*n^n*exp(-n) end:

stirling := proc(n) local A, S;
S := [1/12,1/288,-139/51840,-571/2488320];
A := 1+S[1]/n+S[2]/n^2+S[3]/n^3+S[4]/n^4;
A*sqrt(2*Pi)*n^(n+1/2)*exp(-n) end:


To measure the performance of the formulas we count the number of exact decimal digits.

edd := proc(a,t) evalf(-log[10](abs(1-a/t)),60) end:

for V in [100,1000,10000] do VF := V!:
printf("Stirling  %.5d! %6.1f \n",V,edd(stirling(V),VF));
printf("Nemes-G   %.5d! %6.1f \n",V,edd(nemes(V),VF));
printf("Wehmeier  %.5d! %6.1f \n",V,edd(wehmeier(V),VF));
printf("Gosper    %.5d! %6.1f \n",V,edd(gosper(V),VF));
printf("New       %.5d! %6.1f \n",V,edd(new(V),VF)); od:

 Exact decimal digits 100! 1000! 10000! Stirling 13.1 18.1 23.1 Nemes-G 15.2 21.2 27.2 Wehmeier 15.9 21.9 27.9 Gosper 17.5 23.1 29.1 New 21.5 30.5 39.5

Comparing the exact decimal digits of the Stirling, Gosper and the new formula over a larger range is shown in the plot below.

The peak in the blue edd-function is due to the fact that Gosper's approximation is exact at x = 67.0033148435486248...