Chapitre 5
Diviser pour régner

5.1 Quelques idées sur la complexité

5.1.1 Notion de taille des données

Rappelons que les différentes mémoires sont des éléments essentiels des ordinateurs. On distingue en général deux sortes de mémoires

Toutes les données (nombres, caractères, noms, adresses, etc.) manipulées doivent être stockées dans l’ordinateur. A cet effet, toutes ces mémoires sont divisées en petites cases appelées octets ; ces octets peuvent stocker une suite de 8 chiffres égaux à 0 ou à 1 (des bits), soit 256 éléments différents. Ces octets sont numérotés dans chaque mémoire de 0 jusqu’à la capacité de la mémoire (moins un). L’ordinateur se charge d’effectuer une traduction entre les données entrées et un ensemble d’octets, c’est l’opération de codage. Typiquement

A partir de là, on voit que toute donnée (la liste des élèves d’une classe et leurs notes au cours de l’année, l’annuaire du téléphone) a une taille que l’on peut exprimer en octets. Cette taille dépend de la méthode de codage utilisée (par exemple MS-DOS code les caractères sur un octet alors que Windows les code sur deux octets) mais obéit à quelques règles simples

Nous allons dans la suite de cette section étudier le rapport entre la taille des données et le temps d’exécution d’un traitement concernant ces données. Il nous suffit donc d’étudier le rapport entre le nombre N d’objets élémentaires à traiter et le temps d’exécution.

5.1.2 Types de complexité

Rappelons que si l’on dispose de deux fonctions f et g définies sur N et à valeurs réelles positives, on dit que f(n) = O(g(n)) s’il existe une constante K 0 telle que n ∈ N, f(n) Kg(n).

Considérons une donnée constituée d’un certain nombre d’ objets et un certain traitement concernant ces données. Soit T(N) le temps d’exécution de ce traitement lorsqu’il concerne N objets. Ce temps d’exécution dépend évidemment de la machine utilisée et du codage des données, mais il possède certains invariants : en particulier, si l’on utilise deux machines de performances différentes, les temps T1(N) et T2(N) sont en gros proportionnels (le facteur de proportionnalité mesurant justement la différence de rapidité des deux machines).

Ceci va nous permettre de différencier certains types d’algorithmes, du point de vue de leur complexité, c’est à dire de la rapidité de croissance de T(N) quand N devient grand. On dit qu’un algorithme de traitement de données est à complexité

Ce temps d’exécution dépend bien entendu du nombre d’opérations élémentaires effectuées sur ces données et de la plus ou moins grande rapidité de ces opérations élémentaires. Parmi ces opérations élémentaires on peut citer

Dans l’évaluation de ce temps d’exécution, on sera amené à préciser quelles sont les opérations élémentaires à prendre en considération : si l’on effectue à peu près autant d’additions que de multiplications, on peut négliger les additions qui sont beaucoup plus rapides et ne compter que les multiplications, les échanges en mémoire interne peuvent souvent être négligés vis à vis d’échanges entre la mémoire interne et la mémoire externe, etc.

5.1.3 Comparaison de temps d’exécution

Prenons une machine (très rapide) capable d’accomplir 108 opérations par seconde, et imaginons des algorithmes qui effectuent un traitement donné dans un temps T(N). Nous donnons ci dessous un tableau des temps de traitement que l’on peut espérer pour N éléments suivant la complexité de l’algorithme (où s,m,h,j,a désignent respectivement la seconde, la minute, l’heure, le jour, l’année et u l’âge estimé de l’univers, et où les nombres totalement inimaginables n’ont pas été affichés)










log 2N N N log N N2 N4 2N 4N N!









N = 52.10-8s5.10-8s8.10-8s3.10-7s6.10-6s3.10-7s1.10-5s1.10-6s
N = 102.10-8s1.10-7s2.10-7s1.10-6s1.10-4s1.10-5s1.10-2s4.10-2s
N = 153.10-8s2.10-7s4.10-7s2.10-6s5.10-4s3.10-4s 1.101s 2h
N = 203.10-8s2.10-7s6.10-7s4.10-6s2.10-3s1.10-2s 2h 6.102a
N = 303.10-8s3.10-7s1.10-6s9.10-6s8.10-3s 1.101s 3.102a 6.106u
N = 504.10-8s5.10-7s2.10-6s3.10-5s6.10-2s 1.102j 2.104u 6.1038u
N = 1005.10-8s1.10-6s5.10-6s1.10-4s 1s 2.104u 4.1034u
N = 5006.10-8s5.10-6s3.10-5s3.10-3s1.101m
N = 10007.10-8s1.10-5s7.10-5s1.10-2s 2.100h
N = 50009.10-8s5.10-5s4.10-4s3.10-1s 7.101j
N = 100009.10-8s1.10-4s9.10-4s 1s 3a
N = 500001.10-7s5.10-4s5.10-3s 3.101s 2.103a









On constate que les algorithmes hyperexponentiels deviennent impraticables dès que N dépasse une dizaine, que les algorithmes exponentiels dépassent les limites acceptables dès que N atteint la trentaine. Ces deux types d’algorithmes sont donc incapables de traiter raisonnablement même de petits volumes de données.

On constate également, que pour rechercher un nom parmi 50000 personnes, entre le temps d’exécution d’une recherche linéaire en O(N) et celui d’une recherche dichotomique en O(log 2N), il existe un facteur 1000. De même, pour trier 10000 nombres réels (ce qui est courant en synthèse d’images avec une grille 100 × 100), on gagne un facteur 1000 en utilisant un tri en O(N log 2N) comme le tri fusion, par rapport à un tri en O(N2) comme le tri par insertion : sur la même machine, l’un des tris pourra prendre 4 secondes alors que l’autre tri demandera une heure !

Remarque  Il convient de nuancer un peu la dernière argumentation. Les algorithmes performants sont souvent complexes à écrire et cette complexité peut se ressentir sur les temps d’exécution pour de petits volumes de données. Autrement dit, pour un volume de données de l’ordre de quelques unités, il arrive souvent que les algorithmes bêtes soient plus rapides que les algorithmes intelligents. Entre un algorithme bête en N22 et un algorithme intelligent en 3N log 2N, la différence ne se fera sentir en faveur du second qu’à partir de N = 30. Il peut parfois être nécessaire d’adapter l’algorithme utilisé au volume des données à traiter. N’oublions pas également que les algorithmes sophistiqués utilisent souvent plus de mémoire, ce qui peut être un inconvénient même si celle-ci n’est plus une denrée rare ; il faut savoir établir un compromis entre rapidité et occupation mémoire.

5.2 Principes généraux de ”diviser pour régner”

5.2.1 Aperçus philosophiques

La méthode de programmation généralement désignée par diviser pour régner procède d’une philosophie dont nous allons tenter de donner une idée. Supposons que nous ayions à traiter un problème algorithmique concernant un objet de taille N (cet objet pouvant être lui même un ensemble d’objets). On peut procéder en trois étapes

Cette méthode est bien connue des adeptes de Machiavel : pour détruire un ensemble d’ennemis, les diviser en sous ensembles que l’on écrase un à un, la fusion des résultats ne posant alors guère de difficultés.

Beaucoup d’algorithmes performants se fondent sur ces principes, très souvent avec p = 2. C’est le seul cas que nous examinerons par la suite, mais le lecteur n’aura pas de mal à adapter raisonnements et méthodes à des cas plus complexes. Dans ce cas la méthode devient donc

5.2.2 Evaluations de temps de calcul

Remarque  En classe de Mathématiques Supérieures, on pourra sauter les démonstrations de cette section et se contenter de retenir les énoncés des deux théorèmes ci-dessous. Ces résultats seront repris dans un cadre un peu plus général en classe de Mathématiques Spéciales.

Appelons T(N) le temps nécessaire pour traiter un objet de taille N, P(N) le temps nécessaire pour partitionner l’objet de taille N en 2 sous-objets de taille approximative N∕2 d’une manière adéquate et F(N) le temps nécessaire pour fusionner deux solutions de taille N∕2 en une solution de taille N. On voit que la fonction T(N) va vérifier

T(N ) = P(N )+ 2T(N-)+ F (N )
                  2

Supposons en particulier que N est de la forme 2k et posons t(k) = T(2k). On aura alors

         k               k
t(k) = P (2) + 2t(k- 1)+ F (2)

Nous nous intéresserons tout particulièrement au cas où P et F sont des fonctions polynomiales de N. Quitte à les majorer, on peut alors supposer que P(N) = αNp et F(N) = βNp, auquel cas on obtient une suite t(k) vérifiant

t(k) = 2t(k- 1)+ γ2kp

avec γ = α + β. Posons xk = 2-kt(k). La suite xk vérifie la relation de récurrence

xk = xk-1 + γ2k(p-1)

Si p = 1 (c’est à dire si le temps nécessaire à la partition et à la fusion est proportionnel à la taille N des données ce qui est un cas courant), on obtient xk = x1 + γ(k - 1) ~ γk, donc t(k) = 2kxk ~ γk2k. Autrement dit T(N) ~ γN log 2N lorsque N est de la forme 2k. Si N n’est pas une puissance de 2, choisissons k de telle sorte que 2k-1 < N < 2k. On a alors T(2k-1) T(N) T(2k) où k est égal à la partie entière de log 2N équivalente à log 2N et 2k = 2 2k-1 < 2N, ce qui nous donnera une majoration du type T(N) δN log 2N.

Si p > 1, on a

          ∑k  i(p-1)        2(k+1)(p-1) --22(p-1)   k(p- 1)
xk = x1 +γ   2     = x1 + γ     2p-1 - 1     ~ δ2
          i=2

avec δ = γ2p- 1∕2p-1 - 1. On a alors t(k) = 2kxk ~ δ2kp. On obtient donc T(N) ~ δNp lorsque N est de la forme 2k. Lorsque N n’est pas de la forme 2k, on choisit k de telle sorte que 2k-1 < N < 2k. On a alors

T(N) ≤ T(2k) ≤ δ′2kp ≤ δ′′Np

puisque 2k = 22k-1 < 2N.

Théorème 5.2.1 On suppose que les temps de partition et de fusion sont des O(Np) pour un problème de taille N. Alors le temps de calcul dans une méthode ”diviser pour régner” est

En fait les méthodes diviser pour régner sont parfois un peu plus complexes et ne permettent pas tout à fait de traiter les deux sous problèmes de taille N∕2 de manière indépendante. Comme nous le verrons par la suite, il peut être nécessaire en fait de traiter q sous problèmes de taille N∕2 ce qui conduit, pour le temps de calcul T0(N) à des relations de récurrence du type T0(N) = qT0(N∕2) + P(N) + F(N). Supposons de nouveau que P(N) + F(N) γNp. Alors T0(N) est clairement majoré par la solution de la récurrence T(N) = γT(N∕2) + γNp. Posons alors t(k) = T(2k). Nous avons donc

t(k) = qt(k- 1)+ γ2kp

Soit ω = log 2q et xk = 2-ωktk. On a alors

     -ωk ( ω ω(k-1)        kp)
xk = 2    2 2     xk-1 + γ2

soit encore

             k(p- ω)
xk = xk-1 + γ2

On en déduit que xk = x0 + γ i=1k2i(p-ω).

Si p > ω, soit 2p > q, ceci nous donne

          2(k+1)(p-ω) - 2p-ω
xk = x0 +γ-----p-ω--------= O (2k(p- ω))
              2   - 1

soit encore

T(N ) = T (2k) = t(k) = 2ωkO (2k(p-ω)) = O(2kp) = O (N p)

Par contre, si p < ω, soit 2p < q, la suite xk admet une limite. On a donc xk = O(1), ce qui montre que

T (N) = T(2k) = t(k) = 2ωkO(1) = O (2ωk) = O (Nω )

Enfin, si p = ω, soit 2p = q, on a xk = x0 + γ(k - 1) = O(k). On en déduit que

T(N ) = T (2k) = t(k) = 2ωkO (k) = O(k2ωk) = O (N ω log2N )

Les mêmes méthodes par encadrement que dans le cas q = 2 montre que ces estimations restent valables lorsque N n’est pas de la forme 2k. On a donc

Théorème 5.2.2 On suppose que dans une méthode ”diviser pour régner”, la méthode conduise à résoudre q problèmes de taille N∕2 avec des processus de partition et de fusion en O(Np), si bien que le temps de calcul T(N) vérifie une relation

T (N ) = qT (N ∕2)+ O(N p)

Alors, en posant ω = log 2q

5.2.3 Calcul des puissances d’un entier

Nous allons évaluer le temps de calcul T(N) de la puissance xN d’un entier x en fonction de l’exposant N. Nous avons déjà vu deux méthodes pour calculer cette puissance

Dans les deux méthodes, le temps de calcul est proportionnel à N. Nous allons donc utiliser une méthode diviser pour régner pour diminuer le temps de calcul. Pour cela nous remarquons que xN = (x2)p si N = 2p est pair et xN = x × (x2)p si N = 2p + 1 est impair. Nous obtenons ainsi T(N) = T(N∕2) + α si N est pair et T(N) = T(N - 1∕2) + 2α si N est impair, où α est le temps nécessaire à une multiplication. Dans tous les cas, on a donc T(N) = T(N∕2) + O(1), ce qui d’après le théorème précédent nous donne T(N) = O(log 2N).

La fonction puissance ainsi définie est naturellement récursive et on obtient une fonction Caml


#let rec puiss2 x n =  
       if n = 0 then 1  
       else if n mod 2 = 0  
           then puiss2 (x*x) (n/2)  
           else x*puiss2 (x*x) (n/2);;  
puiss2 : int -> int -> int = <fun>  
#puiss2 2 8;;  
- : int = 256


Programme 5.1: calcul récursif rapide des puissances d’un nombre entier

Une version non récursive est un peu plus difficile à construire. Définissons une suite xn par x0 = x et xn+1 = (xn)2. On a facilement par récurrence n 0, xn = x2n . Décomposons alors n en base 2 sous la forme n = m020 + m121 + ⋅⋅⋅ + mk2k avec mi ∈{0,1}. On a alors xn = x0m0x1m1xkmk. Or on sait calculer les mi : il suffit de définir les ni et mi par récurrence par

n0 = n, mi = ni mod2,  ni+1 = ni div2

où l’on désigne par p div q le quotient de la division entière de p par q. On stocke alors dans trois références d’une part ni, d’autre part xi et enfin x0m0x1m1ximi. Ceci conduit à la fonction suivante :


#let puiss2it x n =  
     let y = ref x and p = ref n and res = ref 1 in  
       while !p<>0 do  
           if (!p mod 2) = 1 then res := !res * !y;  
           y := !y * !y; p := (!p/2)  
      done;  
      !res;;  
puiss2it : int -> int -> int = <fun>  
#puiss2it 2 27;;  
- : int = 134217728


Programme 5.2: calcul itératif rapide des puissances d’un nombre entier

Proposition 5.2.1 La fonction puiss2it calcule xn en un temps O(log 2N)

Démonstration: L’assertion sur le temps est évidente puisque p étant à itération divisé par 2 et étant initialisé à n au départ, le corps de la boucle ne peut être effectué qu’au plus log 2N fois. Ceci montre en même temps que la boucle s’achève bien.

Il nous suffit donc d’exhiber un invariant de la boucle. Or on constate que l’on a l’invariant suivant à la fin de la i-ième itération

y = xi, p = ni et res = x0m0x1m1xi-1mi-1

C’est vrai après la 0-ième itération, c’est à dire à l’initialisation de la boucle puisque y = x = x0, p = n = n0 et res = 1. De plus, si cette condition est vérifiée après la i-ième itération, on a alors p mod 2 = ni mod2 = mi et en distinguant suivant que cet entier vaut 1 ou 0, on obtient après l’affectation de res,

res = xm00xm11 ...xmi-i1-1ymi = xm00xm11 ...xmii--11xmii

Après les affectations de y et p, on a y = xi2 = xi+1 et p = ni mod2 = ni+1, donc la propriété est encore vérifiée après la (i + 1)-ième itération.

A la sortie de la boucle, on a alors ni = 0, donc i = k + 1 et alors res = x0m0x1m1xkmk = xn, ce que l’on voulait montrer.

Les méthodes (récursive ou itérative) qui permettent de calculer la puissance n-ième d’un entier s’adaptent sans aucune difficulté au calcul de la puissance n-ième de tout objet pour lequel cela peut avoir un sens (nombre réel, polynôme, matrice, etc). Revenons par exemple au calcul des nombres de Fibonacci Fn. On vérifie immédiatement que

( Fn+1 )   (  1 1 ) (  Fn   )
   Fn    =    1 0     Fn -1

d’où bien entendu

(       )   (      )  (    )
   Fn+1       1  1  n   F1
    Fn    =   1  0      F0

Le calcul de la puissance n-ième de la matrice pouvant être effectué en un temps O(log 2N), il en est donc de même du calcul de Fn. De la méthode récursive brutale d’ordre 2 à la méthode récursive dichotomique (diviser pour régner), nous sommes donc passés d’un temps de calcul en O(2n) à un temps de calcul en O(log 2n) ; de performances intermédiaires, la méthode itérative ou la méthode récursive d’ordre 1 ont un temps de calcul en O(n).

5.2.4 La multiplication des polynômes

Pour simplifier, nous travaillerons avec des polynômes à coefficients entiers. Si P(X) = i=0naiXi ∈ Z[X], nous stockerons les coefficients dans un tableau de taille n + 1, la valeur de ai étant stockée dans l’élément de numéro i du tableau. Pour des raisons de commodité, nous définirons une fonction qui renvoie le degré d’un polynôme (égal à la longueur du tableau diminuée de 1) et une autre fonction qui initialise un polynôme P(X) = i=0naiXi avec tous les ai égaux à 0. Pour cela nous utiliserons les deux fonctions prédéfinies de Caml : la fonction vect_length qui renvoie la longueur d’un tableau et la fonction make_vect qui renvoie un tableau de longueur donné avec tous les éléments initialisés à une valeur donnée :


#let deg p = (vect_length p)-1;;  
deg : ’a vect -> int = <fun>  
#let init_poly n = make_vect (n+1) 0;;  
init_poly : int -> int vect = <fun>


Programme 5.3: utilitaires sur les polynômes

La somme de deux polynômes est facile à définir à l’aide de la fonction suivante :


#let rec add_poly p q =  
       let m = deg p and n = deg q in  
         if m<n then add_poly q p  
         else  
           begin  
             let somme = init_poly m in  
               for i = 0 to n do  
                  somme.(i)<- p.(i)+q.(i)  
               done;  
               for i = n+1 to m do  
                  somme.(i)<- p.(i)  
               done;  
              somme  
           end;;  
add_poly : int vect -> int vect -> int vect = <fun>  
#add_poly [|1;2;3|]  [|4;5;6;7;8|];;  
- : int vect = [|5; 7; 9; 7; 8|]


Programme 5.4: somme de deux polynômes

Note de programmation  On a ici utilisé une petite astuce de programmation en rendant la fonction récursive. Ceci nous permet, dans le cas où le degré de p est strictement plus petit que le degré de q, d’intervertir facilement les paramètres en se contentant de rappeler la fonction add_poly avec les paramètres q et p, grâce à la commutativité de l’addition. Ceci améliore sensiblement la lisibilité du programme sans ajouter de complexité supplémentaire. En fait ici, la récursivité n’est que de façade.

La structure même des deux boucles indexées successives montre la proposition suivante :

Proposition 5.2.2 La fonction add_poly calcule la somme de deux polynômes p et q en un temps O(max(degp,degq)).

En ce qui concerne le produit de deux polynômes, dans une première méthode, nous utiliserons le schéma de calcul suivant

( m      ) ( n     )    m  ( n         )
  ∑  aXi   (∑  b Xj)  = ∑  (∑  a b Xi+j)
  i=0  i     j=0 j       i=0  j=0 i j

qui est le plus facile à programmer. Il suffit donc de faire parcourir à i l’intervalle [0,m], puis de faire parcourir à j l’intervalle [0,n] en ajoutant le produit aibj au terme de degré i + j du polynôme produit. Ce polynôme produit aura été préalablement initialisé à 0.


#let mult_poly p q =  
     let m = deg p and n = deg q in  
         let prod = init_poly (m+n) in  
             for i = 0 to m do  
               for j = 0 to n do  
                  prod.(i+j)<-prod.(i+j)+p.(i)*q.(j)  
               done  
             done;  
             prod;;  
 
#mult_poly [| 1;2;3;4 |] [| 1;2;1 |];;  
- : int vect = [|1; 4; 8; 12; 11; 4|]


Programme 5.5: algorithme quadratique de calcul du produit

L’imbrication des deux boucles indexées montre que le temps de calcul du produit de deux polynômes de degrés respectifs m et n est proportionnel à (m + 1)(n + 1). Pour deux polynômes de degrés inférieurs ou égaux à N, le temps de calcul par cette méthode est donc un O(N2) (d’où la qualification d’algorithme quadratique).

Prenons maintenant deux polynômes P(X) et Q(X) de degrés strictement inférieurs à N = 2k = 2M et posons P(X) = P1(X) + XMP2(X), Q(X) = Q1(X) + XMQ2(X) (partition). On a donc

P (X)Q (X ) = P1(X )Q1(X) + XM (P1(X)Q2(X )+ P2(X )Q1 (X ))+ X2M Q1 (X )Q2(X )

Il semble qu’il faille calculer 4 produits de polynômes de degré strictement inférieurs à M = N∕2 pour calculer le produit P(X)Q(X). Comme le temps nécessaire à l’addition de polynômes de degré N est un O(N), on trouve un temps de calcul T(N) = 4T(N∕2) + O(N). Avec les notations du paragraphe précédent, on a donc p = 1 et q = 4 ; on a 2p < q et donc T(N) = O(Nlog 24) = O(N2). Nous n’avons donc rien gagné par rapport à la méthode brutale.

En fait un peu d’astuce va nous éviter de calculer un des produits. Il suffit de remarquer que

P1Q2 + P2Q1 = (P1 + P2)(Q1 +Q2 )- P1Q1 - Q1Q2

autrement dit, si nous posons R1 = P1Q1, R2 = (P1 + P2)(Q1 + Q2), R3 = Q1Q2, nous avons (fusion)

PQ = R  + (R - R  - R )XM  + R X2M
       1    2    1   3        3

Nous n’avons plus que trois produits de polynômes de degrés strictement inférieurs à M = N∕2 à effectuer pour calculer le produit P(X)Q(X). Le temps de calcul nécessaire pour effectuer la partition P(X) = P1(X) + XMP2(X), Q(X) = Q1(X) + XMQ2(X) et celui nécessaire pour opérer la fusion

PQ = R1 + (R2 - R1 - R3)XM + R3X2M

étant clairement un O(N), le temps de calcul nécessaire pour calculer par cette méthode le produit de deux polynômes de degrés strictement inférieurs à N = 2k vérifiera

T(N) = 3T(N∕2)+ O (N)

ce qui nous donne grâce au théorème du paragraphe précédent

T(N) = O(N log23)

avec log 23 ~ 1.58, ce qui est un gain non négligeable.

La méthode de calcul ci-dessus, encore appelée méthode de Karatsuba, est clairement récursive, la récursivité portant sur l’entier N = 2k. Pour k = 0, on a N = 1, auquel cas les deux polynômes sont constants, leur produit est alors évident ce qui initialise la récursion.

Les polynômes de degré strictement inférieur à N = 2k seront stockés dans des tableaux de longueur N. Nous utiliserons quelques fonctions utilitaires :


let scinde_poly p =  
   let N = vect_length p in  
               (sub_vect p 0 (N/2),sub_vect p (N/2) (N/2));;  
 
let add_poly p q =  
        let N = vect_length p in  
            let somme = make_vect N 0 in  
                for i = 0 to N-1 do  
                   somme.(i)<- p.(i)+q.(i)  
                done;  
                somme;;  
 
let sub_poly p q =  
        let N = vect_length p in  
            let diff = make_vect N 0 in  
                for i = 0 to N-1 do  
                   diff.(i)<- p.(i)-q.(i)  
                done;  
                diff;;


Programme 5.6: procédures utilitaires sur les polynômes

La fonction fusion reçoit en paramètre trois polynômes r1, r2 et r3 de degrés strictement inférieurs à N et l’entier N = 2M, et renvoie le polynôme R1(X) + R2(X)XM + R3(X)X2M grâce à un calcul évident.


let fusion r1 r2 r3 N =  
          let p = make_vect (2*N) 0 in  
               for i = 0 to N-2 do  
                    p.(i)<-r1.(i)  
               done;  
              for i = 0 to N-2 do  
                    p.(N+i)<- r3.(i)  
              done;  
              for i = 0 to N-2 do  
                   p.(N/2+i)<-p.(N/2+i)+r2.(i)  
              done;  
              p;;


Programme 5.7: fusion de trois polynômes

A partir de là, il suffit de suivre la formulation mathématique pour obtenir la procédure de calcul du produit de deux polynômes :


let rec mult_poly p q =  
       let N = vect_length p in  
          if N = 1 then  
              begin  
                 let res = make_vect 1 (p.(0)*q.(0)) in res  
              end  
          else  
             begin  
                match scinde_poly p with (p1,p2) ->  
                match scinde_poly q with (q1,q2) ->  
                   let r1 = mult_poly p1 q1 and r3 = mult_poly p2 q2 in  
                       let r22 = mult_poly (add_poly p1 p2)  (add_poly q1 q2) in  
                            let r2 = sub_poly (sub_poly r22 r1) r3 in  
                                 fusion r1 r2 r3 N  
              end;;


Programme 5.8: produit rapide de polynômes

Un petit programme nous a permis de calculer le nombre de multiplications de nombres entiers effectuées au cours du calcul de (1 + X)2n par la méthode quadratique et par la méthode diviser pour régner (sachant que la multiplication est de loin l’opération la plus lente parmi toutes celles qui sont effectuées sur les coefficients)




nméthode quadratiquediviser pour régner



3 38 39
4 119 120
5 408 363
6 1497 1092
7 5722 3279



Le gain est peu sensible (voire inexistant) pour des polynômes de petit degré et ne commence à devenir sensible que pour n = 6 (résultat de degré 26 = 64). Pour un résultat de degré 128, on voit que la deuxième fonction est environ deux fois plus rapide que la fonction initiale, ce qui n’est pas négligeable1.

Note de programmation  Notre deuxième fonction n’est évidemment pas universelle. Elle suppose a priori que les polynômes sont tous deux entrés dans des tableaux de taille N = 2k. Si l’on veut en faire une fonction d’usage général, pour multiplier deux polynômes p et q de degrés respectifs m et n, il faut commencer par déterminer le plus petit entier k tel que m < 2k et n < 2k, ensuite recopier nos polynômes dans deux tableaux de taille 2k, appliquer la fonction mult_poly et enfin recopier le tableau obtenu (de taille 2k+1) dans un polynôme de degré m+n en supprimant les 2k+1 -(m+n+1) derniers coefficients qui sont nuls. Le lecteur vérifiera que ceci ne change pas la complexité de la fonction dont le temps d’exécution est donc un O(max(m,n)log 23).

5.2.5 Les opérations sur les grands nombres entiers

Soit b un entier strictement positif appelé base de numération. On sait que tout nombre entier n peut s’écrire de manière unique sous la forme n = akbk + ak-1bk-1 + + a0 avec 0 ai b - 1. Introduisons alors le polynôme Pn(X) = akXk + ak-1Xk-1 + + a0. On a n = Pn(b). On en déduit que si m et n sont deux entiers, on peut écrire mn = Pm(b)Pn(b) = Q(b) où Q(X) est le produit Pm(X)Pn(X). Le polynôme à coefficients entiers Q(X) s’écrit sous la forme Q(X) = αlXl + + α0, seulement il ne vérifie plus la condition 0 αi b - 1. Pour calculer les chiffres du produit mn dans la base b, il nous faut donc renormaliser le polynôme Q(X). Pour cela, il nous suffit de remarquer que si l’on prend i ∈ [0,k] et qi ∈ Z, le polynôme Q(X) et le polynôme Qi(X) = αlXl + + (αi+1 + qi)Xi+1 + (αi - bqi)Xi + + α0, on a Qi(b) = Q(b) ; en effet

                             i+1            i        i+1      i
Qi(X )- Q(X ) =   (αi+1i+1+ qi)X  i +(αi - bqi)X i+-1 αi+1X   - αiX
              =   qiX    - bqiX  = qi(X - b)X
s’annule au point X = b. En prenant en particulier pour qi le quotient de la division euclidienne de αi par b, on a alors 0 αi - bqi b - 1, sans perturber les coefficients α0,i-1. Il nous suffit donc pour renormaliser le polynôme Q de parcourir tous les coefficients de Q(X) à partir du bas en remplaçant successivement αi par ri = αi -bqi et αi+1 par αi+1 + qi. On finit ainsi par aboutir à un polynôme R(X) = βl+λXl+λ + + β0 (de degré supérieur à celui de Q) qui vérifie à la fois 0 βi b- 1 et R(b) = Q(b) = Pm(b)Pn(b) = mn.

L’usage d’une base de numération permet de dépasser les limitations du langage de programmation en termes de taille des entiers manipulés. Il suffit en effet de remplacer l’entier n = akbk + ak-1bk-1 + + a0 par un tableau noté à la Caml [|a0;;ak|], qui coïncide justement avec la notation utilisée pour le polynôme Pn(X) = akXk + ak-1Xk-1 + + a0 vérifiant n = Pn(b). Il nous suffit donc de définir une opération de renormalisation d’un polynôme relativement à la base b pour adapter les procédures de calcul sur les polynômes aux calculs sur les grands nombres entiers stockés sous formes de tableaux en base b.

Note de programmation  La procédure de renormalisation d’un polynôme P par rapport à la base b suit l’algorithme décrit ci-dessus en ce qui concerne tous les coefficients, sauf celui de plus haut degré. Il n’est en effet pas possible en Caml de faire croître le tableau au fur et à mesure des besoins pour tenir compte de l’augmentation du degré. Mais il suffit de remarquer qu’une fois le polynôme P(x) écrit sous la forme P(X) = βlXl + + β0 avec i ∈ [0,l - 1], 0 βi b - 1, si l’écriture de βl en base b est βl = γλbλ + + γ0, alors le polynôme normalisé de P(X) n’est autre que

(γλXλ+ ...+ γ0)Xl+ ...+ β0 = γλXλ+l+...+γ0Xl+ βl-1Xl-1+ ...+ β0

Nous utiliserons donc une fonction table_decimale qui renvoie dans un tableau les chiffres de l’écriture en base b d’un entier n. Pour cela nous stockons au fur et à mesure les chiffres dans une liste qui croît dynamiquement suivant l’algorithme évident par divisions euclidiennes successives : le coefficient de b0 est le reste de la division euclidienne de n par b, puis on remplace n par son quotient par b et on recommence jusqu’à obtenir n = 0. Une fois cette liste construite, on la retourne à l’aide de la fonction prédéfinie rev et on la transforme en tableau à l’aide de la fonction prédéfinie vect_of_list.


#let table_decimales n b =  
     let x = ref n and l = ref [] in  
        while !x<>0 do  
           l := (!x mod b)::!l;  
           x := !x/b  
        done;  
     vect_of_list (rev !l);;  
table_decimales : int -> int -> int vect = <fun>  
#table_decimales 12345678 10;;  
- : int vect = [|8; 7; 6; 5; 4; 3; 2; 1|]


Programme 5.9: calcul des chiffres de n en base b

Il nous suffit maintenant de renormaliser le polynôme suivant la base b par l’algorithme suivant :


#let renormalise p b =  
    let N = vect_length p and pp = p in  
       for i = 0 to N-2 do  
          let r = pp.(i) mod b and q = pp.(i)/b in  
             pp.(i)<-r; pp.(i+1)<-pp.(i+1)+q  
       done;  
       concat_vect  
          (sub_vect pp 0 (N-1))  
          (table_decimales pp.(N-1) b);;  
renormalise : int vect -> int -> int vect = <fun>  
#renormalise [|1;10;3;1245|] 10;;  
- : int vect = [|1; 0; 4; 5; 4; 2; 1|]


Programme 5.10: renormalisation d’un polynôme en base b

Il est alors facile de définir l’addition et la multiplication de deux nombres entiers écrits en base b dans deux tableaux :


#let add_entier m n b =  
      renormalise (add_poly m n) b;;  
#let mult_entier m n b =  
      renormalise (mult_poly m n) b;;  
add_entier : int vect -> int vect -> int -> int vect = <fun>  
mult_entier : int vect -> int vect -> int -> int vect = <fun>


Programme 5.11: opérations sur les entiers longs en base b

A titre d’exemples, nous allons calculer 100 ! en base 10 ; nous utiliserons une petite fonction d’affichage qui se contente de juxtaposer dans le bon ordre (c’est-à-dire l’ordre inverse) les éléments du tableau


#let facto n b =  (* calcule n! en base b *)  
     let m = ref [|1|] in  
        for i = 2 to n do  
          m := mult_entier !m (table_decimales i b) b  
        done;  
     !m;;  
facto : int -> int -> int vect = <fun>  
#let affiche_entier t =  
      for i = (vect_length t)-1 downto 0 do  
         print_int t.(i)  
      done;  
      print_newline();;  
affiche_entier : int vect -> unit = <fun>  
#affiche_entier (facto 100 10);;  
9332621544394415268169923885626670049071\  
5968264381621468592963895217599993229915\  
6089414639761565182862536979208272237582\  
51185210916864000000000000000000000000


Programme 5.12: factorielle en entiers longs

On constate que les polynômes manipulés par cette méthode peuvent être de très grands degrés (environ 160 pour celui décrivant l’écriture de 100! en base 10) et que la méthode diviser pour régner prend tout son intérêt pour de tels polynômes.

5.3 Recherches et tris

5.3.1 Recherche binaire ou dichotomique

La recherche d’éléments dans un ensemble est une opération essentielle en algorithmique. Dans le cas où l’ensemble n’est muni d’aucune structure, la seule possibilité est une recherche exhaustive : parcourir tout les éléments de l’ensemble jusqu’à trouver l’élément choisi. Clairement, si l’ensemble a N éléments, le temps de recherche est un O(N). Le placage d’une structure sur l’ensemble permet souvent de réduire de manière drastique les temps de recherche.

Nous allons montrer que si l’ensemble est un tableau ordonné, le temps de recherche peut être réduit en un O(log 2N) par une méthode de type diviser pour régner. Considérons en effet un ensemble X = {a0,,aN-1} ordonné de telle sorte que a0 a1 aN-1 et soit b un élément de même type que les ai. Nous voulons savoir si b appartient à X et éventuellement calculer un indice i tel que ai = b. Donnons nous pour cela m ∈{0,,N - 1}. Si b am, alors b ∈ X si et seulement si b ∈{a0,,am} ; si par contre b > am, alors b ∈ X si et seulement si b ∈{am+1,,aN-1}.

En choisissant un m de l’ordre de N∕2, on voit que le temps de recherche T(N) d’un élément b dans un ensemble à N éléments vérifiera une récurrence du type T(N) = T(N∕2) + O(1). Il s’agit ici d’une méthode diviser pour régner où la partition s’effectue en temps constant (la recherche du milieu m et la comparaison à am) et où la fusion est triviale (s’effectue en temps nul). On sait d’après un paragraphe précédent que dans ce cas T(N) = O(log 2N).

La recherche dans un tableau trié peut donc se faire de manière récursive (nous avons choisi de renvoyer l’indice i de l’élément ai tel que ai = b si b ∈ X, et - 1 si b∈∕X) ; pour cela nous définissons une procédure récursive cherche qui effectue le gros du travail : elle reçoit en paramètre les indices i et j qui sont les bornes du sous-tableau ai,ai+1,,aj dans lequel doit s’effectuer la recherche.

La procédure de recherche dichotomique se contente de recevoir les paramètres b et X, de créer l’environnement nécessaire à le recherche et d’appeler la procédure cherche avec comme paramètres les limites du tableau initial.


#let recherche_bin b X =  
   let N = vect_length X in  
      let rec cherche i j =  
         if i = j then  
            if X.(i) = b then i else -1  
         else  
            let m = (i+j)/2 in  
               if b <= X.(m) then cherche i m  
                  else cherche (m+1) j  
      in cherche 0 (N-1);;  
recherche_bin : int -> int vect -> int = <fun>  
#recherche_bin 2 [|1;3;5;7;9;11;13|];;  
- : int = -1  
#recherche_bin 5 [|1;3;5;7;9;11;13|];;  
- : int = 2


Programme 5.13: recherche dichotomique récursive

Proposition 5.3.1 La fonction de recherche binaire récursive effectue bien la recherche d’un élément b dans un tableau trié de N éléments en un temps O(log 2N).

Démonstration: l’assertion sur le temps de calcul a déjà été démontrée. Il nous suffit donc de démontrer que la fonction cherche renvoie bien le résultat de la recherche de b dans le sous-tableau ai aj. Nous allons le faire par récurrence sur j -i. C’est clair, si j -i = 0. Si j -i 1, m est la partie entière de i+ j∕2 si bien que i m i+ j∕2 < j.

En conclusion, l’appel cherche 0 (N-1) renvoie bien le résultat de la recherche de b dans {a0,,aN-1} donc dans X.

Note de programmation  La vérification de la relation i m i+j∕2 < j est une étape essentielle de la démonstration. Il serait en effet catastrophique d’avoir un croisement des indices c’est à dire qu’à un moment donné on puisse appeler la fonction cherche i j avec i > j ; c’est une des erreurs fréquentes dans les méthodes de dichotomie.

Il est clair que la récursivité ci-dessus est terminale (la fonction se termine par l’appel récursif). Ceci nous permet de construire facilement une version itérative de la recherche en introduisant


#let recherche_bin b X =  
     let N = vect_length X in  
       let i = ref 0 and j = ref (N-1) in  
          while (!j-!i)>0 do  
             let m = (!i+!j)/2 in  
               if b <= X.(m) then j :=m  
                  else i := m+1  
          done;  
          if X.(!i) = b then !i else -1;;  
recherche_bin : int -> int vect -> int = <fun>


Programme 5.14: recherche dichotomique itérative

Proposition 5.3.2 La fonction de recherche binaire itérative effectue bien la recherche d’un élément b dans un tableau trié de N éléments en un temps O(log 2N).

Démonstration: l’assertion sur le temps de calcul a déjà été démontrée. Il nous faut donc démontrer que la boucle se termine et que la fonction renvoie bien le résultat de la recherche. Pour cela nous allons exhiber un invariant de la boucle booléenne. Posons f(i,j) = -1 si b∕∈{ai,,aj} et f(i,j) = k si k est le plus petit indice tel que b = ak avec i k j. Soit in et jn les valeurs des références i et j à l’entrée dans la n-ième itération. Montrons alors que la condition in < jn et la valeur de f(in,jn) est un invariant de la boucle booléenne à l’entrée dans l’itération. A la première exécution du corps de la boucle, on a i1 = 0, j1 = N - 1 et puisque ce corps est exécuté c’est que 0 < N - 1 ; quant à f(i1,j1) il vaut tout simplement f(0,N - 1) (c’est à dire la valeur recherchée). Au début de la n-ième itération, on sait que in < jn (puisque le corps de la boucle est exécuté) ; on montre comme dans la démonstration de la fonction récursive que in m in + jn∕2 < jn. On a donc in m < m + 1 jn. Si b am, on a alors f(in,jn) = f(in,m) = f(in+1,jn+1) et à l’itération suivante, soit in+1 = jn+1 auquel cas on sortira de la boucle, soit in+1 < jn+1. Si par contre, b > am, on a alors f(in,jn) = f(m + 1,jn) = f(in+1,jn+1) et à l’itération suivante, soit in+1 = jn+1 auquel cas on sortira de la boucle, soit in+1 < jn+1.

La stricte décroissance de la suite d’entiers naturels jn - in montre que la boucle s’achève au bout d’un nombre fini d’itérations. A ce moment on a in = jn et f(in,jn) = f(i0,j0) = f(0,N - 1) puisque cette valeur est un invariant de la boucle. Or f(in,jn) = f(in,in) vaut in si b = ain et -1 sinon. On a donc bien calculé f(0,N - 1), qui était l’entier recherché.

Remarque  Le fait que le tableau soit trié est bien évidemment essentiel. Mais il est aussi essentiel que la structure soit un tableau, ce qui permet d’accéder à l’élément médian am en un temps constant. Si au lieu d’un tableau on travaillait avec une liste, cette recherche ne pourrait se faire qu’en un temps O(N) (puisqu’on ne peut accéder aux éléments que séquentiellement) et la recherche dichotomique perdrait tout son intérêt. C’est toute la différence entre une structure à accès direct (où l’on peut accéder directement à tous les éléments) et une structure à accès séquentiel. Remarquons également qu’il est stupide de commencer par trier un tableau avant d’effectuer une recherche dans ce tableau puisque le tri lui-même demande un temps en O(N log 2N) au minimum ; dans ce cas une banale recherche exhaustive en O(N) est plus performante. Le tri, puis la recherche dichotomique dans un tableau ne peuvent se justifier que si l’on effectue de nombreuses recherches dans ce tableau. Si l’on effectue p recherches dans un tableau à N éléments,

En gros, on peut dire que dès que p devient de l’ordre de log 2N, il peut être intéressant d’effectuer un tri préalable.

5.3.2 Le tri fusion

Supposons qu’un enseignant se trouve devant deux paquets de copies, chacun des paquets étant trié par ordre croissant de notes. Il cherche à réunir les deux paquets en un seul trié par ordre croissant. La solution est évidente : il prend à chaque étape la copie ayant la plus petite note parmi les sommets des deux paquets jusqu’à avoir épuisé l’un des deux.

Autrement dit la fusion de deux tableaux triés par ordre croissant peut se faire grâce à la fonction suivante, où il faut simplement prendre garde à l’épuisement possible d’un des deux tableaux :


#let fusion a b =  
    let m = vect_length a and n = vect_length b in  
      let c = make_vect (m+n) 0 and i = ref 0 and j = ref 0 in  
         for k = 0 to m+n-1 do  
            if (!i<m) then  
              begin  
               if (!j<n)  
                 then  
                   begin  
                    if a.(!i) <= b.(!j)  
                      then begin c.(k)<- a.(!i); i := !i+1 end  
                      else begin c.(k)<- b.(!j); j := !j+1 end  
                   end  
                 else  
                   begin  
                     c.(k)<-a.(!i); i := !i+1  
                   end  
               end  
            else  
               begin  
                 c.(k) <- b.(!j); j := !j+1  
               end  
         done;  
      c;;  
fusion : int vect -> int vect -> int vect = <fun>


Programme 5.15: fusion de deux tableaux triés

Dans le cas où l’on dispose d’une place supplémentaire à la fin des deux tableaux a et b, on peut y disposer une sentinelle dont on est certain qu’elle est supérieure à tous les éléments des deux tableaux. Ceci permet alors de se passer des deux tests de débordement puisque lorsque i = m, on est certain que ai = am bj et lorsque j = n, on est certain que bj = bn ai. Cependant l’usage d’une telle sentinelle n’est pas évident dans le cas général : ceci oblige à trouver un élément plus grand que tous les éléments des deux tableaux et ensuite à réserver une place supplémentaire dans chacun d’eux pour y stocker la sentinelle.

Dans tous les cas, on constate que le temps d’exécution de la fusion de deux tableaux triés de longueurs m et n est un O(m + n) puisque le corps de la boucle indexée est exécuté m + n fois.

Une fois l’idée de la fusion de deux tableaux triés acquise, la procédure de tri par fusion d’un tableau ne pose plus de difficulté. Etant donné un tableau de taille N, on le partage en deux tableaux de taille (approximative) N∕2 que l’on trie récursivement et que l’on fusionne ensuite. Bien entendu les appels récursifs se terminent lorsque le tableau est de taille 1 auquel cas il est tout trié.

La partition du tableau se fait en temps constant et la fusion des deux tableaux en un temps O(N ∕2 + N∕2) = O(N) si bien que le temps T(N) d’exécution du tri par fusion sur un tableau de N éléments vérifie T(N) = 2T(N∕2) + O(N) + O(1) = 2T(N∕2) + O(N). On en déduit que T(N) = O(N log 2N).

Il nous faut utiliser une procédure de fusion de deux sous-tableaux consécutifs d’un même tableau, plutôt que de fusion de deux tableaux indépendants. La procédure de fusion des deux tableaux triés {al,al+1,,am} et {am+1,am+2,,au} peut s’écrire de la même manière que ci-dessus :


#let fusion l m u =  
     let b = make_vect (u-l+1) 0 and i = ref l and j = ref (m+1) in  
        for k = 0 to u-l do  
           if (!i <= m)  
             then  
               begin  
                 if (!j <= u) then  
                   begin  
                     if a.(!i) <= a.(!j)  
                       then begin b.(k)<- a.(!i); i := !i+1 end  
                       else begin b.(k)<- a.(!j); j := !j+1 end  
                   end  
                 else  
                   begin b.(k)<- a.(!i); i := !i+1 end  
               end  
             else begin b.(k)<- a.(!j); j := !j+1 end  
         done;  
       for k = 0 to u-l do  
         a.(l+k)<-b.(k)  
       done;;  
fusion : int vect -> int -> int -> int -> unit = <fun>


Programme 5.16: fusion de deux sous tableaux consécutifs triés

En fait, plutôt que de réinitialiser à chaque fusion un nouveau tableau b, il est plus simple de définir un tableau auxiliaire ce qui conduit simplement à supprimer le b=make_vect (u-l+1) 0 à condition de disposer dans l’environnement de la fonction de fusion d’un tableau auxiliaire b. Ceci conduit à une procédure de tri d’un tableau par fusion sous la forme :


#let tri_fusion a =  
    let N = vect_length a in  
       let b = make_vect N 0 and i = ref 0 and j = ref 0 in  
          let fusion l m u =  
             i := l; j := m+1;  
             for k = 0 to u-l do  
               if (!i <= m)  
                 then  
                   begin  
                     if (!j <= u) then  
                        begin  
                          if a.(!i) <= a.(!j)  
                            then begin b.(k)<- a.(!i); i := !i+1 end  
                            else begin b.(k)<- a.(!j); j := !j+1 end  
                          end  
                     else  
                       begin b.(k)<- a.(!i); i := !i+1 end  
                   end  
                 else begin b.(k)<- a.(!j); j := !j+1 end  
             done;  
           for k = 0 to u-l do  
             a.(l+k)<-b.(k)  
           done (* fin de fusion *)  
       in  
         let rec tri l u =  
           if u>l then  
             begin  
               let m = (l+u)/2 in  
                  tri l m;  
                  tri (m+1) u;  
                  fusion l m u  
             end  
         in tri 0 (N-1);;  
let v = [|9;4;5;6;8;2;1;10;2;10;1;5;8;7;3|] in tri_fusion v; v;;  
tri_fusion : int vect -> unit = <fun>  
#- : int vect = [|1; 1; 2; 2; 3; 4; 5; 5; 6; 7; 8; 8; 9; 10; 10|]


Programme 5.17: tri par fusion d’un tableau

Note de programmation  L’auteur dans une première version avait fait une erreur difficilement décelable en remplaçant l’instruction let b=make_vect N 0 par let b=a ; il oubliait que dans une telle liaison,dans la mesure où un tableau est une structure modifiable en place, le tableau b devenait physiquement égal au tableau a ; toute modification du tableau b (et il y en a) se répercutait sur le tableau a avec toutes les conséquences que l’on peut imaginer. C’est un point à ne pas oublier dans la programmation en Caml.

Proposition 5.3.3 Le tri par fusion trie un tableau de N éléments en un temps O(N log 2N).

Démonstration: Nous allons commencer par montrer que la procédure de fusion fusionne bien deux sous-tableaux triés al al+1 am et am+1 am+2 au en un tableau b0 b1 bl-u qui est ensuite recopié en place dans le tableau a.

L’achèvement de la procédure de fusion est garanti par le fait qu’il s’agit d’une boucle indexée. Il nous suffit donc d’exhiber un invariant de la boucle. Pour cela nous allons utiliser un invariant quadruple à l’entrée dans le corps de la boucle

  1. l i m + 1
  2. m + 1 j u + 1
  3. i + j = l + m + 1 + k
  4. on a
                   { ai ≤ ...≤ am   si i ≤ m
b0 ≤ ...≤ bk- 1 ≤ aj ≤ ... ≤ au  si j ≤ u

    et la famille b0,,bk-1 est à une permutation près la famille

    al,...,ai- 1, am+1,...,aj-1

C’est clair avant la première itération puisqu’alors i = l, j = m + 1 et k = 0. Supposons que ce soit vrai avant la n-ième itération. Les trois premières conditions sont clairement conservées par le corps de la boucle : les deux premières puisque i n’est pas modifié si i = m + 1 et j n’est pas modifié si j = u + 1, la troisième puisque lors de l’exécution du corps de la boucle, soit i, soit j est augmenté de 1 et ce de manière exclusive, et qu’ensuite k est augmenté de 1 pour l’exécution suivante du corps de la boucle, ce qui conserve l’égalité i + j = l + m + 1 + k pour la (n + 1)-ième itération. Si maintenant ai bj alors le corps de la boucle affecte ai à bk et augmente i de 1, si bien que l’on a encore b0 bk-1 bk = ai-1 min(ai,bj) et la quatrième condition est bien vérifiée. De manière similaire, la quatrième condition est bien conservée si bj ai. Avant l’itération d’indice k = u - l, on a donc

  1. l i m + 1
  2. m + 1 j u + 1
  3. i + j = l + m + 1 + u - l = u + m + 1
  4. on a
                     {
b0 ≤ ...≤ bu-l-1 ≤  ai ≤ ...≤ am  si i ≤ m
                   aj ≤ ...≤ au   si j ≤ u

    et la famille b0,,bu-l-1 est à une permutation près la famille

    al,...,ai- 1, am+1,...,aj-1

Comme i + j = u + m + 1 avec i m + 1 et j u + 1, on a soit i = m + 1 et j = u, soit j = u + 1 et i = m. Dans le premier cas, le corps de la boucle se contente d’affecter au à bk et d’augmenter i, donc après la sortie de la boucle on a b0 bk = au et cette famille est à une permutation près la famille

a,...,a  ,a   ,...a
 l     m  m+1     u

ce qui est ce que l’on souhaitait. Le raisonnement est similaire dans le second cas.

Une fois que l’on sait que la procédure de fusion effectue bien son travail, la démonstration de la procédure tri est une banale récurrence sur u-l. Si u = l la procédure ne fait rien et le tableau est tout trié, donc elle le trie bien. Si u-l 1, la procédure trie deux sous tableaux de taille strictement inférieure, ce qu’elle fait bien par hypothèse de récurrence, puis fusionne les deux tableaux, ce qui est correct.

L’assertion sur le temps d’exécution a déjà été justifiée. Ceci achève la démonstration.

5.3.3 Le tri rapide

Il procède des mêmes idées que le tri fusion, en partant d’une autre méthode de partition. Considérons un tableau a0,,aN-1 de N éléments. Effectuons alors une permutation du tableau pour obtenir un nouveau tableau b0,,bN-1 tel que a0 = bm avec

∀i < m, bi ≤ bm  et ∀i > m, bi > bm

Il est clair que si maintenant on trie les sous-tableaux b0,,bm-1 et bm+1,,bN-1, le tableau obtenu sera trié.

Contrairement au tri fusion dans lequel la partition est triviale et la fusion plus délicate, dans le tri rapide, la fusion est inexistante (les sous-tableaux étant triés en place) et la partition cruciale. De plus, contrairement au tri fusion, le tri rapide n’a pas besoin de tableau auxiliaire pour effectuer des recopies, il utilise donc moins de mémoire que le tri fusion. Par contre, on verra que dans certains cas très particuliers, il n’a de rapide que le nom !

Pour effectuer la partition d’un sous-tableau al,,au nous allons effectuer une boucle indexée par i ∈ [l + 1,u] avec un deuxième indice m de manière à maintenir l’invariant suivant

PICT

On initialise la boucle avec i = m = l. Lorsqu’on augmente i de 1, deux cas sont à considérer

Après exécution de la totalité de la boucle, on aura donc la situation

PICT

Il suffit alors d’échanger al et am pour avoir la situation

PICT

qui est exactement ce que l’on souhaitait.

L’implémentation en Caml de cette procédure de partition peut se faire sous la forme suivante qui partitionne le tableau et retourne l’entier m où est situé le pivot :


#let partition a l u =  
    let echange i j =  
        let temp = a.(i) in  
          a.(i)<- a.(j);  
          a.(j)<- temp  
        and m = ref l and pivot = a.(l)  
    in  
      for i = l+1 to u do  
        if a.(i) <= pivot then  
           begin  
             echange i (!m+1);  
             m := !m+1  
           end  
      done;  
      echange l !m;  
      !m;;  
partition : int vect -> int -> int -> int = <fun>  
#let v = [|4;5;8;1;2;4;5;1;2|] in  
    print_int(partition v 0 8); v;;  
5- : int vect = [|2; 1; 2; 4; 1; 4; 5; 8; 5|]


Programme 5.18: partition d’un tableau

Il suffit ensuite d’insérer cette procédure de partition à l’intérieur d’une procédure de tri récursive en tenant compte de ce qu’un sous-tableau qui a au plus un élément est tout trié, et d’appeler cette procédure récursive pour le tableau entier.


let tri_rapide a =  
   let N = vect_length a and m = ref 0 and  
      echange i j =  
       let temp = a.(i) in  
          a.(i)<- a.(j);  
          a.(j)<- temp  
   in  
     let partition l u =  
       m := l;  
       let pivot = a.(l) in  
       for i = l+1 to u do  
         if a.(i) <= pivot then  
           begin  
             echange i (!m+1);  
             m := !m+1  
           end  
       done;  
       echange l !m  
     in  
       let rec tri l u =  
          if u>l then  
            begin  
              partition l u;  
              tri l (!m-1);  
              tri (!m+1) u  
            end  
        in tri 0 (N-1);;  
tri_rapide : int vect -> unit = <fun>  
 
#let v = [|9;1;2;8;6;5;8;9;7;10;1;1;5;4;6;7;8;9|] in tri_rapide v;v;;  
- : int vect = [|1; 1; 1; 2; 4; 5; 5; 6; 6; 7; 7; 8; 8; 8; 9; 9; 9; 10|]


Programme 5.19: tri rapide d’un tableau

Proposition 5.3.4 Le tri rapide trie un tableau de N éléments en un temps O(N log 2N) dans le cas moyen et O(N2) dans le cas le plus défavorable (tableau déjà trié).

Démonstration: La construction même de la procédure de partition à l’aide d’un invariant de boucle montre qu’elle effectue bien son travail, c’est à dire que lors du retour partition l u, l’indice m a été modifié de telle sorte que l m u, i ∈ [l,m-1], ai am et i ∈ [m+1,u], ai > am. Montrons maintenant par récurrence sur u - l que la procédure tri l u trie bien le sous-tableau de l à u. C’est clair si u-l 0, puisque la procédure ne fait rien et que le sous-tableau est tout trié. Si u - l 1, alors la procédure tri l u commence par appeler la procédure de partition. A la sortie de la procédure, on a donc l m u, i ∈ [l,m - 1], ai am et i ∈ [m + 1,u], ai > am ; ensuite la procédure s’appelle récursivement sur les deux sous-tableaux de tailles strictement inférieures al,,am-1 et am+1,,au ; par hypothèse de récurrence la procédure trie bien ces sous-tableaux et à la fin de la procédure on a donc

al ≤ ...≤ am -1(am)am+1 ≤ ...≤ au

ce qui achève la démonstration par récurrence.

Le temps d’exécution de la procédure tri l u dépend essentiellement de la taille des sous-tableaux al,,am-1 et am+1,,au. Si on note t(l,u) le temps d’exécution de la procédure tri l u, comme le temps d’exécution de la procédure de partition est clairement proportionnel à u - l, on aura t(l,u) = t(l,m - 1) + t(m + 1,l) + α(l - u). Le cas le plus défavorable est le cas m = l ou m = u c’est à dire où l’élément al est le plus petit ou le plus grand élément du sous-tableau. Dans ce cas, l’un des sous-tableaux est vide, et l’autre est de taille u - l. Si ceci se produit tout au long du tri, c’est à dire si le tableau initial est déjà trié par ordre croissant ou décroissant, la formule de récurrence concernant le temps T(N) de tri du tableau sera T(N) = T(N - 1) + O(N) et on sait alors que T(N) = 0(N2).

Si par contre, le tableau n’est en aucune façon ordonné et donc parfaitement aléatoire, on peut espérer2 que m est de l’ordre de u+ l∕2, auquel cas les deux sous-tableaux sont de taille moitié. La formule de récurrence devient alors T(N) = 2T(N∕2)+O(N) et on obtient alors T(N) = O(N log 2N).