=Paper=
{{Paper
|id=Vol-2113/paper19
|storemode=property
|title=The density method and permutations with a prescribed descent set
|pdfUrl=https://ceur-ws.org/Vol-2113/paper19.pdf
|volume=Vol-2113
|authors=Philippe Marchal
|dblpUrl=https://dblp.org/rec/conf/gascom/Marchal18
}}
==The density method and permutations with a prescribed descent set==
https://ceur-ws.org/Vol-2113/paper19.pdf
The density method and permutations
with a prescribed descent set
Philippe Marchal
LAGA (UMR CNRS 7539), Université Paris 13
marchal@math.univ-paris13.fr
Abstract
We give a formula to compute the number of permutations with a
prescribed descent set in quadratic time. We give the generating func-
tion of the number of permutations with a periodic descent set. We
introduce the density method algorithm, which generates uniformly dis-
tributed random permutations with a prescribed descent set.
1 Statement of the results
Let σ be permutation of {1, 2, . . . N }. The descent set of σ, D(σ), is the set of integers i ∈ [1, N − 1] such that
σ(i) > σ(i + 1). For instance, D(id) = ∅. If σ is an alternating permutation, that is, if
σ(1) > σ(2) < σ(3) > σ(4) . . .
then D(σ) = [1, N − 1] − 2N.
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:
Theorem 1.1 Fix an integer N ≥ 2 and a set A ⊂ [1, N − 1]. Define by descending induction a sequence of
polynomials (fi , 1 ≤ i ≤ N ) as follows:
fN = 1
Copyright © by the paper’s authors. Copying permitted for private and academic purposes.
In: L. Ferrari, M. Vamvakari (eds.): Proceedings of the GASCom 2018 Workshop, Athens, Greece, 18–20 June 2018, published at
http://ceur-ws.org
179
If i ∈ A, Z x
fi (x) = fi+1 (y) dy
0
If i ∈
/ A, Z 1
fi (x) = fi+1 (y) dy
x
Let QN (A) be the number of permutations of {1, 2, . . . N } with descent set A. Then
Z 1
QN (A) = N ! f1 (y) dy
0
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
X
QN (A) = (−1)|A|−|B| gN (B) (1)
B⊂A
where the function gN is defined as follows. If B = ∅, then gN (B) = 1. Otherwise, there exists k ∈ [1, N − 1]
such that B has the form B = {i1 < i2 < . . . < ik } and in that case
N!
gN (B) =
(N − ik )!(ik − ik−1 )! . . . (i2 − i1 )!i1 !
However, in order to use (1), one needs to compute 2N −|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 fd requires N − d elementary
computations. Hence the number of elementary computations in Theorem 1 is O(N 2 ). Another efficient method
to compute QN is given by Viennot [Vie79].
Theorem 1 also gives access to comparison results. For a permutation σ of [1, N ], say that m ∈ [2, N − 1] is a
local extremum if either σ(m − 1) < σ(m) > σ(m + 1) or σ(m − 1) > σ(m) < σ(m + 1). Then we have
Corollary 1.1 For B ⊂ [2, N − 1], let F (B) be the number of permutations of [1, N ] with set of local extrema
B. Then F : P({2, 3, . . . n − 1}) → N is an increasing function.
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
∞
[
A∞ = (A + np)
n=0
For each integer n ≥ 1, let dn be the number of permutations of {1, . . . , n} with descent set
A∞ ∩ [1, n − 1]
Let ω1 , . . . , ωp be the p-th roots of 1 (resp. -1) if p − |A| is even (resp. odd). For k, l ∈ [1, p], put
∞
X (ωk t)np+l
gk,l (t) =
n=0
(np + l)!
Then there exists a rational function RA ∈ Z(X1 , . . . , Xp(p+1) ), such that
∞
X dn tn
= RA (ω1 , . . . , ωp , g1,1 (t), . . . , gp,p (t))
n=1
n!
180
Theorem 2 generalizes the classical theorem by André [And1879] for alternating permutations:
∞
X dn tn 1 + sin(t)
=
n=0
n! cos(t)
where we agree that d0 = 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 RA 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] among others. Let Y = (Y1 , . . . YN ) be a sequence
of distinct reals in [0, 1]. We can construct from Y a permutation σY as follows. Let k1 ∈ {1, 2, . . . N } be the
integer such that Yk1 is minimal in {Y1 , . . . YN } and put σY (k1 ) = 1. Then, let k2 ∈ {1, 2, . . . N } − {k1 } be the
integer such that Yk2 is minimal in {Y1 , . . . YN } − {Yk1 }, put σY (k2 ) = 2 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
C1{D(Y )=A} (2)
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 (U1 , . . . UN ), uniform on [0, 1], and the functions fi defined in Theorem
1, construct a sequence Y = (Y1 , . . . YN ) as follows:
Y1 is the real in [0, 1] such that
Z Y1 Z 1
f1 (y)dy = U1 f1 (y) dy
0 0
For i ∈ [1, N − 1], Yi+1 is the only solution in [0, 1] of the equation
fi (Yi+1 ) = Ui+1 fi (Yi ) (3)
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
Z Yi+1
fi+1 (y) dy = Ui+1 fi (Yi ) (4)
0
Since the functions fi , fi+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 .
181
Using our algorithm, we need to compute the functions fi , 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 Yi , 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:
R1
Y1 has density f1 / 0 f1 (y) dy.
If i ∈ A, conditionally on Yi , Yi+1 has density 1[0,Yi ] fi+1 /fi (Yi ).
If i ∈
/ A, conditionally on Yi , Yi+1 has density 1[Yi ,1] fi+1 /fi (Yi ).
As a consequence, the density of the sequence Y = (Y1 , . . . , YN ) with respect to the Lebesgue measure on [0, 1]N
is equal to the telescopic product
f1 (Y1 ) f2 (Y2 ) fN (YN )
R1 × × ... × 1{D(Y )=A}
f (y) dy f1 (Y1 )
0 1
fN −1 (YN −1 )
1
= R1 1{D(Y )=A}
f (y) dy
0 1
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 Z Z1 1
dy1 . . . dyN 1{D(Y )=A}
0 0
But this probability is also equal to
QN (A)
N!
Using the formula for the density of Y , we find that
Z 1 Z 1
1
R1 dy1 . . . dyN 1{D(y)=A} = 1
f (y)dy 0
0 1 0
whence Z 1
QN (A) = N ! f1 (y) dy
0
which proves Theorem 1.
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
182
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 fi , f˜i as in Theorem 1.
Now define the polynomials hi , 1 ≤ i ≤ N by hN = 1 and hi (x) = fi (x) if i ∈ A while hi (x) = fi (1 − x) if
i∈/ A. Remark that for i ≤ N − 1, the hi are increasing functions on [0, 1] such that hi (0) = 0. Besides, one can
directly define the hi by descending induction as follows. First, hN = 1 and then, if i ∈
/ B,
Z x
hi (x) = hi+1 (y) dy (5)
0
while if i ∈ B, Z x
hi (x) = hi+1 (1 − y) dy (6)
0
Likewise, define the polynomials h̃i , 1 ≤ i ≤ N by h̃N = 1, h̃i (x) = f˜i (x) if i ∈ Ã and h̃i (x) = f˜i (1 − x) if
i∈/ Ã.
Since B̃ = B ∪ {m}, we have h̃i = hi for i > m. Moreover, using (5) and (6), we get that h̃m (x) > hm (x) for
every x ∈ (0, 1] and by induction, for every i < m and every x ∈ (0, 1], h̃i (x) > hi (x). Integrating this inequality
for i = 1, we get the result.
Remark 3.1 The coefficients of the polynomials fi can be computed by induction. Put
i
(i)
X
fN −i (x) = ak xk
k=0
(0) (i+1)
First, a0 = 1. Next, if N − i ∈ A, then a0 = 0 and for k ≥ 0,
(i)
(i+1) ak
ak+1 =
k+1
On the other hand, if N − i ∈
/ A, then for k ≥ 0,
(i)
(i+1) ak
ak+1 = −
k+1
and
i (i)
(i+1)
X a k
a0 =
k+1
k=0
Therefore, if N − i ∈
/ A,
i (i−j)
(i+1)
X a c
a0 = 0
(−1)|A ∩[N −i+1,N −1]|
j=0
(j + 1)!
By recursion, this leads to (1), which provides an alternative proof of Theorem 1.
4 Proof of Theorem 1.2 (outline)
Define by induction the sequence of polynomials (fn ) by
f0 = 1
If i ∈ A∞ , i ≥ 1 Z x
fi (x) = fi−1 (y)dy
0
If i ∈
/ A∞ , i ≥ 1 Z 1
fi (x) = fi−1 (y)dy
x
183
We claim that Z 1
dN = N ! fN −1 (x)dx
0
Indeed, Theorem 1 tells us that the right-hand side is the number of permutations with length N and descent
set N − DN . Using the symmetries of the problem, we see that this is equal to the number of permutations with
length N and descent set DN .
Consider the bivariate generating function
∞
X
F (x, t) = tn fn (x)
n=0
Then we have
∞ Z 1
X dn
tn = tF (x, t)dx
n=1
n! 0
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,
∂ p F (x, t)
= (−1)µ tp F (x, t)
∂xp
Solutions of this differential equation have the general form
p
X
F (x, t) = ak (t)eωk tx
k=1
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 ak are power
series: X
ak (t) = ak,n tn
n≥0
Therefore we can rewrite Z 1 p X ωk t
X e −1
tF (x, t)dx = ak,n tn
0 ωk
k=1 n≥0
We can also re-express
p
n X
X ω l xl k
fn (x) = ak,n−l
l!
l=0 k=1
and so
p
n X
0
X ω l+1 xl k
fn+1 (x) = ak,n−l
l!
l=0 k=1
It is thus possible to compute the the coefficients ak,n using the fact that
If n ≥ 1 is a descent, fn (0) = 0 whence
p
X
ak,n = 0
k=1
0
and moreover fn+1 = −fn , whence
p
X p
X
ak,n ωkj = ak,n ωkj+1
k=1 k=1
If n ≥ 1 is an ascent, fn (0) = 1 whence
p p
n X
X X ak,n−j ω j k
ak,n = −
j=1 k=1
j!
k=1
184
0
and moreover fn+1 = fn , whence
p
X p
X
ak,n ωkj = − ak,n ωkj+1
k=1 k=1
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
∞ Z 1 p X
p
X dn X eωk t − 1
tn = tF (x, t)dx = hk,l (t) (7)
n=1
n! 0 ωk
k=1 l=1
where, for all k, l ∈ [1, p], the function
∞
X
hk,l (t) = ak,np+l tnp+l
n=0
has a rational expression, with integer coefficients, in terms of the ωk and of the functions
∞
X (ωk t)np+l
gk,l (t) =
n=0
(np + l)!
This proves Theorem 2.
5 Implementation
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 Yi
in the algorithm.
As pointed out in the introduction, generating the variables Yi amounts to solving N equations of the form (4).
Of course, we can only get approximate solutions and since the Yi 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 Yi and its approximation
is smaller than the minimal value of 2|Yi − Yj | over all i, j ∈ [1, nmax].
Since equations of the form (4) are polynomial and since fi0 = ±fi−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.
100 200 300 400 500 600 700 800 900 1000
0.18 0.93 2.21 4.46 6.75 11.36 16.14 23.36 31.73 37.93
0.19 0.97 2.52 4.65 7.36 12.29 17.70 24.66 34.59 41.15
In the next simulations, we have taken p=4 and we only give the integer part of the run time (in seconds).
1200 1400 1600 1800 2000 2200 2400 2600 2800 3000
68 90 119 174 213 269 351 627 678 766
3300 3600 3900 4200 4500 4800 5100 5400 5700 6000
840 1210 1410 1822 2385 3283 5395 5877 6380 9121
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
1 The code is available at http://www.math.univ-paris13.fr/∼marchal/algopermut.m
185
descent profile would be sufficient, since the other polynomials can be obtained by differentiation. However, this
would mean computing these polynomials twice.
Acknowledgements
I am very grateful to Cyril Banderier for his help in the implementation of the algorithm. I also thank Olivier
Bodini and Michael Wallner for useful comments and references. This work is partially supported by the research
project “Combinatoire à Paris”.
References
[And1879] D. André. Développements de sec x et tan x. Comptes Rendus de l’Académie des Sciences, Paris,
88:965-967, 1879.
[BMW18] C. Banderier, P. Marchal, W. Michael. Rectangular Young tableaux with local decreases and the
density method for uniform random generation. Preprint, 2018.
[BHR03] E. A Bender, W. J. Helton, L.B. Richmond. Asymptotics of permutations with nearly periodic patterns
of rises and falls. Electronic Journal of Combinatorics, 10, 2003.
[Bas16] N. Basset. Counting and generating permutations in regular classes. Algorithmica, 76:989-1034, 2016.
[Bon12] M. Bóna. Combinatorics of permutations. Second edition. Series: Discrete Mathematics and its Appli-
cations (Boca Raton). CRC Press, Boca Raton, FL, 2012.
[BRS12] O. Bodini, O. Roussel, M. Soria. Boltzmann samplers for first-order differential specifications. Discrete
Applied Mathematics, 160 (18):2563-2572, 2012.
[BR70] N. G.de Bruijn. Permutations with given ups and downs. Nieuw Archi. Wisk., 3(18):61-65, 1970.
[Car75] L. Carlitz. Generating functions for a special class of permutations. Proceedings of the American
Mathematical Society, 47:251-256, 1975.
[ELR02] R. Ehrenborg, M. Levin, M. A. Readdy. A probabilistic approach to the descent statistic. Journal of
Combinatorial Theory Series A, 98(1):150-162, 2002.
[FS09] P. Flajolet, R. Sedgewick. Analytic combinatorics. Cambridge University Press, Cambridge, 2009.
[LM83] D. J. Leeming, A. MacLeod. Generalized Euler number sequences: asymptotic estimates and congru-
ences. Canadian Journal of Mathematics, 35:526-546, 1983.
[Luc14] J. M. Luck. On the frequencies of patterns of rises and falls. Physica A: Statistical Mechanics and its
Applications, 407:252-275, 2014.
[Mar16] P. Marchal. Rectangular Young tableaux and the Jacobi ensemble. Discrete Mathematics and Theo-
retical Computer Science Proceedings, BC:839-850, 2016.
[MRR10] A. Mendes, J. B. Remmel, A. Riehl. Permutations with k-regular descent patterns. Permutation
patterns. London Mathematical Society Lecture Note Series, 376:259-285, Cambridge University Press,
Cambridge, 2010.
[Niv68] I. Niven. A combinatorial problem of finite sequences. Nieuw Arch. Wiskunde, 3(16):116-123, 1968.
[Vie79] G. Viennot. Permutations ayant une forme donnée. Discrete Mathematics, 26(3):279-284, 1979.
186