# PRF - TP 1 [Sujet](https://moodle.insa-lyon.fr/pluginfile.php/47811/mod_resource/content/3/2018/4tc-prf_tp1-rng.pdf) ###### tags : `PRF` `TP` ## Introduction ![](https://i.imgur.com/D2spkQM.png) ## Why random number generation ? ![](https://i.imgur.com/AreTSoA.png) ## Generating random numbers uniformly in $[0, 1]$ ![](https://i.imgur.com/IKPG3E5.png) ![](https://i.imgur.com/j4wMMOj.png) La propriété cyclique de ces séquences peut poser problème car elle leur retire le côté aléatoire (On a toujours les mêmes valeurs sur lesquelles on boucle et les enchaînement sont toujours les mêmes). ![](https://i.imgur.com/skwguGC.png) Les séquences générées doivent donc être beaucoup plus longues que le nombre de tirages que l'on souhaite réaliser. ### Linear congruential generators ![](https://i.imgur.com/VYkJQCq.png) ```lcg.c``` ```c= /* * Linear Congruential Generator */ #include <stdio.h> #include <stdlib.h> #define LCG_A 3 #define LCG_C 5 #define LCG_M 31 int main(int argc, const char* argv[]) { int seed; /* LCG seed */ int i = 1; /* current iteration */ int x; /* current random value */ /* process commmand line */ if(argc != 2) { fprintf(stderr, "\nSyntax: %s <seed>\n\n", argv[0]); return(0); } seed = atoi(argv[1]); fprintf(stdout, "# LCG, a=%d, c=%d, m=%d, seed=%d\n\n", LCG_A, LCG_C, LCG_M, seed); /* generate 2*LCG_M - 1 random values */ x = seed; while(i < 2*LCG_M - 1) { fprintf(stdout, "%4d\n", x); x = ((LCG_A*x) + LCG_C) % LCG_M; i++; } return(0); } ``` Séquence obtenue : ``` # LCG, a=3, c=5, m=19, seed=12 12 3 14 9 13 6 4 17 18 2 11 0 5 1 8 10 16 15 12 3 14 9 13 6 4 17 18 2 11 0 5 1 8 10 16 ``` ```cycle_length = 18``` ![](https://i.imgur.com/9iDkSGB.png) Si $a=0 \longrightarrow x_i=c\cdot\mod m$ Donc dès la 2ème valeur on va boucler sur $x_i=c\cdot\mod m$ ![](https://i.imgur.com/JuUUzLf.png) ``` # LCG, a=3, c=5, m=31, seed=9 9 1 8 29 30 2 11 7 26 21 6 23 12 10 4 17 25 18 28 27 24 15 19 0 5 20 3 14 16 22 9 1 8 29 30 2 11 7 26 21 6 23 12 10 4 17 25 18 28 27 24 15 19 0 5 20 3 14 16 22 ``` La taille du cycle est toujours égale à $m$ ![](https://i.imgur.com/2W9odVs.png) La taille du cycle est toujours $m$ mais le seed $x_0$ ne fait pas partie du cycle. ![](https://i.imgur.com/xOjvYCi.png) ### How random is a random sequence? ![](https://i.imgur.com/pC2O7sj.png) ```lcg_bis.c``` ```c= /* * Linear Congruential Generator */ #include <stdio.h> #include <stdlib.h> #define LCG_A 7 #define LCG_C 0 #define LCG_M 1000 int main(int argc, const char* argv[]) { int seed; /* LCG seed */ int i = 1; /* current iteration */ int x; /* current random value */ /* process commmand line */ if(argc != 2) { fprintf(stderr, "\nSyntax: %s <seed>\n\n", argv[0]); return(0); } seed = atoi(argv[1]); fprintf(stdout, "# LCG, a=%d, c=%d, m=%d, seed=%d\n\n", LCG_A, LCG_C, LCG_M, seed); /* generate 2*LCG_M - 1 random values */ x = seed; while(i < 2*LCG_M - 1) { fprintf(stdout, "%4f\n",((double)x)/(LCG_M-1)); x = ((LCG_A*x) + LCG_C) % LCG_M; i++; } return(0); } ``` On modifie juste l'affichage des valeurs : ```lcg.c``` ```c= fprintf(stdout, "%4d\n", x); ``` ```lcg_bis.c``` ```c= fprintf(stdout, "%4f\n",((double)x)/(LCG_M-1)); ``` ![](https://i.imgur.com/zsFB7ni.png) ![](https://i.imgur.com/RBXuqB0.png) On ne peut pas savoir si une séquence est meilleure que l'autre sans les avoir générées. ![](https://i.imgur.com/3L3HaS7.png) ```Séquence 1 : cycle_length=30``` ``` # LCG, a=3, c=0, m=31, seed=1 0.033333 0.100000 0.300000 0.900000 0.633333 0.866667 0.533333 0.566667 0.666667 0.966667 0.833333 0.433333 0.266667 0.800000 0.333333 1.000000 0.933333 0.733333 0.133333 0.400000 0.166667 0.500000 0.466667 0.366667 0.066667 0.200000 0.600000 0.766667 0.233333 0.700000 0.033333 0.100000 0.300000 0.900000 0.633333 0.866667 0.533333 0.566667 0.666667 0.966667 0.833333 0.433333 0.266667 0.800000 0.333333 1.000000 0.933333 0.733333 0.133333 0.400000 0.166667 0.500000 0.466667 0.366667 0.066667 0.200000 0.600000 0.766667 0.233333 0.700000 ``` ```Séquence 2 : cycle_length=3``` ``` # LCG, a=5, c=0, m=31, seed=1 0.033333 0.166667 0.833333 0.033333 0.166667 0.833333 0.033333 0.166667 0.833333 0.033333 0.166667 0.833333 0.033333 0.166667 0.833333 0.033333 0.166667 0.833333 0.033333 0.166667 0.833333 0.033333 0.166667 0.833333 0.033333 0.166667 0.833333 0.033333 0.166667 0.833333 0.033333 0.166667 0.833333 0.033333 0.166667 0.833333 0.033333 0.166667 0.833333 0.033333 0.166667 0.833333 0.033333 0.166667 0.833333 0.033333 0.166667 0.833333 0.033333 0.166667 0.833333 0.033333 0.166667 0.833333 0.033333 0.166667 0.833333 0.033333 0.166667 0.833333 ``` ```Séquence 3 : cycle_length=30``` ``` # LCG, a=13, c=0, m=31, seed=1 0.033333 0.433333 0.466667 0.900000 0.333333 0.200000 0.533333 0.733333 0.233333 0.966667 0.166667 0.100000 0.266667 0.366667 0.633333 1.000000 0.600000 0.566667 0.133333 0.700000 0.833333 0.500000 0.300000 0.800000 0.066667 0.866667 0.933333 0.766667 0.666667 0.400000 0.033333 0.433333 0.466667 0.900000 0.333333 0.200000 0.533333 0.733333 0.233333 0.966667 0.166667 0.100000 0.266667 0.366667 0.633333 1.000000 0.600000 0.566667 0.133333 0.700000 0.833333 0.500000 0.300000 0.800000 0.066667 0.866667 0.933333 0.766667 0.666667 0.400000 ``` ```Séquence 4 : cycle_length=20``` ``` # LCG, a=7, c=0, m=1000, seed=1 0.001001 0.007007 0.049049 0.343343 0.401401 0.807808 0.649650 0.543544 0.801802 0.607608 0.249249 0.743744 0.201201 0.407407 0.849850 0.943944 0.601602 0.207207 0.449449 0.143143 ``` Les séquences 1 et 3 sont ***full-period*** (tous les nombres entre $0$ et $m-1$ sont générés) #### Uniformity test ![](https://i.imgur.com/bJFlX6J.png) ```Séquence 1``` ``` prng_pdf 0.050 1.0000 prng_pdf 0.150 1.0000 prng_pdf 0.250 1.0000 prng_pdf 0.350 1.0000 prng_pdf 0.450 1.0000 prng_pdf 0.550 1.0000 prng_pdf 0.650 1.0000 prng_pdf 0.750 1.0000 prng_pdf 0.850 1.0000 prng_pdf 0.950 1.0000 ``` ```Séquence 2``` ``` prng_pdf 0.050 3.3333 prng_pdf 0.150 3.3333 prng_pdf 0.250 0.0000 prng_pdf 0.350 0.0000 prng_pdf 0.450 0.0000 prng_pdf 0.550 0.0000 prng_pdf 0.650 0.0000 prng_pdf 0.750 0.0000 prng_pdf 0.850 3.3333 ``` ```Séquence 3``` ``` prng_pdf 0.050 1.0000 prng_pdf 0.150 1.0000 prng_pdf 0.250 1.0000 prng_pdf 0.350 1.0000 prng_pdf 0.450 1.0000 prng_pdf 0.550 1.0000 prng_pdf 0.650 1.0000 prng_pdf 0.750 1.0000 prng_pdf 0.850 1.0000 prng_pdf 0.950 1.0000 ``` ```Séquence 4``` ``` prng_pdf 0.050 1.5015 prng_pdf 0.150 0.4955 prng_pdf 0.250 1.5015 prng_pdf 0.350 0.5005 prng_pdf 0.450 1.4965 prng_pdf 0.550 0.5005 prng_pdf 0.650 1.5015 prng_pdf 0.750 0.5005 prng_pdf 0.850 1.5015 prng_pdf 0.950 0.5005 ``` On remarque que contrairement aux séquences 1 et 3, les séquences 2 et 4 ne sont pas uniformément distribuées. Cependant, le test ne nous permet pas de déterminer si un des générateurs est meilleur que les autres. ![](https://i.imgur.com/ov9Hjm0.png) Non, les distributions sont les mêmes. ![](https://i.imgur.com/py9JB0q.png) Les séquences full-period contiennent tous les nombres entre $0$ et $m-1$. Si on répète $k$ fois le cycle, on a $k$ fois chaque valeur. L'intervalle $[0,m-1]$ (ou $[0,1]$ si on normalise les valeurs) est donc uniformément distribué. #### 2-dimensional test ![](https://i.imgur.com/kAuGqCX.png) ```2-dimensional.awk``` ```awk BEGIN { prev = -1; } { # skip comment or empty lines if($1 == "#" || $0 == "") next; if(prev != -1) print prev, $1; prev = $1; } ``` Pour chaque ligne, ce script affiche le couple $(x_{n-1},x_n)$ où $n$ est le numéro de la ligne. ![](https://i.imgur.com/4GTS6HS.png) On *plot* les données de ```file.log``` en utilisant les 2 premières colonnes du fichier. Le titre est "". On affiche des points. ![](https://i.imgur.com/2Mh8ZAp.png) ![](https://i.imgur.com/tzfhFk8.png) ```Séquence 2``` ![](https://i.imgur.com/kLzNi4E.png) On a 3 valeurs qui se répètent donc on a 3 points ```Séquence 4``` ![](https://i.imgur.com/xYfd8TH.png) 20 valeurs réparties de manière non uniforme donc 20 points désordonnés ![](https://i.imgur.com/0VqapUl.png) ```Séquence 1``` ![](https://i.imgur.com/53WKr5h.png) ```Séquence 3``` ![](https://i.imgur.com/0fT1oQe.png) On se rend compte que pour un tirage donné de la séquence 1, on a plus de chance de tirer une valeur se situant dans un espace proche (à quelques exceptions près) que pour un tirage donné de la séquence 3. La séquence 3 semble mieux distribuée. #### LCGs of practical use ![](https://i.imgur.com/fEZxzrb.png) ```c= /* * Linear Congruential Generator, * implemented using Schrage's algorithm */ #include <stdio.h> #include <stdlib.h> #define LCG_A 16807 #define LCG_C 0 #define LCG_M 2147483647 int main(int argc, const char* argv[]) { long int seed; /* LCG seed */ long int i = 1; /* current iteration */ long int x; /* current random value */ /* process commmand line */ if(argc != 2) { fprintf(stderr, "\nSyntax: %s <seed>\n\n", argv[0]); return(0); } seed = atoi(argv[1]); fprintf(stdout, "# LCG, a=%d, c=%d, m=%ld, seed=%ld\n\n", LCG_A, LCG_C, (long int)LCG_M, seed); /* define internal variables */ long int q = LCG_M / LCG_A; long int r = LCG_M % LCG_A; long int x_div_q, x_mod_q; /* generate 100000 random values */ x = seed; while(i < 100000) { fprintf(stdout, "%.9f\n", ((double)x)/(LCG_M-1)); x_div_q = x / q; x_mod_q = x % q; x = LCG_A*x_mod_q - r*x_div_q; if(x < 0) x += LCG_M; i++; } return(0); } ``` ```Uniform test``` ``` prng_pdf 0.050 0.9861 prng_pdf 0.150 1.0097 prng_pdf 0.250 0.9790 prng_pdf 0.350 0.9900 prng_pdf 0.450 1.0025 prng_pdf 0.550 1.0053 prng_pdf 0.650 1.0094 prng_pdf 0.750 1.0050 prng_pdf 0.850 0.9992 prng_pdf 0.950 1.0136 ``` ```2D test``` ![](https://i.imgur.com/neBFa46.png) A la vue des 2 tests, on peut dire que RANDU est un bon générateur. Le résultat n'est pas parfait mais on pouvait s'y attendre car on ne génère que 100 000 valeurs (sur un cycle qui en contient $2^{31}$). ![](https://i.imgur.com/79QA5Y6.png) ```3-dimensional.awk``` ```awk BEGIN { prev = -1; prev_bis = -2; } { # skip comment or empty lines if($1 == "#" || $0 == "") next; if(prev_bis != -1 && prev != -1 && prev_bis != -2) print prev_bis, prev, $1; prev_bis = prev; prev = $1; } ``` ![](https://i.imgur.com/7vKJAVC.png) ```3D test``` ![](https://i.imgur.com/YjumKWV.png) ![](https://i.imgur.com/PgGFlOw.png) ```Uniform test``` ``` prng_pdf 0.050 0.9861 prng_pdf 0.150 1.0097 prng_pdf 0.250 0.9790 prng_pdf 0.350 0.9900 prng_pdf 0.450 1.0025 prng_pdf 0.550 1.0053 prng_pdf 0.650 1.0094 prng_pdf 0.750 1.0050 prng_pdf 0.850 0.9992 prng_pdf 0.950 1.0136 ``` ```2D test``` ![](https://i.imgur.com/OjBmlRA.png) ```3D test``` ![](https://i.imgur.com/X61k8qa.png) On remarque que le test 3D pour MINSTD donne une répartition plus uniforme que RANDU. ![](https://i.imgur.com/4Z4fl1o.png)