# Generating terms of the Stirling series

The Stirling series starts as follows: $$n!=\left(\frac{n}{e}\right)^{n}\sqrt{2\pi n} \left\{1+\frac{1}{12n}+\frac{1}{288n^{2}}-\frac{139}{51840 n^{3}}-\frac{571}{2888380 n^{4}}+O\left(n^{-5}\right)\right\}.$$ Is it possible to use Mathematica to explicitly compute more terms?

• Yes: Series[n!, {n, ∞, k}] computes k terms. – Dr. Wolfgang Hintze Jun 24 '15 at 9:49
• – dionys Jun 24 '15 at 11:07

Using Series as suggested by Wolfgang Hintze:

Series[n!, {n, Infinity, 4}]/(Sqrt[2 π n] n^(n Exp[n])) // FullSimplify

(*Exp[(-1 - Log[1/n]) n + SeriesData[n, DirectedInfinity, {}, -1, 1, 1]^5] *
(SeriesData[n, Infinity, {1, 0, 1/12, 0, 1/288, 0, -139/51840, 0, -571/2488320},
0, 10, 2]/n^(E^n*n))*)


Also, Eric W. Weisstein provides some documentation on Stirling series. His notebooks include code for an explicit computation, as well as this approach for Stirling coefficients:

PermutationCyclesD[r_, n_, k_] := 0 /; n <= k r - 1
PermutationCyclesD[r_, n_, 1] := (n - 1)!
PermutationCyclesD[r_, n_, k_] := PermutationCyclesD[r, n, k] =
(n - 1) PermutationCyclesD[r, n - 1, k] +
(-1)^(r - 1) Pochhammer[-(n - 1), r - 1] * PermutationCyclesD[r, n - r, k - 1]

StirlingCoefficient := 1
StirlingCoefficient[n_Integer?Positive] := Module[{k},
Sum[(-1)^k PermutationCyclesD[3, 2 n + 2 k, k] / 2^(n + k)/(n + k)!, {k, 2 n}]]


Computing the first few terms, as given in your question:

Table[StirlingCoefficient[i], {i, 0, 4}]

(*{1, 1/12, 1/288, -139/51840, -571/2488320}*)


Weisstein's method turns out to be a bit complicated when compared to the following procedure, due to Marsaglia and Marsaglia:

SetAttributes[marsagliaA, Listable];
marsagliaA = marsagliaA = 1;
marsagliaA[k_Integer?Positive] := marsagliaA[k] =
(marsagliaA[k - 1] - Sum[j marsagliaA[j] marsagliaA[k - j + 1],
{j, 2, k - 1}])/(k + 1)

m = 5; (* the asymptotic part *)
Sum[(2 k + 1)!! marsagliaA[2 k + 1] n^-k, {k, 0, m - 1}] + O[n, ∞]^m
1 + 1/(12 n) + 1/(288 n^2) - 139/(51840 n^3) - 571/(2488320 n^4) + O[n, ∞]^5


Multiply by the Sqrt[2 π n] (n/E)^n factor and you have the Stirling series.

(Temme gives a similar but less efficient recurrence formula.)