The Number of Inversions and the Major Index of Permutations are Asymptotically Joint-Independently-Normal

Andrew BAXTER^{1} and Doron ZEILBERGER^{1}

Abstract: We use recurrences (alias difference equations) to prove the longstanding conjecture that the two most important permutation statistics, namely the number of inversions and the major index, are asymptotically joint-independently-normal. We even derive more-precise-than-needed asymptotic formulas for the (normalized) mixed moments.

Human Statistics

Human statistics are numerical attributes defined on humans, for example, longevity, height, weight, IQ, and it is well-known, at least empirically, that these are, each separately, asymptotically normal, which means that if you draw a histogram with the statistical data, it would look like a bell-curve. It is also true that they are usually joint-asymptotically-normal, but usually not independently so. But if you compute empirically the correlation matrix, you would get, asymptotically (i.e. for “large” populations) that they are close to being distributed according to a multivariate (generalized) Gaussian exp(−Q(x_{1}, x_{2}, ...)) with Q(x_{1}, x_{2}, …) a certain quadratic form that can be deduced from the correlation matrix.

Permutation Statistics

Let our population be the set of permutations of {1, 2, …, n}. They too, can be assigned numerical attributes, and the great classical combinatorialist, Dominique Foata, (who got his 3rd-cycle doctorate in statistics!) coined the term permutation statistics for them.

The most important permutation statistic is the number of inversions, inv(π), that counts the number of pairs 1 ≤ i < j ≤ n such that π_{i} > π_{j} (and ranges from 0 to n(n − 1)/2). For example, inv(314625) = 5, corresponding to the set of pairs {(1, 2), (1, 5), (3, 5), (4, 5), (4, 6)}. It features in the definition of the determinant, and Netto proved that their probability generating function (the polynomial in q such that its coefficient of q^{i} is the probability that a uniformly-at-random n-permutation has i inversions) is given by

\frac{(1)(1 + q)(1 + q + q^{2}) ⋯ (1 + q + q^{2} + … + q^{n−1})}{n!} = \frac{∏_{i=1}^{n}(1 − q^{i})}{n!(1 − q)^{n}} .

The second most important permutation statistic is the major index, maj(π), that is the sum of the places i, where π_{i} > π_{i+1} For example, maj(314625) = 1 + 4 = 5, because at i = 1 and i = 4 we have descents. Major Percy Alexander MacMahon[M] famously proved that the probability generating function for the major index is also given by that very same formula. In other words the permutation statistics inv and maj are equidistributed. Dominique Foata[Fo] gave a seminal lovely bijective proof that proved the stronger statement that inv and maj are equidistributed also when restricted to permutations ending with a given integer.

1 Department of Mathematics, Rutgers University (New Brunswick), Hill Center-Busch Campus, 110 Frelinghuysen

Rd., Piscataway, NJ 08854-8019, USA. [baxter,zeilberg] at math dot rutgers dot edu , http://www.math.rutgers.edu/~[baxter,zeilberg] . April 5, 2010. Accompanied by the Maple package InvMaj downloadable from http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/invmaj.html .

The work of both authors was supported in part by the USA National Science Foundation.

William Feller ([Fe], p.257) proved that the number of inversions (and hence also the major index) is asymptotically normal in the following sense. Feller easily computed the expectation,

E[inv] = m_{n} = n(n − 1)/4 ,

and the variance,

(σ_{n})^{2} = \frac{2n^{3} + 3n^{2} − 5n}{72} .

If we denote by Z_{n} the centralized and normalized random variable

Z_{n} = \frac{inv − m_{n}}{σ_{n}} ,

then Z_{n} → N in distribution as n → ∞, where N is the Gaussian distribution whose probability density function is e^{−x^{2}/2}/\sqrt{2π}.

A computer-generated proof, which gives much more detail regarding the precise asymptotics, can be

obtained using Zeilberger’s Maple package http://www.math.rutgers.edu/~zeilberg/tokhniot/AsymptoticMoments

that accompanies the article [Z].

Hence both inv and maj are individually asymptotically normal, but what about their interaction, in other words, what can you say about the limit of

\frac{1}{da db}P r(σ_{n}a ≤ inv(π) − m_{n} ≤ σ_{n}(a + da) AND σ_{n}b ≤ maj(π) − m_{n} ≤ σ_{n}(b + db) ) ,

as n → ∞ and da, db → 0?

In this article, we prove that this limit exists and equals (2π)^{−1}e^{−a^{2}/2−b^{2}/2}. In other words, inv and maj are asymptotically joint-independently-normal.

A Brief History

It all started when the great Swedish probabilist Svante Janson (member of the Swedish Academy of Science, that awards the Nobel prizes) asked Donald Knuth (one of the greatest computer scientists of all time, winner of the Turing and Kyoto prizes, among many other honors) about the asymptotic covariance of inv and maj. Neither of these luminaries knew the answer, so Don Knuth asked one of us (DZ). DZ didn’t know the answer either, so he asked his beloved servant, Shalosh B. Ekhad, who immediately ([E]) produced, not just the asymptotics, but the exact answer!. It turned out to be n(n−1)/8. In particular, the correlation, Cov(inv, maj)/σ_{n}^{2} = \frac{\frac{n(n−1)}{8}}{\frac{2n^{3}+3n^{2}−5n}{72}} = \frac{9}{2n} +O(1/n^{2}) tends to zero as n goes to infinity. It followed that for large n, inv and maj are practically uncorrelated.

But there are lots of pairs of random variables that are uncorrelated yet not independent. A convenient way to prove that X_{n} := (inv − m_{n})/σ_{n} and Y_{n} := (maj − m_{n})/σ_{n} are asymptotically independent (we already know that they are both normal) is to use the method of moments, and to prove that the mixed moments

M_{r,s}(n) := E[(X_{n})^{r}(Y_{n})^{s}] ,

tend to the mixed moments of N × N as n → ∞. In other words, for r, s ≥ 1:

lim_{n→∞} M_{2r,2s}(n) = \frac{(2r)!}{2^{r}r!} \frac{(2s)!}{2^{s}s!} ,  (EE)
lim_{n→∞} M_{2r−1,2s}(n) = 0 ,  (OE)
lim_{n→∞} M_{2r,2s−1}(n) = 0 ,  (EO)
lim_{n→∞} M_{2r−1,2s−1}(n) = 0 .  (OO)

Ekhad’s brilliant approach to the Janson-Knuth question merely settled the case r = 1, s = 1 of (OO). Of course, because of symmetry (OE) and (EO) are trivially true (before taking the limits!, i.e. M_{2r,2s−1}(n) ≡ 0 and M_{2r−1,2s}(n) ≡ 0).

One natural approach would be to extend Ekhad’s brilliant answer for M_{1,1}(n) to the general case, and try to derive closed-form expressions for M_{r,s}(n) for larger r and s. Since Ekhad’s proof is so brief, we can cite it here in full.

“Svante Janson asked Don Knuth, who asked me, about the covariance of inv and maj. The answer is \binom{n}{2}/4. To prove it, I asked Shalosh to compute the average of the quantity (inv(π) − E(inv))(maj(π) − E(maj)) over all permutations of a given length n, and it gave me, for n = 1, 2, 3, 4, 5, the values 0, 1/4, 3/4, 3/2, 5/2, respectively. Since we know a priori^{2} that this is a polynomial of degree ≤ 4, this must be it! ”.

Obviously this brute-brute-force approach would be hopeless for deriving polynomial expressions for the moments M_{r,s}(n) for larger r and s. As we will soon see, the degree of the polynomial M_{2r,2s}(n) is 3(r + s), so for example, in order to (rigorously) guess M_{10,10}(n), we would need 31 data points, requiring the computer to examine more than 31! > 0.822 · 10^{34} permutations.

It is perhaps surprising, then, that an inspired, “empirical” (yet fully rigorous), “brute-force” approach does indeed work. The first step is to have a more efficient way to compute the moments M_{r,s}(n), for specific n and specific r and s. We will do this by first designing an efficient way to generate the probability generating function, let’s call it G(n)(p, q), for the pair of statistics (inv, maj). Above we see Netto’s beautiful closed-form expressions for G(n)(1, q) and G(n)(p, 1), but no such closed-form expression seems to exist for the bi-variate generating function, so the best that we can hope for is to find a recurrence scheme.

2 This is the old trick to compute moments of combinatorial ‘statistics’, described nicely [GKP], section 8.2, by

changing the order of summation. It applies equally well to covariance. Rather than actually carrying out the gory details, we observe that this is always a polynomial whose degree is trivial to bound.

A Combinatorial Interlude

Setting aside probability for a few moments, and we focus on a fast algorithm for computing

H(n)(p, q) := ∑_{π∈S_{n}} p^{inv(π)}q^{maj(π)} ,

for n up to, say, 50.

Define the weight of a permutation π to be p^{inv(π)}q^{maj(π)}. Suppose that π ∈ S_{n} ends with i, so we can write π = π′i, where π′ is a permutation of {1, …, i − 1, i + 1, … n}.

When you chop-off i from π you always lose n − i inversions (that is, π′ has n − i fewer inversions than π). The major index, however, decreases by n − 1 if the last letter of π′, let’s call it j, is larger than i, and if j < i the major index does not decrease at all. So writing π = π′′ji, we have

inv(π′′ji) = inv(π′′j) + n − i ,

maj(π′′ji) = \{ \begin{matrix}maj(π′′j), & if j < i ; \\ maj(π′′j) + n − 1, & if j > i. \end{matrix}

Combining, we have

weight(π′′ji) = \{ \begin{matrix}p^{n−i}weight(π′′j), & if j < i ; \\ p^{n−i}q^{n−1}weight(π′′j), & if j > i. \end{matrix}

So in order to compute H(n)(p, q), we need to introduce the more general weight-enumerators of those permutations in S_{n} that end with an i. Let’s call these F (n, i)(p, q). In symbols:

F (n, i)(p, q) := ∑_{{π∈Sn \atop πn=i}} p^{inv(π)}q^{maj(π)} .

It follows that (let’s omit the arguments (p, q) from now on):

F (n, i) = p^{n−i} ∑_{j=1}^{i−1} F (n − 1, j) + p^{n−i}q^{n−1} ∑_{j=i}^{n−1} F (n − 1, j) .  (F ni)

Note that when we chop off the last entry, i, from π = π′′ji, we see π′′j is a permutation of {1, …, i − 1, i + 1, n}. We then “reduce” π′′j to a permutation of {1, …, n − 1} by reducing all the entries larger than i by 1, and so the summation ranges from j = i to j = n − 1 rather than from j = i + 1 to j = n.

Replacing i by i + 1 in the above equation, we have:

F (n, i + 1) = p^{n−i−1} ∑_{j=1}^{i} F (n − 1, j) + p^{n−i−1}q^{n−1} ∑_{j=i+1}^{n−1} F (n − 1, j) .

Subtracting the former equation from p times the latter we get

F (n, i) − pF (n, i + 1) =

p^{n−i} ∑_{j=1}^{i−1} F (n − 1, j) + p^{n−i}q^{n−1} ∑_{j=i}^{n−1} F (n − 1, j)

−p^{n−i} ∑_{j=1}^{i} F (n − 1, j) − p^{n−i}q^{n−1} ∑_{j=i+1}^{n−1} F (n − 1, j)

= p^{n−i} ( ∑_{j=1}^{i−1} F (n − 1, j) − ∑_{j=1}^{i} F (n − 1, j) ) + p^{n−i}q^{n−1} ( ∑_{j=i}^{n−1} F (n − 1, j) − ∑_{j=i+1}^{n−1} F (n − 1, j) )

= −p^{n−i}F (n − 1, i) + p^{n−i}q^{n−1}F (n − 1, i) = p^{n−i}(q^{n−1} − 1)F (n − 1, i) .

Rearranging, we get:

F (n, i) = pF (n, i + 1) + p^{n−i}(q^{n−1} − 1)F (n − 1, i) .  (RecF)

This enables us to compute F (n, i) for i < n starting with F (n, n). We still need to specify F (n, n), and for this we do need the ∑ symbol, namely we use Eq. (F ni) with i = n:

F (n, n) = ∑_{j=1}^{n−1} F (n − 1, j) .  (F nn)

The recurrence (RecF) together with the final condition (F nn), and the trivial initial condition F (1, i) = δ_{i,1}, enables us to efficiently compute F (n, i) for numeric (n, i), for {(n, i)|1 ≤ i ≤ n ≤ N} for any finite N (not too large, but not too small either: e.g., N = 100 is still plausible). In particular, we can compile a table of H(n)(p, q) = F (n + 1, n + 1)(p, q), for n ≤ N − 1.

A crash course in multivariable enumerative probability

Suppose that you have a finite set of objects S and several statistics f_{1}(s), …, f_{r}(s). The multivariable generating function (weight-enumerator under the weight x_{1}^{f_{1}(s)} ⋯ x_{r}^{f_{r}(s)}) is defined to be:

∑_{s∈S} x_{1}^{f_{1}(s)} ⋯ x_{r}^{f_{r}(s)} .

Suppose that you pick an element s ∈ S uniformly at random and you want the multivariable generating function such that the coefficient of x_{1}^{a_{1}} ⋯ x_{r}^{a_{r}} would give you the probability that f_{1}(s) = a_{1}, …, f_{r−1}(s) = a_{r−1}, andf_{r}(s) = a_{r}. It is given by:

P (x_{1}, …, x_{r}) = \frac{1}{|S|} ∑_{s∈S} x_{1}^{f_{1}(s)} ⋯ x_{r}^{f_{r}(s)} .

The expectations, f̄_{1}, …, f̄_{r} are simply

f̄_{i} = (\frac{∂}{∂x_{i}} P ) (1, …, 1) .

The centralized probability generating function, which has mean 0, is

P̃(x_{1}, …, x_{r}) = \frac{P (x_{1}, …, x_{r})}{xf̄_{1} ; 1 ⋯ xf̄_{r} ; r} .

The mixed moments (about the mean) of the statistics f_{1}(s), …, f_{r}(s) are defined by

Mom[a_{1}, …, a_{r}] := \frac{1}{|S|} ∑_{s∈S}(f_{1}(s) − f̄_{1})^{a_{1}} ⋯ (f_{r}(s) − f̄_{r})^{a_{r}} .

Often it is more convenient to consider the mixed factorial moments, using the combinatorial “powers” z^{(r)} := z(z − 1) ⋯ (z − r + 1), better known as the falling-factorials. The mixed factorial moments are defined analogously by

F M[a_{1}, …, a_{r}] = \frac{1}{|S|} ∑_{s∈S}(f_{1}(s) − f̄_{1})^{(a_{1})} ⋯ (f_{r}(s) − f̄_{r})^{(a_{r})} .

Once you know the F M’s for all a_{1}, …, a_{r} ≤ M, you can easily figure out the Mom’s, using Stirling numbers of the second kind (see [Z]). For our present purposes the sets S lie in the family of finite sets S_{n}, and so the leading terms of the F M’s and the Mom’s are the same. Hence for the leading asymptotics it suffices to consider the easier F M’s.

The best way to compute F M[a_{1}, …, a_{r}] for all 0 ≤ a_{1}, …, a_{r} ≤ N for some fixed positive integer N is via the Taylor expansion of P̃(x_{1}, …, x_{r}) around (x_{1}, …, x_{r}) = (1, …, 1), or equivalently, the Maclaurin expansion of P̃(1 + x_{1}, …, 1 + x_{r})

P̃(1 + x_{1}, …, 1 + x_{r}) = ∑_{α_{1},...,α_{r}≥0} \frac{F M(α_{1}, …, α_{r})}{α_{1}! ⋯ α_{r}!}x_{1}^{α_{1}} ⋯ x_{r}^{α_{r}} .

Back to inv-maj

In the present approach, we need to tolerate the more general discrete function F (n, i)(p, q) even though ultimately we are only interested in H(n)(p, q) = F (n + 1, n + 1)(p, q). We will prove the stronger statement that even if you restrict attention to those (n − 1)! permutations that end with a specific i, it is still asymptotically-joint-independently normal.

Since the averages of both inv and maj over the permutations that end in i is n−i+(n−1)(n−2)/4, the centralized probability generating function corresponding to F (n, i)(p, q) is:

G(n, i)(p, q) := \frac{F (n, i)(p, q)}{(n − 1)!(pq)^{n−i+(n−1)(n−2)/4}} .

The recurrence (RecF) becomes

G(n, i) = \frac{1}{q}G(n, i + 1) + \frac{p^{n−i}(q^{n−1} − 1)}{(pq)^{n/2}(n − 1)}G(n − 1, i) ,  (RecG)

and the final condition becomes

G(n, n) = \frac{1}{n − 1} ∑_{j=1}^{n−1}(pq)^{n/2−j}G(n − 1, j) .  (Gnn)

We also need the obvious initial condition G(1, i) = δ_{i,1}.

Guessing Polynomial Expressions for the Factorial Moments

Equipped with these efficient recurrences, our computer computes G(n, i)(p, q) for many values of n and i. Then for each of these it computes the factorial moments F M(r, s)(n, i), for many small numeric (r, s) by computing the (intial terms of the) taylor series centered at (p, q) = (1, 1). Then it fixes numeric values of (r, s) and uses polynomial interpolation to guess explicit polynomial expressions for F M(r, s)(n, i) as polynomials in (n, i), for that pair (r, s). The process is repeated for many pairs (r, s). We note that it is obvious, both from the combinatorics and from the recurrences, that F M(r, s)(n, i) are always polynomials in (n, i) for any fixed numeric r and s.

It would have been nice if we could guess closed-form expressions for F M(r, s)(n, i) for symbolic (r, s), but no such closed-form exists as far as we know, and besides, it is too much to ask for and more than we need. But to prove asymptotic normality we only need the leading terms. Viewing the leading terms, our beloved computer easily conjectures the following expressions, for integers r, s ≥ 0

F M(2r, 2s)(n, i) = \frac{(2r)!}{2^{r}r!} \frac{(2s)!}{2^{s}s!} (\frac{1}{36})^{r+s} n^{3r+3s} + (lower − degree − terms − in − (n, i)) ,

F M(2r, 2s − 1)(n, i) = \frac{(2r)!}{2^{r}r!} \frac{(2s)!}{2^{s}s!} (\frac{1}{36})^{r+s−1} n^{3r+3s−6}[−(s − 1)n^{3} − 6rn^{2}i + 18rni^{2} − 12ri^{3}]

+(lower − degree − terms − in − (n, i)) ,

F M(2r−1, 2s)(n, i) = −\frac{(2r)!}{2^{r}r!} \frac{(2s)!}{2^{s}s!} (\frac{1}{36})^{r+s−1} (r−1)n^{3r+3s−3}+(lower−degree−terms−in−(n, i)) ,

F M(2r−1, 2s−1)(n, i) = \frac{(2r)!}{2^{r}r!} \frac{(2s)!}{2^{s}s!} (\frac{1}{36})^{r+s−1} \frac{9}{2}n^{3r+3s−6}(n−2i)^{2}+(lower−degree−terms−in−(n, i)) .

Nice conjectures but what about proofs?

While we prefer the empirical approach of guessing, an alternative approach to finding many F M(r, s)(n, i)’s is to first use (RecG) and (Gnn). Write G(n, i)(1 + p, 1 + q) as an infinite generic Taylor series around (0, 0), and write-down the implied infinite-order recurrence expressing F M(r, s)(n, i) in terms of F M(r′, s′) with r′ + s′ < r + s. Note that in order to generate F M(r, s) we only need finitely many terms. Of course, as we have already commented, there is no hope for finding a general expression for F M(2r, 2s)(n, i) , F M(2r, 2s − 1)(n, i) , F M(2r − 1, 2s)(n, i) and F M(2r − 1, 2s − 1)(n, i), depending explicitly on r and s as well as on n and i, but to prove by induction on r and s, that the above leading terms are valid, all we need is verify that the leading terms of the implied recurrences for the F M(r, s)’s (easily derivable by hand, although we used the computer) are consistent with the above explicit expressions.

The implication of (RecG) is

F M(r, s)(n, i)−F M(r, s)(n, i+1) = sF M(r, s−1)(n, i+1)−sF M(r, s−1)(n−1, i)+(lower−order−terms) ,

while the implication of Gnn starts with

F M(r, s)(n, n) = \frac{1}{n − 1} ∑_{i=1}^{n−1} F M(r, s)(n − 1, i) − \frac{s}{2(n − 1)} ∑_{i=1}^{n−1}(2i − n)F M(r, s − 1)(n − 1, i)−

\frac{r}{2(n − 1)} ∑_{i=1}^{n−1}(2i − n)F M(r − 1, s)(n − 1, i) + \frac{rs}{4(n − 1)} ∑_{i=1}^{n−1}(2i − n)^{2}F M(r − 1, s − 1)(n − 1, i)

+(lower − order − terms) .

The next step is to refine these recurrences into the four cases where (r, s) are (even, even), (even, odd), (odd, even), and (odd, odd). Once we have these four recurrences, we (or better, our computer) plugs in the above conjectured expressions for the leading terms of F M(r′, s′)(n, i) for r′ + s′ < r + s. When the dust settles, we are left with an expression with leading terms matching those in our conjecture for F M(r, s)(n, i). The proof follows from a routine yet intricate proof by induction on r + s.

The Maple package InvMaj

All the nitty-gritty calculations described above, that constitutes a fully rigorous proof, may be found in the Maple package InvMaj accompanying this article. This package is available from the webpage of the present article:

http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/invmaj.html ,

where the reader can also find some sample input and output. The direct url of the package is:

http://www.math.rutgers.edu/~zeilberg/tokhniot/InvMaj .

La Grande Finale

The special case r = 1, s = 0 and r = 0, s = 1 give

F M(2, 0)(n, i) = \frac{1}{36}n^{3} + O(n^{2})

F M(0, 2)(n, i) = \frac{1}{36}n^{3} + O(n^{2}) .

So

\frac{F M(2r, 2s)}{F M(2, 0)^{r}F M(0, 2)^{s}} = \frac{(2r)!}{2^{r}r!} \frac{(2s)!}{2^{s}s!} + O(1/n) ,

\frac{F M(2r, 2s − 1)}{F M(2, 0)^{r}F M(0, 2)^{s−1/2}} = o(1/n) ,

\frac{F M(2r − 1, 2s)}{F M(2, 0)^{r−1/2}F M(0, 2)^{s}} = o(1/n) ,

\frac{F M(2r − 1, 2s − 1)}{F M(2, 0)^{r−1/2}F M(0, 2)^{s−1/2}} = O(1/n) .

And we see that as n → ∞ these indeed converge to the mixed moments of the famous mixed moments of the bivariate independent normal distribution e^{−a^{2}/2−b^{2}/2}/(2π) .

Encore: A more refined asymptotics for the (Normalized) Mixed Moments

With more effort, we (or rather, our computer) can guess-and-prove the following asymptotics for the case of interest (n + 1, n + 1), i.e. the asymptotic expressions for the (genuine, not factorial) normalized mixed-moments divided by the appropriate powers of the variances (also known as mixed-alpha coefficients), let’s call them α(r, s)(n) . Indeed, according to S. B. Ekhad, we have:

α(2r, 2s)(n) = \frac{(2r)!}{2^{r}r!} \frac{(2s)!}{2^{s}s!} (1 − \frac{9(r^{2} + s^{2} − r − s)}{25} · \frac{1}{n} + O(\frac{1}{n^{2}} ))

α(2r − 1, 2s − 1)(n) = \frac{(2r)!}{2^{r}r!} \frac{(2s)!}{2^{s}s!} (\frac{9}{2n} + (\frac{−81}{50}(r^{2} + s^{2}) + \frac{243}{50}(r + s) − \frac{1773}{100})\frac{1}{n^{2}} + O(\frac{1}{n^{3}} )) .

Of course, by symmetry α(2r, 2s − 1) and α(2r − 1, 2s) are identically (not just asymptotically!) zero.

References

[E] Shalosh B. Ekhad, The joy of brute force: the covariance of the major index and the number of inversions, Personal Journal of S.B.Ekhad and D. Zeilberger, http://www.math.rutgers.edu/ Tilde

zeilberg/pj.html, ca. 1995.

[Fe] William Feller, “An Introduction to Probability Theory and Its Application”, volume 1, three editions. John Wiley and sons. First edition: 1950. Second edition: 1957. Third edition: 1968.

[Fo] Dominique Foata, On the Netto inversion number of a sequence, Proc. Amer. Math. Soc. 19 (1968), 236-240.

[GKP] Ronald Graham, Oren Patashnik, and Donald E. Knuth, Concrete Mathematics: A Foundation for Computer Science, Addison Wesley, Reading, 1989.

[M] Percy A. MacMahon, The indices of permutations, Amer. J. Math. 35 (1913), 281-322.

[Z] Doron Zeilberger, The automatic central limit theorems generator (and much more!), in: “Advances in Combinatorial Mathematics” (in honor of Georgy P. Egorychev), Ilias S. Kotsireas and Eugene V. Zima, eds., Springer, 2009.