Générateurs de nombres aléatoires
TP1 - Générateurs à Congruence Linéaire (GCL)

A. Ridard

Le générateur à congruence linéaire

Question 1 :

  • Programmer une fonction monGCL$(a,b,m)$ qui change l'état de $s$ (variable globale) et retourne le nouvel entier : $as + b$ mod $m$
  • A partir des graines suivantes, générer les entiers avec $a=25$, $b=16$, $m=256$ jusqu'à la première répétition, et préciser la période :
    $s_0=125$, $s_0=96$, $s_0=50$ et $s_0=10$
  • Expliquer pourquoi ce générateur n'est pas satisfaisant
In [ ]:
# Réponse

Attention : la première répétition ne concerne pas toujours $s_0$

Question 2 :

A partir de la graine $s_0=5$, générer les entiers avec $a=6$, $b=2$, $m=2^{24}$ jusqu'à la première répétition, et préciser la période

In [ ]:
# Réponse

Question 3 :

  • Programmer une fonction monGCL_N$(a, b, m, g, N)$ qui retourne la liste des $N$ entiers $s_0, s_1, ..., s_{N-1}$ générés à partir de la graine $s_0 = g$ et la relation de récurrence $s_n = f(s_{n-1})$ avec $f(s) = as + b$ mod $m$
  • Tester la fonction pour retrouver, au passage, le résultat de la question 2
In [ ]:
# Réponse

Remarque : pour initialiser la graine, on pourra utiliser l'horloge de l'ordinateur

In [ ]:
import datetime as dt
date = str(dt.datetime.now())
s = int(date[-4:])

Le générateur multi-récursif

Question 4 :

Programmer une fonction monGMR_N$(a1,a2,a3,m,g,N)$ qui retourne la liste des $N$ triplets $\left(s_0^{(1)}, s_0^{(2)}, s_0^{(3)}\right), \left(s_1^{(1)}, s_1^{(2)}, s_1^{(3)}\right), ..., \left(s_{N-1}^{(1)}, s_{N-1}^{(2)}, s_{N-1}^{(3)}\right)$ générés à partir de la graine $\left(s_0^{(1)}, s_0^{(2)}, s_0^{(3)}\right) = g$ et la relation de récurrence $\left(s_n^{(1)}, s_n^{(2)}, s_n^{(3)}\right) = f\Big(\left(s_{n-1}^{(1)}, s_{n-1}^{(2)}, s_{n-1}^{(3)}\right)\Big) = \left(s_{n-1}^{(2)}, s_{n-1}^{(3)}, a_1s_{n-1}^{(1)} + a_2s_{n-1}^{(2)} + a_3s_{n-1}^{(3)} \mod m\right)$

In [ ]:
# Réponse

Le générateur MRG32k3

Question 5 :

Programmer une fonction monMRG32k3_N$(g,N)$ qui retourne la liste des $N$ réels entre 0 et 1 à partir de la graine $g$

In [ ]:
# Réponse

Qualité d'un générateur de nombres aléatoires

Supposons qu’on génére une suite de nombres réels entre 0 et 1. Évidemment il est difficile de comparer une telle suite déterministe (sauf pour le choix de la graine) avec une suite aléatoire. Pour cela, il faut d’abord choisir un critère de qualité. Mais, ce n’est pas si simple, il y a beaucoup de choses à vérifier pour être aléatoire, et le problème est de trouver des critères simples et pertinents. Un premier critère assez naturel, est de tester l’indépendance de deux termes consécutifs générés. Pour cela, on représente graphiquement les points $(x_k, x_{k+1})$ dans le carré $[0, 1]^2$. Un générateur sera bon (pour ce critère) s’il remplit assez rapidement tout le carré.

Question 6 :

Construire un vecteur $(x_0, x_1, \dots, x_{100})\in [0, 1]^{101}$, pour chacun des générateurs suivants, puis représenter les points $(x_k, x_{k+1})$ :

  • Le GCL de la question 2
  • Le générateur MRG32k3
  • Le générateur de Python

Commenter les résultats obtenus.

In [ ]:
# Réponse

Question 7 :

Considérons le GCL défini par :

  • $S = \{1, 2, \dots, 100\}$
  • $f(s) = as$ mod $101$
  • $g(s) = \displaystyle\frac{s}{101}$

En prenant respectivement les valeurs 51, 7 et 12 pour le paramètre $a$ du générateur, discuter sa qualité selon le critère précédent.

In [ ]:
# Réponse

Pour valider (ou plutôt invalider) la qualité d'un générateur, on peut aussi réaliser des tests statistiques (cf. pages 26 à 28 de cet article) :

  • les tests d'adéquation
  • le test de collision et/ou le test d’espacement des anniversaires

Attaque sur un GCL

Soit $x_0, x_1, ..., x_N$ les entiers générés à partir de la graine $x_0$ et la relation de récurrence $x_n = ax_{n-1} + b$ mod $m$.
Posons $y_k = x_{k}-x_{k-1}$ pour $k\geq 1$.

Cas 1 : $m$ connu

Question 8 :

Dans le cas où $y_1$ est inversible modulo $m$, montrer : $a = y_2y_1^{-1}$ mod $m$ et $b = x_1 - ax_0$ mod $m$.

Réponse :

Question 9 :

  • En sachant que $m = 1023$, déterminer le GCL ayant généré les entiers 97, 188, 235, 293, 604, 596, 412
  • Vérifier le résultat à l'aide de monGCL_N$(a, b, m, g, N)$ de la question 3

Remarque : une fonction permettant de calculer l'inverse modulo $m$ est fournie ci-dessous

In [ ]:
# Fonctions utiles

def pgcde(a, b):
    """ pgcd étendu avec les 2 coefficients de Bézout u et v
        Entrée : a, b entiers
        Sorties : r = pgcd(a,b) et u, v entiers tels que a*u + b*v = r
    """
    r, u, v = a, 1, 0
    rp, up, vp = b, 0, 1
    while rp != 0:
        q = r//rp
        rs, us, vs = r, u, v
        r, u, v = rp, up, vp
        rp, up, vp = (rs - q*rp), (us - q*up), (vs - q*vp)
    return (r, u, v)

#print(pgcde(6, 15))
#print(pgcde(35, 6))

def invModulo(y, m):
    """Sortie : inverse de (l'entier) y modulo (l'entier) m s'il existe, et -1 sinon"""
    res = -1
    r, u, v = pgcde(y, m)
    if r == 1:
        res = u%m # le modulo m permet d'éviter les négatifs
    return res

#print(invModulo(35, 6))
In [ ]:
# Réponse

Cas 2 : $m$ inconnu

Question 10 :

  • Montrer $y_{k+1} = ay_k$ mod $m$
  • En déduire $z_k = y_{k+1}y_{k-1}-y_k^2 = 0$ mod $m$

Réponse :

Remarque : tous les $z_k$ sont donc des multiples de $m$, autrement dit $m$ est un diviseur commun de tous les $z_k$ et "très probablement" le PGCD de tous les $z_k$

Question 11 :

  • A l'aide de la remarque précédente, déterminer le GCL ayant généré les entiers suivants :
    234, 1227, 12158, 2475, 26787, 30101, 12498, 18328, 76, 11400
  • Combien d'entiers suffisent pour déterminer ce GCL ?

Remarque : une fonction permettant de calculer le PGCD de $n$ entiers est fournie ci-dessous

In [ ]:
# Fonctions utiles

def pgcd(a, b):
    """Sortie : pgcd de 2 entiers"""
    while b != 0: 
        a, b = b, a%b
    if a < 0:
        a = -a
    return a

#print(pgcd(-129101980, 171993983))

def pgcdn(lst_n):
    """Sortie : pgcd de n (>=2) entiers"""
    p = pgcd(lst_n[0], lst_n[1])
    for x in lst_n[2:]:
        p = pgcd(p, x)
    return p

#print(pgcdn([-129101980, 171993983, -623162806]))
In [1]:
# Réponse