Using functional equations to enumerate -avoiding permutations
We consider the problem of enumerating permutations with exactly occurrences of the pattern and derive functional equations for this general case as well as for the pattern avoidance () case. The functional equations lead to a new algorithm for enumerating length permutations that avoid . This approach is used to enumerate the -avoiders up to . We also extend those functional equations to account for the number of inversions and derive analogous algorithms.
Let be a sequence of distinct positive integers. We define the reduction of this sequence, denoted by , to be the length permutation that is order-isomorphic to (i.e., if and only if for every and ). Given a (permutation) pattern , we say that a permutation contains the pattern if there exists such that , in which case we call an occurrence of . We also define to be the number of occurrences of pattern in the permutation . For example, if the pattern , the permutation avoids the pattern (so ), whereas the permutation contains two occurrences of (so ).
For a pattern and non-negative integer , we define the set
and also define . For the case, we say that the two patterns and are Wilf-equivalent if for all . Additionally, two patterns that are Wilf-equivalent are said to belong to the same Wilf-equivalence class. Note that the case corresponds with “classical” pattern avoidance, which has been well-studied. Work on the more general problem () has usually been restricted to patterns of length and small . Some recent work include [6, 5, 7, 8, 11, 14, 19, 20].
A little more is known for the pattern avoidance problem, but the problem quickly gets difficult. For each , it is well known that (the Catalan numbers) . For the length patterns, there are three cases (Wilf-equivalence classes) to consider: , , and . The enumeration for the -avoiding permutations was solved by Gessel in . Later, Bóna solved the case for the -avoiding permutations in . The pattern , however, has been notoriously difficult to enumerate.
There is currently no non-recursive formula known for computing , and precise asymptotics are not known either. Marinov and Radoičić developed an approach in  using generating trees and computed for . Another approach (using insertion encoding) was used by Albert et al. in  to compute for (five more terms). Given the difficulty of this pattern, Zeilberger has even conjectured that “Not even God knows ” .
Given the difficulty of exact enumeration, work has also been done on studying how the sequence grows for various patterns. We define the Stanley-Wilf limit of a pattern to be:
Thanks to results by Arratia  and Marcus and Tardos , we know that for each pattern , the limit exists and is finite. For patterns of length three, the Stanley-Wilf limit is known to be . Additionally, Regev  showed that , while Bóna’s result in  gives us . The exact limit for the pattern , however, is still unknown. The best known lower bound is by Albert et al. , who showed that . The best known upper bound has seen some improvements in recent years. The recent “best” upper bound was improved by Claesson, Jelínek, and Steingrímsson in  to . That approach was then refined by Bóna in  to show that .
Additionally, Claesson, Jelínek, and Steingrímsson conjectured that the number of length permutations avoiding with exactly inversions was non-decreasing in (for each fixed ). They show that if this conjecture holds, then . Neither the current lower bound nor this potential new upper bound appear to be “close” to the exact value of . For example, Steingrímsson’s survey paper  considers empirical data and suggests that it may be close to . Some data we consider in this paper suggests a similar story.
The paper is organized in the following manner. In Section 2, we derive functional equations for computing . Furthermore, the approach is specialized to the avoidance case () to derive an algorithm for enumerating the -avoiding permutations. In Section 3, we describe technical details on the algorithm and use it to compute up to , giving us new terms. We use this new data to make some empirical observations. In Section 4, we extend the functional equations to keep track of the number of inversions. We conclude with a few final remarks in Section 5.
2 Functional equations for the pattern
We begin by extending the functional equations approach in  to the pattern . We will first derive functional equations that can be used to compute . The approach will be presented in full detail so that this article is self-contained. We then specialize this approach to the case and develop a new algorithm for enumerating the -avoiding permutations.
2.1 A general approach to
Given a non-negative integer , we define the polynomial (in the variable )
Observe that the coefficient of in is exactly equal to .
In addition to the variable , we introduce catalytic variables with and catalytic variables with . Note that the subscripts of the two sets of catalytic variables range over different quantities. We define the weight of a permutation to be
where it is always assumed that . For example,
In essence, the weight stores information about patterns, patterns (which may become the “” of a occurrence), and patterns (which may become the “” of a occurrence) through the exponents of the variable , the variables , and the variables , respectively.
We will define a multi-variate polynomial on all the previously defined variables. For notational convenience, we first write the variables and the variables as matrices of variables:
where we will disregard the entries below the diagonal in and the entries above the diagonal in .
For each , we now define the multi-variate polynomial:
Observe that is our desired polynomial, where is the matrix of all ’s. For a fixed , our goal is to quickly compute the coefficient of in , which is exactly . We will do this by deriving a functional equation for . This follows readily from the following result:
Let and suppose that . If , then
where is the set of substitutions given by
We assume to be a fixed value. Observe that is equal to the number of occurrences of in plus the number of occurrences of in , where the term corresponding to the “” is greater than .
If we re-insert at the beginning of , we would shift all the terms up by . This gives us the substitutions:
We make a few more observations. First, in , the exponents of and are equal and the exponents of and are equal for each (since ). This gives the substitutions if and if .
Second, the number of patterns that include the first term is the sum of the exponents of for . The substitution changes to (for ).
Third, the number of patterns that include the first term (i.e., the “” term is equal to ) and whose “” term is at least is equal to the sum of the exponent of for . The substitution changes to (for ).
This gives us our substitution set . Finally, the new “” would create new patterns and would require an extra factor of for the weight. ∎
Now, we define the operator on an arbitrary square matrix and to be:
If , then is defined to be the matrix obtained by deleting the -th row and -th column from . In essence, the operator deletes the -th row, merges the -th and -th columns (via term-by-term multiplication), and multiplies this new column by a factor of . It is important to note that while this operator is defined on any matrix, it will only be applied to our “matrix of variables” to get a smaller matrix.
In addition, we define another operator on two square matrices and and to be:
If , then is defined to be the matrix obtained by deleting the -st row and -st column from . In essence, the operator modifies by deleting the -th row, merging the -th column with the -th column (via term-by-term multiplication), and scaling that new column by products of terms from .
Lemma 1 now directly leads to the following:
For the pattern ,
Although all the entries in the matrices are changed for consistency and notational convenience, we will continue to disregard the entries below the diagonal in subsequent matrices and the entries above the diagonal in subsequent matrices . We can apply the same computational tricks shown in [18, 17]. For example, it is not necessary to compute completely symbolically and substitute and at the end. Instead, we can apply functional equation (FE1324) directly to and subsequent terms. We may also use the following lemma from , which is obvious from the definition of the operator :
Let be a square matrix where every row is identical (i.e., the -th row and the -th row are equal for every ). Then, will also be a square matrix with identical rows.
By Lemma 2, repeated applications of the operator to the all ones matrix will still result in a matrix with identical rows. It is therefore sufficient to keep track of only a single row. It is also helpful to note that repeated applications of to the matrix will result in a matrix whose entries are powers of .
While the lemma does not apply to the operator, this still allows us to simplify the polynomial and its functional equation by reducing the number of catalytic variables from variables to variables. Let denote the polynomial where every entry of the matrices and are powers of and every row in is . We get the analogous functional equation:
Repeatedly applying this recurrence to allows us to compute our desired polynomial since . Extracting the coefficient of from this polynomial gives us .
Additionally, for a fixed , the sequence can be computed more quickly by discarding higher powers of (just as in [18, 17]). Although this approach is too memory intensive for larger , for small , this method is still much faster than naive methods that construct the set . The approach has been implemented in the procedure F1324rN(r,N) in the accompanying Maple package F1324.
For example, the Maple call F1324rN(0,19); for the first terms of produces the sequence:
and the Maple call F1324rN(1,17); for the first terms of produces the sequence:
2.2 Specializing to
Unfortunately, the previous algorithm developed for the pattern is very memory intensive, even for . In this subsection, we outline how to extract a simpler recurrence specifically for the pattern avoidance case.
We will specialize for the case beginning at functional equation (FE1324c). Recall that is the polynomial given by where every entry of the matrices and are powers of and every row in is . We had the functional equation
and wanted to compute and extract the coefficient of , which is exactly .111Recall that is the matrix where every entry is .
Since all the variables (from matrix ) and represent powers of , it is sufficient to keep track of powers of through most of the algorithm. This allows us to consider the analogous function , where is an matrix of non-negative integers and each is a non-negative integer. More precisely, is the polynomial , where and are matrices, for every , and every row of is .
In addition, we define the analogous operator on an square matrix (of non-negative integers), a length vector of non-negative integers , and :
If , then is defined to be the matrix obtained by deleting the -st row and -st column from . In essence, the operator modifies by deleting the -th row, merging the -th column with the -th column (via term-by-term addition), and adding partial sums of into the new column.
We now have the functional equation (analogous to Eq. (FE1324c)):
where . Observe that is now our desired polynomial .222We denote the matrix consisting of all zeros by .
Since we are specifically considering the case, we can make additional observations and simplifications. First, we are only interested in the constant term of . As in [18, 17], we only need to keep track of polynomials of the form in intermediate computations, where represents permutations with at least one occurrence of the pattern we are tracking. Because of this, we may consider all matrices/vectors used in to be - matrices/vectors. After every addition (for example, in the term in ), we can take the minimum of the resulting sum and .
Next, observe that only appear in the operator, and in particular, in the partial sums for in Eq. 7. Suppose that some of the are equal to , and let be the smallest number such that . Then,
In particular, the variables are unnecessary, and it is sufficient to keep track of how many ’s there are. We can then consider this slightly simpler function
Finally, observe that whenever in (FE1324e), we can discard the entire term since we are only interested in the constant term of the final polynomial. Observe that if and only if for some . This observation (combined with how the operator “modifies” the matrix ) implies that we only need to keep track of the left-most within each row of . If there are multiple ’s on a row, the left-most is sufficient to force as long as it is not in the -th column. Therefore, we can consider a function of the form
where and for each and is the matrix where the -th row is if and otherwise is .
This approach is implemented in the Maple package F1324 in the procedure AV1324(n). An improved implementation in C++ is provided in the program av1324.cpp and is discussed in greater detail in the next section.
3 Computational details and results
3.1 Algorithmic details for enumerating -avoiders
For notational convenience, let denote the number of 1324-avoiding permutations of length . The values up to were previously computed in  and subsequently listed in A061552 of OEIS (). To extend this list, we have written a C++ implementation of the algorithm described in the previous section. The implementation is in the program av1324.cpp, which is available from the authors’ websites.
The main part of the program is a recursive function which takes as input an integer , an integer , and a vector of integers satisfying (these correspond to the in as , being zero-aligned since this is more natural in C++). We represent as an array of bytes, zero-padded to a fixed maximum length of 32 (sufficient to compute up to ). The output of is an integer, which we represent as a 128-bit unsigned integer (sufficient to compute up to at least , since and all recursive calls to must produce values that are no larger than the output).
We use full memorization to reduce the number of recursive calls that have to be made. The std::map type in the C++ standard library is used to associate vectors with output values , using one such map for each pair . The std::map type implements a self-balancing binary tree with insertion and lookup time where is the size of the cache. The keys are compared lexicographically by casting to 64-bit integers and processing eight bytes at a time.
We compiled the program with GCC 4.3.4 and ran it on the MACH computer at the Johannes Kepler University of Linz, using a 2.66 GHz Intel Xeon E7-8837 CPU with 1024 GiB of memory allocated to the process. The memory limit allowed computing up to . The new terms are:
We computed all consecutively in one run to benefit from already cached function values (computing a single in isolation would not give any significant memory savings). With already computed, the evaluation of took 33 hours and used 920 GiB of memory, and the computation as a whole took 50 hours.
Detailed results from the computation are presented in Table 4 (in the Appendix). Besides the running time and total memory usage, we record the number of cache hits (calls to replaced by cache lookups) and the number of cache misses (calls to that require evaluation). Up to , we also show the number of calls to used by the algorithm with caching disabled, measured in a separate run of the program.
All quantities appear to grow slightly faster than exponentially. Over the measured range, the number of function calls (as well as the size of the cache in memory) is roughly proportional to . If the memory overhead of the implementation were reduced by half (or if we had a computer with twice as much memory), we could thus compute roughly one more entry. With caching disabled, the number of recursive function calls grows much faster, making this version impractical (for , the number of calls is roughly ).
At present, we do not know if the algorithm can be modified to significantly reduce the memory consumption required by full memorization, without significantly increasing the running time. Such a modification might allow computing several more entries in the sequence , particularly if the job could be parallelized.
3.2 Observations on the asymptotics
Since the enumeration problem is solved for the patterns and , it is not hard to derive asymptotic information on the sequences enumerating their respective pattern avoiders. For the pattern , we have
and for the pattern , we have
The convergence occurs fairly quickly, even when observing the first terms. On the other hand, it is not even known if the asymptotics for look like (for constants and ). Shalosh B. Ekhad was kind enough to compute some numerical data for the authors.
Using the first terms, the first terms, and the first terms, the approximate value for and for the patterns and are:
However, when the same guessing is done for the pattern , the convergence is much slower. The values are:
It should be re-emphasized that it is not known whether this sequence fits the asymptotic form of . The empirical asymptotics suggests that either the convergence is much slower than the other length patterns or that the sequence does not have that asymptotic form to begin with (unlike the other patterns). In addition, the values would suggest that the Stanley-Wilf limit is at least . More detailed numerical data on the asymptotics (from Shalosh B. Ekhad) can be found on the authors’ websites.
4 Extending to inversions
In this section, we show how the previous functional equations can be adapted to refine the values by the number of inversions. The number of inversions in a permutation is one of the most commonly studied permutation statistic and in essence, quantifies how “unsorted” a permutation is. Given a permutation , the inversion number of , denoted by , is the number of pairs such that and . An equivalent definition is that , the number of patterns in .
We again consider the pattern . For each , we define the bivariate polynomial
Observe that is exactly from Section 2.
For the pattern , the polynomial can now be “generalized” as
Note that is exactly from Section 2.
We now make an important observation. Given a permutation , suppose that . Then, is equal to the number of inversions in plus the number of terms in that are less than (which is exactly ). For any previously define functional equation, it is enough to insert a factor of into the summation.
We can now quickly derive the modified functional equations for the pattern. The functional equation (FE1324) now becomes:
For the pattern ,
Similarly, the functional equation (FE1324c) for now becomes the analogous:
For the pattern ,
The functional equation (FE1324e) for , which corresponds specifically to the pattern-avoidance case, also becomes the analogous:
For the pattern ,
The enumeration algorithm derived from the functional equation (qFE1324e) has been implemented in the procedure qAv1324r(n,r,q) in the Maple package F1324. It is also worth noting that the same extension can be done for tracking non-inversions by inserting into the summations (instead of ). This has also been implemented in the Maple package F1324 in the procedure pAv1324r(n,r,p).
In , Claesson, Jelínek, and Steingrímsson consider refining the number of -avoiding permutations by the inversion number to prove a new upper bound on the Stanley-Wilf limit. They conjectured that for each fixed , the number of -avoiders with exactly inversions is non-decreasing in . If this conjecture is proven true, their result would improve the upper bound of the growth rate. Using our algorithm, we are able to empirically confirm this conjecture for , and the explicit values can be found on the authors’ websites. We suspect that a more careful analysis of the functional equations and the algorithm may lead to new insights regarding this conjecture.
In this paper, we extended the functional equations approach of [18, 17] to derive functional equations as well as an enumeration algorithm for computing . We were able to specialize this approach to to develop a new enumeration algorithm for computing the number of -avoiding permutations. This new approach was used to compute new terms for the sequence as well as make some empirical observations on the asymptotics. The functional equations were also extended to refine by the number of inversions.
While some of the key contributions of this paper are algorithms, it is important to note that all the intermediate steps were rigorous functional equations. In particular, the specialized functional equations for enumerating -avoiders can be viewed as recursively defined functions on - matrices. We suspect that a more careful analysis of these functions may yield insight into the sequence enumerating the -avoiders and perhaps the Stanley-Wilf limit . In addition, the functional equations for tracking inversions were used to confirm the conjecture by Claesson, Jelínek, and Steingrímsson for up to . We suspect that a more careful analysis of these functional equations may also provide insight toward resolving their conjecture, which would then lower the upper bound on to .
Acknowledgments: We would like to thank Doron Zeilberger and Manuel Kauers for their very helpful comments and suggestions. We would also like to thank Shalosh B. Ekhad for computing numerical data on the asymptotics for us.
-  M. H. Albert, M. Elder, A. Rechnitzer, P. Westcott, and M. Zabrocki. On the Stanley-Wilf limit of 4231-avoiding permutations and a conjecture of Arratia. Adv. in Appl. Math., 36(2):96–105, 2006.
-  Richard Arratia. On the Stanley-Wilf conjecture for the number of permutations avoiding a given pattern. Electron. J. Combin., 6:Note, N1, 4 pp. (electronic), 1999.
-  Miklós Bóna. A new upper bound for -avoiding permutations. arxiv:1207.2379 [math.co], 2012.
-  Miklós Bóna. Exact enumeration of -avoiding permutations: a close link with labeled trees and planar maps. J. Combin. Theory Ser. A, 80(2):257–272, 1997.
-  Miklós Bóna. The number of permutations with exactly -subsequences is -recursive in the size! Adv. in Appl. Math., 18(4):510–522, 1997.
-  Miklós Bóna. Permutations with one or two -subsequences. Discrete Math., 181(1-3):267–274, 1998.
-  Alexander Burstein. A short proof for the number of permutations containing pattern 321 exactly once. Electron. J. Combin., 18(2):Paper 21, 3, 2011.
-  David Callan. A recursive bijective approach to counting permutations containing -letter patterns. arxiv:math/0211380 [math.co], 2002.
-  Anders Claesson, Vít Jelínek, and Einar Steingrímsson. Upper bounds for the Stanley-Wilf limit of 1324 and other layered patterns. J. Combin. Theory Ser. A, 119(8):1680–1691, 2012.
-  Murray Elder and Vince Vatter. Problems and conjectures presented at the third international conference on permutation patterns. arxiv:math/0505504 [math.co], 2005.
-  Markus Fulmek. Enumeration of permutations containing a prescribed number of occurrences of a pattern of length three. Adv. in Appl. Math., 30(4):607–632, 2003.
-  Ira Gessel. Symmetric functions and P-recursiveness. J. Combin. Theory Ser. A, 53(2):257–285, 1990.
-  Donald E. Knuth. The art of computer programming. Vol. 1: Fundamental algorithms. Addison Wesley, Reading, Massachusetts, 1973.
-  Toufik Mansour and Alek Vainshtein. Counting occurrences of 132 in a permutation. Adv. in Appl. Math., 28(2):185–195, 2002.
-  Adam Marcus and Gábor Tardos. Excluded permutation matrices and the Stanley-Wilf conjecture. J. Combin. Theory Ser. A, 107(1):153–160, 2004.
-  Darko Marinov and Radoš Radoičić. Counting 1324-avoiding permutations. Electron. J. Combin., 9(2):Research paper 13, 9 pp. (electronic), 2002/03. Permutation patt