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 ....
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.
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.
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...».
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;
}
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.
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;
}
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 :
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.
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;
}
}
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).
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.
#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 ? ?).
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);
}
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 ?
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 :
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.
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.
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);
}
/* 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.
exp(µ) = -µ*ln(random)
où µ est la moyenne désirée et random suit une distribution uniforme dans l'intervalle ]0,1].
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).
Bonne programmation .....
Pierre Larbier, novembre 1994