Comment "casser" un générateur à congruence linéaire

Le but est de retrouver les paramètres qui ont abouti à une séquence pseudo-aléatoire donnée.

Cas 1: m est connu

Soit la suite obtenue par la formule xn=a·xn-1+b mod m. On a besoin des trois premiers nombres de la suite: x0, x1 et x2. Posons y1=x1-x0 et y2=x2-x1. On a la relation: y2 = a·y1 mod m, d'où l'on tire immédiatement: a = y2y1-1 mod m, où y1-1 mod m est obtenu avec l'algorithme d'Euclide étendu. On obtient b facilement: b = (x1 - a·x0) mod m.

Exemple 1

On connaît le début de la séquence 97, 188, 235, 293, 604, 596, 412. On sait que le modulo est 1023.
  1. On calcule y1 = 188-97 = 91 et y2 = 235-188 = 47.
  2. L'algorithme d'Euclide étendu nous dit que l'inverse modulo 1023 de 91 = 697 = y1-1.
  3. On obtient donc a = 47·697 mod 1023 = 23 et b = (188 - 23·97) mod 1023 = 3.

Exemple 2

On connaît le début de la séquence 99, 234, 270, 75, 705, 873, 645. On sait que le modulo est 1023.
  1. On calcule y1 = 234-99 = 135 et y2 = 270-234 = 36.
  2. L'algorithme d'Euclide étendu nous dit que 99 n'a pas d'inverse modulo 1023!
On a affaire à un mauvais générateur: tous les nombres de la suite sont des multiples de 3.

Cas 2: m n'est pas connu (méthode de Plumstead-Boyar)

Prenons directement un exemple pour expliquer la méthode (on se référera aux références pour une théorie plus approfondie).
Soit le début de la séquence 234, 1227, 12158, 2475, 26787, 30101, 12498, 18328, 76, 11400.

Étape 1: calcul des différences
On calcule d'abord yi = xi-xi-1, pour i=1, ..., n. Le nombre (n) de termes à choisir se détermine ainsi: il faut le plus petit n tel que le pgcd des n premiers termes soit 1 (mais prendre plus de termes n'est pas gênant).

On obtient dans notre exemple y1=993, y2=10931, y3=-9683, y4=24312, y5=3314, y6=-17603, y7=5830, y8=-18252, y9=11324.

Ici on peut prendre n = 2, puisque pgcd(993, 10931) = 1.

Étape 2:
Il faut trouver la solution de l'équation diophantienne c1y1 + c2y2 + ... + cnyn = 1. La fonction Mathematica ExtendedGCD s'en charge (l'algorithme est une variante du théorème d'Euclide étendu)! On calcule ensuite la valeur a' = c1y2 + c2y3 + ... + cnyn+1.

On obtient dans notre exemple c1 = 1365, c2 = -124 et a' = 16'121'507.
Étape 3:
Poser a1 = a', et m1 = . Ensuite, pour j > 1, faire:
  • y'j = aj-1yj-1 mod mj-1
  • mj = pgcd(mj-1, y'j-yj)
  • aj = aj-1 mod mj
On s'arrête quand mj-1=mj

Dans notre exemple:

m1 = a1 = 16'121'507
y'2 = 16'008'656'451 m2 = 16'008'645'520 a2 = 16'121'507
y'3 = 129'092'297 m3 = 129'101'980 a3 = 16'121'507
y'4 = 108'843'519 m4 = 32'767 a4 = 143
y'5 = 3314 m5 = 32'767 a5 = 143 fin

Étape 4:
a = afin, m = mfin, b = x1 - a·x0 mod m.

Dans notre exemple, on trouve finalement: a = 143, m = 32'767, b = 532

Fichier Mathematica

Le fichier téléchargeable ci-dessous permet de casser un générateur à congruence linéaire (m connu ou pas).

casser.nb (Mathematica 4)

Références

Didier Müller, 21.9.03