Generating terms of the Stirling series

Question

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?

Answer 1

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, 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[0] := 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}*)

Answer 2

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[0] = marsagliaA[1] = 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.)