Factorisation par la méthode du crible quadratique

Sujet proposé par Luc Maranget

Luc.Maranget@inria.fr

Important, ce projet est accompagné d'une page de suivi, qui propose en outre des documents complémentaires.


1  Factorisation des grand nombres

La factorisation effective d'entiers toujours plus grands, que l'on pouvait autrefois qualifier d'amour de l'art, voire de hobby (mais aussi d'effort de recherche fondamentale), connaît des applications médiatisées en cryptographie, car la sécurité des systèmes de cryptage modernes semble bien reposer en dernière analyse sur la difficulté de répondre à cette question toute bête : Soit ce nombre n, dont je sais qu'il n'est pas premier, quels sont ses facteurs ?

Dans ce projet nous allons aborder une méthode de factorisation qui a fait date. La méthode du crible quadratique est une introduction aux méthodes générales de factorisation modernes, qui en sont souvent des raffinements. Par ailleurs, le crible quadratique s'appuie sur des concepts raisonnablement simples de la théorie des nombres, ce qui place sa mise en oeuvre à notre portée.

2  Principe

La méthode de factorisation repose sur une idée assez ancienne (1920) du mathématicien Maurice Kraitchick. Soit n un entier que nous souhaitons factoriser. Supposons que, par un moyen quelconque, nous connaissons deux autres entiers X et Y, tels que X2 ≡ Y2 (mod n ) soit vrai, Il existe donc un entier k tel que (X-Y)(X+Y) = k n. Avec un peu de chance, X-Y et n ont des facteurs non triviaux en commun et pgcd(X-Y, n) est un facteur de n. (en cas de malchance on essaie une autre congruence). Cette idée est un raffinement d'une autre, encore plus ancienne, et due à Pierre de Fermat : pour factoriser n on peut essayer de l'exprimer comme la différence de deux carrés n = X2 - Y2 (et donc n = (X-Y)(X+Y)). Pour dire toute la vérité, la méthode de Kraitchick ne fonctionne que si n peut s'écrire comme produit de deux facteurs premiers entre eux, c'est-à-dire si n n'est pas une puissance (et n'est pas premier bien évidemment). Ce n'est pas bien grave.

En ce qui nous concerne, la supériorité de l'idée de Kraitchick tient à ce qu'on peut trouver ses deux fameux entiers X et Y pour un prix plus modique que ceux de Fermat. On peut en effet procéder ainsi :
  1. Définir les petits nombres premiers en considérant les nombres premiers inférieurs à un certain entier B. On note p1, p2, …, pb les petits nombres premiers, que l'on complète par p0 = -1 pour pouvoir factoriser des entiers négatifs.

  2. Chercher les entiers xi tels la décomposition en facteurs premiers de yi = xi2 - n s'écrit uniquement à l'aide de petits nombres premiers. Un tel yi est dit B-régulier. On obtient ainsi des congruences de la forme xi2 ≡ p0α0, i p1α1, i ⋯ pbαb, i (mod n ) (on peut toujours poser αk,i = 0 si pk ne divise pas yi). Pour augmenter la possibilité de trouver des yi B-réguliers, on essaie les yi à priori petits, en prenant les xi voisins de n1/2.

  3. Multiplier membre à membre certaines des équations de l'étape précédente afin d'obtenir une égalité X2 = Y2 modulo n. Il n'y aucun problème pour obtenir le carré X2 (c'est le produit de certains xi). Pour obtenir le carré Y2, on utilise les décompositions en petits facteurs premiers, en sélectionnant les équations de sorte que la somme des αk, i soit paire pour tous les k.
Un exemple simple illustrera la méthode, cherchons à factoriser le nombre n = 13483. Un peu au hasard nous choisissons B = 20. On sélectionne ensuite des valeurs du polynôme P(x) = x2 - n (dit polynôme de Kraitchick) qui sont B-régulières, pour x voisin de ⌊n1/2⌋ = 116. On aura ainsi par exemple, les équations (modulo n) :
1252 = 2 ⋅ 32 ⋅ 7 ⋅ 17,         1272 = 2 ⋅ 33 ⋅ 72,         912 = -1 ⋅ 2 ⋅ 32 ⋅ 172,         1132 = -1 ⋅ 2 ⋅ 3 ⋅ 7 ⋅ 17,         1672 = 2 ⋅ 3 ⋅ 74
(On notera la présence de -1, due à des valeurs de xi inférieures ou égales à ⌊n1/2⌋.)

En multipliant toutes les équations sauf la seconde entre elles on obtient :
(125 ⋅ 91 ⋅ 113 ⋅ 167)2 ≡ (-1)2 ⋅ 24 ⋅ 36 ⋅ 76 ⋅ 174 (mod n )
Soit encore X ≡ 8265 (mod n ) et Y = 22 ⋅ 33 ⋅ 73 ⋅ 172 ≡ 13269 (mod n ). Au final, on calcule pgcd(8265 - 13269, n) qui vaut 139. Et effectivement on a bien n = 139 ⋅ 97.

Pour rendre la méthode de factorisation décrite effective, il faut résoudre deux questions.

3  Le crible quadratique

Le crible quadratique est une genéralisation du crible d'Eratosthène. En effet, pour tout polynôme P(x) et tout entier p nous avons P(x) ≡ P(x + p) (mod p ). Soit, supposant connu xa tel que P(xa) ≡ 0 (mod p ), nous pouvons rapidement conclure que P(xa+p), P(xa+2p), …, P(xa+ℓ p), … sont divisible par p (dans le crible d'Eratosthène, on a P(x) = x et xa = p). En combinant ces parcours de p en p pour divers facteurs premiers, nous allons identifier des P(x) qui sont B-réguliers.

On considère donc une base de facteurs p0, p1, …, pb (p0 = -1, p1 = 2, etc.) en supposant bien évidemment que n n'est divisible par aucun des « vrais » pk (auquel cas nous connaîtrions un facteur de n). Cette base ne regroupe pas tous les nombres premiers inférieurs à B. En effet, il n'y a pas nécessairement de racine modulo p de l'équation P(x) = 0 pour p quelconque. Nous supposons donc que la base de facteurs regroupe les petits nombres premiers tels que P(x) a des racines modulo pk. Dans le cas du polynôme de Kraitchick (qui est quadratique) et pour un pk (k ≥ 1) donné, il y aura une ou deux solutions. Pour p1 = 2, par hypothèse n impair, nous avons une seule solution qui est x ≡ 1 (mod 2 ). Pour tout autre pk premier, il y a deux solutions (car n n'est pas divisible par pk) qui sont d'ailleurs opposées modulo pk, et notées rk et r'k. Par exemple, dans le cas n = 13483, B=20 la base de facteurs se restreint à -1, 2, 3, 7, 17 avec r1 = 1, r2 = 1, r3 = 1 et r4 = 6 (et aussi r'2 = 2, r'3=6, r'4=11).

Étant donnée une taille T et une valeur initiale x0 nous considérons l'intervalle [x0… x0+T[ d'entiers notés x0, x1, … xT-1. Dans le même esprit nous posons yi = P(xi). Soit encore un tableau w d'entiers wi valant initialement zéro. Pour un pk donné nous pouvons détecter tous les yi multiples de pk et enregistrer ce fait dans le tableau w. Il suffit de calculer x0 mod pk et d'en déduire xa et xa', racines de P(x) = 0 modulo pk telles que a et a' sont minimaux (pour k=1 on a a = a'). Ensuite, on parcourt le tableau w de pk en pk à partir de a et de a', en ajoutant le logarithme de pk aux wi ainsi repérés. En pratique, on utilise une approximation entière du logarithme de base 2, qui conduit à des calculs peu coûteux. Une fois ce parcours effectué on applique l'algorithme de factorisation par division des facteurs de la base aux yi pour lesquels wi est suffisamment proche du logarithme de yi. On obtient ainsi une nouvelle équation xi2 = … presqu'à coup sûr.

Dans notre exemple, posons T=32 et criblons à partir de x0 = 117. Le nombre x0 est impair, comme r1, on a donc a=0. Après un premier passage le tableau w est donc
1   0   1   0   1   0   1   0   1   0   1   0   1   0   1   0   1   0   1   0   1   0   1   0   1   0   1   0   1   0   1   0  
Considérons maintenant le passage pour pk = 3. On a x0 mod 3 =0, r2 = 1, r'2 = 2, et donc a = 1 et a' = 2. En approximant log2(3) par 2 et en criblant, il vient :
1   2   3   0   3   2   1   2   3   0   3   2   1   2   3   0   3   2   1   2   3   0   3   2   1   2   3   0   3   2   1   2  
Il reste à effectuer le criblage des multiples de 7 (log2(7) ∼ 3) et de 17 (log2(17) ∼ 4). On obtiendra alors les états successifs du tableau w.
1   5   3   3   3   2   1   2   6   0   6   2   1   2   3   3   3   5   1   2   3   0   6   2   4   2   3   0   3   5   1   5  
1   5   3   3   3   2   1   2   10   0   6   2   1   6   3   3   3   5   1   2   3   0   6   2   4   6   3   0   3   5   5   5  
La case w8 attire l'attention, on a d'une part w8 = 10 et d'autre part log2(y8) ∼ 12 (y8 = 2141). L'entier y8 est sans doute un produit de facteurs de la base (de fait ici on sait que y8 est divisible par 2 ⋅ 3 ⋅ 7 ⋅ 17 et que le quotient est proche de 22), on lance une tentative de factorisation de y8 sur la base de facteurs. Il vient :
1252 ≡ 2 ⋅ 32 ⋅ 7 ⋅ 17 (mod n )
On notera que le décalage entre w8 et log2(y8) peut s'expliquer ici parce que nous n'avons pas pris la peine de cribler selon les puissances des facteurs de la base, mais qu'il provient aussi dans le cas général de l'usage d'une approximation entière du logarithme.

Le crible s'applique en fait à des nombres plus grands que la valeur de n donnée en exemple. Il se révèle alors très efficace pour repérer les nombres B-réguliers, à condition de considérer attentivement certains points.

4  Algèbre linéaire

Après criblage, nous sommes donc en possession de (nombreuses) congruences xi2 ≡ p0α0, i p1α1, i ⋯ pbαb, i (mod n ), que nous cherchons à multiplier afin de produire un carré à droite. Cela revient à trouver un ensemble I d'indices i tel que, pour tous les k, les sommes Σi ∈ I αk,i sont paires. Cela revient de fait à trouver une famille liée dans l'espace vectoriel F2b+1.

À chaque congruence, associons un vecteur zi = (z0i   z1i   ⋯   zb+1i) de F2b+1, avec zki = αk,i mod 2. Nous cherchons en fait l'ensemble I tel que Σi ∈ I zi = 0. Opérationnellement on cherche à démontrer que la matrice Z = (zki) est singulière en exhibant une combinaison linéaire (somme sur F2b+1) de ses lignes qui soit nulle. L'algorithme classique du pivot de Gauss permet d'atteindre ce résultat (garanti dès que la matrice a plus de lignes que de colonnes). Considérons à nouveau notre exemple de la factorisation de 13483. En supprimant la seconde ligne, la matrice des parités des exposants possède 4 lignes de 5 colonnes, nous lui adjoignons une autre matrice T qui est la matrice unité 4 × 4.
Z = æ
ç
ç
ç
è
0 1 0 1 1
1 1 0 0 0
1 1 1 1 1
0 1 1 0 0
ö
÷
÷
÷
ø		 ,          T = æ
ç
ç
ç
è
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
ö
÷
÷
÷
ø
L'opération de base du pivot consiste pour un k donné à annuler les zki sauf pour un i donné, on y parvient par combinaison linéaire (ici somme) des lignes de Z avec une ligne i0 telle zki0 n'est pas nul. Les combinaisons de lignes sont enregistrées dans les lignes de la matrice T. L'élément conservé est dit pivot. Si le nombre de lignes de Z dépasse celui de ses colonnes, on est sûr d'aboutir (à condition de prendre tous les pivots dans des lignes différentes). Ici choisissons d'abord d'éliminer les 1 de la première colonne, en prenant z01 comme pivot. Il vient :
Z0 = æ
ç
ç
ç
è
0 1 0 1 1
1 1 0 0 0
0 0 1 1 1
0 1 1 0 0
ö
÷
÷
÷
ø		          T0 = æ
ç
ç
ç
è
1 0 0 0
0 1 0 0
0 1 1 0
0 0 0 1
ö
÷
÷
÷
ø
La troisième ligne de la matrice T0 est la somme des deuxième et troisième lignes de T : elle nous dit comment obtenir cette troisième ligne de Z0 à partir de Z.

On choisit ensuite le pivot z10 (dans Z0).
Z1 = æ
ç
ç
ç
è
0 1 0 1 1
1 0 0 1 1
0 0 1 1 1
0 0 1 1 1
ö
÷
÷
÷
ø		          T1 = æ
ç
ç
ç
è
1 1 0 0
0 1 0 0
0 1 1 0
1 0 0 1
ö
÷
÷
÷
ø
Prenons enfin le pivot z22 :
Z2 = æ
ç
ç
ç
è
0 1 0 1 1
1 0 0 1 1
0 0 1 1 1
0 0 0 0 0
ö
÷
÷
÷
ø		          T2 = æ
ç
ç
ç
è
1 1 0 0
0 1 0 0
0 1 1 0
1 1 1 1
ö
÷
÷
÷
ø
On voit bien ici que la quatrième ligne de T2 s'obtient comme la somme des deux dernières lignes de T1. Elle nous dit toujours comment la quatrième ligne de Z2 s'obtient en sommant certaines des lignes de la matrice Z de départ (ici toutes). Notons qu'un pivot de Gauss traditionnel échange les lignes des matrices, (afin de systématiquement placer le pivot en zkk), ce qui nous donnerait ici :
Z2 = æ
ç
ç
ç
è
1 0 0 1 1
0 1 0 1 1
0 0 1 1 1
0 0 0 0 0
ö
÷
÷
÷
ø		          T2 = æ
ç
ç
ç
è
0 1 0 0
1 1 0 0
0 1 1 0
1 1 1 1
ö
÷
÷
÷
ø
Les échanges de lignes dans Z sont aussi effectués dans T. Ainsi, rien ne change à la lecture que nous faisons des lignes de T comme exprimant des sommes de lignes de la matrice Z de départ.

L'implémentation efficace du pivot de Gauss dans notre cas particulier repose sur deux remarques.
  1. On peut représenter les lignes des matrices Z et T comme des vecteurs de bits (classe BitSet de java). La sommation est alors l'opération « ou exclusif ».
  2. Les matrices produites par le crible sont creuses, dans le sens qu'elles contiennent bien plus de 0 que de 1. Cette propriété autorise des spécialisations du pivot de Gauss, décrites dans [1] (fourni dans les documents complémentaires).

5  Travail demandé

Réaliser un programme Qs qui prend comme argument un nombre n adapté (n n'a pas de petit facteurs et possède au moins deux facteurs premiers distincts) et produit une factorisation de n. Par exemple : 
# java Qs 350243405507562291174415825999
75576435361433 * 4634293795844903
Il s'agit d'un exemple « facile » dont la factorisation est quasi instantanée. La page de suivi donne d'autres exemples. En particulier, la factorisation du septième nombre de Fermat, F7 = 2128 + 1, exploit mondial dans les années 70 est largement à notre portée (elle peut prendre moins d'une minute).

Une répartition possible du travail dans le binôme est évidemment entre crible et pivot de Gauss. Attention, toutefois, il faut se rendre compte que, pour les tailles de nombres que nous traitons (pas plus de soixante chiffres), la réalisation d'un pivot de Gauss très optimisé semble moins critique que celle d'un crible efficace.

On utilisera donc les classes BigInteger et BitSet. Ces classes fournissent toutes les opérations de base dont vous aurez besoin à l'exception des « racines carrées » entières et modulo p. Je fournis (sur la page de suivi) le code nécessaire, afin que vous puissiez vous concentrer sur la programmation de l'algorithme (plus que sur la théorie des nombres, certes fascinante). L'objectif à atteindre est la factorisation de nombres d'environ soixante chiffres en quelques heures.

6  Prolongements possibles

Vous pouvez par exemple :
  1. Proposer un programme de factorisation pratique qui s'applique aux entiers quelconques, on intégrera par exemple des tests de primalité, l'algorithme de division voire l'algorithme dit ρ de Pollard. Il faudra aussi détecter les puissances. Un bon point de départ pour aborder les algorithmes nécessaires est certainement [2] (à la bibliothèque).
  2. Accélérer significativement votre crible. Vous pouvez commencer par déterminer si il vaut la peine ou pas de cribler selon les puissances des facteurs de la base, puis raffiner l'algorithme. Les méthodes sont nombreuses (polynômes multiples, distribution du criblage, etc.).

Références

[1]
C. Pomerance   J. W. Smith. Reduction of huge, sparse matrices over a finite field via created catastrophes. Experimental Math. 1 (1992), 90-94, 1992.

[2]
D. E. Knuth. The Art of Computer Programming, Volume 2: Seminumerical Algorithm, third Edition. Addison-Wesley, 1997.

[3]
C. Pomerance. A tale of two sieves. The Notices of the Amer. Math. Soc. 43 (1996), 1473-1485, 1996.
Ce document a été traduit de LATEX par HEVEA.