/* Ensemble de routines de gestion de nombres pseudo-al‚atoires **************************************************************************** Auteur : Pierre Larbier D‚velopp‚ le : 3.4.1994 DerniŠre MAJ : 3.4.1994 Toutes les routines de ce fichier sont librement utilisables (si vous les employez dans un programme commercial, envoyez-moi un petit mot ... ‡a fait plaisir). Elles ont ‚t‚ test‚es mais c'est sans aucune garantie ... ***************************************************************************/ #include #include #include #include #include #include /* Prototypes des fonctions de ce fichier */ /******************************************/ /* Impl‚mentation du Standard Minimal de Park & Miller */ long alea_stdmin(void); /* G‚n‚rateur de Mitchell & Moore sur 32 bits */ unsigned long alea_Mitchell_Moore(void); /* G‚n‚rateur polynomial x^15+x^1+1 */ char alea_bit(void); /* G‚n‚rateur uniforme, sortie dans l'intervalle [0,1[ */ #define alea_uniforme() (double)alea_Mitchell_Moore()/4294967296. /* G‚n‚rateur normal */ double alea_normal(double moyenne , double ecart_type); /* G‚n‚rateur exponentiel */ double alea_exp(double moyenne); /* Fixe un nouveau point de d‚part "al‚atoire" pour tous les g‚n‚rateurs */ void alea_randomize(void); /* Choisi un nouveau point de d‚part pour tous les g‚n‚rateurs */ void alea_srand(unsigned long depart); /* Fixe la borne sup‚rieure pour la fonction alea_random() */ void alea_fixe_borne(unsigned long max); /* G‚nŠre un nombre entier inf‚rieur … une borne fix‚e */ unsigned long alea_random(void); /* Variables globales utilis‚es par les divers g‚n‚rateurs , elles sont initialis‚es correctement */ long alea_var = 933757703L; unsigned long alea_tabMM[64] = { 2002705692L , 1963366013L , 58860689L , 1429122403L , 1803119173L, 1882197794L , 1664203448L , 1440332008L , 1224389472L , 1125550350L , 2088769674L , 1036733409L , 1843576952L , 1103773348L , 1154917050L , 1733657764L , 527917052L , 1446947207L , 736889421L , 362306498L , 1169172641L , 809207237L , 332095808L , 224246503L , 77175436L , 7430064L , 323034122L , 395828838L , 1938425507L , 1790571159L , 1441123902L , 1648850048L , 1093775848L , 630659016L , 1654283967L, 79855660L , 2104281892L , 1905060048L , 1510533613L , 2134242504L , 794408887L , 724330410L , 1883889674L , 34859550L , 1768904866L , 220473794L , 1093764683L , 443008861L , 324122678L , 1511320354L , 324612962L , 1161588954L , 51715001L , 1590628419L , 1815400277L , 2132282610L , 66725134L , 462863404L , 1159461594L , 804397480L , 1098888495L , 659571265L , 103665041L , 689106370L }; short alea_ptMM=55, alea_pt55MM=0 , alea_pt24MM=31; short alea_polyreg = 1; double alea_sav; short alea_etat = 0; unsigned long alea_borne = 4294967295L; short alea_dec = 0; unsigned long alea_init = 0; /************************************************************************* Proc‚dure : alea_stdmin Auteur / Date : Pierre Larbier 3/4/1994 Fonction : Produit un nombre pseudo-al‚atoire entier en utilisant le g‚n‚rateur de Park & Miller : x(n+1) = x(n)*16807 mod (2^31-1) Les nombres produits sont des long dans l'intervalle [1,2^31-2]. La p‚riode du g‚n‚rateur est 2^31-2. L'‚tat du g‚n‚rateur est maintenu dans la variable globale alea_var (qui ne doit jamais valoir 0). L'algorithme de calcul ne n‚cessite ni division ni multiplication par un nombre de plus de 16 bits. La routine est donc facilement codable en asm sur la plupart des machines. Entr‚es : aucune Sortie : le nombre pseudo-al‚atoire produit. *************************************************************************/ long alea_stdmin(void) { long haut , bas , inter; 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; } /************************************************************************* Proc‚dure : alea_Mitchell_Moore Auteur / Date : Pierre Larbier 3/4/1994 Fonction : Produit un nombre pseudo-al‚atoire entier en utilisant le g‚n‚rateur de Mitchell & Moore : x(n) = [x(n-24) + x(n-55)] mod 2^32 Les nombres produits sont des long dans l'intervalle [0,2^32-1]. La p‚riode du g‚n‚rateur est un multiple de 2^55-1. L'‚tat du g‚n‚rateur est maintenu dans le tableau alea_tabMM[64] ainsi que dans les variables alea_ptMM, alea_pt24MM et alea_pt55MM. Il ne m'a pas sembl‚ utile de rajouter un m‚lange suppl‚mentaire des donn‚es produites comme certains le conseillent (mais sans en donner de justification). Afin d'‚viter un test, on utilise un tableau de 64 entr‚es au lieu des 55 qu'on trouve habituellement. ATTENTION : on utilise le fait que le cast implicite du r‚sultat de l'addition en unsigned long r‚alise le modulo 2^32 (je n'ai jamais vu de compilateur C pour lequel ce ne soit pas vrai). Entr‚es : aucune Sortie : le nombre pseudo-al‚atoire produit. *************************************************************************/ unsigned long alea_Mitchell_Moore(void) { unsigned long x; 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; } /************************************************************************* Proc‚dure : alea_bit Auteur / Date : Pierre Larbier 3/4/1994 Fonction : Produit un bit al‚atoire en utilisant le polynome g‚n‚rateur : x^15 + x + 1. La p‚riode du g‚n‚rateur est de 2^15-1. L'‚tat du polynome est maintenu dans la variable alea_polyreg (qui ne doit jamais valoir 0). Entr‚es : aucune Sortie : un char valant 1 ou 0. *************************************************************************/ char alea_bit(void) { short 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; } } /************************************************************************* Proc‚dure : alea_normal Auteur / Date : Pierre Larbier 3/4/1994 Fonction : Produit un nombre al‚atoire flottant suivant une distribution normale avec l'algorithme de Box-Muller-Marsaglia (voir D.Knuth, "Seminumerical Algorithms" p. 117) Le g‚n‚rateur uniforme servant de base est le Mitchell-Moore. Afin de limiter les calculs, on conserve le r‚sultat interm‚diaire dans la variable globale alea_sav et la position du g‚n‚rateur dans alea_etat. Entr‚es : moyenne : la moyenne des nombres d‚sir‚e. ecart_type : l'‚cart_type des nombres d‚sir‚. Sortie : le nombre pseudo-aleatoire. *************************************************************************/ double alea_normal(double moyenne , double ecart_type) { double v1 , v2 , x; if (alea_etat == 0) { do { v1 = 2.*alea_uniforme() - 1.; v2 = 2.*alea_uniforme() - 1.; x = v1*v1 + v2*v2; } while ((x >= 1.) || (x == 0.)); x = sqrt(-2.*log(x)/x); alea_sav = x*v1; x *= v2; } else x = alea_sav; alea_etat ^= 1; return (x*ecart_type + moyenne); } /************************************************************************* Proc‚dure : alea_exp Auteur / Date : Pierre Larbier 3/4/1994 Fonction : Produit un nombre al‚atoire flottant suivant une distribution exponentielle. Le g‚n‚rateur uniforme servant de base est le Mitchell-Moore. Entr‚es : moyenne : la moyenne des nombres d‚sir‚e. Sortie : le nombre pseudo-aleatoire. *************************************************************************/ double alea_exp(double moyenne) { double x; do { x = alea_uniforme(); } while (x == 0.); return (-moyenne*log(x)); } /************************************************************************* Proc‚dure : alea_randomize Auteur / Date : Pierre Larbier 3/4/1994 Fonction : Fixe un nouveau point de d‚part "al‚atoire" pour tous les g‚n‚rateurs de ce fichier. ATTENTION : cette proc‚dure n'est pas portable. En l'‚tat elle ne fonctionne QUE sur un compatible PC. Pour choisir le point de d‚part, on utilise le contenu du timer 0 du 8254. L'horloge qui l'anime est … 1.193167MHz donc sa valeur … l'instant du premier appel de la routine est "al‚atoire". Aux appels suivants, on fait tourner un g‚n‚rateur Park & Miller pour produire les nouvelles graines des g‚n‚rateurs. On peut donc produire 2^31-2 s‚quences diff‚rentes. Entr‚es : aucune Sortie : aucune *************************************************************************/ 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) */ outp(0x43,0); asm pushf asm cli alea_init = (unsigned long)time(NULL); lb = (unsigned long)inp(0x40); hb = (unsigned long)inp(0x40); asm popf alea_init = ((alea_init << 16) | (hb << 8) | lb) & 0x7fffffffL; } alea_var = alea_init; alea_init = alea_stdmin(); alea_srand(~alea_init); } /************************************************************************* Proc‚dure : alea_srand Auteur / Date : Pierre Larbier 3/4/1994 Fonction : Fixe un nouveau point de d‚part pour tous les g‚n‚rateurs de ce fichier. La graine est fournie par l'utilisateur. Le g‚n‚rateur Mitchell & Moore est initialis‚ par le Park & Miller. On effectue pas le test pour s'assurer que les 55 premiers nombres ne sont pas tous pairs parce que j'ai v‚rifi‚ que le Standard Minimal ne peut pas produire 55 nombres pairs de rang. On peut donc produire 2^31-2 s‚quences diff‚rentes. Si la graine fournie vaut 0, on prend 1 … la place pour ne pas "coincer" les g‚n‚rateurs. Entr‚es : depart : point de d‚part de tous les g‚n‚rateurs : les 15 bits de poids faible sont pris pour le g‚n‚rateur polynomial les 31 bits de poids faible sont pris pour le Park & Miller. le Mitchell & Moore est initialis‚ par le Park & Miller. le Box-Muller-Marsaglia est remis … 0. Sortie : aucune *************************************************************************/ void alea_srand(unsigned long depart) { unsigned long x; short i; /* G‚n‚rateur polynomial */ x = depart & 0x7fff; if (x == 0) x = 1; alea_polyreg = (short)(x); /* Standard Minimal */ x = depart & 0x7fffffffL; if (x == 0) x = 1; alea_var = (short)(x); /* Etat du Box-Muller-Marsaglia */ alea_etat = 0; /* Etat du Mitchell-Moore */ alea_ptMM = 55; alea_pt55MM = 0; alea_pt24MM = 31; for (i=0 ; i<64 ; i++) alea_tabMM[i] = alea_stdmin(); /* Je ne fait pas le test pour v‚rifier que les 55 premiers nombres ne sont pas tous pairs car dans le cas du Standard Minimal cela ne peut pas arriver (j'ai fait le test syst‚matique). */ alea_var = (short)(x); } /************************************************************************* Proc‚dure : alea_fixe_borne Auteur / Date : Pierre Larbier 3/4/1994 Fonction : Fixe la borne sup‚rieure des nombres fournis par alea_random(). Si la borne vaut 0, on prend 2^32-1 … la place. C'est ‚galement la valeur par d‚faut avant le premier appel de cette proc‚dure. Entr‚es : max : valeur maximale que pourront prendre les nombres fournis par alea_random(). Sortie : aucune *************************************************************************/ void alea_fixe_borne(unsigned long max) { if (max == 0) max = 4294967295L; alea_borne = max; alea_dec = 0; while ((max & 0x80000000L) == 0) { alea_dec++; max = max << 1; } } /************************************************************************* Proc‚dure : alea_random Auteur / Date : Pierre Larbier 3/4/1994 Fonction : d‚livre un nombre entier compris entre 0 et la borne sup‚rieure initialis‚e dans alea_fixe_borne() inclus. Cette proc‚dure utilise la m‚thode de r‚jection au lieu de faire une division. Si la borne sup‚rieure n'est pas initialis‚e avec alea_fixe_borne(), la valeur max par d‚faut est 2^32-1. Entr‚es : aucune Sortie : le nombre pseudo-al‚atoire. *************************************************************************/ unsigned long alea_random(void) { unsigned long x; do { x = alea_Mitchell_Moore() >> alea_dec; } while (x > alea_borne); return x; } /************************************************************************* Petit programme de d‚monstration de l'usage de quelques routines de ce fichier. *************************************************************************/ void main(void) { long i , j , k , d1 , d2 , d3 ,max , min , des[216]; double x , y; /* O— on ‚crit la petite pr‚sentation */ printf("\n\n"); printf("Exemple d'utilisation des routines de g‚n‚ration de nombres pseudo-al‚atoires\n"); printf("*****************************************************************************\n"); printf(" Attention : les tests effectu‚s ne correspondent en rien … ceux qui sont\n"); printf(" faits pour valider un g‚n‚rateur pseudo-al‚atoire\n\n"); /* L'aiguille de Buffon */ printf("On va tester le g‚n‚rateur uniforme en simulant 100 000 lanc‚s d'aiguilles\n"); printf(" de Buffon. Le r‚sultat obtenu doit s'approcher de PI. \n"); for (j=0 , i=0 ; i<100000L ; i++) { x = sin(alea_uniforme()*M_PI_2); y = 2*alea_uniforme(); if ( ((x+y) >= 2) || ((y-x) <= 0) ) j++; } printf("R‚sultat de la 'simulation' : %.5f\n\n",2*100000./(double)j); printf("Tapez sur une touche pour continuer\n"); getch(); printf("\n\n"); /* Les d‚s sont-ils pip‚s ?? */ printf("A pr‚sent, on simule 33 333 lanc‚s de 3 d‚s.\n"); for (i=0 ; i<216 ; i++) des[i] = 0; alea_randomize(); alea_fixe_borne(5); for (i=0 , j=0 ; i<33333L ; i++) { j = alea_random(); j += 6*alea_random(); j += 36*alea_random(); des[j]++; } j = 33334L; k = -1; for (i=0 ; i<216 ; i++) { if (des[i] > k) { max = i; k = des[i]; } if (des[i] < j) { min = i; j = des[i]; } } d1 = min/36; d2 = (min-d1*36)/6; d3 = min%6; printf("On a obtenu l'occurence minimale de %ld lanc‚s avec la combinaison %ld-%ld-%ld\n",j, 1+d1 , 1+d2 , 1+d3); d1 = max/36; d2 = (max-d1*36)/6; d3 = max%6; printf("On a obtenu l'occurence maximale de %ld lanc‚s avec la combinaison %ld-%ld-%ld\n",k, 1+d1 , 1+d2 , 1+d3); printf("La moyenne ‚tant de 154.3, un tel ‚cart vous inquiŠte sur autant de lanc‚s ?\n"); printf("On se fait un petit Chi2 pour voir sur les 33 333 tirages pr‚c‚dents?\n"); for(i=0 , x=0 ; i<216 ; i++) x += (double)des[i]*(double)des[i]; x = 216.*x/33333. - 33333.; printf("Le g‚n‚rateur est suspect si le Chi2 obtenu : %.2f n'est pas\n",x); printf("dans l'intervalle 182.1 - 250.1 (intervalle de confiance 5%% - 95%%).\n\n"); printf("P.S : comme on a fait un randomize() avant les lanc‚s , je ne peux pas\n"); printf(" pr‚voir le Chi2 calcul‚ 4 lignes plus haut. Si jamais il sort ou tangente\n"); printf(" les 2 bornes de l'intervalle de confiance, pas de panique : ‡a peut\n"); printf(" arriver (1 fois sur 10 en moyenne).\n\n"); printf("Tapez sur une touche pour continuer\n"); getch(); printf("\n\n"); /* Ecart-type du g‚n‚rateur uniforme sur 100 000 tirages */ for (x = 0, i = 0 ; i<100000L; i++) { y = alea_uniforme() - 0.5; x += y*y; } printf("Sur cent mille tirages, on a obtenu une variance\n"); printf(" de 1/%.5f avec le g‚n‚rateur uniforme (id‚alement 1/12)\n\n",100000./x); /* Moyenne du g‚n‚rateur exponentiel sur 100 000 tirages */ for (x = 0, i = 0 ; i<100000L; i++) x += alea_exp(5.); printf("Sur cent mille tirages, on a demand‚ une moyenne de 5. au\n"); printf(" g‚n‚rateur … distribution exponentielle. On a mesur‚ : %.5f\n\n",x/100000.); /* Ecart-type du g‚n‚rateur normal sur 100 000 tirages */ for (x = 0, i = 0 ; i<100000L; i++) { y = alea_normal(0 , 5.); x += y*y; } printf("Sur cent mille tirages, on a demand‚ un ‚cart-type de 5. au\n"); printf(" g‚n‚rateur … distribution normale. On a mesur‚ : %.5f\n\n",sqrt(x/100000.)); printf("\nEt c'est fini\n"); }