User:Peter Luschny/FactorialFunction

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

$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

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

The first few coefficients a_n 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+{\cfrac {c_{1}}{n+{\cfrac {c_{2}}{n+{\cfrac {c_{3}}{n+\ddots }}}}}}}}$

The first few coefficients c_n 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 W_n 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 H_n 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 G_n (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,
(h(k-1)/k-add((h(k-j)*h(j))/(j+1),j=1..k-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...

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)