Le tri linéaire multi-usage

16 décembre 2006
c radix sort

Radix Sort

Je veux avoir une routine de tri qui ne soit pas aussi lente qu'un Heap Sort ou autre QuickSort, et pour cela utiliser une méthode digne de ce qu'on apprend quand on programme en assembleur 68000 avec un nombre ridicule de cycles par frames : le radix !

Le plus facile pour comprendre la méthode est de la voir à l'oeuvre. Prenons la liste de nombres suivante qui se trouve en désordre :

{ 12, 203, 3, 17, 801, 76, 45 }

Si on trie cette liste par rapport aux chiffres de unités on obtient :

{ 801, 12, 203, 3, 45, 76, 17 }

Si on trie cette liste intermédiaire par rapport aux dizaines :

{ 801, 203, 3, 12, 17, 45, 76 }

On voit ici que les cinq derniers chiffres sont triés. Passons aux centaines à présent :

{ 3, 12, 17, 45, 76, 203, 801 }

Et la liste est triée !

Le but est donc d'implémenter cette méthode non pas sur des chiffres en base 10 mais en base 2 bien évidemment. Pour réaliser ces passes successives on peut constater que l'on peut exécuter des passes de préparation sur la liste de nombres à trier (une passe par radix) : une passe pour compter les chiffres et une autre pour en déduire des segments. Je m'explique... Pour la liste de départ on obtiendrait les compteurs suivant :

{ 0, 1, 1, 2, 0, 1, 1, 1, 0, 0, 0 } pour les unités
{ 3, 2, 0, 0, 1, 0, 0, 1, 0, 0, 0 } pour les dizaines
    { 5, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0 } pour les centaines

Ces tableaux contiennent, à l'indice i, le nombre d'éléments qui possèdent le chiffre i à la position 0 (unités), 1 (dizaines) et 2 (centaines). La seule originalité de l'algorithme est de constater que, de ces listes de compteurs, on peut déduire des listes d'offsets qui indiquent le départ des segments des nombres qui possèdent le radix considéré :

{ 0, 0, 1, 2, 4, 4, 5, 6, 7, 7, 7 } pour les unités
{ 0, 3, 5, 5, 5, 6, 6, 6, 7, 7, 7 } pour les dizaines
    { 0, 5, 5, 6, 6, 6, 6, 6, 6, 6, 7 } pour les centaines

A partir de ces listes d'offsets on peut exécuter les n passes de tri avec une inner loop simpliste, que l'on peut illustrer en revenant à la liste de départ de nombres à trier, en n'exécutant que la passe des unités :

Liste d'offset pour la passe considérée : { 0, 0, 1, 2, 4, 4, 5, 6, 7, 7, 7 }
            Résultat de la passe courante : { X, X, X, X, X, X, X }
Nombre '12', radix=2, offset=1 que l'on incrémente dans sa liste d'offset
Liste d'offset pour la passe considérée : { 0, 0, 2, 2, 4, 4, 5, 6, 7, 7, 7 }
    Résultat de la passe courante : { X, 12, X, X, X, X, X }
Nombre '203', radix=3, offset=2
Liste d'offset pour la passe considérée : { 0, 0, 2, 3, 4, 4, 5, 6, 7, 7, 7 }
    Résultat de la passe courante : { X, 12, 203, X, X, X, X }
Nombre '3', radix=3, offset=3
Liste d'offset pour la passe considérée : { 0, 0, 2, 4, 4, 4, 5, 6, 7, 7, 7 }
    Résultat de la passe courante : { X, 12, 203, 3, X, X, X }
Nombre '17', radix=7, offset=6
Liste d'offset pour la passe considérée : { 0, 0, 2, 4, 4, 4, 5, 7, 7, 7, 7 }
    Résultat de la passe courante : { X, 12, 203, 3, X, X, 17 }
Nombre '801', radix=1, offset=0
Liste d'offset pour la passe considérée : { 0, 1, 2, 4, 4, 4, 5, 7, 7, 7, 7 }
    Résultat de la passe courante : { 1, 12, 203, 3, X, X, 17 }
Nombre '76', radix=6, offset=5
Liste d'offset pour la passe considérée : { 0, 0, 2, 4, 4, 4, 6, 7, 7, 7, 7 }
    Résultat de la passe courante : { 1, 12, 203, 3, X, 76, 17 }
Nombre '45', radix=5, offset=4
Liste d'offset pour la passe considérée : { 0, 0, 2, 4, 5, 4, 6, 7, 7, 7, 7 }
    Résultat de la passe courante : { 1, 12, 203, 3, 45, 76, 17 }

Et la passe est terminée : les nombres sont effectivement triés selon leurs unités. En première implémentation on obtient le code suivant (qui ne fonctionne pas encore pour les nombres négatifs) :

void Radix_sort(int* values, int nbValues, bool verbose)
{
    int i, j;
    //On établit les décomptes des chiffres en fonctions de leurs positions
    int**   Counters = new int*[4];
    for(i=0; i<4; i++)      for(Counters[i] = new int[256], j=0; j<256; j++)        Counters[i][j] = 0;
    //On procède aux décomptes d'occurence des chiffres pour toutes les valeurs à trier
    int Value;
    for(i=0; i<nbValues; i++)   for(Value = values[i], j=0; j<4; j++, Value>>=8)    Counters[j][Value&0xFF]++;
    //Arrivé ici, on a par exemple : Counters[2][4] représente le nombre de valeurs à trier qui ont le chiffre '4' en deuxième position
    if (verbose)    for(i=0; i<4; i++)      Radix_dump(Counters[i], 256);
    //On cumule les compteurs précédents pour qu'ils représentent les offsets des valeurs 'triées' à l'intérieur d'une liste de tri
    int CurOffset;
    for(j=0; j<4; j++)
        for(i=0, CurOffset = 0; i<256; i++)     {
        int CurCounter = Counters[j][i];
        Counters[j][i] = CurOffset;
        CurOffset += CurCounter;
    }
    //Arrivé ici Counters[2][4], par exemple, représente l'offset de départ auquel on peut placer les valeurs qui ont 4 en deuxième position.
    if (verbose)    for(i=0; i<4; i++)      Radix_dump(Counters[i], 256);
    //On va trier sur deux listes en parallèle
    int*    List1 = new int[nbValues];
    int*    List2 = new int[nbValues];
    int*    SrcList = values;
    int*    DestList = List1;
    for(j=0; j<4; j++)      {
        int*    Offsets = Counters[j];
        //Voilà la ligne immonde !
        for(i=0; i<nbValues; i++)       *(DestList + Offsets[(SrcList[i]>>(j*8))&0xFF]++) = SrcList[i];
        if (SrcList == values)      SrcList = List2;
        int*    TmpList = SrcList;
        SrcList = DestList;
        DestList = TmpList;
    }
    Radix_dump(DestList, nbValues);
    //Releases des listes temporaires
    ReleaseArray(List1)
    ReleaseArray(List2)
    //Release des listes internes du radix
    for(i=0; i<4; i++)      ReleaseArray(Counters[i])
    ReleaseArray(Counters)
}
On peut voir que le réel code de tri ne prend qu'un quart de la source.

La ligne de code immonde de la boucle de tri interne est l'équivalente des lignes ci-dessous, plus lisibles :

int Value = SrcList[i];
int Digit = (Value>>(j*4))&0xF;
int*    Ptr = DestList + Offsets[Digit];
*Ptr = Value;
Offsets[Digit]++;

Pour tester cette routine, appelons-là avec le code suivant :

#define     NBVALUES    12
int     Src[NBVALUES] = { 170, 13, 75, 45, 90, 1985, 24, 802, 2, 66, 3, 76 };
    Radix_sort(Src, NBVALUES, false);

Une analyse rapide nous donne, pour N éléments à trier :

Occupation mémoire : 2*N+128+8 int (pour info, l'opération de tri de 1000 éléments entraîne l'allocation, en chiffrage-4, de 2136 int = 8544 octets = 8.34 ko

Temps machine : Init1{4*256*1W} + Init2{N*(1R+1W)} + Init3{4*256*(1R+1W)} + BoucleInterne{4*(1R+N*(2R+2W))} = 1028R+2048W + N*(9R+9W)

On est donc bien en O(N) pour la mémoire ET pour le temps machine.

Ã?a n'est pas encore fini !

Ceci étant fait, il reste encore beaucoup de petites choses à faire et à arranger :

voir comment trier des clés mais garder leur correspondance s'il y en a une gérer des doublons de valeurs à trier : ils sont gérés nativement, sans avoir rien à faire ! gérer les valeurs négatives : les valeurs négatives étaient déjà séparées des valeurs positives dans le tri de base. Il a juste fallu faire une recopie à la fin... J'aurais préféré faire l'équivalent d'un buffer cyclique mais ça ne marche pas des masses. gérer les flottants faire un reverse sort une liste Counters[x] qui ne contient qu'une seule valeur différente de zéro peut entraîner une optim ! ne peut-on faire du in-place ? Du in-place avec une seule liste à priori non, mais avec aucune autre liste que celles d'entrées OUI ! Le seul souci étant que la liste avec les valeurs a trier est trashée après l'appel... faire du radix un service qui s'initialise en créant ses tableaux de compteurs une seule fois. Transformer les fonctions en une classe fournissant un service de tri thread-safe. 
void Radix_dump(int* toDump, int nb, int itemSize)
{
    int i, j;
    for(i=0; i<nb*itemSize; i+=itemSize)        {
        for(j=0; j<itemSize; j++)   {
            if (j)      printf("-");
            printf("%d", toDump[i+j]);
        }
        if (i != (nb-1)*itemSize)       printf(", ");
    }
    printf("\n");
}

//Le tableau de compteurs, utilisé pour chaque appel à la fonction de tri
int**   Counters = NULL;

//Allocation mémoire des listes internes
bool Radix_init()
{
    if (!Counters)      {
        Counters = new int*[4];
        memset(Counters, 0, 4*sizeof(int*));
    }
    int i;
    for(i=0; i<4; i++)      {
        if (!Counters[i])   Counters[i] = new int[256];
        memset(Counters[i], 0, 256*sizeof(int));
    }
    return true;
}

//Release des listes internes du radix
bool Radix_release()
{
    int i;
    for(i=0; i<4; i++)      ReleaseArray(Counters[i])
    ReleaseArray(Counters)
}

//Trie par la méthode "Radix" une liste de n-tuples <clé>-<donnée>.
//Les données peuvent avoir n fois la taille d'une clé (<itemSize>=n) ou être carrément absent (<itemSize>=0).
//La liste triée figure dans <result>, la liste à trier, <values>, est modifiée (attention !).
bool    Radix_sort(int* values, int nbValues, int itemSize, int* result)
{
    int i, j, Value;
    if (!Radix_init())      return false;
    //On procède aux décomptes d'occurence des chiffres pour toutes les valeurs à trier
    int NbNegativeValues = 0;
    for(i=0; i<nbValues*itemSize; i+=itemSize)      {
        Value = values[i];
        if (Value < 0)      NbNegativeValues++;
        for(j=0; j<4; j++, Value>>=8)   Counters[j][Value&0xFF]++;
    }
    //On cumule les compteurs pour qu'ils représentent les offsets des valeurs 'triées' à l'intérieur d'une liste de tri
    int CurOffset;
    for(j=0; j<4; j++)
        for(i=0, CurOffset = 0; i<256; i++)     {
        int CurCounter = Counters[j][i];
        Counters[j][i] = CurOffset;
        CurOffset += CurCounter*itemSize;
    }
    //On va trier en parallèle sur les deux listes fournies en entrée : values et result
    int*    SrcList = values;
    int*    DestList = result;
    for(j=0; j<4; j++)      {
        int*    Offsets = Counters[j];
        int*    SrcPtr = SrcList;
        for(i=0; i<nbValues*itemSize; i+=itemSize)      {
            int     ValueIndex = ((*SrcPtr)>>(j*8))&0xFF;
            int*    DestPtr = DestList + Offsets[ValueIndex];
            while (SrcPtr < SrcList+i+itemSize)     *DestPtr++ = *SrcPtr++;
            Offsets[ValueIndex] += itemSize;
        }
        int*    TmpList = SrcList;  SrcList = DestList; DestList = TmpList;     //Swap
    }
    //On réordonne les deux segments : valeurs négatives puis valeurs positives
    SrcList = values+(nbValues-NbNegativeValues)*itemSize;
    DestList = result;
    for(i=0; i<NbNegativeValues*itemSize; i++)          *DestList++ = *SrcList++;
    SrcList = values;
    for(i=0; i<(nbValues-NbNegativeValues)*itemSize; i++)   *DestList++ = *SrcList++;
    //
    return true;
}
Un code un peu plus stable et sûr.

Tri inversé

Pour que la même routine puisse fournir le tri ascendant (comme ici) ou descendant il suffit de changer la partie qui détermine les offsets à partir des compteurs d'occurence de radix. Comme suit :

    //reverse est un bool en entrée de la fonction de tri.
    int iLimit = reverse ? -1 : 256;
    int iStep = reverse ? -1 : 1;
    for(j=0; j<4; j++)      {
        int CurOffset = 0;
        i = reverse ? 255 : 0;
        while (i != iLimit)     {
            int CurCounter = Counters[j][i];
            Counters[j][i] = CurOffset;
            CurOffset += CurCounter*itemSize;
            i += iStep;
        }
    }

Les nombres flottants

La gestion des nombres flottants se fait naturellement avec un tri qui ordonne les éléments en fonction de radix : ce qui marche naturellement sur des entiers fonctionne également sur des flottants IEEE 754, avec le bout de code suivant tout de même :

bool    Radix_sortf(float* values, int nbValues, int itemSize, float* result, bool reverse)
{
    //Première passe : on trie toutes les valeurs, les négatives seront inversées.
    if (!Radix_sort((int*) values, nbValues, itemSize, (int*) result, reverse))     return false;
    //Deuxième passe : on ne trie que les négatives pour les inverser.
    bool    Ret = Radix_sort((int*) result, Radix_NbNegativeValues, itemSize, (int*) values, !reverse);
    int i;
    for(i=0; i<Radix_NbNegativeValues; i++)     result[i] = values[i];
    return Ret;
}

Pourquoi ce bout de code ? Parce que sans rien changer, les flottants positifs sont gérés correctement mais les négatifs sont inversés, d'où...

Il n'y a pas de cas pathologique (même NaN et inf sont correctement gérés).

Utilisation du Radix Sort

Module serveur :: requêtes géographiques

Accueil