# How good is this approximation of $n$ factorial?

For fixed $n\geq 1$, consider the sum $$S_n:=\sum_{k=1}^\infty \frac{k^n}{e^k}.$$ If you compute this sum for some small values of $n$, you will see that it is remarkably close to $n!$. The fact that it is somehat close to $n!$ is not too surprising because $S_n$ "should" be similar to the following integral $$n! = \int_0^\infty \frac{x^n}{e^x}dx.$$ However, the precision of the approximation seems to be surprisingly good (check some small values on Wolfram Alpha and see for yourself!).

Question: Is there a good way to show that $S_n$ is bounded above by $(1+\varepsilon) n!$ for every $n\geq1$ where $\varepsilon$ is some explicit positive constant which is "close" to $0$?

After checking some small values of $n$, it appears that $n=3$ might be the worst value for an upper bound of this type. Note that, for many values of $n$ (perhaps for most $n$), it is actually the case that $S_n<n!$.

If you think of $S_n$ as a Riemann sum approximating the integral, then the $k$th summand is approximating the area under the curve between $x=k-1$ and $x=k$ by the value at the right endpoint of this interval (namely, by $\frac{k^n}{e^k}$). The function $\frac{x^n}{e^x}$ is increasing for $0\leq x<n$ and decreasing for $x>n$. Therefore, for $k\leq n$, the $k$th summand is overestimating the area and for $k\geq n+1$ it is underestimating it. Magically, it seems that the overestimating terms and the underestimating terms cancel out beautifully and give us a result which is really close to $n!$. Does anyone know whether this phenomenon occurs for all $n$ and whether there is a good reason why this would be occur?

By the way, there is a connection with Eulerian numbers here, if anyone is interested. It turns out that $$S_n = \frac{e\cdot \sum_{m=0}^{n-1}A(n,m)e^m}{(e-1)^{n+1}}$$ where $A(n,m)$ is the number of permutations of $1,\dots,n$ with exactly $m$ "ascents." Interestingly, $A(n,0)+A(n,1)+\cdots +A(n,n-1)=n!$ (since every permutation has some number of ascents) but this doesn't seem to be terribly helpful since there are these factors of $e$ and $e-1$ floating around everywhere.

For any $\alpha > 0$ we have that $\sum_{k\geq 1}e^{-\alpha k}=\frac{1}{e^{\alpha}-1}$, hence by differentiating both sides $n$ times with respect to $\alpha$ we get that $$(-1)^n\sum_{k\geq 1} k^n e^{-\alpha k} = \frac{d^n}{d\alpha^n}\,\frac{1}{e^\alpha-1}=\frac{d^n}{d\alpha^n}\left[\frac{1}{\alpha}+\sum_{m\geq 1}\frac{B_m}{m!}\alpha^{m-1}\right] \tag{1}$$ hence: $$(-1)^n \sum_{k\geq 1} k^n e^{-\alpha k} = \frac{n!(-1)^n}{\alpha^{n+1}}+\sum_{m\geq n+1}\frac{B_m}{m!}\cdot\frac{(m-1)!}{(m-n-1)!}\alpha^{m-n-1}\tag{2}$$ and by evaluating at $\alpha=1$ it comes the magic: $$\sum_{k\geq 1}\frac{k^n}{e^k} = \color{red}{n!}+(-1)^n \sum_{m\geq n+1}\frac{B_m}{m(m-n-1)!}\tag{3}$$ The behaviour of Bernoulli numbers (I am going to talk about Bernoulli numbers with even index, since every Bernoulli number with odd index is zero, with the exception of $B_1$) is quite erratic: till $B_{12}$ they are all less than one in absolute value, then their absolute value starts growing pretty fast: $|B_{2 n}| \sim 4 \sqrt{\pi n} \left(\frac{n}{ \pi e} \right)^{2n}$ for large values of $n$. Since the remainder series in $(3)$ converges pretty fast and the first Bernoulli numbers are essentially negligible, for small values of $n$ (namely $n\leq 16$) $S_n$ and $n!$ are very close, as conjectured.

• Thanks for this! I had not come across the Bernoulli numbers, and so the second equality in (1) looks rather foreign to me. I will have to do some reading. – Jon Noel May 31 '17 at 23:27
• A question: Assuming that Bernoulli numbers are all positive, your final formula seems to suggest that $S_n$ should be slightly smaller than $n!$ when $n$ is odd and slightly larger than $n!$ when $n$ is even. However, for $n=3$ and $n=7$ the sum is larger than $n!$ (and for $n=5$ it is smaller). – Jon Noel May 31 '17 at 23:29
• @JonNoel: they are not positive: Bernoulli numbers with odd index are zero and Bernoulli numbers with even index have alternating signs. Everything is pretty well-known and easy to find, starting from Wikipedia. – Jack D'Aurizio May 31 '17 at 23:31
• @JonNoel note the Bernoulli numbers are of an alternating sorts in themselves, so it kinda cancels. – Simply Beautiful Art May 31 '17 at 23:31
• Ah, I have just read that $B_n$ is negative if $n$ is divisible by $4$. So, this partially explains some of the strange behaviour. – Jon Noel May 31 '17 at 23:32

By the Euler-Maclaurin formula,

$$S_n\approx n!-\frac1{2e}$$

You can get better approximations by taking more terms.

• But considering the function is zero at both endpoints, is there another term apart from the remainder? – Chappers May 31 '17 at 23:05
• @Chappers That's what I was wondering... but I wasn't sure. – Simply Beautiful Art May 31 '17 at 23:08
• The essence of the EML formula lies in Bernoulli numbers and polynomials, so why not employ them directly instead of relying on a uncertain approximation (uncertain due to the remainder term, of course)? – Jack D'Aurizio May 31 '17 at 23:10
• @JackD'Aurizio Thinking the remainder is basically your answer. – Simply Beautiful Art May 31 '17 at 23:12
• @SimplyBeautifulArt: that's right :D – Jack D'Aurizio May 31 '17 at 23:12

$\newcommand{\bbx}{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,} \newcommand{\braces}{\left\lbrace\,{#1}\,\right\rbrace} \newcommand{\bracks}{\left\lbrack\,{#1}\,\right\rbrack} \newcommand{\dd}{\mathrm{d}} \newcommand{\ds}{\displaystyle{#1}} \newcommand{\expo}{\,\mathrm{e}^{#1}\,} \newcommand{\ic}{\mathrm{i}} \newcommand{\mc}{\mathcal{#1}} \newcommand{\mrm}{\mathrm{#1}} \newcommand{\pars}{\left(\,{#1}\,\right)} \newcommand{\partiald}[]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\root}[]{\,\sqrt[#1]{\,{#2}\,}\,} \newcommand{\totald}[]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}} \newcommand{\verts}{\left\vert\,{#1}\,\right\vert}$ \begin{align} \left.\sum_{k = 1}^{\infty}{k^{n} \over \expo{k}}\right\vert_{\ n\ \geq\ 0} & = -\delta_{n0} + \sum_{k = 0}^{\infty}{k^{n} \over \expo{k}} = -\delta_{n0} + \int_{0}^{\infty}{k^{n} \over \expo{k}}\,\dd k + \left.{1 \over 2}{k^{n} \over \expo{k}}\right\vert_{\ k\ =\ 0} -2 \int_{0}^{\infty}{\Im\pars{\bracks{\ic x}^{n}/\expo{\ic x}} \over \expo{2\pi x} - 1}\,\dd x \\[5mm] & = -\,{1 \over 2}\,\delta_{n0} + n! - 2\left\{\begin{array}{lcl} \ds{-\pars{-1}^{n/2}\int_{0}^{\infty}{x^{n}\sin\pars{x} \over \expo{2\pi x} - 1}\,\dd x} & \mbox{if} & \ds{n}\ \mbox{is}\ even \\[3mm] \ds{\phantom{-\,}\pars{-1}^{\pars{n - 1}/2}\int_{0}^{\infty}{x^{n}\cos\pars{x} \over \expo{2\pi x} - 1}\,\dd x} & \mbox{if} & \ds{n}\ \mbox{is}\ odd \end{array}\right. \end{align}

$$\left.\vphantom{\large A}n!\,\right\vert_{\ n\ \geq\ 0} = \sum_{k = 1}^{\infty}{k^{n} \over \expo{k}} + {1 \over 2}\,\delta_{n0} + 2\left\{\begin{array}{lcl} \ds{-\pars{-1}^{n/2}\int_{0}^{\infty}{x^{n}\sin\pars{x} \over \expo{2\pi x} - 1}\,\dd x} & \mbox{if} & \ds{n}\ \mbox{is}\ even \\[3mm] \ds{\phantom{-\,}\pars{-1}^{\pars{n - 1}/2}\int_{0}^{\infty}{x^{n}\cos\pars{x} \over \expo{2\pi x} - 1}\,\dd x} & \mbox{if} & \ds{n}\ \mbox{is}\ odd \end{array}\right.$$

Relative Error plot: Absolute Error plot: It's worsening as $\ds{n > 16}$. Such behaviour is fully discussed in @Jack D'Aurizio answer. • What is $\mathfrak{I}$? – Jon Noel Jun 1 '17 at 8:11
• $\Im z$ means "Imaginary part of $z$". – Felix Marin Jun 1 '17 at 21:50
• Your 'remainder' can be calculated explicitly in terms of $\coth$ derivatives...so why not writing it down? – tired Jun 2 '17 at 20:32
• @tired I guess it leads to the $\texttt{@Jack D'Aurizio}$ answer via the relation between $\left(\mathrm{e}^{2\pi x} - 1\right)^{-1}$ and the $\mathsf{Bernoulli\ numbers}$. So, I decided the best idea was to give $\texttt{@Jack D'Aurizio}$ answer as a reference. – Felix Marin Jun 2 '17 at 21:50