Counting Parabolic Double Cosets in Symmetric Groups

Billey, Konvalinka, Petersen, Solfstra, and Tenner recently presented a method for counting parabolic double cosets in Coxeter groups, and used it to compute $p_n$, the number of parabolic double cosets in $S_n$, for $n\leq13$. In this paper, we derive a new formula for $p_n$ and an efficient polynomial time algorithm for evaluating this formula. We use these results to compute $p_n$ for $n\leq5000$ and to prove an asymptotic formula for $p_n$ that was conjectured by Billey et al.


Introduction
For a Coxeter system (W, S), each subset I ⊆ S generates a subgroup W I of W . These subgroups of W are called parabolic subgroups. A parabolic double coset in the Coxeter system (W, S) is a double coset of the form W I wW J for w ∈ W and I, J ⊆ S. Parabolic double cosets have several properties that make them interesting objects to study. For example, the parabolic double cosets of a finite Coxeter group (W, S) are rank-symmetric intervals in the Bruhat order on W [3], and the poset of presentations of parabolic double cosets gives a boolean complex that retains many properties of the Coxeter complex [6].
The problem of counting parabolic double cosets is considered in [1]. Counting parabolic double cosets is hard because a given parabolic double coset C might arise from multiple choices of I and J. To avoid this problem, Billey, Konvalinka, Petersen, Slofstra, and Tenner take the approach of counting lex-minimal presentations [1]. The lex-minimal presentation of a parabolic double coset C is the unique presentation C = W I wW J such that w is the minimal element of C in the Bruhat order and such that (|I| , |J|) is lexicographically minimal among all presentations of C. Theorem 4.26 of [1] gives a formula for the number of parabolic double cosets with a given minimal element w ∈ W , for any Coxeter group W . In the case of the symmetric group, W = S n , summing this formula over all n! elements of S n will compute p n , the total number of distinct parabolic double cosets in S n . The sequence {p n } may be found at [5, A260700]. The values of p n for n ≤ 13 have been computed using this method [1]. However, the summation over all n! elements of S n limits the number of terms of {p n } that can be computed with this approach.
In this paper we restrict our attention to counting parabolic double cosets in S n . Two ways of depicting parabolic double cosets in S n are the balls-in-boxes model considered in [1] and two-way contingency tables, which are matrices of nonnegative integers with nonzero row sums and column sums. Diaconis and Gangolli use the balls-in-boxes model to give a bijection between parabolic double cosets in the double quotient W I \S n /W J = {W I wW J : w ∈ S n } and two-way contingency tables with prescribed row sums and column sums [2]. They attribute the idea behind the bijection to Nantel Bergeron. This bijection is also used in [6], where it provides an alternate construction of the two-sided Coxeter complex of S n in terms of two-way contingency tables.
There are three main results of this paper. The first result is the explicit formula where the values f n,k = n l=k l! n−k l−k are generalizations of the Fubini numbers. The second result is an efficient polynomial time algorithm to evaluate this formula, which is used to compute the values of p n for n ≤ 5000. These values of p n are summarized in the appendix. The third result is an asymptotic formula for p n that proves the following conjecture. We conclude this introduction by outlining the two key ideas behind formula (1). The first idea is to associate each parabolic double coset with a canonical two-way contingency table. In particular, there is a bijection between parabolic double cosets in S n and two-way contingency tables with sum n that are "maximal" (as defined at the end of section 2). This reduces the problem of computing p n to the problem of counting maximal two-way contingency tables with sum n. Two-way contingency tables have both a restriction on the rows (nonzero row sums) and a restriction on the columns (nonzero column sums). These two restrictions are entangled since any nonzero entry in the table both makes its row sum nonzero and its column sum nonzero. The second idea is to break this entanglement by transforming the problem from counting two-way contingency tables to counting pairs of weak orders (a weak order is a binary relation that is transitive and has no incomparable pairs of elements).
This transformation is applied in sections 3.2 and 3.3 to the problem of counting maximal two-way contingency tables with sum n, as well as to the simpler problem of counting all two-way contingency tables with sum n. More generally, this transformation is applicable whenever counting two-way contingency tables with a restriction on the locations of nonzero entries, but not on their specific values.

Background
In this section, we give an overview of parabolic double cosets in S n , and describe two ways of depicting presentations of parabolic double cosets in S n . We will use s i to denote the ith adjacent transposition in S n (the transposition that swaps i and i + 1).
Definition 2.1. A parabolic subgroup W I in the symmetric group S n is a subgroup W I = I = s : s ∈ I of S n generated by a collection I ⊆ {s 1 , . . . , s n−1 } of adjacent transpositions in S n . Definition 2.2. A parabolic double coset C in the symmetric group S n is a double coset of the form C = W I wW J for parabolic subgroups W I and W J in S n and an element w ∈ S n . The triple (I, w, J) is called a presentation of C. The total number of distinct parabolic double cosets in S n is denoted by p n . The sequence {p n } may be found at [5, A260700].
A parabolic double coset C will usually have many presentations. For example, if (I, w, J) is a presentation of C, then so is (I, w , J) for any w ∈ C. As the next example shows, there could also be multiple choices for the collections I and J. In other words, C might be an element of more than one double quotient For each element w ∈ S n , the ball-in-boxes picture of w is given by placing n balls in an n × n grid of boxes in the layout of the permutation matrix of w. A ball is placed in column a and row b if and only if w(a) = b. We label the columns left-to-right and the rows bottom-to-top, as in the Cartesian coordinate system. Left-multiplication acts by permuting the rows of the grid and right-multiplication acts by permuting the columns of the grid.
Consider a presentation (I, w, J) of a parabolic double coset C in S n . The group W I acts on the balls-inboxes picture of w by permuting certain rows of the grid. The group W J acts on the balls-in-boxes picture of w by permuting certain columns of the grid. For each adjacent transposition s i / ∈ I, we draw a horizontal wall between row i and row i + 1. For each adjacent transposition s j / ∈ J, we draw a vertical wall between column j and column j + 1. Then W I acts on the balls-in-boxes picture of w by permuting rows of the grid within the horizontal walls, and W J acts on the balls-in-boxes picture of w by permuting columns of the grid within the vertical walls. The resulting picture is called the balls-in-boxes picture with walls for the presentation (I, w, J) [1].
Example 2.4 (Example 3.5 in [1]). Let w = 3512467 ∈ S 7 , let I = {s 1 , s 3 , s 4 }, let J = {s 2 , s 3 , s 5 , s 6 }, and let C = W I wW J . Figure 1 depicts the balls-in-boxes picture with walls for the presentation (I, w, J) of C. The walls naturally divide the n × n grid into larger rectangular cells. Notice that the number of balls in each cell remains unchanged under the actions of W I and W J . Thus the number of balls in a given cell is constant over all elements of C. By counting the number of balls in each cell, we obtain a matrix of nonnegative integers such that the sum of all of the entries is n and such that each row sum and column sum is strictly positive.
Definition 2.5. A two-way contingency table is a matrix of nonnegative integers such that each row sum and column sum is strictly positive. Let p n,0 denote the number of two-way contingency tables with sum n. The sequence {p n,0 } may be found at [5, A120733].
Thus, each presentation (I, w, J) of C gives rise to a two-way contingency table with sum n. Conversely, given a two-way contingency table with sum n, it is possible to recover the collections I and J, as well as the parabolic double coset C ∈ W I \S n /W J [2, Section 3B]. As a consequence, two-way contingency tables with sum n are in bijection with triples (C, I, J) consisting of two collections I, J ⊆ {s 1 , . . . , s n−1 } and a parabolic double coset C ∈ W I \S n /W J [6]. Another way to put this is that p n,0 = I J |W I \S n /W J |, because the right hand side counts triples (C, I, J). Even though a given parabolic double coset might have multiple distinct two-way contingency tables arising from multiple choices for the collections I and J, the following proposition shows that it is always possible to choose I and J to be simultaneously maximal, resulting in a canonical two-way contingency table associated to each parabolic double coset.
Combining this with the earlier inclusion W I wW J ⊆ C gives the equality W I wW J = C. Thus, (I, w, J) is a presentation of C. Finally, (I, w, J) is maximal since we already showed that if (I , w , J ) is any presentation of C, then I ⊆ I and J ⊆ J.
For an alternative proof of Proposition 2.7, see Proposition 3.7 in [1]. Proposition 2.7 gives a bijection  We remark that the maximal two-way contingency table associated to a given parabolic double coset has the smallest possible dimensions and the largest possible entries.
Putting this all together gives a bijection parabolic double cosets C ∈ Sn

←→
Maximal two-way contingency tables with sum n . Corollary 2.9. The sequence {p n,0 } counts the number of two-way contingency tables with sum n. The sequence {p n } counts the number of maximal two-way contingency tables with sum n.
We will now turn our attention to counting maximal two-way contingency tables with sum n.
3 Computation of p n

Characterization of Maximality
The following proposition describes how to determine whether the conditions s i C = C and Cs i = C are satisfied from the balls-in-boxes picture with walls. We will use the term cell to mean one of the rectangular regions formed by the walls of the balls-in-boxes picture with walls.
Proposition 3.1. Let C = W I wW J be a parabolic double coset in S n , and let P be the balls-in-boxes picture with walls for the presentation (I, w, J). Proof. Let s i ∈ {s 1 , . . . , s n−1 } \ I. First suppose that the row of cells of P directly below the s i -wall and the row of cells of P directly above the s i -wall each have only one nonempty cell, the two of which are vertically adjacent. Let A be the nonempty cell of P directly below the s i -wall and let B be the nonempty cell of P directly above the s i -wall. Since A and B are each the only nonempty cells in their row of cells, the ball in row i must lie in A and the ball in row i + 1 must lie in B. Then the ball in column w −1 (i) must lie in A and the ball in column w −1 (i + 1) must lie in B. Since the cells A and B are vertically adjacent, w −1 (i) and w −1 (i + 1) lie in the same column of cells. Then the transposition w −1 s i w that swaps w −1 (i) and w −1 (i + 1) is an element of W J . Taking w −1 s i w ∈ W J and multiplying on the left by w gives s i w ∈ wW J ⊆ C. Since w ∈ C was arbitrary, this shows that For the converse, suppose that P has a nonempty cell A directly below the s i -wall and a nonempty cell B directly above the s i -wall such that A and B are not vertically adjacent. This is just the negation of the original assumption on P . Since A and B are nonempty, we can multiply w on the left by an element x ∈ W I to make A contain the ball in row i and to make B contain the ball in row i + 1. Then s i xw has fewer balls in cells A and B than xw. However, the number of balls in a given cell is constant over all elements of C.
Since xw ∈ C, this forces s i xw ∈ C. In particular, s i C = C. This proves the first statement. The proof of the second statement is similar.
By definition, a two-way contingency table is maximal precisely when neither of the "if and only if"s in Proposition 3.1 are ever satisfied. This gives a characterization of maximal two-way contingency tables.
Corollary 3.2. Let T be a two-way contingency table. Then T is maximal if and only if T satisfies both of the following two conditions: 1. There does not exist a pair of vertically adjacent nonzero entries, each of which is the only nonzero entry in its row.
2. There does not exist a pair of horizontally adjacent nonzero entries, each of which is the only nonzero entry in its column.

Sequence Transformation
Note that the conditions in Corollary 3.2 only depend on the configuration of nonzero entries of T . Because of this, it will be helpful to ignore the specific values of the nonzero entries of T .

Definition 3.4.
A pattern is a two-way contingency table consisting of 0's and 1's.
1. Let P m,0 denote the collection of all patterns with sum m.
2. Let P m denote the collection of all patterns with sum m that satisfy the conditions of Corollary 3.2.
3. For a two-way contingency table T , we define the pattern associated to T to be the two-way contingency table formed by replacing all nonzero entries of T with 1's.
For each pattern P with sum m, there are exactly n−1 n−m two-way contingency tables with sum n whose associated pattern is P . This is because two-way contingency tables with sum n whose associated pattern is P are in bijection with ways to surjectively distribute n identical balls among m distinguishable urns.
Then combining Corollary 2.9 with Corollary 3.2 gives the formulas Here we use the conventions −1 0 = 1 and also k−1 k = 0 for k ≥ 1, since n−1 n−m is counting ways to surjectively distribute n identical balls among m distinguishable urns. We now invoke the identity where · · denotes unsigned Stirling numbers of the first kind and · · denotes Stirling numbers of the second kind (see the exercise on Lah numbers in [8]). If we define Proposition 3.5. The sequence {q n,0 } counts the number of functions from n distinguishable balls to a rectangular matrix of boxes, such that each row and column is nonempty. The sequence {q n } counts the number of such functions that satisfy the following two conditions: 1. There does not exist a pair of vertically adjacent nonempty boxes, each of which is the only nonempty box in its row.
2. There does not exist a pair of horizontally adjacent nonempty boxes, each of which is the only nonempty box in its column.
Proof. Relabeling indices (from q k,0 to q n,0 and from q k to q n ) gives Then the proposition follows from the fact that k! n k counts the number of surjective functions from n distinguishable balls to the k nonempty boxes prescribed by the pattern in P k,0 or P k .

Pairs of Weak Orders
Definition 3.6. A weak order on a set S is a binary relation ≤ on S that is transitive and has no incomparable pairs of elements (so for all s, t ∈ S, s ≤ t or t ≤ s). The total number of weak orders on {1, . . . , n} is denoted by f n . The sequence {f n } is known as the Fubini numbers or the ordered Bell numbers, and may be found at [5, A000670].
Given a weak order on a set S, it is possible for two elements s, t ∈ S to be tied with each other, in the sense that both s ≤ t and t ≤ s. The "tied" relation is an equivalence relation and partitions the set S into tied subsets. Furthermore, the weak order on S gives a total order on this partition. In other words, we obtain an ordered set-partition S = S 1 ∪ · · · ∪ S k . Conversely, an ordered set-partition S = S 1 ∪ · · · ∪ S k determines a weak order on S by setting s ≤ t whenever s ∈ S i and t ∈ S j and i ≤ j. Thus, there is a bijection between weak orders on S and ordered set-partitions on S. In what follows, we will freely translate back and forth between weak orders and ordered set-partitions.
Similarly, a right consecutive embedding is a consecutive pair (i, i + 1) r such that T i ∪ T i+1 ⊆ S j for some (necessarily unique) S j . The phrase consecutive embedding will refer to either a left consecutive embedding or a right consecutive embedding.
Since weak orders can be viewed as ordered set-partitions, it makes sense to talk about consecutive embeddings in a pair of weak orders.  Of these, the five in the first row have no consecutive embeddings, whereas the first two in the second row each have a right consecutive embedding, and the last two in the second row each have a left consecutive embedding. Thus, q 2 = 5. Then (1 + 9) = 5 and

Inclusion-Exclusion
Recall that we have reduced the computation of the sequences {p n,0 } and {p n } to the computation of the sequences {q n,0 } and {q n }. The formula q n,0 = f 2 n from Proposition 3.8 gives a fast polynomial-time algorithm for computing the sequence {q n,0 }. We now restrict our attention to computing the sequence {q n }. This subsection will give an inclusion-exclusion formula for q n . 3. Let q n,k count the number of pairs of weak orders on {1, . . . , n} with k distinguished consecutive embeddings.
The notation q n,k does not conflict with the notation q n,0 because pairs of weak orders with 0 distinguished consecutive embeddings are in bijection with (ordinary) pairs of weak orders. The first few values of q n,k are given in Table 1. There are two main observations from this table that we will prove below. First, q n,k = 0 for k > n (Corollary 3.16). Second, n k=0 (−1) k q n,k = q n (Lemma 3.17). n q n q n,0 q n,1 q n,2 q n,3 q n,4 0  Table 1: Small Values of q n and q n,k We will need the notion of a chain in a pair of weak orders with distinguished consecutive embeddings.  Let (i, i + 1) l be a left consecutive embedding and let (i , i + 1) r be a right consecutive embedding. Then S i ∪ S i+1 ⊆ T j for some (necessarily unique) T j and T i ∪ T i +1 ⊆ S j for some (necessarily unique) S j . If these two consecutive embeddings have an element of {1, . . . , n} in common, then T j must be either T i or T i +1 , and S j must be either S i or S i+1 . Then which is a contradiction.
We can now determine q n,k for k ≥ n − 1.
Proof. Consider a pair of weak orders on {1, . . . , n} with k ≥ 1 distinguished consecutive embeddings. Let c ≥ 1 be the number of chains. Note that c + k counts the number of parts that are contained in the chains, because each chain has one more part than distinguished consecutive embedding. Note that these c + k parts have no elements in common, since any two parts within a chain have no elements in common, and Lemma 3.14 states that distinct chains have no elements in common. This gives the inequality c + k ≤ n.
If k ≥ n, then this is impossible, which shows that q n,k = 0. Now suppose that k = n − 1. Then c = 1, so there is one chain consisting of n − 1 overlapping distinguished consecutive embeddings. Thus, one side (either left or right) has n parts in one of the n! possible permutations, and the other side has just one part consisting of the unordered set {1, . . . , n}. This gives 2 · n! possibilities. Proof. If k > n, then k ≥ 1 so Proposition 3.15 gives q n,k = 0.
We now give the inclusion-exclusion formula for q n in terms of q n,k . Proof. If a pair of weak orders on {1, . . . , n} has m total consecutive embeddings, then m ≤ n by Corollary 3.16, so that pair of weak orders is counted m k=0 (−1) k m k times by the sum. This is equal to 0 if m ≥ 1 and equal to 1 if m = 0. Thus, the sum counts pairs of weak orders with no consecutive embeddings.

The Final Count
It remains to count pairs of weak orders on {1, . . . , n} with k distinguished consecutive embeddings. To deal with the fact that consecutive embeddings can overlap, we will work with chains. The first step of our formula will be to sum over the number of chains, which we will call c, and the total number of elements from {1, . . . , n} to put inside these chains, which we will call t. Keep in mind Lemma 3.14, which states that distinct chains have no elements in common, even if they are on different sides. In other words, each of the t elements will end up in exactly one of the c chains.
There are n t ways to choose which t elements from {1, . . . , n} to put inside the c chains. In the proof of Proposition 3.15 we saw that the c chains contain c + k total parts, so t ≥ c + k. Then there are (c + k)! t c+k ways to partition these t elements into an ordered sequence of c + k parts. Now we cut the ordered sequence of c + k parts into c chains, each of which must consist at least 2 parts. Since the parts are already ordered, this just requires choosing how many parts to put into each chain. The number of ways to do this equals the number of ways to distribute c + k identical balls among c distinguishable urns, each of which must get at least 2 balls. This is the same as the number of ways to surjectively distribute k identical balls among c distinguishable urns. There are k−1 k−c such ways. At this point, we have an ordered sequence of c chains, along with n − t non-chain elements. We can minimize the distinction between the c chains and the n − t elements by collapsing each chain into a single element, as depicted in Figure 3. For example, consider the top right chain in Figure 3. Both the left and right side of this chain get collapsed to a single element (which is labeled c 2 in Figure 3). The key difference between the left c 2 and the right c 2 is that the right c 2 must be in a part by itself. Thus, the circling in the collapsed form indicates which side of the collapsed chain must be in a part by itself. We will need a generalization of the Fubini numbers that allows for the possibility of elements that are required to be placed in a part by themselves. For 0 ≤ k ≤ n, let f n,k count the number of weak orders on {1, . . . , n} where each of the elements 1, . . . , k must be placed in a part by itself. These generalized Fubini numbers have the formula f n,k = n l=k l! n−k l−k , because there are n−k l−k ways to partition the set {k + 1, . . . , n} into l − k parts, and l! ways to order these l − k parts along with the remaining k elements {1, . . . , k}. The ordinary Fubini numbers are the special case when k = 0.
Returning to the situation prior to Figure 3, we now sum over the number of chains on the left side, which we will call j. This forces the number of chains on the right side to be c − j. We already have an ordering on the c chains, so we will take the first j chains to be the chains on the left side, and the last c − j chains to be the chains on the right side. After collapsing each of the c chains into a single element on each side, we have n − t + c total elements to distribute. On the left side, j of these elements are marked as having to be placed in a part by themselves. On the right side, c − j of these elements are marked as having to be placed in a part by themselves. Then there are f n−t+c,j f n−t+c,c−j ways to produce a valid pair of weak orders. However, since we already have an ordering for the c chains, we must divide by j! and (c − j)!. This gives the formula q n,k = We should point out that the second sum is zero for c > n−k. Plugging this into Lemma 3.17 and rearranging the summations gives the formula This proves our formula for p n (equation 1 from the introduction).
where the values f n,k = n l=k l! n−k l−k are generalizations of the Fubini numbers.

Equation 2
has been written so that there are several independent components inside the formula. The next subsection shows how we can precompute these components to speed up computation.

The Algorithm
We now give an algorithm to compute the values p 1 , . . . , p N in O(N 2 ) memory and O(N 3 ) time. The statement "O(N 2 ) memory and O(N 3 ) time" assumes that storing a number requires constant memory and that arithmetic operations require constant time. This assumption is not true in practice since the computation involves numbers that grow super-exponentially. However, the number of digits of these numbers grows polynomially. Then storing a number requires polynomial memory and arithmetic operations require polynomial time. Thus, even with these practical considerations, the algorithm uses polynomial memory and polynomial time. Empirically, the runtime for the algorithm grows like N 4.75 .
Here is the algorithm: 1. Compute the values n k , n k , n k , for 0 ≤ k ≤ n ≤ N using the standard recurrences.

2.
Compute the values f n,k = n l=k l! n − k l − k for 0 ≤ k ≤ n ≤ N .

Compute the values g
In this calculation, note that −1 0 = 1 and also k−1 k = 0 for k ≥ 1. We will give a faster and simpler algorithm for h t,c in the next subsection.  (2) for q n so that the outside sum is a sum over the values of n−t+c, which avoids expensive recomputation of the values g n−t+c,c . Version 2 is limited by time, rather than by memory.

Combinatorial Interpretation of h t,c
We conclude this section by giving a combinatorial interpretation of h t,c . This will give a fast and simple algorithm for h t,c , as well as a bound on h t,c that will be used when proving asymptotic formulas. Proof. First, we rewrite h t,c as Recall that the quantity (2) counts the number of ways to choose an ordered set-partition of {1, . . . , t} into c + k parts, which are further grouped into c chains of at least two parts each. Then the quantity counts the number of ways to choose an ordered set-partition of {1, . . . , t} into k parts, which are further grouped into c consecutive blocks of at least two parts each. Here's example with c = 3, k = 7, and t = 10: 3 5 4 8,10 2 1,7, 9 6 Considering the c blocks gives an ordered set-partition of {1, . . . , t} into c parts. If we let Ω denote the collection of all ordered set-partitions of {1, . . . , t} into c parts, then we can write where a π counts the number of ways that π ∈ Ω arises as the c blocks of an ordered set-partition of {1, . . . , t} into k parts, which are further grouped into c consecutive blocks of at least two parts each. If π is the ordered set-partition {1, . . . , t} = π 1 ∪ · · · ∪ π c , then we can write a π as Thus, Putting this all together gives We now invoke the identity n k=2 (−1) k k! n k = 2 if n is even, 0 if n is odd which follows from setting x = −1 in the identity where (x) k = x(x − 1) · · · (x − k + 1) denotes the falling factorial. Thus, if each π i has even cardinality, then π contributes (−1) c 2 c to h t,c , and otherwise π contributes nothing to h t,c . The result follows. The combinatorial interpretation of h t,c given in Proposition 3.19 also gives a recurrence for h t,c . It will be convenient to introduce an auxiliary sequence. Definition 3.21. We define integers T (n, k) for 0 ≤ k ≤ n by • T (n, 0) = 0 for n ≥ 1, • T (n, n) = 1 for n ≥ 0, The values T (n, k) maybe found at [5, A036969].
The significance of Proposition 3.22 is that it reduces the computation of h t,c in step 4 of the algorithm from cubic time to quadratic time. Let S(n, k) = (−1) k h 2n,k /2 k , which by Proposition 3.19 counts the number of ordered set-partitions of {1, . . . , 2n} into k parts of even cardinality. The proof of Proposition 3.22 boils down to showing S(n, k) satisfies the recurrence S(n, k) = k 2 S(n − 1, k) + k(2k − 1)S(n − 1, k − 1). This recurrence is known [5, A241171]. We include a proof of Proposition 3.22 for completeness.
The first property follows from the observation that if n ≥ 1, then there are no ordered set-partitions of {1, . . . , 2n} into 0 parts. For the second property, note that an ordered set-partition of {1, . . . , 2n} into n parts of even cardinality consists of pairing up {1, . . . , 2n} into an ordered sequence of n unordered pairs. There are (2n)!/2 n such pairings. It remains to show the third property (the recurrence). Consider the following two ways of constructing an ordered set-partition of {1, . . . , 2n} into k parts of even cardinality: • Start with an ordered set-partition {1, . . . , 2n − 2} = S 1 ∪ · · · ∪ S k into k parts of even cardinality, and choose two (possibly equal) indices i, j ∈ {1, . . . , k}. Consider the smallest element of the union S i ∪ S j and flip its position (if it was in S i , then move it to S j , and vice versa). Finally, add 2n − 1 to S i and add 2n to S j .
• Start with an ordered set-partition {1, . . . , 2n − 2} = S 1 ∪ · · · ∪ S k−1 into k − 1 parts of even cardinality, and choose an index k 0 ∈ {1, . . . , k}. Inserting the empty set at position k 0 and relabeling gives {1, . . . , 2n − 2} = S 1 ∪ · · · ∪ S k , which would be an ordered set-partition except that S k0 is empty. Now choose two (possibly equal) indices i, j ∈ {1, . . . , k}, at least one of which is equal to k 0 . Consider the smallest element of the union S i ∪ S j and flip its position (if i = j = k 0 , then do nothing). Finally, add 2n − 1 to S i and add 2n to S j .
We claim that every ordered set-partition of {1, . . . , 2n} into k parts of even cardinality arises uniquely from one of these two constructions. This is because the process can be reversed: • Start with an ordered set-partition of {1, . . . , 2n} = S 1 ∪ · · · ∪ S k into k parts of even cardinality. Let S i be the set containing 2n − 1 and let S j be the set containing 2n. Remove 2n − 1 from S i and remove 2n from S j . Consider the smallest element of the union S i ∪ S j and flip its position (if S i = S j = ∅, then do nothing). If S i and S j are still nonempty, then we are in the first case. If either S i or S j is now empty, then we are in the second case.
In the first case, there are S(n − 1, k) ways of choosing the initial ordered set-partition, and k 2 ways of choosing the indices i, j. In the second case, there are S(n − 1, k − 1) ways of choosing the initial ordered set-partition, k ways of choosing the index k 0 , and 2k − 1 ways of choosing the indices i, j. This proves the recurrence S(n, k) = k 2 S(n − 1, k) + k(2k − 1)S(n − 1, k − 1).

Asymptotics of f n,k
We first give a formula for the generalized Fubini numbers f n,k in terms of the ordinary Fubini numbers f n .
Proposition 4.1. If k ≥ 0 and n ≥ 0 are integers, then where each n i is a nonnegative integer.
Proof. Recall that f n+k,k counts the number of weak orders on {1, . . . , n + k} where each of the elements 1, . . . , k must be placed in a part by itself. We first choose the permutation of the elements 1, . . . , k. This leaves k + 1 regions between the elements 1, . . . , k where we must place the remaining n elements k + 1, . . . , n + k. We then choose how many elements n i should go into each of the k + 1 regions. There are n! n0!···n k ! ways to choose how to distribute the n elements k + 1, . . . , n + k into these k + 1 regions. Finally, there are f ni choices for the weak order on the n i elements within each of the k + 1 regions.
We remark that Proposition 4.1 gives a generating function identity It is possible to prove Lemma 4.3 below using complex-analytic generating function techniques by looking at the poles of k!/(2 − e x ) k+1 . However, we will prove Lemma 4.3 using more direct real-analytic techniques. We will need the following technical lemma.
n=0 be a sequence of real numbers and let k ≥ 0 be an integer. If a n → 1, then 1 n+k k n0+···+n k =n all ni≥0 a n0 · · · a n k → 1 as n → ∞.
Proof. We first remark that n+k k is the number of terms in the sum, since the number of solutions of n 0 + · · · + n k = n equals the number of ways to distribute n identical balls into k + 1 distinguishable urns. Let ε > 0. Then there exists an m ≥ 0 such that a n ∈ [1 − ε, 1 + ε] for n ≥ m. Let C = max |a n |. For n ≥ (k + 1)m, we will split up the sum as n0+···+n k =n a n0 · · · a n k = n0+···+n k =n all ni≥m a n0 · · · a n k + n0+···+n k =n some ni<m a n0 · · · a n k .
The number of terms of the first sum equals the number of solutions to n 0 +· · ·+n k = n with all n i ≥ m, which equals the number of solutions to n 0 + · · · + n k = n − (k + 1)m. Then the first sum has n−(k+1)m+k k terms (by the same argument as at the start of the proof), each of which lies in the interval The second sum has n+k k − n−(k+1)m+k k terms, each of which lies in the interval [−C k+1 , C k+1 ]. Thus, a n0 · · · a n k and n0+···+n k =n a n0 · · · a n k ≤ (1 + ε) k+1 n − (k + 1)m + k k + C k+1 n + k k − n − (k + 1)m + k k .

Asymptotics of p n
We will need the following technical lemma.
Lemma 4.6. For each fixed nonnegative integer k, we have n n − k ∼ n 2k 2 k k! as n → ∞. We also have n n − k ≤ n 2k 2 k k! for all n ≥ k.
We can now prove the asymptotic formula for p n (Theorem 1.2 from the introduction).   Exponentiating and rearranging terms gives This suggests the following conjecture.

Congruence Conjecture
Theorem 2 of [7] states that if p is prime and n ≥ m, then where ϕ is Euler's totient function. In other words, the Fubini numbers are periodic modulo p m . Squaring both sides gives the congruence q n+ϕ(p m ),0 ≡ q n,0 (mod p m ) for p prime and n ≥ m. Replacing q n,0 with q n suggests the following conjecture, which is supported by the computed values of q n for n ≤ 5000.
Conjecture 5.2. If p is prime and n ≥ m, then q n+ϕ(p m ) ≡ q n (mod p m ).
Unfortunately, this conjecture does not appear to imply anything about the behavior of the sequence {p n } modulo p m , because of the division by n! when converting from {q n } to {p n }.

Acknowledgements
Many thanks to Sara Billey for proposing the problem and for helpful discussions, and to the WXML (Washington Experimental Mathematics Lab) for providing the project that led to this paper. Also thanks to Sara Billey and the anonymous referees for their comments and suggestions on this paper.

Web Resources
The basic and memory-optimized implementations of the algorithm described in section 3.6 and the values of p n for n ≤ 5000 can be found at https://github.com/tb65536/ParabolicDoubleCosets