Recherche dichotomique dans une liste triée

Où l'on fait une recherche efficace, et l'on prouve que ce que l'on fait est correct…

IPT SupNSI 1èreAlgo

Dans cet article, nous nous intéressons à l'algorithme de recherche dichotomique dans une liste triée. Nous présentons l'algorithme de base, quelques variantes en comparant leurs vitesses, et parlons preuve de programme.

Présentation de l'algorithme

L'idée de base est simple :

  • On dispose d'un tableau t trié de valeurs, et on cherche à déterminer si une valeur v est présente dans le tableau.

  • Pour cela, on procède par dichotomie :

    • On regarde l'élément du milieu du tableau et on le compare à v.
    • S'ils sont égaux, on a gagné, sinon, on poursuit la recherche dans la première ou la second moitié du tableau.
  • Si à la fin, on a réduit la portion du tableau au point qu'il ne contient plus aucun élément, la recherche s'arrête sur un echec : v n'est pas présent dans t.

Illustrons cela :

29

3

1

9

2

13

3

15

4

21

5

27

6

30

7

35

8

39

9

40

10

46

11

52

12

Voilà pour la théorie, passons à la pratique. Et là, cela se complique. C'est un problème connu que cet algorithme est délicat à programmer sans erreur (voir la page Wikipedia en anglais à ce sujet et plus particulièrement ).

Although the basic idea of binary search is comparatively straightforward, the details can be surprisingly tricky, and many good programmers have done it wrong the first few times they tried.
— D. E. Knuth, The Art of Computer Programming

Pour bien comprendre le comportement de l'algorithme et donc le programmer correctement, il est indispensable d'avoir au moins une petite idée de la correction et de la terminaison de l'algorithme. Pour le premier, la notion d'invariant de boucle permet de formaliser cela.

Invariants de boucle, petit rappel

Exposons brièvement quelques idées sur les invariants de boucle. On peut les voir comme la version « preuve de programmes » de l'hypothèse de récurrence. Dans une démonstration par récurrence, on a typiquement une hypothèse de récurrence Pn{`P_n`} que l'on utilise dans deux étapes :

  • lors de l'initialisation, pour montrer qu'à un certain rang n0{`n_0`} de départ, l'hypothèse de récurrence est bien vérifiée,
  • lors de l'hérédité, pour montrer que que si l'hypothèse de récurrence est vrai à un certain rang n{`n`} (supérieur ou égal à n0{`n_0`}), elle est aussi vraie au rang suivant n+1{`n + 1`} :
Pn    Pn+1{`P_n \\implies P_{n + 1}`}

On peut alors en déduire que l'hypothèse de récurrence est vrai pour tout rang supérieur à n0{`n_0`} :

nn0, Pn{`\\forall \\, n \\geq n_0,\\ P_n`}

Un invariant de boucle procède de la même idée :

  • Si une propriété est vrai en entrant dans une boucle,
  • et si elle est préservée par le corps de la boucle (autrement dit si elle est vérifiée au début, elle l'est encore après exécution du corps de la boucle),

alors elle vérifiée à chaque entrée de boucle… et aussi en sortant de la boucle (on parle ici d'une boucle while).

La version classique de l'algorithme

On peut écrire l'algorithme ainsi :

def dichotomie(t, v):
    a = 0
    b = len(t) - 1
    while a <= b:
        m = (a + b) // 2
        if t[m] == v:
            # on a trouvé v
            return True
        elif t[m] < v:
            a = m + 1
        else:
            b = m - 1
    # on a a > b
    return False

Les variables a et b représentent les bornes entre lesquels on recherche v dans t. Autrement dit, avec les notations de Python, on cherche v dans t[a:b + 1]. Comment exprimer cela proprement en termes d'invariants ?

Preuve de correction

Ici, notre invariant est :

vt    vt[a:b+1]&#123;`\\mathtt v \\in \\mathtt t \\implies \\mathtt v \\in \\mathtt{t[a:b + 1]&#125;`}

Autrement dit, si v est présente dans t, alors c'est nécessairement entre les indices a et b (inclus).

Pour prouver qu'il s'agit bien d'un invariant de boucle, on va supposer que notre propriété est vérifiée au début du corps de la boucle, et on va prouver que c'est encore le cas à la fin.

Après avoir défini m, examinons les tests successifs. Tout d'abord, si l'on a bien t[m] == v, alors on renvoie True ce qui est correct, ayant effectivement m dans t. Dans ce cas, l'invariant ne sert à rien puisque l'on a effectivement trouvé l'élément recherché.

Maintenant, deux cas sont possibles :

  • si t[m] < v, puisque t est trié, la valeur v ne peut pas être présente dans t[a:m + 1]. Ainsi, si v in t[a:b + 1], alors nécessairement v in t[m + 1:b + 1]. Ainsi, en posant a = m + 1, on a v in t[a:b + 1].
  • Sinon, sachant de plus que t[m] != v, cela signifie que t[m] > v et donc que si v in t[a:b + 1], alors v in t[a:m]. Ainsi, en posant b = m - 1, on a v in t[a:b + 1].

Dans tous les cas, à la fin de l'exécution de la boucle (si l'on y est toujours), la propriété

vt    vt[a:b+1]&#123;`\\mathtt v \\in \\mathtt t \\implies \\mathtt v \\in \\mathtt{t[a:b + 1]&#125;`}

est encore vérifiée. C'est donc bien un invariant de la boucle.

Ce que l'on vient de prouver correspond à la phase d'hérédité d'une démonstration par récurrence. Pour terminer ce type de démonstration, il faut aussi une phase d'initialisation. Ainsi, si notre hypothèse de récurrence est vérifiée pour un certain rang, et qu'elle est préservée en passant d'un rang au suivant, alors elle est vérifiée pour tous les rangs suivants.

De même, vérifions que notre invariant est vérifié avant d'entrer dans la boucle. À ce moment, on a a == 0 et b == len(t) - 1 et donc t[a:b + 1] == t[0:len(t)] == t. Notre invariant correspond alors à la tautologie

vt    vt,&#123;`\\mathtt v \\in \\mathtt t \\implies \\mathtt v \\in \\mathtt{t&#125;,`}

il est donc bien vérifié.

Pour résumer, notre propriété est vérifiée en entrée de boucle et c'est un invariant de la boucle. En conséquent, elle reste vérifiée en sortant de la boucle (quelque soit le nombre d'itérations effectué). De plus, à ce moment, on a a > b et donc t[a:b + 1] = [] et on peut donc affirmer que v in t => v in [].

Comme v in [] est nécessairement faux, puisque la liste vide ne contient aucune valeur, on en déduit que v not in t : on peut conclure qu'en sortie de boucle, v n'est pas présent dans la liste.

Comme annoncé précédemment, notre invariant de boucle ne sert pas à montrer la correction en cas de succés (puisque l'on a effectivement trouvé un indice m tel que t[m] == v) mais en cas d'échec : bien que l'on n'a pas inspecté toutes les valeurs de t, on a prouvé que v n'est pas présent dans t. On a cherché aux bons endroits et puisque l'on n'a pas trouvé cette valeur, cela signifie qu'elle n'y est pas.

Terminaison

Insistons sur le fait que la présence d'un invariant de boucle ne présume en rien de la terminaison de l'algorithme ; il ne prouve pas que l'on va sortir de la boucle. Par contre, si l'on sort de la boucle, alors le résultat est correct.

Pour la terminaison, prouvons que si au début d'une itération, a et b sont tels que ba2n2&#123;`b - a \\leq 2^n - 2`&#125; pour un certain entier n&#123;`n`&#125;, alors à la fin de l'itération, on a

ba2n12&#123;`b - a \\leq 2^{n - 1&#125; - 2`}

En effet, on a

ma=ba2etbm=ba2&#123;`m - a = \\Bigl\\lfloor \\dfrac {b - a&#125; 2 \\Bigr\\rfloor \\qquad \\mathrm{et} \\qquad b - m = \\Bigl\\lceil \\dfrac {b - a} 2 \\Bigr\\rceil`}

et, par croissance de   &#123;`\\lfloor \\ \\cdot \\ \\rfloor`&#125; et   &#123;`\\lceil \\ \\cdot \\ \\rceil`&#125;, si ba2n2&#123;`b - a \\leq 2^n - 2`&#125;, alors ma&#123;`m - a`&#125; et bm&#123;`b - m`&#125; sont tous les deux inférieurs ou égaux à 2n11&#123;`2^{n - 1&#125; - 1`}. Ainsi, dans tous les cas, après l'affectation a = m + 1 ou b = m - 1, on a bien

ba2n12&#123;`b - a \\leq 2^{n - 1&#125; - 2`}

L'exposant de 2&#123;`2`&#125; dans la majoration décroit donc de 1&#123;`1`&#125; à chaque itération, etaprès un nombre suffisamment grand d'itérations, l'écart b - a va donc devenir négatif, ce qui est incompatible avec la condition de boucle b >= a.

Ainsi, on est sûr de sortir de la boucle, l'algorithme termine.

Analyse de complexité

La preuve précédente de terminaison nous permet de déterminer la complexité au pire des cas de l'algorithme. En effet, le corps de la boucle est en O(1)&#123;`O(1)`&#125; et si 12n2&#123;`\\ell - 1 \\leq 2^n - 2`&#125; (avec &#123;`\\ell`&#125; la longueur de t), alors on effectue au maximum n itérations successives avec

n=log2(+1)=log2()+1&#123;`n = \\bigl\\lceil \\log_2(\\ell + 1) \\bigr\\rceil = \\bigl\\lfloor \\log_2(\\ell) + 1\\bigr\\rfloor`&#125;

La complexité au pire des cas est donc en O(log2)&#123;`O(\\log_2 \\ell)`&#125;.

Une autre version

Dans la version « classique » de l'algorithme, on effectue deux tests par itération, puisque l'on effectue le test d'égalité à chaque fois. On peut accélerer l'algorithme en ne faisant qu'un unique test d'égalité, en sortie de boucle. Elle décrite ici et peut s'écrire ainsi en Python (avec un peu de refactoring, puisque la variable a correspond à L du code original tandis que b correspond à R + 1) :

def dicho2(t, v):
    a = 0
    b = len(t)
    if b == 0:
        # il faut traîter le cas
        # où la liste est vide
        return False
    while b > a + 1:
        m = (a + b) // 2
        if t[m] > v:
            b = m
        else:
            a = m
    return t[a] == v

On pourra noter que sur la page Wikipedia indiquée plus haut, le test pour savoir si la liste est vide ne figure pas. L'algorithme est donc faux. Cela va dans le sens de la remarque initiale : la recherche par dichotomie, sous son apparente simplicité, est plus délicate à programmer que l'on croit.

Cette fois-ci, l'invariant de boucle est vt    vt[a:b]&#123;`\\mathtt v \\in \\mathtt t \\implies \\mathtt v \\in \\mathtt{t[a:b]&#125;`} et les preuves de correction et de terminaison sont similaires à la version classique (voire plus faciles).

Concernant la complexité, dans dicho2, on ne fait qu'un test par itération au lieu de 2 dans dicho mais globalement, dicho2 fait plus d'itérations que dicho. En particulier, en cas de succés, dans dicho on s'arrête tout de suite alors que dans dicho2, la boucle continue jusqu'à n'avoir plus qu'un seul élément dans l'encadrement.

Néanmoins, dicho2 est en général entre 15&#123;`15`&#125; et 20%&#123;`20\\%`&#125; plus rapide, comme l'illustre le test suivant :

>>> import timeit
>>> import random
>>> l1 = [random.randint(1, 1000000) for _ in range(10000)]
>>> l1.sort()
>>> l2 = [x + 1 for x in range(1000000)]
>>> timeit.timeit("[dicho(l1, x) for x in l2]",
                  "from __main__ import dicho, l1, l2",
                  number=1)
3.9517973890001485
>>> timeit.timeit("[dicho2(l1, x) for x in l2]",
                  "from __main__ import dicho2, l1, l2",
                  number=1)
3.2702691959998447

Si il y a beaucoup de succés et peu d'échecs, dicho2 reste plus rapide que dicho. Pour illustrer cela, considérons le cas extrême où l1 et l2 contiennent toutes les valeurs d'un intervalle donné d'entiers.

>>> l1 = [x + 1 for x in range(1000000)]
>>> l2 = l1
>>> timeit.timeit("[dicho1(l1, x) for x in l2]",
                  "from __main__ import dicho, l1, l2",
                  number=1)
5.1099315159999605
>>> timeit.timeit("[dicho2(l1, x) for x in l2]",
                  "from __main__ import dicho2, l1, l2",
                  number=1)
4.5058743929998855

La raison est simple : si l'on suppose que la recherche des différents éléments du tableau est équiprobable (on ne s'intéresse ici qu'aux succés, pour lesquels le nombre d'itérations varie sensiblement), alors dans dicho, une valeur va être accessible après une itération, deux valeurs au bout de deux itérations, quatres valeurs au bout de trois itérations et, plus généralement, 2k1&#123;`2^{k - 1&#125;`} valeurs au bout de k&#123;`k`&#125; itérations. Si le tableau contient 2n1&#123;`2^n - 1`&#125; éléments (ce qui correspond à un arbre binaire complet, ou parfait selon la terminologie utilisée), le nombre moyen d'itérations en cas de succés est donc :

12n1k=1nk2k1=n2n2n+12n1n1&#123;`\\frac 1 {2^n - 1&#125; \\sum_{k = 1}^n k 2^{k - 1} = \\frac {n 2^n - 2^n + 1}{2^n - 1} \\geq n - 1`}

Le gain de s'arrêter plus tôt en cas de succès ne fait gagner en moyenne qu'au plus une itération pour dicho par rapport à dicho2, ce qui ne compense pas la différence du nombre de tests par itération.

Notons que l'on peut encore accélérer la fonction dicho2. La cible est le + 1 du test de boucle b > a + 1. En effet, après calcul de m, cela revient à m != a (à vous de voir pourquoi). On peut donc écrire :

def dicho2bis(t, v):
    a = 0
    b = len(t)
    if b == 0:
      return False
    while True:
        m = (a + b) // 2
        if a == m:
          return t[a] == v
        if t[m] > v:
            b = m
        else:
            a = m

Une petite comparaison de temps d'exécution illustre le gain de temps :

>>> l1 = [random.randint(1, 200000) for _ in range(100000)]
>>> l1.sort()
>>> l2 = [x + 1 for x in range(200000)]
>>> timeit.timeit("[dicho1(l1, x) for x in l2]",
                  "from __main__ import dicho1, l1, l2",
                  number=10)
8.873610426000141
>>> timeit.timeit("[dicho2(l1, x) for x in l2]",
                  "from __main__ import dicho2, l1, l2",
                  number=10)
7.439924989999781
>>> timeit.timeit("[dicho2bis(l1, x) for x in l2]",
                  "from __main__ import dicho2bis, l1, l2",
                  number=10)
6.409183945000223

Le diable est dans les détails

Comme indiqué au début de l'article, l'implémentation de la recherche dichotomique est très délicate. L'exemple avec l'erreur dans la présentation par Wikipedia de l'algorithme amélioré en est un bon exemple.

Pour finir, en dépit de tout le soin apporté aux preuves de correction de boucle et de terminaison, les fonctions que j'ai présentées sont incorrectes… dans la plupart des langages de programmation. Bien sûr, la grande majorité des cas, elles retourneront le bon résultat mais, dans l'absolu, il peut arriver que ce ne soit pas le cas.

Le problème vient de l'addition des indices a et b pour en trouver calculer la moyenne. En général, a et b sont des pointeurs, des adresses mémoires. Et leur somme peut dépasser la taille de représentation des adresses. En effet, ce n'est pas parce que a et b sont représentables en 64 bits que a + b l'est aussi. Il peut y avoir un dépassement. Dans ce cas (qui peut arriver en pratique), le résultat obtenu en effectuant (a + b) // 2 ne correspondra pas à la moyenne attendue.

Cette erreur est restée pendant plusieurs décennies dans divers bibliothèques (en Java, notamment) avant d'être détectée, comme raconté ici en anglais.

Pour éviter cela, une version rigoureuse est

m = a + (b - a) // 2

En Python, cependant, on n'a pas ce problème (en Python3 du moins), car les entiers peuvent être arbitrairement grands. De là à dire que c'est une bonne chose…