1. A challenging Email
On January 22, 2020 I received the following email by Doron Zeilberger:
Dear Carsten,
I (and Shalosh) just posted a paper
https://arxiv.org/abs/2001.06839with a challenge to you (see the middle of page 4)
Can you (and Sigma) extend theorem 5 of that paper
to the general case with k absentminded passengers?
These expressions should be expressions in n,k,
the harmonic numbers H_{n1}, H_{k1} and their
generalizations (the partial sum , to n1, k1, respectively of
zeta(n) and zeta(k)).
If you and Sigma can do the fourth moment, and derive
the asymptotic in n (with a fixed but arbitrary k), I will
donate $100$ to the OEIS in your honor.
...
Best wishes,
Doron
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 OnLine 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 introduced in the following.
2. The underlying problem and symbolic summation
The combinatorial problem of absentminded passengers can be introduced as follows.
Definition 2.1.
Consider a plane with seats and suppose that passengers enter the plane stepwise taking their seats. In addition, suppose that the first passengers are absentminded, i.e., they lost their seat ticket and take a seat (uniformally) at random. The remaining not absentminded (but shy) passengers take their dedicated seats (as given in the plane ticket) if it is still free; otherwise, they choose (uniformally) at random one of the still available free seats. For , let
be the probability that exactly
passengers sit in the wrong seat and letbe 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
and its variance is  
The situation of one absentminded passenger () has been considered in [17, 5] and has been explored further in [7] for the general case . Among other fascinating results, closed forms for and have been obtained in [7] by skillful combinatorial arguments. More precisely, the definite sum representations
(2.1)  
and  
(2.2)  
have been derived and simplified to
(2.3)  
and  
(2.4)  
in terms of the harmonic numbers
of order (often they are also denoted by with the special case ). In the following we will prefer to write such expressions in terms of the modified harmonic numbers
yielding, e.g., the more compact expressions
(2.5)  
and  
(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 carried out. Here I want to point out that such (usually painful) classical manipulations are meanwhile obsolete. For instance, by loading in the computer algebra packages
In[1]:=
Sigma  A summation package by Carsten Schneider © RISCJKU
In[2]:=
HarmonicSums by Jakob Ablinger © RISCJKU
In[3]:=
EvaluateMultiSums by Carsten Scneider © RISCJKU
into Mathematica the multisum expression for the variance
In[4]:=
can be simplified to a closed form within seconds by executing the command^{1}^{1}1Within the Sigma package (and later the HarmonicSums package) is denoted by S[o,n].
In[5]:=
Out[5]=
Remark. Internally, the command EvaluateMultiSums of the package EvaluateMultiSums.m [14] uses the summation paradigms of telescoping, creative telescoping and recurrence solving of the summation package Sigma.m [12]; 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, 15, 16]. In addition, the calculations are supported by special function algorithms of the package HarmonicSums.m [2, 3]. We note that the user is completely dispensed from carrying out the summation tools explicitly. However, all the calculation steps are accomplished by proof certificates (based on the creative telescoping paradigm introduced in [18]). Thus if necessary, a rigorous correctness proof can be extracted.
3. The generating function approach for large moments
As elaborated in details in [6] the expectation, the variance and higher moments can be elegantly described with the generating function approach. Namely, consider the generating function
of the probability introduced in Definition 2.1. Then the expectation and variance (see also Definition 2.1) can be straightforwardly connected to the generating function via
and
More generally, define
and
(3.1) 
for ; note that . Then one can define the th moment by
(3.2) 
In the above email but also in [6] the task was assigned to compute (besides the known moments and ) further moments (at least for ).
The combinatorial approach of [7] to derive multisum expressions of for larger 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 for the special case and . 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]:
(3.3) 
here denotes the Pochhammer symbol. As a consequence, one can calculate straightforwardly any value of and thus of for particularly chosen and . E.g., we can compute
In[6]:=
In[7]:=
Out[7]=
In[8]:=
Out[8]=
In[9]:=
Out[9]=
In[10]:=
Out[10]=
and activating the formula (3.2) for , and yields the third moment
In[11]:=
Out[11]=
But even more is possible with the hypergeometric sum representation (3.3). As already promoted in [6] one can derive easily the closed form expressions of for symbolic and . For this task we observe that
(3.4) 
i.e., for any at most summands contribute and the remaining summands vanish with the evaluation . Moreover the differentiation of the arising building blocks in (3.4) can be carried out easily. For instance if we differentiate twice w.r.t. by using the Mathematicacommand D and the Sigmacommands ToSigma and CollectProdSum we get
In[12]:=
Out[12]=
In this way, one can compute and rediscovers (2.5). Similarly, one gets, e.g.,
(3.5) 
with
Note that in the special case the factor does not contribute; this comes from the fact that in the summation of (3.4) the lower bound should be not negative.
In a nutshell, one can calculate sufficiently many and using the formula (3.2) one gets the desired moments . E.g., with given in (2.5) and (3.5) one can reproduce as given in (2.6). Similarly, one obtains, e.g.,
and  
For instance, specializing to gives
for one absentminded passenger (as derived in [6, Theorem 5]) and
for two absentminded passengers.
Following this strategy we succeeded in computing the moments up to 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 [19] enabled us to compute even more moments. Namely, the calculation of
(3.6) 
turns out to be much faster (and less memory consuming) than the calculation of (3.4). Given the moments of this socalled exponential moments, one gets back the ordinary moments by using the following formula from [19]:
(3.7) 
where denote the Stirling Numbers of the Second kind. In summary, we carried out the following steps:

Calculate the exponential moments with the formula (3.6).

Calculate the ordinary moments with the formula (3.7).

Finally, calculate the desired moments with the formula (3.2).
Using this approach, we succeeded in calculating generously^{2}^{2}2Originally only was required in the challenge. the moments up to in a recent amount of time and memory. The moments up to order 16 are online available and the link can be found in [11]. The corresponding timings^{3}^{3}3The calculations have been distributed (parallelized) neatly using about 10 Mathematica subkernels. for the different steps, the total time and the size of the moments are summarized in the table of Figure 1.
l  time for step 1  time for step 2  time for step 3  total time  size of 

1  0.16 s  0.016 s  0.96 s  1.14 s  88 B 
2  1.86 s  0.025 s  0.068 s  1.95 s  1536 B 
3  3.06 s  0.067 s  0.182 s  3.307 s  6 KB 
4  4.61 s  0.15 s  0.48 s  5.24 s  25 KB 
5  6.66 s  0.33 s  1.18 s  8.18 s  73 KB 
6  9.55 s  0.75 s  2.81 s  13.12 s  210 KB 
7  13.81 s  1.70 s  6.51 s  22.02 s  531 KB 
8  21.53 s  3.73 s  15.14 s  40.40 s  1.3 MB 
9  35.79 s  8.03 s  34.50 s  78.33 s  2.9 MB 
10  63 s  18 s  75 s  156 s  6.4 MB 
11  115 s  37 s  160 s  313 s  13 MB 
12  222 s  77 s  335 s  634 s  27 MB 
13  454 s  233 s  678 s  1364 s  52 MB 
14  1101 s  579 s  1344 s  3024 s  98 MB 
15  2559 s  1063 s  2611 s  6233 s  180 MB 
16  5249 s  2380 s  4995 s  12625 s  326 MB 
17  11510 s  4238 s  9521 s  25270 s  573 MB 
18  22357 s  8807 s  17669 s  48834 s  993 MB 
19  48597 s  17843 s  32300 s  98740 s  1.7 GB 
20  95457 s  30467 s  59384 s  185309 s  2.8 GB 
21  180954s  57621s  126000s  364576s  4.7 GB 
We note that for the closed form expression can be given in terms of the harmonic numbers and (resp. and ) with . We note further that the sequences generated by the harmonic numbers (resp. ) are algebraically independent among the rational sequences, i.e., the representation of in terms of the harmonic numbers is optimal. Interesting enough, the algebraic independence can be shown with the help of the summation paradigm of parameterized (creative) telescoping; see [13, 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:
Finally, to get the donation in your honor you have to complete the challenge. Find asymptotic expressions for $m_l(n)$ for at least l=7 (it would be nice to also have l=8), and prove (hopefully automatically) that
equals 0 if l is odd and (2*l1)(2*l3)...1 if l is even. This would be a partial elementary reproof of the asymptotic normality proved for arbitrary r using "fancy" probability of the asymptotic normality, in the HenzeLast paper.
Given the calculations of we proceed as follows. For we load in the computed closed form expression
In[13]:=
and use the HarmonicSums command to calculate the first terms of the expansion in :
In[14]:=
Out[14]=
where is Euler’s constant, , and denotes the RiemannZeta function. Summarizing, we have calculated
(4.1) 
Analogously, we obtain, e.g., the expansion of :
In total we calculated these expansions of up to with the following timings
Here we parallelized the calculations on 10 Mathematica kernels. E.g., we needed hours to obtain the expansion for ; however summing up all the timings of the used kernels we needed in total days of CPU time to tackle the case . The expansions up to order 16 are online available and the corresponding link can be found in [11].
To fully win Doron Zeilberger’s challenge, we also dealt with the limit
Looking at (4.1) the leading term in the expansion of equals . In addition the leading term in is also . Thus
Similarly, given the other computed expansions we get
and
As a consequence we can confirm that the first values agree with the double factorial of odd numbers (sequence A001147 in OEIS), i.e.,
holds for all .
5. Conclusion
When I received Doron Zeilberger’s email, the first calculations were straightforward – except that I was (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
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 largescale QCDcalculations 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.References
 [1] J. Ablinger, A. Behring, J. Blümlein, A. De Freitas, A. von Manteuffel, and C. Schneider. Calculating Three Loop Ladder and VTopologies for Massive Operator Matrix Elements by Computer Algebra. Comput. Phys. Comm., 202:33–112, 2016. arXiv:1509.08324 [hepph].
 [2] J. Ablinger, J. Blümlein, and C. Schneider. Harmonic sums and polylogarithms generated by cyclotomic polynomials. J. Math. Phys., 52(10):1–52, 2011. [arXiv:1007.0375 [hepph]].
 [3] J. Ablinger, J. Blümlein, and C. Schneider. Analytic and algorithmic aspects of generalized harmonic sums and polylogarithms. J. Math. Phys., 54(8):1–74, 2013. arXiv:1302.0378 [mathph].
 [4] J. Ablinger and C. Schneider. Algebraic independence of sequences generated by (cyclotomic) harmonic sums. Annals of Combinatorics, 22(2):213–244, 2018. arXiv:1510.03692 [cs.SC].
 [5] B. Bolobás. The Art of Mathematics: Coffee Time in Memphis. Cambridge University Press, 2006.
 [6] S.B. Ekhad and D. Zeilberger. The absentminded passengers problem via computer algebra. The Personal Journal of Shalosh B. Ekhad and Doron Zeilberger, 2020. arXiv:2001.06839 [math.CO].
 [7] N. Henze and G. Last. Absentminded passengers. Amer. Math. Monthly, 126(10):867–875, 2019. arXiv:1809.10192 [math.PR].
 [8] M. Karr. Summation in finite terms. J. ACM, 28:305–350, 1981.
 [9] M. Petkovšek, H. S. Wilf, and D. Zeilberger. . A K Peters, Wellesley, MA, 1996.
 [10] H. Prodinger, C. Schneider, and S. Wagner. Unfair permutations. Europ. J. Comb., 32:1282–1298, 2011.
 [11] C. Schneider. Moments and their expansions up to order 16. March 3, 2020, online available at: https://www.risc.jku.at/people/cschneid/data/AMPassenger.tar.gz.
 [12] C. Schneider. Symbolic summation assists combinatorics. Sém. Lothar. Combin., 56:1–36, 2007. Article B56b.
 [13] C. Schneider. Parameterized telescoping proves algebraic independence of sums. Ann. Comb., 14:533–552, 2010. [arXiv:0808.2596]; for a preliminary version see FPSAC 2007.
 [14] C. Schneider. Simplifying multiple sums in difference fields. In C. Schneider and J. Blümlein, editors, Computer Algebra in Quantum Field Theory: Integration, Summation and Special Functions, Texts and Monographs in Symbolic Computation, pages 325–360. Springer, 2013. arXiv:1304.4134 [cs.SC].
 [15] C. Schneider. A difference ring theory for symbolic summation. J. Symb. Comput., 72:82–127, 2016. arXiv:1408.2776 [cs.SC].
 [16] C. Schneider. Summation Theory II: Characterizations of extensions and algorithmic aspects. J. Symb. Comput., 80(3):616–664, 2017. arXiv:1603.04285 [cs.SC].
 [17] P. Winkler. Mathematical Puzzles: A Connoisseur’s Collection. A K Peter, 2004.
 [18] D. Zeilberger. The method of creative telescoping. J. Symbolic Comput., 11:195–204, 1991.

[19]
D. Zeilberger.
The automatic central limit theorems generator (and much more!).
In Advances in combinatorial mathematics, pages 165–174. Springer, Berlin, 2009.
Comments
There are no comments yet.