User Tools

Site Tools


start:oldies:random

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 :

  1. la vitesse : il faut que le calcul du nombre pseudo-aléatoire suivant soit rapide. Il n'est pas rare de devoir générer des millions de nombres et on ne peut pas passer sa vie à attendre.
  2. la simplicité : ça c'est pour le confort du programmeur, une méthode très compliquée sera très difficile à programmer et surtout à tester. En général, on dit même en forme de boutade que plus l'algorithme est compliqué moins la séquence sera aléatoire.
  3. la méthode ne doit pas souffrir de faille grave : on en reparlera plus loin mais l'histoire des générateurs pseudo-aléatoire est pleine d'algorithmes qui se « coincent » lorsqu'ils arrivent sur un nombre particulier. Ce qui est très très gênant car c'est souvent source de plantages qui semblent eux aléatoires et qui sont par conséquent sont très difficiles à mettre en évidence.
  4. les nombres produits ne doivent pas faire apparaître de suite logique, quelle que soit la façon de les regarder. Ce critère est le plus difficile à quantifier car il dépend fortement de l'application. Par exemple, imaginons qu'on veuille utiliser un générateur pour simuler un lancé de dés. On va vouloir absolument que tous les nombres de 1 à 6 aient exactement les mêmes chances d'apparaître. Mais si c'est pour jouer au 421, on va vouloir en plus que toutes les suites de deux et trois nombres consécutifs aient aussi exactement les mêmes chances d'arriver. Si la séquence 421 a moins de chances de se produire que 666, certains ne vont pas hésiter à prétendre que les dés sont pipés, et c'est le début des ennuis…

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.

Générateurs congruentiels linéaires

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

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 :

$x_{n+1} = (a * x_n + 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.

Et ça marche bien ?

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

$x_{n+1} = (25 * x_n + 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 :

  • b et c doivent être premiers entre eux
  • a-1 doit être un multiple de p, pour tout p nombre premier diviseur de c
  • a-1 doit être un multiple de 4 si c est un multiple de 4.
  • si c est une puissance de 2, le bit de poids faible des nombres produits vaut altérnativement 0 et 1 (ce n'est d'ailleurs pas le seul cas où cela se produit).

(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 :

  • a mod 8 doit être égal à 3 ou 5
  • la période du générateur est 2^(n-2) et tous les bits de poids faible des nombres produits sont identiques (il faut donc décaler le nombre produit de 1 vers la droite).

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

$x_{n+1} = (31415821 * x_n + 1) \mod 100000000$

(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:

  • RANDU défini par $x_{n+1} = 65539 * x_n \mod 2^{31}$. Implanté sur les IBM n+1 n System/370 et appellé par des noms d'oiseaux par tout ceux qui ont eu besoin de l'utiliser sérieusement.
  • Le générateur de Turbo Pascal défini par $x_{n+1} = (129 * x_n + 907633385) \mod 2^{32}$. On sait dans 1 cas sur 129 que $x_{n+1} > x_n$, ce qui est un biais énorme pour une application scientifique. En plus comme 129 est de la forme $2^n + 1$, les nombres produits ne sont pas si « aléatoires que ça » [cf. Knuth p.23].
  • Le générateur d'UNIX (dont le comité ANSI C a heureusement recommandé de ne conserver que les 16 bits de poids fort): $x_{n+1} = (1103515245 * x_n + 12345) \mod 2^{31}$. Même la documentation du Berkeley 4.2 admet que les bits de poids faible des nombres produits ne sont pas vraiment aléatoires.

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…».

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.

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 :

$x_{n+1} = (16807 * x_n \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 $x_{n+1}$ tel que: $16807 * x_n = (2^{31} - 1) * q + x_{n+1}$

Pour ce faire, on évite la division et on calcule r' et q' tels que $16807 * x_n = 2^{31} * q' + r'$. Par conséquent, $2^{31}*(q-q') = q + r' - x$.

  • Si on ne s'est pas trompé sur le quotient, $q=q'$ donc le reste recherché vaut $q'+r'$.
  • Si on s'est trompé, c'est forcément de 1 par excès, et le reste dans ce cas vaut $q'+r'-1- 2^31$. Or comme le nombre attendu ne peut pas être nul (sinon le générateur se coince, voir plus haut) ni négatif (le reste d'une division euclidienne est forcément positif), cela signifie que le calcul de q'+r' a donné un résultat supérieur ou égal à $2^31$. Le voilà notre critère pour savoir si on s'est trompé sur le quotient.

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;
}

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:

$x_{n+1} = ( 137 * x_n + 187) \mod 2^8$ une suite cyclique de 256 nombres, est-ce encore du hasard ?
$x_{n+1} = ( 1664525 * x_n + 1013904223) \mod 2^{32}$ générateur de Knuth & Lewis
$x_{n+1} = 69069 * x \mod 2^{32}$ générateur de Marsaglia au multiplicande remarquable, utilisé sur le VAX
$x_{n+1} = ( 31167285 * x_n + 1) \mod 2^{48}$ générateur de Lavaux & Jenssens
$x_{n+1} = ( 6364136223846793005 * x_n + 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 :

$x_{n+1} = (22 695 477 * x_n + 1) \mod 2^{32}$

et on sort $x_{n+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.

D'autres générateurs

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.

Générateur additif

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

┌───────────────────────────────┐
│  x  = ( x     + x     ) mod m │
│   n      n-24    n-55         │
└───────────────────────────────┘

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 :

┌──────────────────────────────────┐
│  x  = ( x     XOR  x     ) mod m │
│   n      n-24       n-55         │
└──────────────────────────────────┘

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 :

┌───────────────────────────────┐
│  x  = ( x     - x     ) mod m │
│   n      n-24    n-55         │
└───────────────────────────────┘

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 x  = ( x    + x     ) mod 2^16 a une période de 1533.
                n      n-5    n-55

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;
    }

Améliorer un générateur pseudo-aléatoire

Pour remplir complètement une page écran 320×200 de points aléatoires, j'avais choisi la suite s = (5 * s + 1) mod 65536 n+1 n en rejetant les valeurs obtenues supérieures à 64000 et en calculant les coordonnées des points x = s \ 320 et n+1 y = s mod 320. n+1 (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.

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 1) { alea_polyreg |= 0; return 0; } else { alea_polyreg |= 1; return 1; } }

Initialiser un générateur pseudo-aléatoire

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).

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.

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 : sand2) ; où time() renvoie l'heure système en secondes (! ! !) depuis le 1er janvier 1970.

On combine donc là deux maladresses : * srand() qui ne sait initialiser que 65536 séquences différentes * l'appel à time() qui fait qu'on ne sait initialiser une nouvelle séquence que toutes les secondes.

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 ? ?).

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 : * si la fonction randomize() a déjà été appelée, utiliser comme graine du générateur le dernier nombre aléatoire produit. * si la fonction randomize() n'a pas encore été appelée, utiliser les compteurs « hard » de l'ordinateur pour constituer une graine de taille convenable (32 bits généralement)

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 = 3) & 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); }

UTILISATION D'UN GENERATEUR PSEUDO-ALEATOIRE

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. si le nombre produit par rand() est égal à 0 ou 2^32-2, au résultat on aura j = 0. si le nombre produit par rand() est égal à 1 ou 2^32-1, au résultat on aura j = 1. dans tous les autres cas, j = rand(). 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 ?

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é : * Mélanger n éléments d'un ensemble revient à choisir un arrangement parmi les n! qui sont possibles. Si il s'agit d'un jeu de 52 cartes, il y a 80.7E66 arrangements possibles. * Or avec nos générateurs pseudo-aléatoires, on ne sait produire qu'un nombre multiple de la période de séquences différentes. Et c'est souvent petit devant le nombre d'arrangements possible. Par exemple le Standard Minimal avec sa période de 2^31-2 (2.14E9) aura bien du mal avec les 2.6E35 arrangements possibles d'un jeu de 32 cartes. * De ce fait, on ne saura pas produire n'importe quel arrangement. Il y en aura qui seront impossibles. Lesquels ? Ce n'est pas facile à déterminer a priori. Si les mains à 4 as ou avec une quinte flush sont assez fréquentes, il ne fait pas de doute qu'elles sortiront. Mais que les 4 joueurs de belote se retrouvent avec chacun un carré lors d'une distribution, c'est déjà moins sûr.

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.

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.

Distribution normale

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

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 4); 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); }

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*5); }

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.

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].

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

1)
x == 0) || (x == 0x4001
2)
unsigned) time(NULL
3)
alea_init « 16) | (hb « 8) | lb
4)
alea_sav >= 1.) || (alea_sav==0.
5)
double)total/268435456.-6.
This website uses cookies. By using the website, you agree with storing cookies on your computer. Also you acknowledge that you have read and understand our Privacy Policy. If you do not agree leave the website.More information about cookies
start/oldies/random.txt · Last modified: 2017/10/16 14:05 by stephane