OFFSET
2,3
COMMENTS
From Natalia L. Skirrow, Sep 04 2025: (Start)
Occurs in the (formally valid but everywhere-divergent) Laurent series for the error of a variation of Stirling's approximation, from applying the Euler-Maclaurin formula: n! ~ sqrt(tau) * ((n+1/2)/e)^(n+1/2) * exp( Sum_{k>=1} -q(k+1)/(k*n^k) ).
This formula is sometimes attributed to Burnside, but was known to Stirling, who found it before his original approximation.
sqrt(tau) * ((n+1/2)/e)^(n+1/2) = exp(Integral_{x=n-1/2..n+1/2} log(x!) dx), via Raabe's integral; Stirling's original comes from integrating {x=n-1..n} and correcting for the offset by multiplying by sqrt(n).
The log-convexity of Gamma can thus be used to show that Stirling's nominal approximation is a lower bound and his original approximation is an upper bound.
The leading term of -1/(24*n) in this error is -1/2 times that of Stirling's approximation the formulation in the Borwein & Corless link explains this as a consequence of the midpoint rule being asymptotically twice as precise as the trapezoid rule.
This sequence is not explicitly mentioned in Borwein and Corless, but comes out of Theorem 2, eq. (57) via q(n+1) = -n/(-2)^n*Sum_{o=1..k} (2^o-1)*B_{o+1}*C(k-1,o-1)/(o*(o+1)) which can be proven with umbral calculus, using the identity B^n = (1+B)^n (where B^j is B_j).
(End)
LINKS
Natalia L. Skirrow, Table of n, a(n) for n = 2..1024
Jonathan M. Borwein and Robert M. Corless, Gamma and Factorial in the Monthly, Amer. Math. Monthly, 125(5) (2018), 400-424; arXiv version, arXiv:1703.05349 [math.HO], 2017.
E. Chlebus, A recursive scheme for improving the original rate of convergence to the Euler-Mascheroni constant, Amer. Math. Monthly, 118 (2011), 268-274.
Peter Luschny, Approximations to the Factorial Function
Natalia L. Skirrow, Stirling's approximations
FORMULA
From Natalia L. Skirrow, Sep 04 2025: (Start)
With q(n) = a(n)/A189049(n) (when written in lowest terms),
q(n) = (1/(-2)^n - B(n))/n, where B(n) are the Bernoulli numbers.
E.g.f. for q(n): Ei(-x/2) + log(2/(1-e^-x)) - gamma.
E.g.f. for q(n+1): f(x) = (1/x - 1/(2*sinh(x/2))) / exp(x/2).
The series comes out of a Laplace transform representation, via Watson's lemma;
H_n - log(n+1/2) - gamma = Integral_{x=0..oo} f(x)/exp(n*x) dx.
(End)
EXAMPLE
1/(24n^2) - 1/(24n^3) + 23/(960*n^4) - 1/(160n^5) - 11/(8064*n^6) - 1/(896n^7) + 143/(30720*n^8) + ...
MAPLE
a := n -> numer(Zeta(1-n) + (-2)^(-n)/n): seq(a(n), n = 2..32); # Peter Luschny, Sep 09 2025
MATHEMATICA
s = Sum[1/k, {k, 1, n}] - Log[n*(1 + 1/(2*n))] - EulerGamma; CoefficientList[ Series[s, {n, Infinity, 25}], 1/n][[3 ;; -1]] // Numerator (* Jean-François Alcover, Sep 12 2013 *)
Table[(1/(-2)^n-BernoulliB[n])/n, {n, 2, 27}]//Numerator (* Natalia L. Skirrow, Sep 04 2025 *)
CROSSREFS
KEYWORD
sign,frac
AUTHOR
N. J. A. Sloane, Apr 16 2011
EXTENSIONS
Corrected and extended by Jean-François Alcover, Sep 12 2013
STATUS
approved