1 A Challenging Email

On January 22, 2020 I received the following email by Doron Zeilberger:

figure a

When I received this email, I was thrilled: First, if one gets such an email – in particular from Doron Zeilberger, one is automatically eager to solve it. And second, I highly appreciate “The On-Line Encyclopedia of Integer Sequences” (OEIS) at http://www.oeis.org and supporting it by a donation of Doron Zeilberger gave an extra strong motivation.

Summarizing, the above email provoked various heavy calculations by means of computer algebra that will be described in the following.

2 The Underlying Problem and Symbolic Summation

The combinatorial problem of absent-minded passengers can be introduced as follows.

Definition 2.1

Consider a plane with \(n\ge 2\) seats and suppose that n passengers enter the plane step-wise taking their seats. In addition, suppose that the first \(k\ge 1\) passengers are absent-minded, i.e., they lost their seat ticket and take a seat uniformly at random. The remaining \(n-k\) not absent-minded (but shy) passengers take their dedicated seats (as given in the plane ticket) if it is still free; otherwise, they choose uniformly at random one of the still available free seats. For \(0\le r\le n\), let \(p_{n,k,r}\) be the probability that exactly r passengers sit in the wrong seat and let \(X_n\) be the random variable for “the number of passengers sitting in the wrong seat”. Then the expected value for the passengers sitting in the wrong seat is

$$\begin{aligned} E(X_n)&=\sum _{r=1}^{n}r\,p_{n,k,r} \end{aligned}$$

and its variance is

$$\begin{aligned} V(X_n)&=E(X_n^2)-E(X_n)^2=\sum _{r=1}^{n}r^2\,p_{n,k,r}-E(X_n)^2. \end{aligned}$$

The situation of one absent-minded passenger (\(k=1\)) has been considered in [5, 16] and has been explored further in [7] for the general case \(k\ge 1\). Among other fascinating results, closed forms for \(E(X_n)\) and \(V(X_n)\) have been obtained in [7] by skillful combinatorial arguments. More precisely, the definite sum representations

$$\begin{aligned} E(X_n)&=\frac{k (n-1)}{n} + \sum _{i=1}^{-k +n } \frac{k}{1 -i +n } \end{aligned}$$
(2.1)

and

$$\begin{aligned} V(X_n)=&\frac{k (n-1)}{n^2}+\sum _{i=1}^{-k +n } \frac{(1 -i -k +n ) \big ( 1 -\frac{1 -i -k +n }{1 -i +n } \big )}{1 -i +n}\nonumber \\&+2 \left( \frac{(k-1) k}{2 (n-1) n^2} + \sum _{i=1}^k \sum _{j=1}^{-k +n } \frac{\frac{1 -j -k +n }{-j +n } -\frac{1 -j -k +n }{1 -j +n } }{n} \right) \end{aligned}$$
(2.2)

have been derived and simplified to

$$\begin{aligned} E(X_n)&=\frac{k (n-1)}{n}-k S_1(k)+k S_1(n) \end{aligned}$$
(2.3)

and

$$\begin{aligned} V(X_n)=&\frac{ 2 k -k^2 -2 n -2 k n +2 k^2 n +2 n^2 -k n^2 }{(n-1) n^2}\nonumber \\&-\frac{k (2+n) S_1(k)}{n} +\frac{k (2+n) S_1(n)}{n} +k^2 S_2(k) -k^2 S_2(n) \end{aligned}$$
(2.4)

in terms of the harmonic numbers

$$\begin{aligned} S_o(n)=\sum _{i=1}^n\frac{1}{i^o} \end{aligned}$$

of order \(o\ge 1\) (often they are also denoted by \(H^{(o)}_n=S_o(n)\) with the special case \(H_n=S_1(n)\)). In the following we will prefer to write such expressions in terms of the modified harmonic numbers

$$\begin{aligned} {\bar{S}}_o(n)=S_o(n-1)=\sum _{i=1}^{n-1}\frac{1}{i^o} \end{aligned}$$

yielding, e.g., the more compact expressions

$$\begin{aligned} E(X_n)&=-1 +k -k {\bar{S}}_1(k) +k {\bar{S}}_1(n) \end{aligned}$$
(2.5)

and

$$\begin{aligned} V(X_n)&=\frac{(k-1) k}{(n-1) n} -\frac{k (2+n) {\bar{S}}_1(k)}{n} +\frac{k (2+n) {\bar{S}}_1(n)}{n} +k^2 {\bar{S}}_2(k) -k^2 {\bar{S}}_2(n). \end{aligned}$$
(2.6)

While the simplification from (2.1) to (2.3) is straightforward, more work has to be carried out to derive (2.4) from (2.2). In [7] further details are suppressed how these simplifications have been obtained. Here I want to point out that such (usually painful) classical manipulations are meanwhile obsolete. For instance, by loading in the computer algebra packages

figure b

into Mathematica the multi-sum expression for the variance

figure c

can be simplified to a closed form within seconds by executing the commandFootnote 1

figure d

Remark

Internally, the command EvaluateMultiSums of the package EvaluateMultiSums.m [13] uses the summation paradigms of telescoping, creative telescoping and recurrence solving of the summation package Sigma.m [11]; the underlying algorithms generalize the hypergeometric case [9] to the class of indefinite nested sums and products in the setting of difference fields and rings [8, 14, 15]. In addition, the calculations are supported by special function algorithms of the package HarmonicSums.m [2, 3]. We note that the user is completely freed from applying the summation tools explicitly. However, all the calculation steps are equipped with proof certificates (based on the creative telescoping paradigm introduced in [17]). Thus if necessary, a rigorous correctness proof can be extracted.

3 The Generating Function Approach for Large Moments

As elaborated in detail in [6] the expectation, the variance and higher moments can be elegantly described with the generating function approach. Namely, consider the generating function

$$\begin{aligned} f_n^{(k)}(w)=\sum _{r=0}^{n}p_{n,k,r}\,w^r \end{aligned}$$

of the probability \(p_{n,k,r}\) introduced in Definition 2.1. Then the expectation and variance (see also Definition 2.1) can be straightforwardly connected to the generating function via

$$\begin{aligned} E(X_n)=\frac{d}{dw}f_n^{(k)}(w)|_{w=1} \end{aligned}$$

and

$$\begin{aligned} V(X_n)=\big (w\tfrac{d}{dw}\big )^2f_n^{(k)}(w)|_{w=1}-E(X_n)^2. \end{aligned}$$

More generally, define

$$\begin{aligned} M_0(n,k)=1 \end{aligned}$$

and

$$\begin{aligned} M_l(n,k)=\big (w\tfrac{d}{dw}\big )^lf_{n}^{(k)}(w)|_{w=1} \end{aligned}$$
(3.1)

for \(l\ge 1\); note that \(M_1(n,k)=E(X_n)\). Then one can define the l-th moment by

$$\begin{aligned} m_l(n,k)=\sum _{i=0}^l\left( {\begin{array}{c}l\\ i\end{array}}\right) (-1)^{l-i}M_i(n,k)M_{1}^{l-i}(n,k). \end{aligned}$$
(3.2)

In the above email but also in [6] the task was assigned to compute (besides the known moments \(m_1(n,k)=0\) and \(m_2(n,k)=V(X_n)\)) further moments \(m_l(n,k)\) (at least for \(l=3,4\)).

The combinatorial approach of [7] to derive multi-sum expressions of \(m_l(n,k)\) for larger l seems hopeless. Even though it would have been pure fun to challenge Sigma.m, similarly as in in [10], with more complicated sums than (2.3). Another exciting approach has been carried out in [6, Theorem 5] by guessing the moments \(m_l(n,k)\) for the special case \(k=1\) and \(l=2,\dots ,8\). As it turns out (and proposed by Doron Zeilberger in his email), this attempt can be pushed forward by the following powerful formula given in [6, Theorem 3]:

$$\begin{aligned} \begin{aligned} f_n^{(k)}(w)&=\frac{1}{n!}\sum _{r=0}^kr!\left( {\begin{array}{c}k\\ r\end{array}}\right) w^r(1-w)^{k-r} \prod _{i=0}^{n-k-1}(r\,w+i+1)\\&=\frac{1}{n!}\sum _{r=0}^k(1+r\,w)r!\left( {\begin{array}{c}k\\ r\end{array}}\right) w^r(1-w)^{k-r}(2+r\,w)_{n-k-1}; \end{aligned} \end{aligned}$$
(3.3)

here \((x)_k=x(x+1)\dots (x+k-1)\) denotes the Pochhammer symbol. As a consequence, one can calculate straightforwardly any value of \(M_l(n,k)\) and thus of \(m_l(n,k)\) for particularly chosen n and k. E.g., we can compute

figure e

and activating the formula (3.2) for \(n=10\), \(k=2\) and \(l=3\) yields the third moment

figure f

But even more is possible with the hypergeometric sum representation (3.3). As already proposed in [6] one can derive easily the closed-form expressions of \(M_l(n,k)\) for symbolic n and k. For this task we observe that

$$\begin{aligned} M_l(n,k)=&\big (w\tfrac{d}{dw}\big )^lf_{n}^{(k)}(w)|_{w=1}\nonumber \\ =&\big (w\tfrac{d}{dw}\big )^l\left[ \frac{1}{n!}\sum _{r=0}^k(1+r\,w)r! \left( {\begin{array}{c}k\\ r\end{array}}\right) w^r(1-w)^{k-r}(2+r\,w)_{n-k-1}\right] \Big |_{w=1}\nonumber \\ =&\big (w\tfrac{d}{dw}\big )^l\left[ \frac{1}{n!} \sum _{r=\max (0,k-l)}^k(1+r\,w)r!\left( {\begin{array}{c}k\\ r\end{array}}\right) w^r(1-w)^{k-r}(2+r\,w)_{n-k-1}\right] \Big |_{w=1}, \end{aligned}$$
(3.4)

i.e., for any \(l\in {\mathbb Z}_{\ge 0}\) at most \(l+1\) summands contribute and the remaining summands vanish with the evaluation \(w=1\). Moreover the differentiation of the arising building blocks in (3.4) can be carried out easily. For instance if we differentiate \((2+r\,w)_{n-k-1}\) twice w.r.t. w by using the Mathematica-command D and the Sigma-commands ToSigma and CollectProdSum we get

figure g

In this way, one can compute \(M_1(n,k)(=E(X_n))\) and rediscovers (2.5). Similarly, one gets, e.g.,

$$\begin{aligned} \begin{aligned} M_2(n,k)=&(k-1)^2 +\frac{\big ( -2 k +k n -2 k^2 n \big ) {\bar{S}}_1(k)}{n} +k^2 \big ( {\bar{S}}_1(k)\big )^2\\&+\frac{\big ( 2 k -k n +2 k^2 n \big ) {\bar{S}}_1(n)}{n} -2 k^2 {\bar{S}}_1(k) {\bar{S}}_1(n) +k^2 \big ( {\bar{S}}_1(n)\big )^2 +k^2 {\bar{S}}_2(k) -k^2 {\bar{S}}_2(n) \end{aligned} \end{aligned}$$
(3.5)

and

$$\begin{aligned} \begin{aligned} M_3(n,k)=&2 k^3 {\bar{S}}_3(n)-2 k^3 {\bar{S}}_3(k)-\frac{k \left( 6 k-n^2-5 n\right) {\bar{S}}_1(n)}{(n-1) n}-\frac{3 k {\bar{S}}_1(k)^2}{n}-\frac{3 k {\bar{S}}_1(n)^2}{n}\\&+{\bar{S}}_1(k) \left( \frac{6 k {\bar{S}}_1(n)}{n}+\frac{k \left( 6 k-n^2-5 n\right) }{(n-1) n}\right) +\frac{3 k (k n+2 k-1) {\bar{S}}_2(k)}{n}\\&-\frac{3 k (k n+2 k-1) {\bar{S}}_2(n)}{n}-\frac{(k-2) (k-1) k}{(n-2) (n-1) n}+\frac{3 (k-2) k (n-k)}{(n-1)^2 n}\delta (k-2) \end{aligned} \end{aligned}$$
(3.6)

with

$$\begin{aligned} \delta (x)={\left\{ \begin{array}{ll} 1 &{}\text { if } x\ge 0\\ 0&{}\text { if } x<0.\end{array}\right. } \end{aligned}$$

Note that in \(M_3(n,k)\) the factor \(\frac{3 (k-2) k (n-k)}{(n-1)^2 n}\) does not contribute for \(k=1,2\) (\(k=2\) would vanish anyway); this comes from the fact that in the summation of (3.4) the lower bound \(k-l\) with \(l=3\) should be non-negative.

In a nutshell, one can calculate sufficiently many \(M_l(n,k)\) and using the formula (3.2) one gets the desired moments \(m_l(n,k)\). With \(M_1(n,k)=E(X_n)\) given in (2.5), (3.5) and (3.6) one can reproduce \(m_2(n,k)\) as given in (2.6) and can compute in addition

$$\begin{aligned} m_3(n,k)=&-\frac{(k-2) (k-1) k}{(n-2) (n-1) n}+2 k^3 {\bar{S}}_3(n) +\frac{\big ( 6 k^2 -5 k n -k n^2 \big ) {\bar{S}}_1(k)}{(n-1) n} -\frac{3 k \big ({\bar{S}}_1(k)\big )^2}{n}\\&+\frac{\big ( -6 k^2 +5 k n +k n^2 \big ) {\bar{S}}_1(n)}{(n-1) n} +\frac{6 k {\bar{S}}_1(k) {\bar{S}}_1(n)}{n} -\frac{3 k \big ( {\bar{S}}_1(n)\big )^2}{n}+\frac{3 \big ( -k +2 k^2 +k^2 n \big ) {\bar{S}}_2(k)}{n}\\&-2 k^3 {\bar{S}}_3(k) -\frac{3 \big ( -k +2 k^2 +k^2 n \big ) {\bar{S}}_2(n)}{n} +\frac{3 (k-2) k (n-k)}{(n-1)^2 n}\delta (k-2). \end{aligned}$$

Similarly, one obtains, e.g.,

$$\begin{aligned} m_4(n,k)=&\frac{(k-3) (k-2) (k-1) k}{(n-3) (n-2) (n-1) n} +\frac{ 6 k -12 k^2 +30 k^3 -6 k^4 +18 k n -29 k^2 n -7 k^2 n^2 }{(n-1) n} {\bar{S}}_2(n)\\&+\frac{{\bar{S}}_2(k)}{(n-1) n} \big ( -6 k +12 k^2 -30 k^3 +6 k^4 -18 k n +29 k^2 n +7 k^2 n^2 \big )\\&+\frac{{\bar{S}}_1(k)}{(n-2) (n-1) n} \big ( -52 k +36 k^2 -12 k^3 +40 k n -11 k n^2 -k n^3 \big )\\&+\frac{{\bar{S}}_1(n)}{(n-2) (n-1) n} \big ( 52 k -36 k^2 +12 k^3 -40 k n +11 k n^2 +k n^3 \big )\\&+\frac{3 \big ( -2 k +4 k^2 -6 k n +3 k^2 n +k^2 n^2 \big ) \big ({\bar{S}}_1(k)\big )^2}{(n-1) n} -\frac{4 k \big ( {\bar{S}}_1(k)\big )^3}{n}+3 k^4 \big ( {\bar{S}}_2(n)\big )^2\\&-\frac{6 \big ( -2 k +4 k^2 -6 k n +3 k^2 n +k^2 n^2 \big ) {\bar{S}}_1(k) {\bar{S}}_1(n)}{(n-1) n} +\frac{12 k \big ( {\bar{S}}_1(k)\big )^2 {\bar{S}}_1(n)}{n}\\&+\frac{3 \big ( -2 k +4 k^2 -6 k n +3 k^2 n +k^2 n^2 \big ) \big ({\bar{S}}_1(n)\big )^2}{(n-1) n} -\frac{12 k {\bar{S}}_1(k) \big ( {\bar{S}}_1(n)\big )^2}{n} +\frac{4 k \big ({\bar{S}}_1(n)\big )^3}{n}\\&-\frac{6 \big ( 2 k -4 k^2 +2 k^3 +k^3 n \big ) {\bar{S}}_1(k) {\bar{S}}_2(k)}{n} +\frac{6 \big ( 2 k -4 k^2 +2 k^3 +k^3 n \big ) {\bar{S}}_1(n) {\bar{S}}_2(k)}{n}\\&-\frac{4 \big ( 2 k -6 k^2 +6 k^3 +3 k^3 n \big ) {\bar{S}}_3(k)}{n} +6 k^4 {\bar{S}}_4(k) +\frac{6 \big ( 2 k -4 k^2 +2 k^3 +k^3 n \big ) {\bar{S}}_1(k) {\bar{S}}_2(n)}{n}\\&-\frac{6 \big ( 2 k -4 k^2 +2 k^3 +k^3 n \big ) {\bar{S}}_1(n) {\bar{S}}_2(n)}{n} -6 k^4 {\bar{S}}_2(k) {\bar{S}}_2(n)\\&+\frac{4 \big ( 2 k -6 k^2 +6 k^3 +3 k^3 n \big ) {\bar{S}}_3(n)}{n} -6 k^4 {\bar{S}}_4(n)-\frac{2}{(n-2)^2 (n-1)^2 n} \big ( -4 k +36 k^2\\&-30 k^3 +6 k^4 -24 k n +k^2 n +15 k^3 n -4 k^4 n+16 k n^2 -15 k^2 n^2 +3 k^3 n^2 \big )\delta (k-3)\\&+\Big [ -\frac{6 {\bar{S}}_1(k)}{(n-1)^2 n} \big ( -2 k -5 k^2 +3 k^3 +10 k n -7 k^2 n +k^3 n \big )\\&+\frac{6 {\bar{S}}_1(n)}{(n-1)^2 n} \big ( -2 k -5 k^2 +3 k^3 +10 k n -7 k^2 n +k^3 n \big )+\frac{1}{(n-1)^3 n} \big ( -k +37 k^2\\&-42 k^3 +12 k^4 -34 k n +52 k^2 n -18 k^3 n -13 k n^2 +7 k^2 n^2 \big ) \Big ]\delta (k-2). \end{aligned}$$

For instance, specializing \(m_4(n,k)\) to \(k=1\) gives

$$\begin{aligned} m_4(n,1)=&\big ( \frac{14+n}{n} -6 {\bar{S}}_2({n}) \big ) {\bar{S}}_1({n}) +\frac{3 (n-2) {\bar{S}}_1({n})^2}{n} +\frac{4 {\bar{S}}_1({n})^3}{n}\\&-\frac{(18+7 n) {\bar{S}}_2({n})}{n}+3 {\bar{S}}_2({n})^2 +\frac{4 (2+3 n) {\bar{S}}_3({n})}{n} -6 {\bar{S}}_4({n}) \end{aligned}$$

for one absent-minded passenger (as derived in [6, Theorem 5]) and

$$\begin{aligned} m_4(n,2)=&\big ( \frac{2 \big ( -74+13 n+13 n^2\big )}{(n-1) n} -\frac{24 (1+2 n) {\bar{S}}_2({n})}{n} \big ) {\bar{S}}_1({n})+48 {\bar{S}}_2({n})^2\\&+\frac{12 \big ( 5-2 n+n^2\big ) {\bar{S}}_1({n})^2}{(n-1) n} +\frac{8 {\bar{S}}_1({n})^3}{n} -\frac{4 (21+19 n) {\bar{S}}_2({n})}{n}\\&+\frac{16 (7+6 n) {\bar{S}}_3({n})}{n} -96 {\bar{S}}_4({n})+\frac{2 \big ( 51-45 n+19 n^2\big )}{(n-1) n} \end{aligned}$$

for two absent-minded passengers.

Following this strategy we succeeded in computing the moments \(m_l(n,k)\) up to \(l=15\) on a machine with 12 cores and 1.5TB memory, but failed to proceed due to time and space limitations. Luckily, another great trick from [18] enabled us to compute even more moments. Namely, the calculation of

$$\begin{aligned} \begin{aligned} {\bar{M}}_l(n,k)=&\big (\tfrac{d}{dw}\big )^lf_{n}^{(k)}(w)|_{w=1}\\ =&\big (\tfrac{d}{dw}\big )^l\left[ \frac{1}{n!}\sum _{r=\max (0,k-l)}^k (1+r\,w)r!\left( {\begin{array}{c}k\\ r\end{array}}\right) w^r(1-w)^{k-r}(2+r\,w)_{n-k-1}\right] \Big |_{w=1} \end{aligned} \end{aligned}$$
(3.7)

turns out to be much faster (and less memory consuming) than the calculation of (3.4). Given these so-called exponential moments, one gets back the ordinary moments \(M_l(n,k)\) by using the following formula from [18]:

$$\begin{aligned} M_l(n,k)=\sum _{r=1}^lS(l,r){\bar{M}}_l(n,k) \end{aligned}$$
(3.8)

where S(l, r) denotes the Stirling numbers of the second kind. In summary, we carried out the following steps:

  1. Step 1:

    Calculate the exponential moments \({\bar{M}}_l(n,k)\) with the formula (3.7).

  2. Step 2:

    Calculate the ordinary moments \(M_l(n,k)\) with the formula (3.8).

  3. Step 3:

    Finally, calculate the desired moments \(m_l(n,k)\) with the formula (3.2).

Using this approach, we succeeded in calculating generouslyFootnote 2 the moments \(m_l(n,k)\) up to \(l=21\) in a decent amount of time and memory. The moments up to order 16 are available online at

https://www.risc.jku.at/people/cschneid/data/AMPassenger.tar.gz.

Furthermore, a Mathematica notebook is provided at

https://www.risc.jku.at/people/cschneid/data/AMPassenger.nb

that illustrates the main calculations (e.g., up to \(l=7\)). For higher moments the calculations have been distributed (parallelized) neatly using about 10 Mathematica subkernels. Since the code is rather technical and not easily digestible, it is suppressed within the attached Mathematica notebook. The corresponding timings for the various steps, the total time and the size of the moments are summarized in the Table 1).

Table 1 Timings and size of the calculated moments \(m_l(n,k)\)
Full size table

We note that for \(l\in {\mathbb Z}_{\ge 0}\) a closed-form expression for \(m_l(n,k)\) can be given in terms of the harmonic numbers \(S_r(n)\) and \(S_r(k)\) (resp. \({\bar{S}}_r(n)\) and \({\bar{S}}_r(k)\)) with \(1\le r\le l\). We note further that the sequences generated by the harmonic numbers \(\{S_r(n)\mid r\ge 1\}\) (resp. \(\{{\bar{S}}_r(n)\mid r\ge 1\}\)) are algebraically independent among the rational sequences, i.e., the representation of \(m_l(n,k)\) in terms of the harmonic numbers is optimal. Interestingly enough, the algebraic independence can be shown with the help of the summation paradigm of parameterized (creative) telescoping; see [12, Example 6.3]. For more general classes of harmonic sums we refer also to [4].

4 Asymptotic Expansions

Two days after Doron Zeilberger’s first email, I got a slight update of his challenge:

figure h
$$\begin{aligned} \mathtt \lim _{n\rightarrow \infty } \frac{m_l(n)}{m_2(n)^{l/2}} \end{aligned}$$
figure i

Given the calculations of \(m_l(n,k)\) we proceed as follows. For \(l=2\) we enter the computed closed-form expression

figure j

and use the HarmonicSums command to calculate the first terms of the expansion in n:

figure k

where \(\gamma \approx 0.577216\) is Euler’s constant, \(S[2,\infty ]=S_2(\infty )=\zeta (2)=\frac{\pi ^2}{6}\), and \(\zeta (z)=\sum _{i=1}^{\infty }\frac{1}{i^z}\) denotes the Riemann-Zeta function. Summarizing, we have calculated

$$\begin{aligned} \begin{aligned} m_2(n,k)=&\frac{7 (k-1) k}{6 n^3} +\frac{k (-25+18 k)}{12 n^2} +\frac{k (-1+2 k)}{2 n} -\big ( k +\frac{2 k}{n} \big ) S_1({k})\\&+k^2 S_2({k})-k^2 \zeta (2) +\big ( k +\frac{2 k}{n} \big ) (\gamma +\log (n))+O(\tfrac{1}{n^4}). \end{aligned} \end{aligned}$$
(4.1)

Analogously, we obtain, e.g., the expansion of \(m_3(n,k)\):

$$\begin{aligned} m_3(n,k)&=-\frac{k \big ( 31-38 k+8 k^2\big )}{4 n^3} -\frac{k \big ( 73-90 k+12 k^2\big )}{12 n^2}+\frac{\frac{1}{2} k (-1+6 k) -3 k (-1+2 k) \zeta (2) }{n}\\&+(\gamma +\log (n)) \big ( k -\frac{k (-13+12 k)}{2 n^3} -\frac{3 k (-3+2 k)}{n^2} +\frac{6 k}{n} +\frac{6 k S_1({k})}{n} \big )\\&+\big ( -k +\frac{k (-13+12 k)}{2 n^3} +\frac{3 k (-3+2 k)}{n^2} -\frac{6 k}{n} \big ) S_1({k})\\&-\frac{3 k S_1({k})^2}{n} +\big ( 3 k^2 +\frac{3 k (-1+2 k)}{n} \big ) S_2({k}) -2 k^3 S_3({k}) -3 k^2 \zeta (2) +2 k^3 \zeta (3)\\&-\frac{3 k (\gamma +\log (n))^2}{n} +\big ( -\frac{3 (k-2)^2 k}{n^3} +\frac{3 (k-2) k}{n^2} \big ) \delta (k-2)+O(\tfrac{1}{n^4}). \end{aligned}$$

In total we calculated these expansions of \(m_l(n,k)\) up to \(l\le 18\) with the following timings

l

5

6

7

8

9

10

11

12

13

14

15

16

17

18

time

1.2s

3.6s

9.6s

25

65s

167s

405s

986s

2419s

1.6h

3.8h

9.8h

22.7h

53.4h.

Here we parallelized the calculations on 10 Mathematica kernels. E.g., we needed 53.4 hours to obtain the expansion for \(m_{18}(n,k)\); however summing up all the timings of the used kernels we needed in total 16 days of CPU time to tackle the case \(l=18\). The expansions up to order 16 are available online and the corresponding link can be found on page 8.

To fully win Doron Zeilberger’s challenge, we also dealt with the limit

$$\begin{aligned} c_l(k)=\lim _{n\rightarrow \infty } \frac{m_l(n,k)}{m_2(n,k)^{l/2}},\quad \quad l\ge 2. \end{aligned}$$

Looking at (4.1) the leading term in the expansion of \(m_2(n,k)\) equals \(\log (n)\,k\). In addition the leading term in \(m_3(n,k)\) is also \(\log (n)\,k\). Thus

$$\begin{aligned} c_3(k)=\lim _{n\rightarrow \infty } \frac{m_3(n,k)}{m_3(n,k)^{3/2}} =\lim _{n\rightarrow \infty }\frac{1}{\sqrt{k\log (n)}}=0. \end{aligned}$$

Similarly, given the other computed expansions we get

$$\begin{aligned} c_{2l+1}=0,\quad 0\le l\le 8 \end{aligned}$$

and

l

2

4

6

8

10

12

14

16

18

\(c_l(k)\)

1

3

15

105

945

10395

135135

2027025

34459425.

As a consequence we can confirm that the first values \(c_{2l}(k)\) agree with the double factorial of odd numbers (sequence A001147 in OEIS), i.e.,

$$\begin{aligned} c_{2l}(k)=(2l-1)!!=\prod _{i=1}^{l}(2 i - 1)=\frac{(2 l - 1)!}{(l - 1)!2^{l-1}} \end{aligned}$$

holds for all \(1\le l\le 9\).

5 Conclusion

When I received Doron Zeilberger’s email, the first calculations were straightforward – except that I was not (and I am still not) expert in probability theory and thus made some annoying errors in the beginning. However, pushing the calculations further to the \(21^{st}\) moment and computing its asymptotic expansions was a challenge. Here the computer algebra packages developed at RISC (Sigma.m and EvaluateMultiSums.m by myself and HarmonicSums.m by Jakob Ablinger) were of great help to calculate these huge expressions. I hope that this little note advertises how existing computer algebra tools in general and in particular those of the combinatorics group at RISC can push forward interesting (combinatorial) research topics. For instance, these packages have been heavily used in the last years to carry out large-scale QCD-calculations in the research field of elementary particle physics; see, e.g., [1] and references therein. Last but not least, I am very grateful to Doron Zeilberger who challenged me with these problems and initiated this fun project.