Générateur de nombres aléatoires

La génération de nombres aléatoires

La physique numérique fait une consommation immodérée de simulations en tous genres. Et parmi ces simulations, beaucoup traitent de phénomènes aléatoires : désintégration nucléaire, arrivée de clients dans une file d'attente, percolation, marche au hasard, etc. Il est donc nécessaire de disposer d'un outil qui permette à un programme de générer des nombres aléatoires. Pas simple lorsqu'on sait qu'un ordinateur et un programme sont des objets en principe déterministes !
Autant le dire tout de suite, il est impossible de créer avec un algorithme une suite de nombres parfaitement aléatoire. Je ne vais pas en développer ici les raisons et je vous renvoie à la littérature spécialisée. Seuls des générateurs physiques peuvent générer du pur aléatoire, ceux basés sur les désintégrations nucléaires ou le bruit thermique par exemple. On se contentera donc de générer des suites pseudo-aléatoires avec un ordinateur.
Une suite de nombres pseudo-aléatoire est une suite périodique, de période aussi grande que possible, dont les nombres, à l'intérieur de chaque période, sont aussi peu corrélés entre eux que possible (i.e. le plus indépendants les uns de autres). Dans la suite du texte, je parlerai par paresse de "suite aléatoire" : il faut comprendre bien sur "suite pseudo-aléatoire".
Deux physiciens, Park et Miller ("Random number generators: good ones are hard to find" in Communications of the ACM 31,1192-201) ont définis en 1988 les trois qualités d'un générateur de nombres aléatoires:

la suite de nombres produite doit être de période très grande. Classiquement, sur une CPU de 32 bits, la période choisie est de 2^31 - 1. Sur une CPU de 64 bits, la période choisie est de 2^63 - 1. Ce n'est pas infini mais c'est pas mal quand même....
la corrélation entre les nombres doit être la plus faible possible. Il existe un moyen pas très compliqué de calculer cette corrélation (voir par exemple le §1.2 de la troisième partie du livre "Physique numérique" de P.Depondt) mais à notre niveau, je préfère un moyen plus ludique: on plot les nombres aléatoires sur un plan xy et on cherche visuellement un motif quelconque qui traduirait un certain ordre dans la suite aléatoire. C'est ce que je fais dans le programme exemple ci-dessous.
le générateur doit être aussi rapide que possible. Dans certaines applications, comme l'intégration par la méthode de Monte-Carlo, les appels au générateur de nombres aléatoires sont massifs. Il ne s'agit donc pas que celui-ci mette une seconde à délivrer un nombre ! Il faut dire que ce point est moins sensible aujourd'hui qu'en 1988 ! Avec une station quadruple core à 3 GHz, on a pas vraiment de problème...

Reste maintenant à trouver un algorithme susceptible de répondre à ces trois critères. Il en existe un certain nombre... Ceux que ça intéresse peuvent consulter par exemple "Numerical recipes" en C ou Fortran, au chapitre 7 dans les deux ouvrages.
Dans le cadre de cette expérimentation, j'utilise un algorithme de génération d'une distribution uniforme simple de mise en oeuvre et très efficace, que je décris au paragraphe suivant. Il s'agit de l'algorithme de Lehmer (ou schéma de congruence linéaire) qui date de 1951 et a été publié pour la première fois dans "Mathematical Models in large scale computing unit in Proceedings of the Second symposium on Large scale digital computing machinery . En passant, allez voir l'histoire de Dick Lehmer, elle est intéressante....

Distribution uniforme de nombres aléatoires et algorithme de Lehmer

Je vous renvoie à cette page de TangenteX pour établir la définition d'une distribution uniforme de variables aléatoires. Je vais rappeler ici quelques unes de ses propriétés qui nous intéresseront pour tester notre algorithme.
Pourquoi avoir choisi une distribution uniforme ? Parce que c'est la plus simple à analyser et à mettre en oeuvre, parce qu'à partir de celle-ci on peut construire d'autres distributions, mais aussi et surtout parce que c'est la plus courante dans les phénomènes physiques.
Nous cherchons donc un algorithme qui génère des nombres aléatoires selon une distribution uniforme dans le domaine [0,1[. Pourquoi ce domaine ? Par simplicité et parce que c'est suffisant. Les propriétés de la distribution uniforme sur ce domaine m'intéressent particulièrement. Sa moyenne est égale à :

\begin{align} <x> = \int_{0}^{1}xdx \end{align}

et sa variance est:

\begin{align} \sigma^2 = \int_{0}^{1}x^2dx - <x>^2 \end{align}

En calculant ces valeurs, on obtient une moyenne de 0.5 et un écart type de 0.289. Ces valeurs sont caractéristiques de la distribution uniforme sur le domaine choisi.
Voyons maintenant le moyen d'obtenir la génération d'une suite pseudo-aléatoire de distribution uniforme. Comme dit plus haut, j'utilise le schéma de congruence linéaire qui est on ne peut plus simple, regardez :

x _i+1 = (ax_i + b) mod c

Ceux qui voudront comprendre le fonctionnement de l'algo peuvent se plonger dans l'article de Lehmer ou dans "Computational physics" de Landau §10. Les autres me feront confiance.
Evidemment, vous aurez compris que le point sensible de l'algorithme tient dans le choix des paramètres a,b et c. Le paramètre c est le paramètre de normalisation, qui permet de restreindre les valeurs au domaine [0,1[. Il vaut donc 2^31 - 1 pour une machine 32 bits (ou 2^63 - 1 pour une 64 bits).
Le choix des paramètres a et b est plus large. J'ai choisi les valeurs les plus communes dans la littérature (celles de Park and Miller): a = 7^5 et b = 0. Vous pouvez modifier ces valeurs pour voir ce qui se passe, en matière de "randomness" sur la distribution uniforme.

Comment vérifier le caractère aléatoire de la distribution uniforme

Je calcule la moyenne et l'écart type de la distribution pour un nombre de tirages donné et je les compare aux valeurs théoriques. Les plus courageux pourront améliorer le programme en calculant la corrélation entre les nombres générés et vérifier qu'elle est aussi petite que possible (voir la référence ci-dessus).
La méthode que je préfère, même si ce n'est pas la plus précise, c'est la méthode visuelle. A partir de la série générée, je trace chaque point sur un plan en choisissant pour coordonnées du point (x_i, x_i+10). Le i+10 est arbitraire, ça pourrait être i+20 ou i+50 sans problème ! Je fais plusieurs tracés en faisant varier le nombre de tirages. Et je cherche dans les figures obtenues des motifs, des zones non homogènes ou des tâches particulières. Si la corrélation entre nombres est très faible, avec un nombre de tirages grand (10^4 ou 10^5), j'arrive à une couverture quasi-uniforme du domaine.

Le programme C et les résultats

J'ai implémenté cet algorithme en C. J'utilise Gnuplot pour le tracé xy de la distribution. Voici le source :


              //**************************************************************************

              // Générateur de nombres aléatoires 

              // Dominique Lefebvre - Octobre 2009 V1.0 

              // TangenteX 

              // 

              // http://www.tangenteX.com"

//**************************************************************************

              #include <cstdkib>

              #include <iostream>

              #include <math.h>

              #include <time.h>

              using namespace std;

              // déclaration des constantes 

              #define NBTIRAGE 10000 v // déclaration des variables globales 

              long seed;v
              //**************************************************************

              // Génération de nombres aléatoires distribués uniformément entre
              0 et 1 

              // Selon la formule de Lehmer et l'algorithme et les paramètres de
              

              // Parker and Miller (1988): 

              // x(i+1) = a*x(i) mod c 

              // Sa période est fixée par c, soit 2^31 - 1 

//**************************************************************************

              double Random(void) 

              { 

                const double a = 16807.0; // 7^5 

                const double c = 2147483647.0; // 2^31 - 1 

                double x,y; 

                x = a*seed;

                seed = (long)(fmod(x,c)); 

                y = seed/c; // normalisation par c 

                return (y); 

              } 

//***************************************************************************

              // Génération de la constante d'initialisation du générateur
              Random 

              // à partir de la date courante selon Anderson (1990) 

//***************************************************************************

              long GenerationSeed(void) 

              { 

              struct tm *temps; 

              time_t maintenant; 

              int t0,t1,t2,t3,t4,t5; 

              long lseed; 

              // récupération de l'heure actuelle 

              time(&maintenant); 

              // enregistrement de l'heure actuelle 

              temps = localtime(&maintenant);

              // conversion en structure tm 

              t0 = temps->tm_sec; 

              t1 = temps->tm_min; 

              t2 = temps->tm_hour; 

              t3 = temps->tm_mday; 

              t4 = temps->tm_mon; 

              t5 = temps->tm_year; 

              // calcul du seed d'après l'algo d'Anderson (seed varie entre 0 et
              2^31-1 

              lseed = t5 + 70*(t4 + 12*(t3 + 31*(t2 + 23*(t1 + 59*t0)))); 

              // s'assurer que seed est impair 

              if ((lseed%2) == 0) lseed = lseed-1; 

              return lseed; 

              } 

//***************************************************************************

              // Affichage des points aléatoires dans un carré [0..1][0..1] en
              utilisant 

              // GNUPlot 

              // Le fichier Random.p contient les commandes Gnuplot de tracé 

              // Les données sont contenues dans le fichier Donnees.dat 

//***************************************************************************

              void PlotSerie(void) 

              { 

                FILE *pipe; 

                // ouverture d'un pipe vers Gnuplot et lancement du tracé 

                pipe = popen("gnuplot -persist Random.p","w"); v // fermeture du
              pipe 

                pclose(pipe);v return; v } 

              //***************************************************************

              // Programme principal 

              //***************************************************************

              int main(int argc, char *argv[]) 

              { 

              FILE *Donnees; 

              char Buffer[20]; 

              int i; 

              double Serie[NBTIRAGE],Somme, Moyenne, EcartType;

              // initialisation du générateur 

              printf("Generateur de nombres aleatoires selon une distribution
              uniforme" 

              " - Dominique Lefebvre Octobre 2009\n\n"); 

              seed = GenerationSeed();

              // génération d'une série aléatoire de distribution uniforme 

              for(i=0;i<NBTIRAGE;i++) Serie[i] = Random(); 

              // Vérification du caractère aléatoire de la série 

              // pour une distribution aléatoire uniforme entre 0 et 1, la
              moyenne de 

              // la distribution doit être de 0.5 et son écart 0,289 (1/12) 

              Somme = 0.; 

              for(i=0;i<NBTIRAGE;i++) Somme = Somme + Serie[i]; 

              Moyenne = Somme/NBTIRAGE; 

              Somme = 0.; 

              for(i=0;i<NBTIRAGE;i++) Somme = Somme + (Serie[i] -
              Moyenne)*(Serie[i] - Moyenne); v EcartType = sqrt(Somme/NBTIRAGE);
              

              printf("Moyenne : %.5f Ecart type: %.5f\n\n",Moyenne, EcartType);
              

              printf("Pour une distribution aleatoire uniforme entre 0 et 1, " 

              "pour un nombre infini de tirages:\n"); 

              printf("- la moyenne vaut 0.5\n"); 

              printf("- l'ecart-type vaut 1/12 soit 0.289\n"); 

              // stockage des données dans le fichier Donnees.dat pour affichage
              v // ce fichier est dans le même répertoire que le programme
              RandomGenerator 

              Donnees = fopen("Donnees.dat","w"); 

              for(i=0;i<NBTIRAGE-10;i++) fprintf(Donnees,"%.5f
              %.5f\n",Serie[i],Serie[i+10]); 

              fclose(Donnees);

              // Visualisation de la distribution aléatoire uniforme 

              PlotSerie(); 

              // Fin du programme 

              return EXIT_SUCCESS;

              }

Les points particuliers de ce programme :

le nombre de tirages, c'est à dire le nombre d'éléments de la série, est fixé par la constante NBTIRAGE (original, non!).
la routine GenerationSeed : elle génére le germe (seed) nécessaire à l'algo de Lehmer. La valeur de seed doit être différente avant chaque construction de distribution. J'obtiens cette garantie en construisant mon seed à partir de la date/heure de l'instant de lancement. Cette idée est due à S.L Anderson (1990), un camarade de SIAM (Society for Industrial and Applied Mathematics) dans son article "Random number generators on vector supercomputers and other advanced architectures" in "SIAM Review 32, 221-51".
La communication avec GnuPlot se fait par un pipe, technique que je n'avais pas encore présenté. Il faut bien sur que le programme wgnuplot.exe soit installé sur votre machine.
Rien d'autre de particulier, programme très simple qui doit tourner du premier coup.

Le code source du programme RandomGenerator.cpp est disponible dans la bibliothèque de codes de TangenteX.

L'estimation de Pi

Pour s'amuser, voici un petit programme, tiré du premier, qui permet de calculer une valeur estimée de pi à partir d'une série de nombres aléatoires. Je vous laisse réfléchir à la méthode utilisée, qui est très commune dans certaines techniques d'intégration (la célèbre méthode de Monte-Carlo par exemple).


              //*************************************************************************

              // Calcul du nombre Pi par tirages aléatoires 

              // Dominique Lefebvre - Octobre 2009 V1.0 

              // TangenteX 

              // 

              // http://www.tangenteX.com

//**************************************************************************

              #include <cstdkib>

              #include <iostream>

              #include <math.h>

              #include <time.h>

              using namespace std;

              // déclaration des constantes 

              #define NBTIRAGE 10000 // déclaration des variables globales 

              long seed;v
              //***************************************************************************

              // Génération de nombres aléatoires distribués uniformément entre
              0 et 1 

              // Selon la formule de Lehmer et l'algorithme et les paramètres de
              

              // Parker and Miller (1988): 

              // x(i+1) = a*x(i) mod c 

              // Sa période est fixée par c, soit 2^31 - 1 

//******************************************************************************

              double Random(void) 

              { 

                const double a = 16807.0; // 7^5 

                const double c = 2147483647.0; // 2^31 - 1 

                double x,y; 

                x = a*seed;

                seed = (long)(fmod(x,c)); 

                y = seed/c; // normalisation par c 

                return (y); 

              } 

//******************************************************************************

              // Génération de la constante d'initialisation du générateur
              Random 

              // à partir de la date courante selon Anderson (1990) 

//******************************************************************************

              long GenerationSeed(void) 

              { 

              struct tm *temps; 

              time_t maintenant; 

              int t0,t1,t2,t3,t4,t5; 

              long lseed; 

              // récupération de l'heure actuelle 

              time(&maintenant); 

              // enregistrement de l'heure actuelle 

              temps = localtime(&maintenant);

              // conversion en structure tm 

              t0 = temps->tm_sec; 

              t1 = temps->tm_min; 

              t2 = temps->tm_hour; 

              t3 = temps->tm_mday; 

              t4 = temps->tm_mon; 

              t5 = temps->tm_year; 

              > // calcul du seed d'après l'algo d'Anderson (seed varie entre
              0 et 2^31-1 

              lseed = t5 + 70*(t4 + 12*(t3 + 31*(t2 + 23*(t1 + 59*t0)))); 

              // s'assurer que seed est impair 

              if ((lseed%2) == 0) lseed = lseed-1; 

              return lseed; 

              } 

//******************************************************************************

              // Programme principal 

//******************************************************************************

              int main(int argc, char *argv[]) 

              {

              int i, pi4; 

              double x,y,pi;

              // initialisation du générateur 

              printf("Calcul du nombre PI par tirages aleatoires" 

              " - Dominique Lefebvre Octobre 2009\n\n"); 

              seed = GenerationSeed(); 

              // calcul de PI 

              pi4 = 0; 

              for(i=0;i<NBTIRAGE;i++) 

              { 

              x = Random(); 

              y = Random(); v if ((x*x + y*y) < 1) pi4++; 

              } 

              pi = 4.0*pi4/NBTIRAGE; 

              printf("Valeur estimee de PI: %.6f pour %i
              tirages\n",pi,NBTIRAGE); 

              // Fin du programme 

              return EXIT_SUCCESS;

Faites tourner ce programme, vous constaterez que si le nombre de tirages est grand, la valeur de pi obtenue n'est pas complètement idiote.
Le source CalculPi.cpp est disponible dans la bibliothèque de codes de TangenteX.