This site is supported by donations to The OEIS Foundation.

User:Peter Luschny/FactorialFunction

From OeisWiki

Jump to: navigation, search

Approximations to the
Factorial Function.

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

Concerned with sequences:
Formula Numerator Denominator
Stirling A001163A001164
De Moivre A182935A144618
Stieltjes A005146A005147
Lanczos A090674A090675
Nemes A181855A181856
NemesG A182912A182913
Wehmeier A182916A182917
Gosper A182919A182920
/ New / A182914A182915
/ New / A277000A277001
/ New / A277002A277003



The factorial of n is defined as

  n! = \int_0^{\infty} \ x^n \ e^{-x}.

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

 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.

 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.

 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.

 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.

 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


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)

\tau(n) = n \times \cfrac{c_0}{n+

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;
  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:

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;
  W[k] := solve(%=0,a); print(`coeff: `,k,W[k]);
od end:

With the same calling conventions as above we get

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;
  H[k] := solve(%=0,a); print(`coeff: `,k,H[k]);
od end:

With the same calling conventions as above we get

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;
  G[k] := solve(%=0,a); print(`coeff: `,k,G[k]);
od end:

With the same calling conventions as above we get

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
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...

See also

[1] More about factorial formulas can be found on my homepage.

[2] (Added Sep 23 2016) The most comprehensive study of approximations of the factorial (and the Gamma function) to the present day is: W. Wang, Unified approaches to the approximations of the gamma function, J. Number Theory (2016)

Personal tools