σ(1) > σ(2) < σ(3) > σ(4) . . .

The aim of this paper is to study permutations of {1, 2, . . . N } with a prescribed descent set A ⊂ [1, N − 1]. We first give a method to compute this number efficiently, namely, using O(N 2 ) elementary computations, while the classical formula based on an inclusion-exclusion argument uses an exponential number of elementary computations. Next, given a periodic pattern of ascents and descents, we give a complete expression of the exponential generating function of permutations whose descent set follows this pattern. Finally, we introduce an algorithm generating at random, with uniform probability, permutations with a given descent set.

These results make use of a generic method, which we call the density method, and which applies in a more general framework, namely the study of (finite) posets. Generally speaking, if P os is a finite poset, its order polytope is the subset of [0, 1] P os consisting of all functions f : P os → [0, 1] satisfying, for all x, y ∈ P os, x ≤ y =⇒ f (x) ≤ f (y). The density method consists in expressing the marginal densities of the uniform measure on the order polytope by computing iterated integrals. This is not always possible for all types of posets but in the context of this paper, we shall see that this is relevant. An interest of this method is that it provides an explicit algorithm of random generation with polynomial cost. Other uses of this method can be found in [Mar16,BMW18].

Our first theorem is the following:

The reason why we use descending induction, rather than regular induction, will appear in Theorem 3. An easy inclusion-exclusion argument (see for instance [Bon12], Chapter 1) gives

where the function g N is defined as follows.

However, in order to use (1), one needs to compute 2 N −|A| terms. By comparison, Theorem 1 only requires the computations of the coefficients of N polynomials of degree 0 to N − 1, which means computing O(N 2 ) coefficients. Moreover, each coefficient can be computed using only one elementary computation, except for the constant coefficients. It is easily seen that calculating the constant coefficient of f d requires N − d elementary computations. Hence the number of elementary computations in Theorem 1 is O(N 2 ). Another efficient method to compute Q N is given by Viennot [Vie79].

Theorem 1 also gives access to comparison results. For a permutation σ of Of course, such results are more difficult to derive using alternating sums such as (1). The fact that the maximum of F is achieved for alternating permutations was first observed by Niven [Niv68], see also [BR70]. Our next result deals with the case of a periodic descent set.

Theorem 1.2 Let p ≥ 2 be an integer, fix a set A ⊂ {1, 2, . . . p} and put

For each integer n ≥ 1, let d n be the number of permutations of {1, . . . , n} with descent set

Then there exists a rational function R A ∈ Z(X 1 , . . . , X p(p+1) ), such that

Theorem 2 generalizes the classical theorem by André [And1879] for alternating permutations:

where we agree that d 0 = 1. More recent results were established by Carlitz (see [Car75] for instance) and Mendes-Remmel-Riehl [MRR10] in the case µ = 1, and later by Luck [Luc14] and Basset [Bas16]. We shall see how to recover these results from our method in Section 3.

For the first-order asymptotics, see [LM83] and [BHR03]. For every set A, one can compute R A and derive the asymptotics from the zeros of the denominator, using the classical tools of analytic combinatorics. See for instance the analysis of alternating permutations as Example IV.35 in [FS09]. However, our approach does not provide the general form of the denominator.

Finally, we introduce an algorithm generating uniform random permutations with a given descent set. We begin by a simple remark, which was already used in [ELR02]

and so on. To recover σ Y from Y , one can use a sorting algorithm.

One can define the descent set D(Y ) of a sequence of reals Y as for a permutation and clearly, D(Y ) = D(σ Y ). Moreover, if Y is chosen according to the Lebesgue measure on [0, 1] N , then σ Y is uniformly distributed over all permutations of {1, 2, . . . N }. As a consequence, if Y is chosen uniformly at random among all sequences with descent set A, that is, if its density with respect to the Lebesgue measure on [0, 1] N is

for some constant C > 0, then σ Y is uniformly distributed over all permutations with descent set A. Therefore, we want to find an algorithm constructing a random sequence with descent set A.

Algorithm Using iid random variables (U 1 , . . . U N ), uniform on [0, 1], and the functions f i defined in Theorem 1, construct a sequence Y = (Y 1 , . . . Y N ) as follows:

Finally, recover the permutation σ Y from Y by sorting.

Theorem 1.3 The algorithm described above yields a random permutation of {1, 2, . . . N }, uniformly distributed over all permutations with descent set A.

To see that the algorithm is well-defined, suppose for instance that i ∈ A. Then (3) becomes

Since the functions f i , f i+1 are nonnegative on the interval [0, 1], (4) clearly has a unique solution on [0, 1].

We shall not study in detail the complexity of the algorithm. Nevertheless, let us make a few comments. First, remark that a simple rejection algorithm would have exponential complexity. Indeed, if we draw a permutation at random, the most likely profile is the alternating case, and the probability for a random permutation of length N to be alternating is ∼ 4 π (2/π) N (see Example IV.35 in [FS09]). Therefore, the average number of times one has to draw a permutation before finding one with the prescribed descent set is at least c(π/2) N .

Using our algorithm, we need to compute the functions f i , which can be done using O(N 2 ) elementary computations, as already seen. The last step, namely sorting the sequence Y to deduce σ Y , has an average complexity O(N log N ) if one uses randomized quicksort. What is more intricate is to evaluate the time needed to generate the variables Y i , which amounts to solving N equations of the form (4). We shall not discuss this from a theoretical point of view. However, some experimental results are given at the end of the paper. Typically, our algorithm can generate permutations of length a few thousands.

In the particular case of alternating permutations, other algorithms with a better complexity can be found in [BRS12].

We now proceed to the proof of the results. We first prove Theorem 3 in Section 2 and deduce Theorem 1 and Corollary 1 in Section 3. In Section 4, we give an outline of the proof of Theorem 2. In the last section , we comment on the implementation in Maple of our algorithm.

2 Proof Theorem 1.3 First, observe that there are obvious symmetries in the problem. If one replaces σ(i) with N + 1 − σ(i) for each i, then A is replaced with its complement. If σ(i) is replaced with σ(N + 1 − i) for each i, then A is replaced with N − A.

Let us prove Theorem 1.3. One can re-formulate the algorithm, using the notion of conditional density of a random variable, as follows:

As a consequence, the density of the sequence Y = (Y 1 , . . . , Y N ) with respect to the Lebesgue measure on [0, 1] N is equal to the telescopic product

which is exactly the form required in (2). Therefore, σ Y is uniformly distributed over all permutations with descent set A. This proves Theorem 3.

3 Proof of Theorem 1.1

Let us deduce Theorem 1. If a permutation of {1, . . . n} is drawn uniformly at random, the probability that its descent set is A is

But this probability is also equal to Q N (A) N ! Using the formula for the density of Y , we find that

Finally, let us prove Corollary 1. It suffices to prove that if B, B are two subsets of [2, N −1] with B = B ∪{m}, m / ∈ B, then the number of permutations with set of extrema B is smaller than the number of permutations with set of extrema B. Using the symmetries of the problem, it suffices to count permutations with set of extrema B, B ending with a descent. If we impose this condition, let A, Ã be the corresponding descent sets. Using A, Ã we define by descending induction the polynomials f i , fi as in Theorem 1.

Remark that for i ≤ N − 1, the h i are increasing functions on [0, 1] such that h i (0) = 0. Besides, one can directly define the h i by descending induction as follows. First, h N = 1 and then, if i / ∈ B,

Likewise, define the polynomials hi , 1

Since B = B ∪ {m}, we have hi = h i for i > m. Moreover, using ( 5) and ( 6), we get that hm (x) > h m (x) for every x ∈ (0, 1] and by induction, for every i < m and every x ∈ (0, 1], hi (x) > h i (x). Integrating this inequality for i = 1, we get the result.

Remark 3.1 The coefficients of the polynomials f i can be computed by induction. Put

By recursion, this leads to (1), which provides an alternative proof of Theorem 1.

Define by induction the sequence of polynomials (f n ) by

We claim that

Indeed, Theorem 1 tells us that the right-hand side is the number of permutations with length N and descent set N − D N . Using the symmetries of the problem, we see that this is equal to the number of permutations with length N and descent set D N . Consider the bivariate generating function

Then we have

We want to prove that the right-hand side of this equality has the form given by Theorem 2. Let µ be the number of ascents in a period, µ = p − |A|. By differentiating,

Solutions of this differential equation have the general form

where the ω k are the p-th roots of 1 if µ is even and the p-th roots of −1 if µ is odd, and where the a k are power series: All these equations can be summarized in a matrix system with integer coefficients. We skip the details of the resolution of this system (this is a bit long but only uses elementary linear algebra) but eventually we find

where, for all k, l ∈ [1, p], the function

a k,np+l t np+l has a rational expression, with integer coefficients, in terms of the ω k and of the functions

This proves Theorem 2.

We have implemented our algorithm in Maple1 . The length of the permutation is nmax. The random variables U [n], which are drawn independently at random, represent the descent profile: U [n] = 0 or 1 according as whether nmax − n is a descent or an ascent. The reals V [n] in the code correspond to the random variables Y i in the algorithm. As pointed out in the introduction, generating the variables Y i amounts to solving N equations of the form (4). Of course, we can only get approximate solutions and since the Y i are computed recursively, there is a propagation of errors. We shall not discuss this point from a theoretical point of view here. Observe, however, that these errors do not affect the permutation we generate as long as, for each i, the difference between Y i and its approximation is smaller than the minimal value of 2|Y i − Y j | over all i, j ∈ [1, nmax].

Since equations of the form (4) are polynomial and since f i = ±f i−1 , we can directly implement Newton's method in the algorithm. We stop the iterations in Newton's method as soon as two successive approximations of the solution are at distance less that 10 −p , where p is a parameter that can be tuned. We give the time needed to generate permutations of variable length and for various values of p. For a permutation of length 2000, we have run the algorithm for the values p = 4 and p = 5 and we observe empirically that this yields the same permutation.

Here are some tables giving the run time of the algorithm for permutations of various lengths. In the first table, we have also let the parameter p vary. The first line is the length of the permutation, the second line the run time (in seconds) for p=4 and the third line the run time for p=5. Empirically, the run time is more than quadratic. This is due to the fact that we have to perform computations with polynomials of large degree. Also, the quantity of memory is quite large if we want to record all the coefficients of all the polynomials. Of course, only keeping in memory the polynomial of highest degree and the descent profile would be sufficient, since the other polynomials can be obtained by differentiation. However, this would mean computing these polynomials twice.