Combinatoire 1

Les permutations

Pour tout-un-chacun, une permutation de $n$ objets (différents) est une disposition de ces objets dans un certain ordre. Par exemple, $\Box\Diamond\triangle$, $\heartsuit\spadesuit\clubsuit\diamondsuit$,...

Pour nous les objets en question seront des symboles (des lettres) et les permutations pourront être vues comme des mots. On utilisera en général comme alphabet l'ensemble des entiers positifs (en prenant garde à ne pas confondre le mot 123 (un, deux, trois) avec l'entier cent vingt trois.

Les permutations de $123$ sont donc $123,132,213,231,312,321$.

Une permutation peut aussi être vue comme une fonction $\sigma:\ \{1,2,\ldots,n\}\rightarrow \{1,2,\ldots,n\}$. Par exemple, $\sigma=24513$ est la fonction telle que $\sigma(1)=2,\sigma(2)=4$, etc., qu'on réprésente aussi par le tableau à deux lignes $\sigma={12345\choose 24513}$.

Ce sont les fonctions bijectives de cet ensemble dans lui même.

Dénombrement

Il est connu depuis la nuit des temps que le nombre de permutations de $n$ objets est égal à la factorielle $$n! = n\times(n-1)\times(n-2)\times\cdots\times 2\times 1.$$

En effet, on a $n$ choix pour le premier objet, $(n-1)$ pour le second, $(n-2)$ pour le troisième etc.

Ces nombres grandissent très vite :

In [1]:
def factorielle(n): return reduce(lambda x,y:x*y, xrange(1,n+1))
for n in range(1,25): print factorielle(n)
1
2
6
24
120
720
5040
40320
362880
3628800
39916800
479001600
6227020800
87178291200
1307674368000
20922789888000
355687428096000
6402373705728000
121645100408832000
2432902008176640000
51090942171709440000
1124000727777607680000
25852016738884976640000
620448401733239439360000

Avec les machines actuelles et un bon algorithme, le temps pour parcourir les permutations de 10 est d'environ une seconde. Pour parcourir celles de 20, il faudrait donc (en années)

In [2]:
reduce(lambda x,y:x*y, range(11,21))/(3600.*24*365)
Out[2]:
21259.594520547944

Génération

Le plus efficace pour engendrer les permutations en Python est d'utiliser itertools.permutations (codée en C) :

In [3]:
from time import time

from itertools import permutations
ll=range(10)
a=time();print len(list(permutations(ll))); print time()-a
3628800
1.11357498169

Cependant, tout n'est pas prévu dans itertools, et de nombreuses structures combinatoires peuvent être codées à la volée en utilisant une définition récursive simple. Il faut s'entraîner à le faire. Pour les permutations en ordre lexicographique, ça donnerait quelque chose du genre :

In [4]:
# On préférera retourner des tuples qui pourront servir de clefs dans des dictionnaires
def listperms(ll):
    if len(ll)<2: return [tuple(ll)]
    return [ (ll[i],)+x for i in range(len(ll)) for x in listperms(ll[:i]+ll[i+1:]) ]
In [5]:
print listperms('123')
print listperms((1,2,3))
print listperms([1,2,3])
[('1', '2', '3'), ('1', '3', '2'), ('2', '1', '3'), ('2', '3', '1'), ('3', '1', '2'), ('3', '2', '1')]
[(1, 2, 3), (1, 3, 2), (2, 1, 3), (2, 3, 1), (3, 1, 2), (3, 2, 1)]
[(1, 2, 3), (1, 3, 2), (2, 1, 3), (2, 3, 1), (3, 1, 2), (3, 2, 1)]

Évidemment, ce n'est pas très efficace. Mais ça permet d'avoir rapidement des exemples corrects pour de petites valeurs de $n$, lesquels pourront permettre de tester la validité des résultats d'algorithmes plus sophistiqués.

In [6]:
a = time();print len(listperms(ll));print time()-a
3628800
17.2320098877

Permutations et algorithmes de tri

Pour un informaticien, une permutation est une liste en désordre, et il va donc s'employer à la trier.

La première méthode qui vient à l'esprit pour trier avec un ordinateur rudimentaire est le tri par bulles (bubble sort en anglais).

Pour trier une permutation $\sigma = 43127865$, on la parcourt de gauche à droite avec un pointeur $i$. Si on recontre une descente, c'est à dire une valeur de $i$ pour laquelle $\sigma[i]>\sigma[i+1]$, on échange $\sigma[i]$ et $\sigma[i+1]$, et on décrémente $i$ (pour tester si on a fait apparaître une nouvelle descente), sinon on incrémente $i$. Dans l'exemple qui suit, la colonne de droite indique les descentes rencontrées.

\begin{align} 43127865 &\quad 1\\ 34127865 &\quad 2\\ 31427865 &\quad 1\\ 13427865 &\quad 3\\ 13247865 &\quad 2\\ 12347865 &\quad 6\\ 12347685 &\quad 5\\ 12346785 &\quad 7\\ 12346758 &\quad 6\\ 12346578 &\quad 5\\ 12345678 \end{align}

L'élimination d'une descente peut en créer une autre, par exemple $312\rightarrow 132$. Ce que l'algorithme élimine à chaque échange, c'est une inversion.

Une inversion d'une permutation $\sigma$ est un couple $(i,j)$ avec $i<j$ et $\sigma(i)>\sigma(j)$. Par exemple les inversions de $\sigma=43127865$ sont $\{(1,2),(1,3),(1,4),(2,3),(2,4),(5,7),(5,8),(6,7),(6,8),(7,8)\}$.

Les descentes sont donc des inversions spéciales $(i,i+1)$. Chaque fois qu'on échange $\sigma(i)$ et $\sigma(i+1)$, on élimine l'inversion $(i,i+1)$. Dans les autres inversions, les rôles de $i$ et $i+1$ sont échangés, ce qui ne change pas leur nombre. On élimine donc exactement une inversion à chaque échange, ce qui prouve que l'algorithme se termine.

On peut ainsi prédire le nombre d'échanges que le tri effectuera en calculant le nombre d'inversions de la permutation.

Plutôt que d'en faire la liste, il est commode d'écrire en dessous de chaque $\sigma(i)$ le nombre d'inversions dont il est responsable, c'est à dire le nombre de valeurs plus petites sur sa droite : \begin{align} \sigma = &43127865\\ c= &32002210 \end{align} Le nombre d'inversions ${\rm inv}(\sigma)$ est la somme de ces valeurs, ici $3+2+0+0+2+2+1+0=10$.

Le tuple $c(\sigma)=(3,2,0,0,2,2,1,0)$ est appelé code de Lehmer de la permutation. Il permet de la reconstituer sans ambigüité. En effet, $c_1=3$ nous dit que $\sigma(1)$ est la 4ème lettre de $12345678$, donc 4. Ensuite, $c_2=2$ dit que $\sigma_2$ est la 3ème lettre de $1235678$,(on a éliminé 4) donc 3, et ainsi de suite.

Par ailleurs, il est clair que les composantes du code doivent vérifier $0\le c_i\le n-i$. Le code établit donc une bijection entre les permutations et les $n!$ tuples vérifiant ces conditions.

In [7]:
def lehmer_code(p):
    return [len([x for x in p[i+1:] if x<p[i]]) for i in range(len(p))] 

lehmer_code('43127865')
Out[7]:
[3, 2, 0, 0, 2, 2, 1, 0]
In [8]:
def perm_from_code(c):
        res = []; n = len(c); ll = range(1,n+1)
        while ll:
            res.append(ll[c[0]]);ll.remove(ll[c[0]]); c=c[1:]
        return tuple(res)
In [9]:
perm_from_code((3,2,0,0,2,2,1,0))
Out[9]:
(4, 3, 1, 2, 7, 8, 6, 5)

La structure de groupe

Ler permutations, vues comme bijections d'un ensemble dans lui-même, peuvent être composées : $$\tau\circ\sigma(i) =\tau(\sigma(i)).$$

Cette opération est associative, et possède comme élément neutre la permutation identité $$id_n=123\cdots n.$$

In [10]:
def compose(q,p): return tuple([q[p[i]-1] for i in range(len(p))])
In [11]:
p =(2,4,5,1,3)
q =(3,5,2,1,4)
compose(q,p)
Out[11]:
(5, 1, 4, 3, 2)

De plus, toute permutation est inversible. Pour calculer l'inverse de $\sigma={12345\choose24513}$, il suffit d'échanger les deux lignes du tableau : $$\sigma^{-1}={24513\choose 12345}={12345\choose41523}$$ (on peut voir ces tableaux comme des suites de colonnes, qui représentent la même permutation quelque soit leur ordre.)

In [12]:
def inverse(p): return tuple([p.index(j)+1 for j in range(1,len(p)+1)])
In [13]:
q = inverse(p); q
Out[13]:
(4, 1, 5, 2, 3)
In [14]:
compose(q,p)
Out[14]:
(1, 2, 3, 4, 5)

Les permutations d'un ensemble $A$ forment donc un groupe, appelé groupe symétrique de $A$, et noté ${\mathfrak S}(A)$. Pour le choix standard $A=\{1,2,\ldots,n\}$, on le note ${\mathfrak S}_n$.

Action à droite sur les mots

En plus de les composer , on peut faire agir les permutations de ${\mathfrak S}_n$ sur les mots de longueur $n$ de n'importe quel alphabet $A$. Cette action s'écrit à droite : si $w=a_1a_2\cdots a_n\in A^n$, $$w\cdot\sigma := a_{\sigma(1)}a_{\sigma(2)}\cdots a_{\sigma(n)}.$$

Par exemple, $auxb\cdot 3241 = xuba$.

L'intérêt de cette notation vient du fait que si $w$ est lui-même le mot représentant une permutation $\tau$, alors $w\cdot\sigma =\tau_{\sigma(1)}\tau_{\sigma(2)}\cdots \tau_{\sigma(n)}=\tau\circ\sigma$.

On aura ainsi $(w\cdot\tau)\cdot \sigma = w\cdot (\tau\circ\sigma)$, ce qui est moralement satisfaisant.

On pourra omettre les symboles $\circ$ ou $\cdot$ et écrire simplement $w\sigma$ ou $\tau\sigma$ dans les contextes ou on ne risquera pas de confondre avec la concaténation des mots.

On peut donc voir $\tau\sigma$ comme représentant l'action de $\tau$ sur les valeurs de $\sigma$, ou l'action de $\sigma$ sur les positions de $\tau$.

Transpositions

Une transposition est une permutation qui échange deux valeurs $i<j$ : $\tau_{i,j}$, notée aussi $(i,j)$ est définie par $\tau(i)=j$, $\tau(j)=i$ et $\tau(k)=k$ pour $k\not=i,j$.

Il y a donc ${n\choose 2}$ transpositions dans ${\mathfrak S}_n$.

Les transpositions élémentaires $s_i=(i,i+1)$ sont celles qui échangent deux valeurs consécutives.

Tri par bulles, et structure du groupe symétrique

Et que fait le tri par bulles ? Il trie la permutation en lui appliquant à droite une suite de transpositions élémentaires ! Autrement dit, il exprime la permutation inverse $\sigma^{-1}$ comme un produit de transpositions élémentaires : $$43127865\cdot s_1s_2s_1s_3s_2s_6s_5s_7s_6s_5 = 12345678$$ ce qui nous donne en renversant l'odre des $s_i$, une expression de $\sigma$ comme produit de transpositions élémentaires : $$\sigma = s_5s_6s_7s_5s_6s_2s_3s_1s_2s_1.$$

La preuve de terminaison du tri par bulles se traduit donc par le

Théorème

Toute permutation peut se décomposer en produit de transpositions élémentaires. Le nombre minimal de transpositions à utiliser est égal au nombre d'inversions de la permutation.

Une decomposition de longueur minimale est dite réduite. C'est donc ce que calcule le tri par bulles, qui pour cet usage, est en fait optimal (le tri est un effet de bord).

La décomposition réduite n'est en général pas unique. Par exemple, $321=s_1s_2s_1 = s_2s_1s_2$, et $2143=s_1s_3=s_3s_1$.

On peut prouver que toutes les relations satisfaites par les $s_i$ sont des conséquences de

  • $s_i^2=id$,
  • $s_is_j=s_js_i$ si $|i-j|\ge 2$,
  • $s_is_{i+1}s_i = s_{i+1}s_i s_{i+1}$.

Les deux dernières relations sont appelées relations de tresses : elles sont en effet vérifiées si on interprète $s_i$ comme l'opération consistant à croiser le $i$-ème brin d'une tresse par dessus le $i+1$ème.

In [15]:
def bubble_sort(ll):
    mm = list(ll)
    n = len(mm);  rd = []; done = False
    while not done:
        done = True
        for i in range(n-1):
            if mm[i]>mm[i+1]:
                mm[i],mm[i+1] = mm[i+1],mm[i]
                rd.append(i+1)
                done = False
    return mm, list(reversed(rd))

bubble_sort('43127865')
Out[15]:
(['1', '2', '3', '4', '5', '6', '7', '8'], [5, 6, 5, 2, 1, 7, 6, 3, 2, 1])
In [16]:
 def reduced_dec(p):
        return bubble_sort(p)[1]
reduced_dec('43127865')
Out[16]:
[5, 6, 5, 2, 1, 7, 6, 3, 2, 1]
In [17]:
def perm_from_reduced_dec(rd,n):
        ll = range(1,n+1)
        for i in rd: ll[i-1],ll[i] = ll[i],ll[i-1]
        return tuple(ll)
perm_from_reduced_dec([5, 6, 5, 2, 1, 7, 6, 3, 2, 1],8)
Out[17]:
(4, 3, 1, 2, 7, 8, 6, 5)

Statistiques sur les permutations

Le nombre moyen d'échanges effectués par le tri par bulles sur une liste aléatoire de longueur $n$ est donc le nombre moyen d'inversions d'une permutations de ${\mathfrak S}_n$. Pour le calculer, il faudrait donc savoir compter les permutations par nombre d'inversions.

Le code de Lehmer peut fair ça pour nous. Si $c$ est le code de $\sigma$, le monôme $x_1^{c_1}x_2^{c_2}\cdots x_n^{c_n}$ deviendra $q^{{\rm inv}(\sigma)}$ si l'on pose $x_1=x_2=\ldots = x_n=q$. Or on peut donner un expression simple de la somme de tous ces monômes. On sait que $c_k$ prend toutes les valeurs possibles entre 0 et $k-1$. Donc, $$\sum_{\sigma\in{\mathfrak S}_n}x_1^{c_1(\sigma)}x_2^{c_2(\sigma)}\cdots x_n^{c_n(\sigma)} =(1+x_1+x_1^2+\cdots+x_1^{n-1})(1+x_2+\cdots+x_2^{n-2})\cdots(1+x_{n-1})(1)$$.

Par exemple, les codes des permutations de ${\mathfrak S}_3$ sont $000,100,0100,200,110,210$ et $$x_1^{0}x_2^{0}x_3^{0} + x_1^{1}x_2^{0}x_3^{0} + x_1^{0}x_2^{1}x_3^{0} + x_1^{2}x_2^{0}x_3^{0} + x_1^{1}x_2^{1}x_3^{0} + x_1^{2}x_2^{1}x_3^{0} = (1+x_1+x_1^2)(1+x_2)$$

Si l'on pose tous les $x_i=q$ dans cette expression, on obtient un polynôme appelé $q$-factorielle : $$[n]_q! := (1+q+q^2+\cdots+q^{n-1})(1+q+\cdots+q^{n-2})\dots(1+q)1$$ qui devient évidemment $n!$ pour $q=1$.

C'est le produit des $q$-entiers $$[k]_q:=1+q+\cdots+q^{k-1} = \frac{1-q^k}{1-q}$$ pour $k$ de 1 à $n$.

In [18]:
from sympy import *
In [19]:
var('q')
Out[19]:
q
In [20]:
def qfact(n): 
    return expand(reduce(lambda x,y: x*y, [ sum([q**i for i in range(k)]) for k in range(1,n+1)]))
In [21]:
qfact(3)
Out[21]:
q**3 + 2*q**2 + 2*q + 1
In [22]:
qfact(4)
Out[22]:
q**6 + 3*q**5 + 5*q**4 + 6*q**3 + 5*q**2 + 3*q + 1

Au vu de ce resultat, le nombre moyen d'inversions d'une permutation de ${\mathfrak S}_4$ s'obtiendra en calculant $$1\times 6 + 3\times 5 + 5\times 4 + 6\times 3 + 5\times 2 + 3\times 1 + 1\times 0$$ et en divisant le résultat par $4!$.

Cela revient à évaluer la dérivée du polynôme $P_n(q)=[n]_q!$ en $q=1$ et à diviser par $n!=P_n(1)$. C'est donc $$\frac{P'_n(1)}{P_n(1)} = \left.\frac{d}{dq}\right|_{q=1}\ln P_n(q).$$

Comme $P_n(q)$ a le bon goût d'être un produit d'expressions simples, son logarithme est la somme des logarithmes desdites expressions :

$$\ln P_n(q) = \sum_{k=1}^n\left[\ln(1-q^k)-\ln(1-q)\right]$$ de sorte que sa dérivée est $$\frac{d}{dq}\ln P_n(q) = \sum_{k=1}^n\left[\frac{-kq^{k-1}}{1-q^k}+\frac1{1-q}\right]$$ On ne peut pas poser directement $q=1$ dans cette expression, mais on peut trouver sa valeur en posant $q=1-t$ et en faisant un développement limité en $t$: $1-q=t$, $1-q^k=1-(1-t)^k=1-(1-kt+{k\choose 2}t^2+O(t^3))=kt(1-\frac12(k-1)t+O(t^2))$, et $-kq^{k-1}=-k(1-(k-1)t+O(t^2))$, ce qui nous donne en réduisant au même dénominateur $$ \frac{-kq^{k-1}}{1-q^k}+\frac1{1-q}=\frac{k-1}{2}+O(t)$$ et donc finalement $$\frac{P'_n(1)}{P_n(1)}= \sum_{k=1}^n\frac{k-1}{2}= \frac{n(n-1)}{4}.$$

La complexité du tri par bulles est donc $O(n^2)$. Il fallait tout ça pour en arriver là.

Faire une statistique sur les permutations, c'est compter le nombre de permutations possédant une propriété donnée. Cette propriété ${\rm stat}{\sigma}$ est généralement un nombre (par exemple le nombre d'inversions), ou un ensemble de nombres (les descentes). La fonction $f_n(x)=|\{\sigma\in{\mathfrak S}_n|{\rm stat}(\sigma)=x\}|$ est appelée distribution de la statistique. On la représente souvent par sa fonction génératrice, un polynôme en une ou plusieurs variables $X$ $$ F_n(X) = \sum_{\sigma\in{\mathfrak S}_n}X^{{\rm stat}(\sigma)}.$$

Pour le nomre d'inversions ${\rm inv}(\sigma)$, $F_n(q)=[n]_q!$.

In [23]:
import matplotlib.pyplot as plt
%matplotlib inline
In [24]:
P10 = qfact(10)
F10 = P10.as_poly().all_coeffs()
plt.plot(F10)
Out[24]:
[<matplotlib.lines.Line2D at 0x21325490>]

Une autre statistique naturelle est le nombre de descentes ${\rm des}(\sigma)$. On commence par calculer la liste des descentes :

In [25]:
def descents(p):
    return [i+1 for i in range(len(p)-1) if p[i]>p[i+1]]
In [26]:
descents((4,6,2,5,1,7,8,3))
Out[26]:
[2, 4, 7]
In [27]:
def des_num(p): return len(descents(p))

Étudions maintenant sa distribution :

In [28]:
def des_distr(n):
    f = {}
    for p in permutations(range(1,n+1)):
        d = des_num(p)
        if d in f: f[d]+=1
        else: f[d]=1
    return [x[1] for x in sorted(f.items())]
In [29]:
for n in range(10): print des_distr(n)
[1]
[1]
[1, 1]
[1, 4, 1]
[1, 11, 11, 1]
[1, 26, 66, 26, 1]
[1, 57, 302, 302, 57, 1]
[1, 120, 1191, 2416, 1191, 120, 1]
[1, 247, 4293, 15619, 15619, 4293, 247, 1]
[1, 502, 14608, 88234, 156190, 88234, 14608, 502, 1]

Il n'est pas évident de trouver une formule pour ces nombres. Bien que l'étude des descentes des permutations ne remonte qu'au début du XXème siècle, il se trouve qu'Euler avait recontré ces mêmes nombres dès 1755 en étudiant un problème qui n'avait rien à voir avec les permutations : le calcul de la somme des séries $${\mathcal A}_n(t)=\sum_{k\ge 0} k^n t^k.$$ Il obtient la formule $${\mathcal A}_n(t)=\frac{A_n(t)}{(1-t)^{n+1}}$$ où $$A_n(t)=\sum_{k=0}^{n-1}A(n,k)t^{k+1}$$ et les nombres $A(n,k)$, calculés au moyen d'une récurrence trouvée par Euler, comptent apparemment le nombre de permutations de ${\mathfrak S}_n$ ayant exactement $k$ descentes.

Jusqu'à une période récente, on ne savait pas prouver cette coincidence autrement que par récurrence. Et surtout, pour trouver la série génératrice des polynômes $A_n(t)$ à partir de la récurrence, il fallait résoudre une équation différentielle compliquée, ce qu'on ne faisait qu'en parachutant la solution, miraculeusement connue d'avance grâce au travail d'Euler.

On verra dans la suite une méthode moderne permettant de la trouver sans tricher.

In [30]:
var('t')
def A(n):
    aa = des_distr(n)
    return sum([aa[i]*t**(i+1) for i in range(len(aa))])
 

A(4)
Out[30]:
t**4 + 11*t**3 + 11*t**2 + t

Vérifions le résultat d'Euler :

In [31]:
for n in range(1,6): print series(A(n)/(1-t)**(n+1),t,0,8)
t + 2*t**2 + 3*t**3 + 4*t**4 + 5*t**5 + 6*t**6 + 7*t**7 + O(t**8)
t + 4*t**2 + 9*t**3 + 16*t**4 + 25*t**5 + 36*t**6 + 49*t**7 + O(t**8)
t + 8*t**2 + 27*t**3 + 64*t**4 + 125*t**5 + 216*t**6 + 343*t**7 + O(t**8)
t + 16*t**2 + 81*t**3 + 256*t**4 + 625*t**5 + 1296*t**6 + 2401*t**7 + O(t**8)
t + 32*t**2 + 243*t**3 + 1024*t**4 + 3125*t**5 + 7776*t**6 + 16807*t**7 + O(t**8)

L'indice du major

Puisque les descentes d'une permutation de ${\mathfrak S}_n$ sont des entiers compris entre 1 et $n-1$, leur somme $${\rm maj}(\sigma) = \sum_{d\in {\rm Des}(\sigma)} d$$ est comprise entre 0 et $1+2+\cdots+n-1=\frac12 n(n-1)$, donc dans le même intervalle que le nombre d'inversions.

Une découverte surprenante du Major Percy Alexander MacMahon est que cette statistique, qu'il appelait greatest index d'une permutation, devenue indice du Major en français, puis major index en anglais histoire de faire un jeu de mots, et enfin indice majeur, possède la même distribution que le nombre d'inversions :

In [32]:
def maj(p): return sum(descents(p))

def gen_maj(n): 
    g = {}
    for p in permutations(range(1,n+1)):
        m = maj(p)
        if m in g: g[m]+=1
        else: g[m]=1
    return sum([g[m]*q**m for m in g])
In [33]:
gen_maj(4)
Out[33]:
q**6 + 3*q**5 + 5*q**4 + 6*q**3 + 5*q**2 + 3*q + 1
In [34]:
qfact(4)
Out[34]:
q**6 + 3*q**5 + 5*q**4 + 6*q**3 + 5*q**2 + 3*q + 1
In [34]: