<< Retour à la liste <<

PRODUIRE DES NOMBRES « AU HASARD »

 

Avant le préambule :

Dans cet article, je me garde bien de définir rigoureusement le hasard ainsi que de démontrer les résultats avancés. Ceux que le sujet intéresse se reporteront utilement aux références citées.

Les sources et exemples de cet article sont en {C} parce que c'est ce que je « parle » le mieux. Ils se trouvent dans le fichier HASARD.C qui accompagne cet article.

Une petite précision mathématique, les générateurs que nous allons voir sont périodiques, c'est à dire que les séquences produites se répètent. Le nombre de valeurs générées entre deux répétitions de la séquence est la période du générateur (je refuse le barbarisme « longueur de la période »).

Voilà, bein on y va ....

 

Le préambule :

Produire des nombres sans suite apparente est une nécessité dans de nombreux domaines de l'informatique. C'est le cas des jeux, du traitement du signal, de la synthèse d'images pour ne citer que ceux là.

On ne dispose malheureusement pas dans un ordinateur de vraie source de hasard. Tout est synchronisé sur des horloges très stables et les calculs sont fait avec une régularité et une fiabilité que les horlogers suisses ont enviée jusqu'à en mettre plein les montres.

Alors comme le matériel est d'un secours bien pauvre dans ce domaine, c'est vers le logiciel qu'il faut se tourner pour produire des séquences de nombres qui ont toutes les apparences du hasard. Mais un logiciel, quelqu'il soit, ne peut prendre les informations nécessaires pour générer le nouveau nombre que dans ce qu'il a déjà calculé (si on lui donne accès à une base de données externe, autant enregistrer une vraie séquence aléatoire produite par un phénomène physique quelconque...). C'est donc bien mal parti pour faire du hasard tout ça.

On appelle traditionnellement une telle suite une séquence pseudo-aléatoire (Pseudo Random Sequence for Jack AllGood); pseudo pour bien signifier que de toute façon, ce n'est qu'un pis-aller, on fait comme on peut.

Mais qu'est-ce qu'on attend en fait réellement d'un générateur de nombres aléatoires ?
La réponse à cette question conditionne la méthode à choisir, car dans ce domaine comme ailleurs, tout est affaire de compromis. En général on privilégie quatre critères :

On dit qu'un générateur pseudo-aléatoire est acceptable s'il a passé avec succès toute une batterie de tests de statistiques généraux (ces tests dépassent très largement la vérification simple que le générateur ne se « coince » pas; au passage ce n'est pas forcément évident à vérifier comme on le verra plus loin).

On va passer en revue quelques-uns des algorithmes plus couramment utilisés.

 

I. GENERATEURS CONGRUENTIELS LINEAIRES

(barbare comme nom n'est-ce pas ?)

I.1 C'est quoi ?

Il s'agit de l'algorithme le plus utilisé pour produire des nombres aléatoires depuis qu'il a été inventé en 1948 par D.H Lehmer . C'est la suite :

xn+1 = ( a * xn + b) mod c

Il est difficile de croire qu'une formule si simple puisse produire des nombres suffisamment au hasard et pourtant ....

Si on désire produire toujours la même séquence (ce qui est pratique à des fins de tests), on rentre toujours la même valeur de x0 (la graine du générateur).
Si on préfère plutôt que la séquence soit toujours différente, on initialise x0 avec une grandeur toujours différente, l'heure système par exemple (ce que fait la fonction randomize() du C).

Dans tous les cas, les nombres de la suite sont compris entre 0 et c-1.

 

I.2 Et ça marche bien ?

Imaginons à titre d'exemple qu'on choisisse a,b et c tels que :

xn+1 = ( 25 * xn + 16) mod 256

Pour x0 = 12, les nombres produits sont : 60 , 236 , 28 , 204 , 252 , 172 ....
Damned, ils sont tous pairs !!

Pour x0 = 11, les nombres produits sont : 35, 123 , 19 , 235 , 3 , 91 , 243 ..

Aille, ils sont tous impair ! ! !

Pour x0 = 10, les nombres produits sont : 10 ,10 ,10 ,10 ....
Bon sang, voilà qu'il est bloqué maintenant ! ! ! !

La formule est simple mais le choix des trois paramètres ne doit pas être fait à la légère.

Pour ne pas retrouver les défauts graves de notre petit exemple, il faut regarder cette suite d'un peu plus près. Si on choisi b non nul , il est toujours possible d'obtenir une période de longueur c (donc pas de risque de bloquage puisqu'on va retrouver tous les nombres entre 0 et c-1 dans la suite). D. Knuth fait la démonstration des critères que doivent remplir a,b et c pour cela :

(on avait pas respecté le premier critère dans notre petit exemple, vous avez vu les dégâts ?)

Si b=0, la situation est plus compliquée parcequ'il reste toujours la possibilité de bloquage (si x0 = 0) donc la période du générateur ne peut être au mieux que c-1.
L'analyse dans le cas général est assez difficile, pour c = 2^n avec n > 3 on peut dire que :

Bon, pour notre second essai, soyons plus prudent et respectons scrupuleusement les critères précédents :

xn+1 = ( 31415821 * xn + 1) mod 100 000 000

(citée par R. Sedgewick)
Alors là, on a fait attention : 31 415 821 - 1 = 31 415 820 est bien un multiple de 2 et de 5 (les deux seuls nombres premiers qui divisent 100 000 000). C'est aussi un multiple de 4 (car 100 000 000 est lui-même un multiple de 4). Et enfin 1 et 100 000 000 sont bien premiers entre-eux. Par conséquent, on est certain qu'il n'y aura pas de problème de coinçage ou de série de nombres tous pairs ou impairs. Tous les nombres entre 0 et 99 999 999 vont sortir et chacun une fois tous les 100 000 000 de tirages.

Regardons ce que ça donne, par exemple pour x0 = 0. Les nombres produits sont:
1 , 31415822, 40519863, 62952524 , 25482205 , 90965306, 70506227, 6817368 , 12779129, 29199910, 45776111, 9252132, 22780373, 20481234, 81203115, ....
Vous ne remarquez rien ?
Regardez bien ...
Le dernier chiffre ...
C'est du hasard ça ? ?

Avoir la certitude que tous les nombres vont sortir ne suffit pas à garantir qu'ils seront suffisament aléatoires !

Dans un article fameux , Park & Miller passent en revue de nombreux générateurs couramment utilisés et qui la plupart du temps présentent des défauts souvent assez graves. Citons en vrac :

Je laisse la conclusion de tout ceci à Robert Sedgewick : « C'est une règle générale, les générateurs de nombres aléatoires sont délicats et doivent être traités avec respect. Il est difficile d'être certain qu'un générateur donné est bon sans investir d'énormes efforts dans des tests statistiques variés...».

 

I.3 Mais si on en a a besoin, on fait comment ?

Dans l'immense majorité des cas, il faut absolument éviter de chercher à concevoir un générateur pseudo-aléatoire soi-même. Les exemples précédents sont les pièges les plus simples auxquels on s'expose en choisissant à la légère les paramètres du générateur de Lehmer. Mais il y en a bien d'autres et de plus vicieux.
Alors, à moins d'être un mathématicien chevronné, il est plus sage d'utiliser les générateurs qui ont fait leurs preuves.

 

I.3.1 Le Standard Minimal

Afin d'éviter que n'importe qui propose un générateur pseudo- aléatoire sans savoir très bien ce qu'il fait (comme dans les exemples dramatiques précédents), Park & Miller ont proposé un générateur standard convenablement testé et dont le code est portable sur toutes les machines. Ils l'ont appelé le Standard Minimal.

Il est défini par :

xn+1 = 16807 * xn mod (2^31- 1)

(cela correspond à un cas où b est nul donc il faut absolument éviter de partir avec x0 = 0).
La période de ce générateur est 2^31 - 2, ce qui est suffisamment long pour la plupart des applications.
On peut lui reprocher d'avoir un terme multiplicatif un peu petit, mais à part pour les cas pointus, ce n'est pas très grave. Enfin, ce générateur a passé avec succès toute une batterie de tests depuis son invention par Lewis et al. en 1969 même si il a quelques difficultés avec le test du Chi2.

Par code portable, Park & Miller signifient que pourvu qu'une machine sache représenter les nombres dans la gamme -2^31 à 2^31- 1, on peut l'implanter sans risque d'overflow (la gestion du débordement de la multiplication est parfois différente d'un processeur à l'autre et se baser sur ses caractéristiques pour faire tourner un programme est un sport de haute volée). Ils ont utilisé un algorithme dû à Schrage pour en réaliser leur programme en Pascal. En C, cela donne ceci :

    /* Implémentation du Standard Minimal de Park and Miller utilisant
            l'algorithme de Schrage.

            Le nombre aléatoire est maintenu dans la variable globale :
                    alea.
            C'est un long compris entre 1 et 2^31-2 inclus. */

    long alea = 1 ;
    long alea_stdmin(void)
    {
    const long a=16807, m=2147483647L, q=127773L, r=2836 ;
    ldiv_t x ;
    long test ;

            x = ldiv(alea , q) ;            /* la division */
            test = a*x.rem - r*x.quot ;     /* les 2 multiplications */
            if (test >= 0)
                    alea = test ;
            else
                    alea = test + m ;
            return alea ;
    }

Ici, les long font 32 bits, par conséquent, on ne déborde à aucun moment du calcul.
La génération d'un nouveau nombre aléatoire nécessite principalement une division et deux multiplications.

Pour montrer que l'arithmétique qu'on apprend à l'école sert à quelque chose, les lignes qui suivent sont une « démonstration » de ce qu'on peut faire pour développer un algorithme efficace. Je l'ai mis parce que je suis plutôt content de moi sur ce coup là. Enfin, je doute fort que personne ne l'ai trouvé avant moi.... Si vous y êtes réfractaire, on se retrouve quelques lignes plus bas.

Le programmeur en assembleur sur une machine peu puissante sera certainement désolé de constater que cette implémentation nécessite une division par un nombre de 17 bits (car q = 127 773). Le 68000 par exemple ne sait effectuer que des divisions 32 bits par 16 bits. En plus, la division est extrêmement lente sur ce processeur. Et puis le calcul de la variable « test » nécessite de savoir réaliser la multiplication par un nombre de 17 bits. C'est également généralement génant.

Heureusement, on peut réaliser un code aussi portable sans division et sans multiplication de nombres de 17 bits :
L'idée est que 16807*xn \ (2^31 - 1) = 16807*xn \ 2^31 à 1 près. Or comme la division entière par 2^31 peut être réalisée par un décalage plutôt que par une division, il est très rentable de diviser par 2^31 puis de faire ensuite l'ajustement dans le cas (rare en plus) où on se trompe de 1 sur le quotient :

on cherche xn+1 tel que : 16807 * xn = (2^31- 1)*q + xn+1

Pour ce faire, on évite la division et on calcule r' et q' tels que 16807 * xn = 2^31*q' + r'

Par conséquent, 2^31*(q-q') = q + r' - xn+1

Il reste un dernier détail, q'+r' ne peut pas être strictement supérieur à 2^32 car q' <16807 et r'<2^31. Donc il suffit de regarder le bit n°32 de q'+r' pour savoir qu'on s'est trompé de 1 sur le quotient. Enfin, la soustraction de 2^31 se fait simplement en mettant ce dernier bit à 0.

Le codage en C qui découle des considérations précédentes est le suivant :

    long alea_var = 1;
    long alea_stdmin(void)
    {
    long haut , bas , inter , q;

            bas = 16807 * (alea_var & 0xffffL);
            haut = 16807 * (alea_var >> 16);
            inter = (bas >> 16) + haut;
            bas = ((bas & 0xffff) | ((inter & 0x7fff)<<16)) +
                                     (inter >> 15);
            if ((bas & 0x80000000L) != 0)
                    bas = (bas + 1) & 0x7fffffffL;

            alea_var = bas;
            return bas;
    }

 

I.3.2 Il y en a d'autres ?

D. Knuth recense dans « Seminumerical Algorithms » un certain nombre de générateurs congruentiels linéaires de qualité. Si il n'est pas possible d'utiliser le Standard Minimal, on pourra essayer par exemple :

xn+1 = ( 137 * xn + 187) mod 2^8
(une suite cyclique de 256 nombres, est-ce encore du hasard ?)

xn+1 = ( 1664525 * xn + 1013904223) mod 2^32
(générateur de Knuth & Lewis)

xn+1 = 69069 * xn mod 2^32
(générateur de Marsaglia au multiplicande remarquable, utilisé sur le VAX)

xn+1 = ( 31167285 * xn + 1) mod 2^48
(générateur de Lavaux & Jenssens)

xn+1 = ( 6 364 136 223 846 793 005 * xn + 1) mod 2^64
(générateur de Haynes)

Dans tous ces cas, comme le modulant est une puissance de 2, les bits de poids faible des nombres produits ne sont pas eux très « aléatoires ». Il est par conséquent préférable de ne conserver que les poids forts (16 bits par exemple), comme le fait le générateur de Borland C++ par exemple :

xn+1 = (22 695 477 * xn + 1) mod 2^32
et on sort xn+1 >> 16.

Ce faisant, la période du générateur reste la même mais on ne va produire que des nombres compris entre 0 et 65535, chacun se retrouvant 65536 fois dans la suite.

 

II. D'AUTRES GENERATEURS

Les générateurs de Lehmer que nous venons de voir nécessitent de calculer une multiplication (et parfois une division pour faire le modulo). Il est possible de produire des nombres aléatoires sans ces opérations pénalisantes en temps de calcul.

 

II.1 Générateur additif

Mitchell et Moore ont proposé le générateur suivant :

xn = ( xn-24 + xn-55 ) mod m

Si m est une puissance de 2, ce ne générateur ne nécessite qu'une addition.
Il a en plus l'avantage d'avoir une période extrêmement longue : un multiple 2^55 - 1. De ce fait, il n'est pas possible avec les moyens de calculs actuels de faire tourner le générateur de manière à voir une période complète. Donc la vérification pratique que ce générateur ne se « coince » pas est impossible. Mais par bonheur, on peut le démontrer.

Malheureusement, son étude théorique est particulièrement ardue, c'est pourquoi on dispose de peu de preuves théoriques de ses performances statistiques.
Cependant, tous les tests effectués depuis son invention en 1958 ont toujours été concluants.

Pour l'implanter dans un composant programmable, il est plus pratique d'utiliser un générateur très semblable dû à Lewis et Payne :

xn = ( xn-24 XOR xn-55 ) mod m

Si m est de la forme 2^n, on se retrouve en fait avec n polynômes générateurs en parallèle chacun avec son point de départ différent. De ce fait, il faut faire attention à ce que le même bit ne se trouve pas à 0 dans les 55 nombres lors de l'initialisation.

D. Knuth propose lui le générateur :

xn = ( xn-24 - xn-55 ) mod m

En entiers, il ne présente pas d'interêt par rapport au générateur de Mitchell & Moore, mais il peut être implémenté avec des flottants sans modifications.
Cela évite la division de normalisation dans ce cas (on cherche à produire des nombres entre 0.0 et 1.0).

Les nombres 24 et 55 ne sont pas choisis au hasard et on serait très mal intentionné de les changer inconsidérément. Par exemple la suite xn = ( xn-5 + xn-55) mod 2^16 a une période de 1533.

La programmation de tels algorithmes est simple comme bonjour. Il suffit d'utiliser une liste cyclique d'au moins 55 éléments. Pour éviter de tester l'arrivée des pointeurs à 55, la routine C suivante travaille sur une liste de 64 éléments. L'incrémentation modulo 64 est alors faite simplement avec un masquage. Il n'est d'ailleurs pas du tout évident que ce soit plus performant qu'un test, mais bon ...
Pour initialiser le tableau, on pourra utiliser un algorithme plus lent (le Standard Minimal par exemple), en prenant bien garde à ce que les 55 nombres ne soient pas tous pairs, sans quoi le bit de poids faible restera à 0 (la probabilité que cela se produise est très faible mais on ne sait jamais, un petit test ne coute pas bien cher).

    unsigned long alea_Mitchell_Moore(void)
    {
    unsigned long x;
         /* On utilise le fait que le cast implicite en long du
         résultat l'addition réalise le modulo 2^32 désiré */
         alea_tabMM[alea_ptMM] = alea_tabMM[alea_pt24MM] +
                                 alea_tabMM[alea_pt55MM];
         x = alea_tabMM[alea_ptMM];
         alea_ptMM   = (alea_ptMM + 1) & 63;
         alea_pt24MM = (alea_pt24MM + 1) & 63;
         alea_pt55MM = (alea_pt55MM + 1) & 63;

         return x;
    }

 

II.2 Améliorer un générateur pseudo-aléatoire

Pour remplir complètement une page écran 320x200 de points aléatoires, j'avais choisi la suite sn+1 = (5 * sn + 1) mod 65536

en rejetant les valeurs obtenues supérieures à 64000 et en calculant les coordonnées des points x = sn+1 \ 320 et y = sn+1 mod 320.

(en réalité, le calcul des coordonnées n'est pas nécessaire : sn est directement l'adresse du pixel dans la mémoire écran sur ma machine).

Tout l'écran se remplissait bien car les critères donnés plus haut sont respectés et ça allait très vite parcequ'une multiplication par 5 se code avec une addition et un décalage. Mais, au résultat ça n'avait pas l'air vraiment « très aléatoire ». On croyait voir des motifs et des petits segments verticaux, rien de bien grave malgré tout.

Il y avait pourtant une méthode simple pour se sortir de ce mauvais pas : le mélange des nombres produits.
On part d'une séquence pseudo-aléatoire xn. Ces nombres sont inférieurs à m (un générateur de Lehmer par exemple). On s'en sert pour remplir un tableau tab[] de k éléments de 0 à k-1.

La procédure pour générer le nombre pseudo-aléatoire suivant est alors :

  1. calcul du xn par l'algorithme de son choix.
  2. calcul de j = (k * xn) \ m (si k et m sont des puissances de 2, tout ça n'est que des décalages)
  3. le nombre produit sera tab[j], et on rempli tab[j] avec le xn qu 'on a utilisé pour calculer j.

Cet algorithme en apparence tout bête est de C. Bays et S.D.Durham. La période du générateur résultant de ce mélange est la même que celle d'origine.

L'amélioration du caractère aléatoire de la séquence en utilisant un tel mélange est considérable, et ce pour un coût de calcul ridiculement faible.
Par exemple il fait disparaitre les défaillances du Standard Minimal de Park & Miller lors du test du Chi2.
Ahhh si j'avais su quand je remplissais mon écran....

Cette idée de mélanger les nombres donne naissance à toute une classe d'excellents générateurs pseudo-aléatoire dont les périodes sont souvent extrêment longues.

 

II.3 Produire des bits au hasard

Il est parfois nécessaire de produire des bits aléatoires isolés. C'est dommage de devoir utiliser un générateur de mots pour n'en utiliser qu'un seul bit non ? D'autant plus que l'indépendance statistique des bits d'un mot d'un générateur pseudo-aléatoire n'est pas toujours évidente (et à l'inverse, il est très déconseillé de constituer des mots à partir de bits isolés aléatoires).

La méthode courante consiste à utiliser les polynômes primitifs modulo 2 (ça aussi c'est un nom barbare n'est-ce pas ?). Par exemple, le polynôme générateur x^16+x^12+x^5+x^0 (le polynôme CCITT des CRC 16 bits) peut être réalisé avec le montage suivant :

                           Registre à décalage de 16 bits
                    +--+--+--+--+--+--+-+-+-+-+-+-+-+-+-+-+
                    |15|14|13|12|11|10|9|8|7|6|5|4|3|2|1|0|<--+
                    +--+--+--+--+--+--+-+-+-+-+-+-+-+-+-+-+   |
                      |          |               |            |
                      |          |  +-----+      |  +-----+   |
                      |          +->| XOR |      +->| XOR |---+
                      +------------>|     |-------->|     |
                                    +-----+         +-----+

Ici, on réalise le ou exclusif des bits 15,11 et 4. Puis on décale le mot de 1 bit à gauche. Le nouveau bit 0 sera le résultat du ou exclusif. Le choix de faire un décalage à gauche plutôt qu'à droite est tout à fait arbitraire.

La séquence produite (on parle souvent de PRBS ou Pseudo-Random Bit Sequence) a une période de 2^16-1 bits et il faut prendre garde à ne pas initialiser le registre à décalage avec la valeur 0 ,sans quoi ça reste « coincé » à 0. Le bit aléatoire recherché peut être n'importe lequel des 16 bits du registre.

D'une manière générale, la période d'un générateur de ce type est 2^n-1 où n est le plus haut degré du polynôme.

La programmation d'un polynôme générateur n'est pas très pratique. Quand on a le choix, il est préférable pour des raisons de rapidité d'utiliser ceux qui n'ont que trois termes (un seul XOR) comme x^15+x^1+x^0 , c'est celui dont le code C est le suivant :

    short alea_polyreg = 1;
    char alea_bit(void)
    {
    int x;
            x = alea_polyreg & 0x4001;
            alea_polyreg = (alea_polyreg << 1) & 0x7fff;
            if ((x == 0) || (x == 0x4001)) {
                    alea_polyreg |= 0;
                    return 0;
            }
            else {
                    alea_polyreg |= 1;
                    return 1;
            }
    }

 

III INITIALISER UN GENERATEUR PSEUDO-ALEATOIRE

Ce n'est pas tout de savoir produire des nombres qui paraissent sans suite logique, il faut encore initialiser convenablement le générateur.

En C, on dispose de deux fonctions à cette fin : srand() et randomize(). Et pour chacune d'elles, on peut trouver à y redire (il n'y a pas de « yaka » dans le domaine des générateurs de nombres pseudo-aléatoires ... et comme dans d'autres domaines du calcul scientifique, les concepteurs du C ont été ... disons... moyens).

 

III.1 srand(unsigned)

Cette fonction initialise le générateur avec une valeur qui lui est fournie de manière à pouvoir générer toujours la même séquence, ce qui est bien pratique lorsque l'on teste une routine. Le défaut est écrit dans le prototype même de la fonction : on doit lui passer un unsigned int.
Dans de nombreuses implémentations du C, ce type correspond à un nombre de 16 bits. Cela signifie, qu'on ne peut produire que 65536 séquences différentes possibles, et c'est fort peu ....
On serait par exemple bien mal intentionné d'utiliser srand() pour initialiser une séquence utilisée en cryptage, car même si la période du générateur est très longue, on ne fait qu'une bouchée de passer en revue toutes les possibilités....

C'est d'autant plus regrettable qu'il n'y a pas spécialement de problème à initialiser le générateur au moins avec un long (sur 32 bits) dans le cas du C ANSI car c'est avec cette taille de mots qu'il travaille.

 

III.2 randomize(void)

Cette fonction initialise le générateur avec l'heure du système. Elle est donc censée produire de nouvelles séquences à chaque appel. Par exemple, le help du Borland C++ 3.1 indique :
« randomize initialise le générateur interne de nombres aléatoires avec une valeur aléatoire ». Et l'exemple d'utilisation fourni est le suivant :

    #include <stdlib.h>
    #include <stdio.h>
    #include <time.h>

    int main(void)
    {
    int i ;
            randomize() ;
            printf("Dix nombres aléatoires entre 0 et 99\n\n") ;
            for (i=0 ; i<10 ; i++)
                    printf("%d\n",rand()%100);
            return 0 ;
    }

Vu comme ça, ça a l'air pas mal, effectivement ce court programme affiche bien une séquence de dix nombres aléatoires. Et si on change la position de l'appel à randomize() ?

    #include <stdlib.h>
    #include <stdio.h>
    #include <time.h>

    int main(void)
    {
    int i ;
            printf("Dix nombres aléatoires entre 0 et 99 ? ? ?\n\n") ;
            for (i=0 ; i<10 ; i++) {
                    randomize() ;
                    printf("%d\n",rand()%100);
            }
            return 0 ;
    }

Là, catastrophe ! ! dix fois le même nombre sur mon écran ! ! !
Où est l'initialisation aléatoire de ma séquence ? ? ? ?
Et non, ça ne marche pas.... en fouillant dans les sources fournis (il n'y a pas besoin d'aller bien loin, visiblement c'est pareil pour la plupart des compilateurs C), la fonction randomize() est en fait une macro :
sand((unsigned) time(NULL)) ;
où time() renvoie l'heure système en secondes (! ! !) depuis le 1er janvier 1970.

On combine donc là deux maladresses :

De ce fait, plusieurs appels à randomize() dans un même programme sont à proscrire absolument à moins qu'on soit certain que le temps entre deux appels consécutifs est supérieur à une seconde (pas vraiment portable comme recommandation non ? ?).

III.3 Faire « propre »

Pour initialiser le générateur de façon à peu près aléatoire en évitant les pièges décrits précédement, on peut procéder de la façon suivante :

La routine en C sur PC suivante réalise cela. Elle n'est malheureusement pas portable. Pour l'adapter sur d'autres machines il faut rechercher un ou plusieurs compteurs rapides de manière à constituer la graine de taille convenable. Il serait possible de se contenter de la fonction time() car le nombre de secondes depuis le 01/01/1970 tient sur 30 bits mais j'ai voulu faire le puriste ici.

    alea_init = 0;
    void alea_randomize(void)
    {
    unsigned long lb , hb;

            while (alea_init == 0) {
                    /* les 16 bits de poids faible sont donnés par le
                    contenu du timer 0, les 16 bits de poids fort par
                    le temps en seconde depuis le 1/1/1970 (ou 1968 ça
                    dépend des compilateurs) */

                    /* On ne sera pas interrompu pendant la lecture du
                    timer */
                    asm pushf
                    asm cli
                    outp(0x43,0);
                    /* on laisse un peu de temps au 8254 pour figer le
                    compteur */
                    alea_init = (unsigned long)time(NULL);
                    /* puis on le lit */
                    lb = (unsigned long)inp(0x40);
                    hb = (unsigned long)inp(0x40);
                    asm popf
                    alea_init = ((alea_init << 16) | (hb << 8) | lb))
                                                     & 0x7fffffff;
            }

            alea_var = alea_init;
            /* on bidouille la graine obtenue */
            alea_init = alea_stdmin();
            /* Puis on appelle notre propre routine srand() avec sa
            graine sur 32 bits */
            alea_srand(~alea_init);
    }

 

IV UTILISATION D'UN GENERATEUR PSEUDO-ALEATOIRE

IV.1 Produire un nombre entier entre deux bornes

Les algorithmes qu'on a vu produisent des nombres pseudo-aléatoires entre deux bornes très éloignées (pour le Standard Minimal, l'intervalle est [1,2^31-2]). Mais dans la pratique, on a plutôt besoin d'un intervalle réduit. Celui qui veut simuler un lancer de dé doit se limiter à un intervalle [1,6] par exemple. Comme généralement les bits de poids faibles des générateurs ne sont pas particulièrement aléatoires, la ligne C suivante effectuant un tirage entre 0 et n est vivement déconseillée :

j = rand() % n;

Si n est grand devant la valeur maximale du générateur, elle peut en plus introduire un biais considérable. Imaginons par exemple que notre générateur de base produise des nombres avec une distribution uniforme entre 0 et 2^32-1 inclus et qu'on demande un n = 2^32-2.

Conclusion, 0 et 1 ont deux fois plus de chances de sortir que les autres nombres. C'est dommage, même si cet exemple est un cas extrême malgré tout assez peu génant il faut le reconnaître.

Pour vous amuser, regardez ce qui se passe pour n = 6 si la période du générateur est de 8 (il sort tous les nombres entre 0 et 7). Les dés sont franchement pipés ....

Il est donc préférable de procéder comme suit :

j = (int)(n*rand()/(RAND_MAX + 1));

(RAND_MAX est la valeur maximale que peut produire le générateur rand(), la valeur minimale étant 0.)

Si RAND_MAX + 1 est une puissance de 2, il est possible de réduire le calcul à des multiplications entières.

Heureuseusement il est possible de faire la même chose sans multiplication ni division et en restant avec des entiers. C'est l'objet du petit programme C suivant :

    /*      Tirage d'un nombre entre 1 et 6 sans division.
            Le générateur utilisé est alea_Mitchell_Moore(). Le tableau
            tabMM[] aura été convenablement initialisé au préalable.
            On utilise une méthode de réjection */
    char lance_de(void)
    {
    char x = 8;
            while (x>5)
                    x = (char)(alea_Mitchell_Moore() >> 29);

            return x+1;
    }

Ici, notre générateur de base produit des nombres entre 0 et 2^32-1. On ne garde que les trois bits de poids fort grace à un décalage, cela constitue donc un tirage entre 0 et 7. Puis on élimine tous les tirages qui sortent de l'intervalle désiré [0,5].

Dans le pire des cas, la probabilité que l'on fasse un tirage « pour rien » par cette méthode de réjection est de 0.5 (donc 2 tirages en moyenne pour produire un nombre dans l'intervalle désiré). Ici, pour le lancé de dé, elle est de 0.25. L'inconvénient le plus pénalisant est que l'on ne peut pas prévoir exactement le temps nécessaire à la génération du nombre suivant. Mais est-ce si grave ?

 

IV.2 Mélanger un ensemble

On vient de choisir un nombre au hasard parmi n. Une autre application importante des générateurs pseudo-aléatoire est de mélanger n nombres. C'est le cas par exemple lorsqu'on désire simuler un jeu de cartes : il faut battre le jeu avant de le distribuer.

Et là, notre ordinateur a une petite difficulté :

Si la période du générateur est suffisament longue, ce n'est pas si grave car de toutes façons, il faudra beaucoup de temps (voire trop pour notre espérance de vie) pour passer en revue assez d'arrangements qui permettraient de se rendre compte du défaut. Donc pas de danger qu'un des joueurs sorte son colt en hurlant que le croupier (le pauvre ordinateur) triche.
Enfin, cela plaide quant même en faveur d'un générateur à période la plus longue possible.

Ceci dit, comment faire pour mélanger au mieux notre ensemble ? On part d'un tableau tab[] de n nombres à mélanger et on utilise une fonction random(x) qui renvoie un entier pseudo-aléatoire inférieur à x (ce que nous avons fait au paragraphe précédent).
Puis :

  1. j = n
  2. echanger tab[j] et tab[random(j)]
  3. j = j - 1
  4. retourner en 2 tant que j > 0

Voilà, c'est l'algorithme de Moses & Oakford.

Quittons un peu le domaine des entiers...

Si la notion de densité de probabilité ou de fonction de répartition vous est étrangère, c'est que vous n'avez vraisemblablement pas besoin de lire les lignes qui suivent. Alors on se retrouve à la fin de l'article.... A tout de suite.

 

IV.3 Distribution uniforme en flottant

Dans les applications scientifiques, on a besoin la plupart du temps d'un générateur produisant un nombre entre 0 et 1 avec une distribution uniforme. C'est exactement ce que fait la ligne C suivante :

random = (float)rand()/(RAND_MAX + 1);

Attention en double précision car si le RAND_MAX+1 vaut 2^32, les 21 bits de poids faible de la mantisse seront tous à 0.

 

IV.4 Distribution normale

Pour obtenir une distribution normale, il y a deux méthodes principales :

 

IV.4.1 l'algorithme de Box-Muller:

La distribution normale de moyenne µ et d'écart-type s est obtenue en appliquant la formule :

norm(µ,s) = µ + s*sqrt(-2*ln(random1))*cos(2*PI*random2)

où random1 et random2 suivent une distribution uniforme dans l'intervalle ]0,1]. C'est cette implémentation que Mathcad 3.1 recommande.

On peut diviser le temps de calcul quasiment par 2 en réutilisant les calculs déjà effectués, comme dans la procédure suivante :

    /*      Générateur normal. On utilise l'algorithme de
            Box-Muller-Marsaglia tel que décrit par D. Knuth p.117.
            Le générateur uniforme servant de base est le
            Mitchell-Moore */
    double alea_sav;
    double alea_var;
    char alea_etat = 0;
    double alea_normal_BMM(double moyenne , double ecart_type)
    {
    double v1 , v2;
         if (alea_etat == 0) {
            do {
                    v1=2.*(double)alea_Mitchell_Moore()/4294967296.-1.;
                    v2=2.*(double)alea_Mitchell_Moore()/4294967296.-1.;
                    alea_sav = v1*v1 + v2*v2;
            } while ((alea_sav >= 1.) || (alea_sav==0.));
            alea_sav = sqrt(-2.*log(alea_sav)/alea_sav);
            alea_var = v2;
        }
        else
            v1 = alea_var;
        alea_etat ^= 1;
        return (alea_sav*v1*ecart_type + moyenne);
    }

 

IV.4.2 En utilisant le limite centrale :

Si les calculs flottants nécessaires dans l'algorithme de Box-Muller sont génants (lorsqu'on on ne dispose pas d'un coprocesseur arithmétique par exemple), une approximation utilisant le théorème de la limite centrale est digne d'être examinée : La variance des nombres issus d'un générateur uniforme compris entre -0.5 et 0.5 vaut 1/12. En sommant 12 nombres de ce type, on obtient une distribution qui ressemble fort à une gaussienne avec une variance de 1 et de moyenne nulle. C'est ce que fait la procédure suivante :

    /*      Générateur normal. On utilise une méthode approximative
            reposant sur le théorème de la limite centrale. A utiliser
            avec précautions !!
            Le générateur uniforme servant de base est le
            Mitchell-Moore */
    double alea_normal_approx(double moyenne , double ecart_type)
    {
    int i;
    unsigned long total = 0;
            for(i=0 ; i<12 ; i++)
                    total += (alea_Mitchell_Moore() >> 4);

            return (moyenne+ecart_type*((double)total/268435456.-6.));
    }

Il faut malgré tout prendre garde au fait que cette méthode est une approximation d'autant plus mauvaise qu'on s'éloigne de la moyenne: la probabilité de produire un nombre à plus de 6Õ est rigoureusement nulle (contre 6E-8 pour Box-Muller). Multiplier le nombre de nombres à sommer améliore la qualité de l'approximation, mais naturellement au détriment de la vitesse. Le choix de 12 nombres n'est en fait qu'un joli compromis.

 

IV.5 Distribution exponentielle

Une autre distribution très courante (en gestion des files d'attentes par exemple) est la distribution exponentielle. Elle peut être obtenue en appliquant la formule :

exp(µ) = -µ*ln(random)

où µ est la moyenne désirée et random suit une distribution uniforme dans l'intervalle ]0,1].

 

V LA FIN

Et après toutes ces considérations, le source joint HASARD.C présente les différentes procédures qui me semblent les meilleurs compromis sur mon PC 486DX33.
Mis à part la fonction randomize() qui utilise un registre du timer interne du PC, elles sont écrites de manière à être portables. En retirant la procédure main() qui teste les différentes fonctions, cela peut constituer une librairie de base.

Voilà, on a effleuré le sujet des générateurs pseudo-aléatoires. Il y en a bien d'autres, par exemple ceux basés sur la conjecture de Syracuse, l'exponentielle modulaire ou encore les automates cellulaires; mais on a vu les plus courants et les plus faciles à programmer.

Il reste beaucoup de choses à dire. Pour ceux qui veulent plus de détails, il faut absolument lire :

Knuth, D.E. « The Art of Computer Programming : Seminumerical Algorithms », Addison Wesley,1981 (INDISPENSABLE, c'est une bible. Il faut tout lire, même les exercices et leurs corrections)

Park S.K., Miller K.W., « Random Number Generators: Good Ones Are Hard To Find » Comm. of the ACM, Vol 31, 10, Oct.1988, 1192-1201

Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P, « Numerical Recipes in C, Second Edition », Cambridge University Press, 1992. (un autre must à lire absolument !!!, disponible aussi en Pascal et FORTRAN)

Naudin P., Quitté C., « Algorithmique Algébrique », Masson, 1992. (c'est en français, pas très pratique (exemples en ADA !!!) et lisible par les amateurs avertis des groupes cycliques).

 

Après la conclusion :

J'ai fait un effort, mais il est bien possible que j'ai commis des erreurs dans cet article. Si c'est le cas, n'hésitez pas à me le signaler. C'est pas que ça fasse plaisir; mais un erratum, ça vaut mieux que de laisser trainer des bêtises.

Bonne programmation .....

 

Pierre Larbier, novembre 1994

<< Retour à la liste <<