In How Many Ways Can You Reassemble Several Russian Dolls?

Doron ZEILBERGER^{1}

Dedicated to my favorite identical twins: Thotsaporn and Thotsaphon THANATIPANONDA

Preface

In the Fall of 1981, I had the pleasure and honor of living next-door to the eminent combinatorialist Joel Spencer, who was then a visiting professor at the Weizmann Institute of Science in Israel. Joel was always glad to talk to me, provided that I spoke in Hebrew. Joel is also a great punner, and was very proud when he made his first pun in Hebrew, during a volleyball game with other math faculty and students, when he said “kadur sheli” that could mean “my ball” but also “Lee [Segal]’s ball”. One day I got as a present a set of “Russian Dolls” (a.k.a. as Matryoshkas or Babushkas), which is a nested set of dolls. Trying to test Joel’s knowledge of enumeration (after all he is an expert in “Hungarian”, rather than enumerative, combinatorics), I asked him:

In how many ways can one reassemble an n-nested Russian Doll?,

and he immediately replied: This is a Stirling question, it rings a Bell.

Several Russian Dolls

Almost thirty years later, my brilliant student, Thotsaporn “Aek” Thanatipanonda asked me what happens if you have several, say r, identical Russian Dolls. Thanks to Sloane, MathSciNet, and Google scholar, we quickly found out that the case r = 2 goes back to Comtet[C] and is featured in Sloane’s A020554. See also [Ba] and [R], and for an insightful species treatment see [L] and [P]. The case r = 3 was nicely handled by Ed Bender[Be], while the general case was given its coup de grâce by John Devitt and David Jackson[DJ], who used a very ingenious generatingfunctionology approach.

As Herb Wilf[W1] famously said, an “answer” to an enumeration question is an efficient algorithm to generate many terms in the enumerating sequence. Explicit formulas are just one such way, and often not the most efficient one! In this article, I will describe a “calculus” approach, using differential operators, that has the following advantages.

1. It is somewhat faster than using the Devitt-Jackson[DJ] complicated exponential generating function.

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

Rd., Piscataway, NJ 08854-8019, USA. zeilberg at math dot rutgers dot edu , http://www.math.rutgers.edu/~zeilberg . Sept. 16, 2009. Revised: Dec. 16, 2009. Accompanied by Maple package BABUSHKAS downloadable from http://www.math.rutgers.edu/~zeilberg/tokhniot/BABUSHKAS . Exclusively published in the Personal Journal of Ekhad and Zeilberger http://www.math.rutgers.edu/~zeilberg/pj.html and http://arxiv.org . Supported in part by the USA National Science Foundation.

2. It is so flexible that it can handle non-identical Russian Dolls. Given a_{1} single-births, a_{2} pairs of (identical) twins, a_{3} sets of (identical) triplets, a_{4} sets of (identical) quadruplets, …, a_{k} sets of (identical) k-tuplets, …, (and you can’t tell identical twins etc. apart form each other), in how many ways can you partition them into (not necessarily distinct) sets? This is useful if the school principal has to assign them into (non-ordered) classes such that no identical-looking children would be in the same class (or else their teacher won’t be able to tell them apart, and the children can play tricks on her).

3. It can be used in conjunction with Wilf’s [Wi2] celebrated methodology for random selection of combinatorial objects to design quick algorithms for selecting, uniformly at random, a mutli-set set-partition of any given multiset.

4. It beautifully illustrates MacMahon’s lovely method of “differential operators” that transcribes combinatorial operations into differential operators. MacMahon was very fond of it, and he nicely described it in the entry Combinatorial Analysis of the eleventh edition of Encyclopeadia Britannica [M].

5. It beautifully illustrates my favorite methodology of rigorous experimental mathematics. You teach the computer how to do the combinatorics, it derives, all by itself, the (symbolic) differential operator (that for r > 3 would be too complicated to derive by hand, let-alone apply, even for MacMahon), and then the computer goes on and uses it to crank-out as many terms as desired in the enumerating sequence.

6. It beautifully illustrates the notion of catalytic variables. These are variables corresponding to quantities that we may not care about, but are nevertheless needed in order to facilitate the enumeration. At the end of the day, we set them all equal to 1.

7. Last but not least, it enables me to contribute six new sequences to Sloane, the (beginnings of) the enumerating sequence for the number of ways of reassembling r identical n-nested Russian Dolls for 3 ≤ r ≤ 8. So far, only r = 1 (the Bell numbers, A000110) and r = 2 (Comtet’s sequence A020554) are present there.

The Evolution Differential operator

A Baby Example: r = 1

For the sake of pedagogy, let’s first treat the classical case of one Russian Doll.

Suppose that you have a set partition of {1, 2, …, n − 1} with k sets. How can we accommodate a new-comer n? Either she is shy (or anti-social), and decides to form her own new set {n}, or she is outgoing and would like to join an existing set, for which she has k choices. If we call the “number of sets” the “state”, then in the former case the new comer, n, caused the set-partition to move to state k + 1, while in each of the k latter cases, it stayed in state k. If we give state k the weight z^{k}, then each and every set-partition of state k (with weight z^{k}) gives rise to the “evolution”

z^{k} → z^{k+1} + kz^{k} .

In other words

z^{k} → (z + z \frac{d}{dz})z^{k} .

This is true for each and every monomial z^{k}, and by linearity for each polynomial. So if P_{n}(z) is the sum of the weights of all set-partitions of {1, …, n}, then we have the differential-recurrence equation:

P_{n}(z) = D_{1}P_{n−1}(z) ,

where D_{1} is the differential operator

D_{1}f(z) := (z + z \frac{d}{dz})f(z) .

The initial condition is P_{0}(z) = 1. If we are only interested in the total number of set partitions of {1, 2, …, n}, then at the end of the day we plug-in z = 1, getting

B_{n} = P_{n}(1) .

This gives a quick way to crank-out a table of the first one thousand (or whatever) Bell numbers, that is memory efficient. Once we are at day n, and know P_{n}(z), we can let the computer forget about P_{n−1}(z) (and even about B_{n−1}, once it is printed out). Surprisingly, this turned out to be much more efficient (for large n and using Maple) than using

B_{n} = ∑_{i=0}^{n−1} \binom{n − 1}{i}B_{i} ,

or

∑_{n=0}^{∞} \frac{B_{n}}{n!} z^{n} = e^{e^{z}−1} .

A Toddler Example: r = 2

Suppose that we have n−1 pairs of identical twins, already arranged into sets. The only restriction is that no two identical twins can be in the same set, but sets can be repeated. So we have a multiset of sets, whose union, as a multiset, is 1^{2}2^{2} … (n − 1)^{2}.

Of course, no set can be repeated more than twice (why?). Let there be a sets, A_{1}, …, A_{a} that show up once, and let there be b sets B_{1}, …, B_{b} that show up twice. We say that the state of this arrangement is (a, b), and its weight is z_{1}^{a}z_{2}^{b}.

The school principal has to place the newly-arrived pair of twins n and n, but he can’t put them in the same set. There are seven cases.

Case 1: Create two new singleton sets {n}^{2}. There is only one way of doing it, and the new state is (a, b + 1), with weight z_{1}^{a}z_{2}^{b+1}. The corresponding operation on monomials is

f(z_{1}, z_{2}) → z_{2}f(z_{1}, z_{2}) .

Case 2: Create one new singleton set {n}, and place the other n into one of the existing A′_{i}s. There are a ways of doing it, and the new state is (a+1, b), with weight z_{1}^{a+1}z_{2}^{b}. The corresponding operation on monomials is

f(z_{1}, z_{2}) → z_{1} (z_{1} \frac{d}{dz_{1}} ) f(z_{1}, z_{2}) .

Case 3: Create one new singleton set {n}, and place the other n into one of the existing B′_{i}s. There are b ways of doing it, and the new state is (a + 3, b − 1), since one of the B’s became two As, with the help of n. The weight of the new state is z_{1}^{a+3}z_{2}^{b−1}. The corresponding operation on monomials is

f(z_{1}, z_{2}) → z_{1}^{3}z_{2}^{−1} (z_{2} \frac{d}{dz_{2}} ) f(z_{1}, z_{2}) .

Case 4: Place the two twins n and n into two (different, of course) A_{i}’s. There are \binom{a}{2} ways of doing it, and the new state remains (a, b), with weight z_{1}^{a}z_{2}^{b}. The corresponding operation on monomials is

f(z_{1}, z_{2}) → \frac{1}{2} (z_{1}^{2} \frac{d^{2}}{d^{2}z_{1}} ) f(z_{1}, z_{2}) .

Case 5: Place one of the new twins (n or n) into one of the A_{i}’s, and the other one into one of the B_{i}’s. There are ab ways of doing it, and the new state is (a + 2, b − 1), since one of the doubletons B’s, let’s call it B_{i}, was lost, and it became the two distinct sets B_{i} and B_{i} ∪ {n}. The new weight is z_{1}^{a+2}z_{2}^{b−1}. The corresponding operation on monomials is

f(z_{1}, z_{2}) → z_{1}^{2}z_{2}^{−1} (z_{1} \frac{d}{z_{1}} ) (z_{2} \frac{d}{dz_{2}} ) f(z_{1}, z_{2}) .

Case 6: Place the twins (n and n) into two different B_{i}’s. There are \binom{b}{2} ways of doing it, and the new state is (a + 4, b − 2), since two of the doubletons B’s, let’s call them B_{i} and B_{j}, were lost, and they became the four distinct new sets B_{i}, B_{i} ∪ {n}, B_{j} and B_{j} ∪ {n}. The new weight is z_{1}^{a+4}z_{2}^{b−2}. The corresponding operation on monomials is

f(z_{1}, z_{2}) → z_{1}^{4}z_{2}^{−2} \frac{1}{2} (z_{2}^{2} \frac{d^{2}}{dz_{2}^{2}} ) f(z_{1}, z_{2}) .

Case 7: Place both of the new twins (n and n) into each of the two copies of the same B_{i}. There are b ways of doing it, and the new state is the same, (a, b), since single sets stay single and doubletons stay doubletons. The new weight is z_{1}^{a}z_{2}^{b}. The corresponding operation on monomials is

f(z_{1}, z_{2}) → (z_{2} \frac{d}{dz_{2}} ) f(z_{1}, z_{2}) .

Combing, we have just proved:

Fact: Let P_{n}^{(2)}(z_{1}, z_{2}) be the sum of the weights of all multiset set-partitions of the multiset {1^{2} … n^{2}}, with the weight being z_{1} to the power the number of sets that show up once times z_{2} to the power the number of sets that show up twice. Let D_{2} be the partial-differential operator (where D_{1} := \frac{d}{dz_{1}} , D_{2} := \frac{d}{dz_{2}} ),

D_{2} := z_{2}D_{2} + 1/2 z_{1}^{4}D_{2}^{2} + z_{1}^{3}D_{1}D_{2} + 1/2 z_{1}^{2}D_{1}^{2} + z_{1}^{3}D_{2} + z_{1}^{2}D_{1} + z_{2} ,

then

P_{n}^{(2)}(z_{1}, z_{2}) = D_{2}P_{n−1}^{(2)}(z_{1}, z_{2}) .

The Comtet numbers are P_{n}^{(2)}(1, 1).

The general case

Suppose that we already have a multiset set-partition of {1^{r} … (n−1)^{r}}, with a_{1} sets that show-up once, a_{2} sets that show-up twice, …, a_{r} sets that show-up r times. The state of this particular multiset set-partion is (a_{1}, a_{2}, …, a_{r}), and its weight is z_{1}^{a_{1}} … z_{r}^{a_{r}}.

We have to place the r identical new comers n^{r}.

We must make the following decisions

1. How many of them would start their own singleton sets, say, c_{0}

2. For the remaining r − c_{0} new members n, for i = 1, 2, …, r, how many of them, let’s call it c_{i}, would be placed in sets that show up i times.

After these decisions we have a vector of non-negative integers [c_{0}, c_{1}, …, c_{r}], such that c_{0} + c_{1} + … + c_{r} = r.

Once we decided that c_{i} of the new n’s would go to sets that show up i times, we have to decide, amongst those c_{i} siblings which sets should be asked to invite them. These sets can be all different, but they could all be different copies of the same set (if c_{i} ≤ i). This naturally leads to an integer partition λ_{i} = 1^{m_{1}}2^{m_{2}}3^{m_{3}} … i^{m_{i}} (written in multiplicity notation).

We have to place m_{1} of these siblings such that each of them goes to different sets. These make m_{1} of the formerly a_{i}-repeated sets become (a_{i} − 1)-repeated sets and creates m_{1} new singleton sets. So a_{1} → a_{1} + m_{1}, and a_{i−1} → a_{i−1} + m_{1} and a_{i} → a_{i} − m_{1}. In terms of differential operators it is

f(z_{1}, …, z_{r}) → (1/m_{1}!)(z_{1}z_{i−1}D_{i})^{m_{1}}f(z_{1}, …, z_{r}) .

Similarly for 2^{m_{2}} we have

f(z_{1}, …, z_{r}) → (1/m_{2}!)(z_{2}z_{i−2}D_{i})^{m_{2}}f(z_{1}, …, z_{r}) ,

and so on.

So every possible scenario of placing the new identical r siblings of the n family corresponds to an r+1 tuple:

T = [c_{0}, λ_{1}, …, λ_{r}] ,

where c_{0} is an integer, λ_{1}, λ_{2}, …, λ_{r} are integer partitions such that the largest part of λ_{i} is ≤ i, and

c_{0} + |λ_{1}| + |λ_{2}| + … + |λ_{r}| = r .

For each such scenario corresponds the “monomial” operator

P[T] := z_{c_{0}} ∏_{i=1}^{r} Q[λ_{i}] ,

where, writing λ_{i} = 1^{m_{1}}2^{m_{2}}3^{m_{3}} … i^{m_{i}} (m_{1}, m_{2} are now local variables, i.e. they are different, of course, for each λ_{i}),

Q[λ_{i}] = ∏_{j=1}^{i} \frac{1}{m_{j}!}(z_{j}z_{i−j}D_{i})^{m_{j}} .

Here D_{i} := \frac{d}{dz_{i}} and z_{0} := 1.

Finally, we can write down the evolution operator D_{r}:

D_{r} := ∑_{T scenario} P[T] ,

and we have the

Theorem: Let P_{n}^{(r)}(z_{1}, z_{2}, …, z_{r}) be the sum of the weights of all multiset set-partions of the multiset {1^{r} … n^{r}}, with the weight being z_{1} to the power the number of sets that show up once times z_{2} to the power the number of sets that show up twice times … times z_{r} to the power the number of sets that show up r times, then

P_{n}^{(r)}(z_{1}, z_{2}, …, z_{r}) = D_{r}P_{n−1}^{(r)}(z_{1}, z_{2}, …, z_{r}) .

The number of such multi-set set partitions is of course P_{n}^{(r)}(1, 1, …, 1).

Non-Identical Russian Dolls

The operators D_{r} can be combined to yield the

Main Theorem: The number of ways of partitioning into sets a multiset consisting of m_{1} elements that appear once, m_{2} elements that appear twice, …, m_{r} elements that appear r times, or equivalently, the multiset

1 … m_{1}(m_{1} + 1)^{2} … (m_{1} + m_{2})^{2} … (m_{1} + … + m_{r−1} + 1)^{r} … (m_{1} + … + m_{r−1} + m_{r})^{r}

is computed as follows. First compute the polynomial in z_{1}, …, z_{r}:

P (z_{1}, z_{2}, …, z_{r}) = (∏_{i=1}^{r} D_{i}^{m_{i}}) (1) ,

and then plug-in z_{1} = 1, …, z_{r} = 1.

Random Generation

Going back to the Stirling-Bell case, r = 1, the “differential recurrence” P_{n}(z) = (z + z \frac{d}{dz})P_{n−1}(z) is equivalent to the famous recurrence

S(n, k) = S(n − 1, k − 1) + kS(n − 1, k) .

This can be used, according to Wilf[W2], to generate uniformly at random, a set-partion with k sets as follows.

First pre-compute a table of S(n, k) using the recurrence. Now roll a loaded coin with probability of Heads being S(n−1, k−1)/S(n, k) and probability of Tails being kS(n−1, k)/S(n, k). If it lands Heads, recursively generate a random set-partitions of {1, 2, …, n − 1} with k − 1 sets, and adjoin the singleton {n} to it, otherwise generate recursively a random set-partition of {1, 2, …, n − 1} with k sets, and then roll a fair k-sided die, and accordingly decide which of the k members of the set-partition should invite n to join it.

If we want a (uniformly) random set partition, then decide on the number of sets k, by rolling a loaded n-faced die with probabilities of it landing k equalling S(n, k)/B_{n}, and then proceed as before.

The differential-recurrence of the Theorem yields to a partial recurrence for the quantity, let’s call it S^{(r)}(n; a_{1}, …, a_{r}) for the number of multiset set-partitions of 1^{r} … n^{r} with a_{1} sets that show up once, …, a_{r} sets that show up r times. Using this the computer (all by itself!) can use the Wilf Methodology to create a random-generation algorithm. The programming details are a bit daunting, so we leave it as a challenge to the reader.

The Maple package BABUSHKAS

Everything here (except for the random-generation, for which we only have the simple r = 1 case) is implemented in the Maple package BABUSHKAS available, via a link, from the webpage of this

article: http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/babushkas.html
or directly from: http://www.math.rutgers.edu/~zeilberg/tokhniot/BABUSHKAS . That webpage also con-

tains sample input and output, including the sequences for 1 ≤ r ≤ 8.

The main procedure is SeqBrn(r,n) that uses the present approach to generate the first n terms of the enumerating sequence for the number of ways of reassembling r identical Russian Dolls. SeqBrnDJ(r,n) does the same thing using the Devitt-Jackson approach. We are glad to report that they agree! As yet another check, we have the program SSP, that actually constructs the set of all multi-set set partitions of any given multiset, and emables checking, for small values, with the naive count. Procedure B(L) handles the case of non-identical Russian Dolls, or equivalently, an arbitrary multiset, using the Main Theorem. SeqCrn and SeqCrnDJ handle set-partitions of the multiset 1^{r} … n^{r}, in other words, each set can only show up once. This is simply P_{n}^{(r)}(1, 0, …, 0), in the above notation. Full details are available on-line by typing ezra();. The sequences Even though much more data is available in the above-mentioned webpage, and these sequences will soon be submitted to Sloane, let us cite the first ten terms for r = 1, 2, 3, 4. r = 1 : 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975 (the Bell Numbers). r = 2 : 1, 3, 16, 139, 1750, 29388, 624889, 16255738, 504717929, 18353177160 (the Comtet numbers) r = 3 : 1, 4, 39, 862, 35775, 2406208, 238773109, 32867762616, 6009498859909, 1412846181645855 r = 4 : 1, 5, 81, 4079, 507549, 127126912, 55643064708, 38715666455777, 40095856807088486, 58901884724160709571. References

[Ba] G. Baroti, Calcul des nombres de bicouvrements et de birevêtements d’un ensemble fini, employant la méthode fonctionnelle de Rota, in: Combinatorial theory and its Applications I (ed. P. Erdős et al.), North Holland, Amsterdam, 1970), 93-103.

[Be] E.A. Bender, Partitions of multisets, Discrete Mathematics 9(1974), 301-312.

[C] L. Comtet, Birecouvrements et birevetements d’un ensemble fini, Studia Sci. Math. Hungar. 3(1968), 137-152.

[DJ] J.S. Devitt and D.M. Jackson, The enumeration of covers of a finite set, J. London Math. Soc.(2) 25 (1982), 1-6.

[L] G. Labelle, Counting enriched multigraphs, Discrete Math. 217(2000), 237-248.

[M] P.A. MacMahon, Combinatory Analysis, in: Encyclopeadia Britantica, 11th edition (available on-line) vol. 6, 752-758 (1910).

[P] G. Paquin, Dénombrement de multigraphes enrichis, Mémoire, Math. Dept., Univ. Québec à Montréal, 2004.

[R] J. Reilly, “Bicoverings of a Finite Set” by Generating Function Methods, J. Comb. Theory, Series A 28 (1980), 219-225.

[W1] H.S. Wilf, What is an answer?, Amer. Math. Monthly 89 (1982), 289-292.

[W2] H. S. Wilf, A unified setting for sequencing, ranking, and selection algorithms for combinatorial objects, Advances in Mathematics 24(1977), 281-291.