Comprendre la méthode de factorisation du Crible Quadratique

En 2017 dans mon mémento consacré à la cryptologie (voir https://laurent-ott.developpez.com/tutoriels/programmation-excel-vba-tome-4/) j'écrivais au sujet des cryptosystèmes asymétriques (et donc du célèbre chiffrement RSA), que « c'est uniquement la difficulté, actuelle, à factoriser un grand nombre composé de plusieurs centaines de chiffres en deux nombres premiers qui assure la confidentialité de la clé privée. Car un petit nombre de 28 chiffres peut être décomposé en moins d'une seconde sur les sites Internet (http://www.dcode.fr/decomposition-nombres-premiers) ».

Un an après, la curiosité m'a poussé à vouloir en apprendre plus sur les méthodes utilisées pour décomposer un nombre de 28 chiffres en moins d'une seconde. Car la méthode basique, qui consiste à diviser un entier par les nombres premiers pour retrouver ses facteurs, demande un temps de traitement « déraisonnable » pour un nombre de cette taille. Alors même si les serveurs sont extrêmement rapides, il doit bien se cacher derrière cette vitesse de traitement impressionnante quelques algorithmes bigrement efficaces !

Naïvement je pensais trouver sur Internet la réponse à mes interrogations ainsi qu'un code source livré sur un plateau, mais j'ai vite déchanté : les thèses des mathématiciens étant hors de portée pour un néophyte comme moi, je n'ai pas trouvé grand-chose d'explicite et encore moins un algorithme complet.

Alors il m'a fallu un long travail d'investigation et de regroupement d'informations, de la persévérance et un peu de chance, pour mener à bien ce projet et arriver à factoriser une clé de 28 chiffres en une poignée de secondes à l'aide de l'algorithme dit « Crible Quadratique », développé en 1984 par Carl Pomerance, et qui reste encore aujourd'hui une référence incontournable.

Le fruit de ce travail, j'ai voulu le partager, c'est pourquoi je vous présente aujourd'hui cette documentation destinée à vous expliquer le plus simplement possible comment fonctionne cette méthode de factorisation réputée être l'une des plus efficaces.

Vous trouverez dans le fichier joint le code source en VBA, facilement adaptable dans votre langage de programmation préféré.

Entendons-nous bien, il ne s'agit pas de la version « officielle », je ne suis pas un professionnel du sujet, mais simplement d'une version personnelle conçue sur la base de mes différentes lectures.

Je suis bien conscient, d'une part, que ce code en VBA ne vous permettra pas de casser une clé RSA de 1024 bits (309 chiffres) d'un message secret avec votre ordinateur bureautique, car si la factorisation d'un nombre de 28 chiffres ne prend que quelques secondes, pour un nombre seulement deux fois plus grand il vous faudra compter cette fois en heures, et l'on atteint alors les limites des capacités de stockage de la matrice des données. Une « petite » clé d'une centaine de chiffres est donc inaccessible en VBA, aujourd'hui.

De toute façon cette documentation n'a pas pour sujet l'attaque des clés RSA.

Et d'autre part, certains langages de programmation intégrant une fonction de factorisation très performante, il semble donc inutile dans ce cas de perdre son temps à en écrire une autre.

Le but ici est uniquement de comprendre comment marche le crible quadratique et cette documentation s'adresse aux curieux.

Si vous n'avez jamais entendu parler du crible quadratique ou plus généralement des différentes méthodes de factorisation, avant de commencer votre lecture je vous invite à consulter ces deux sites : l'un est une très bonne présentation du sujet « Algorithmes de factorisation à l'envi » ; l'autre a l'avantage de présenter des exemples concrets et m'a été très utile « MT10. Mathématiques pour la cryptographie Partie 4 Factorisation : algorithme du crible quadratique (1) Walter SCHÖN ».

En reprenant l'exemple du dernier lien mentionné dans l'introduction pour le présenter sous un angle différent et peut-être vous permettre de mieux le comprendre, voici en résumé le principe de base de la méthode de factorisation imaginée par John D. Dixon en 1981 en s'inspirant des travaux de Maurice Kraitchik des années 1920, puis les améliorations successives qui ont conduit au crible quadratique :

Soit 2041, l'entier N à factoriser.

En partant de sa racine carrée arrondie à l'entier supérieur, X = 46, on calcule X² – N, soit 46² – 2041 = 75, et sa décomposition en nombres premiers : 75 = 3 x 5 x 5.

On répète ces calculs après avoir incrémenté X, soit 47² – 2041 = 168 = 2 x 2 x 2 x 3 x 7.

Avec X = 48² – 2014 = 263, la factorisation donne 263 car c'est un nombre premier.

Et ainsi de suite jusqu'à 51² – 2041 = 560 = 2 x 2 x 2 x 2 x 5 x 7.

Ce qui donne le tableau suivant où la factorisation de 75 = 3 x 5 x 5 = 3¹ x 5², donc 1 dans la colonne « x 3 » et 2 dans la colonne « x 5 » ;

autre exemple pour la deuxième ligne : 168 = 2³ x 3¹ x 7¹ ;

même chose pour les autres nombres :

Pour trouver dans cette matrice les « relations » qui forment la solution, nous transformons la factorisation obtenue ci-dessus en modulo 2 (les nombres pairs valent 0 et les impaires 1) comme l'ont proposé Michael A. Morrison et John Brillhart :

L'addition des colonnes des lignes de la matrice où X vaut 46, 47, 49, 51 en modulo 2 donne 0 (en bleu les cellules des colonnes « x 2, x 3, x 5, x 7 » font toutes une addition modulo 2 égale à zéro), soit la solution car en notant :

X = (46 x 47 x 49 x 51) modulo 2041 = 311 ;

et Y = Racine(75 x 168 x 360 x 560)  modulo 2041 = 1416 ;

on obtient X – Y modulo 2041 = (311 – 1416) modulo 2041 = –1105 soit 936 (gestion des modulos négatifs) ;

et X + Y modulo 2041 = 1727.

Le calcul du PGCD(X – Y, N) c'est-à-dire, PGCD(936,  2041)  = 13.

Le calcul du PGCD(X + Y, N) c'est-à-dire, PGCD(1727,  2041)  = 157.

Et donc 13 x 157 = 2041.

Tout simplement !

Nous verrons plus loin comment déterminer quelles lignes de la matrice il faut additionner pour trouver une solution, en utilisant le pivot de Gauss, concentrons-nous pour le moment sur la factorisation des nombres calculés par X² – N.

Pour factoriser un nombre il faut le diviser par les nombres premiers qui donnent un quotient entier jusqu'à ce que le quotient vaille 1.

Par exemple pour factoriser 459 (de la 5^e ligne du tableau), il faut diviser par 2, soit 459/2 = 229,5, le quotient n'est pas un entier donc 2 n'est pas un diviseur et l'on passe à 3 : 459/3 = 153 puis 153/3 = 51 puis 51/3 = 17 enfin 17/3 = 5,6 on mémorise donc trois fois que 3 est un diviseur et que le reste vaut 17.

On continue en tentant de diviser 17 par les nombres premiers suivants, c'est-à-dire 5, 11, 13, qui ne sont pas diviseurs et enfin 17/17 = 1, on mémorise une fois que 17 est diviseur, et comme le quotient vaut 1 on sait que la factorisation est terminée (voir la 5^e ligne du 1^er tableau ci-dessus).

Cela nécessite donc de nombreuses opérations, surtout pour factoriser 263.

Les mathématiciens ont par conséquent développé des astuces pour limiter le nombre de divisions à effectuer, très gourmandes en temps de traitement, et ainsi améliorer l'efficacité de cette méthode.

La première amélioration permettant de réduire considérablement les temps de traitement consiste à limiter la liste des premiers que l'on va utiliser pour les factorisations, quitte à faire l'impasse sur certaines relations qui seraient fructueuses.

Une formule que les mathématiciens vous expliqueront mieux que moi détermine le « seuil de friabilité », qui permet d'optimiser la taille de cette liste, dite « base des facteurs », car si elle est trop grande les relations seront plus faciles à trouver mais il y en aura trop, inversement si elle est trop petite, il y aura moins de relations mais elles seront plus difficiles à trouver :

• Seuil de friabilité = Entier(Exp((1/2) x Racine(Log(N) x Log(Log(N)))))

Ou exprimé autrement quand N vaut 2041 :

• Entier(2,718281828459 Puissance(0,5 x Racine(Log(2041) x Log(Log(2041))))) = 7.

Ce qui signifie que la base des facteurs se limitera aux nombres premiers inférieurs ou égaux à 7, soit 2, 3, 5, 7.

Effectivement, dans l'exemple présenté ci-dessus, ces quatre diviseurs suffisent à générer assez de relations fructueuses pour trouver une solution. Et la tentative de factorisation de 263 se limitera à quatre divisions infructueuses avant d'abandonner.

La deuxième amélioration consiste à supprimer de la base des facteurs calculée ci-avant ceux qui n'ont aucune chance de diviser X² – N.

Le test d'Adrien-Marie Legendre renvoie 1 dans la formule suivante si un nombre premier (supérieur à 2) noté p peut servir à factoriser N :

• Résultat = (N modulo p) Puissance ((p-1)/2) modulo p

Dans notre cas, les premiers 2, 3, 5, 7 sont conservés. Un test sur 11 renvoie que ce premier ne divise pas X² – N. Mais de toute façon ce nombre n'était pas dans la base des facteurs. Un exemple sur un nombre plus grand montrera l'intérêt de cette amélioration.

Enfin, Carl Pomerance proposa en 1984 une amélioration décisive en évitant de nombreuses divisions inutiles. Pour cela il utilise la technique du crible, d'où le nom du « crible quadratique ».

La méthode s'inspire du crible d'Ératosthène : imaginez une grille numérotée de 1 à 100. Partez du nombre 2 et rayez les nombres suivants en avançant de 2 en 2 : c'est-à-dire les nombres 4, 6, 8, 10, 12, etc. Puis partez du nombre 3 et rayez les nombres suivants en avançant de 3 en 3 : les nombres 6, 9, 12, 15, 18, etc. Partez du nombre 4 et rayez les nombres suivants en avançant de 4 en 4 : 8, 12, 16, etc.

En répétant la manipulation jusqu'au nombre 7, vous obtenez cette grille où comme par magie les nombres qui ne sont pas rayés (ici en orange) donnent la liste des nombres premiers jusqu'à 100 (hormis le 1 qui n'est pas un nombre premier).

Il n'est même pas nécessaire de connaître ses tables de multiplication.

Vous l'avez deviné, le crible quadratique de Carl Pomerance fait à peu près la même chose pour identifier rapidement les X² – N divisibles et ceux qui ne le sont pas.

Voici un exemple pour N = 48206621 :

Avec les deux améliorations précitées la liste initiale des premiers sera limitée à 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31 et seront éliminés les nombres 3, 13, 17, et 29 par le test de Legendre. Ne reste donc que 2, 5, 7, 11, 19, 23, et 31 dans la base des facteurs.

En partant de la racine carrée de N, 6944, et en ignorant les deux premiers nombres de la liste (nous verrons pourquoi plus loin) nous allons réaliser un crible pour chaque nombre premier, en commençant donc par le nombre 7, noté p ci-dessous.

La première chose à faire est de calculer ses deux « racines » a et b, afin de déterminer le décalage par rapport à l'origine du crible.

• Cela peut être calculé par la méthode de la force brute :

  • Soit T = N modulo p, c'est-à-dire T = 48206621 modulo 7 = 1.
  • Pour i = 1, tant que la racine de (p x i + T) n'est pas un entier alors incrémenter i.
  • Donc a = Racine(p x i + 1) soit a = Racine(7 x 5 + 1) = 6.
  • Et b = p – a, soit b = 7 – 6 = 1.

• Ou par l'algorithme de Tonelli-Shanks, qui est plus rapide (vous trouverez le code source en VBA dans le fichier joint).

Une fois ces racines définies, le crible peut être généré :

• le crible de la racine a part de l'origine (6944) avec un décalage de : p – (origine modulo p) + a, soit 7 – (6944 modulo 7) + 6 = 13. Comme 13 est supérieur ou égal à 7, alors on retient (13 modulo 7) = 6 ;
• le crible de la racine b part de 6944 avec un décalage de : p – (origine modulo p) + b = 8. Comme 8 est supérieur ou égal à 7, alors on retient (8 modulo 7) = 1.

En répétant ces opérations sur les autres nombres de la base des facteurs, vous obtiendrez ce résultat représenté sur une feuille de calcul pour une meilleure compréhension :

La racine a du nombre premier 7 a un décalage de 6 (cellule E9) et crible donc l'élément 6 (cellule M9), puis 13 (cellule T9), 20 (cellule AA9), etc. avec un pas de 7 en 7.

La racine b du nombre premier 7 a un décalage de 1 (cellule E10) et crible donc l'élément 1 (cellule H10), puis 8 (cellule O10), 15 (cellule V10), etc. avec un pas de 7 en 7.

Idem pour les autres nombres premiers avec un pas de p en p.

À quoi ça sert ?

J'ai reporté en ligne 6 la valeur de X² – N, vous constaterez que les cases du crible en orange sont des diviseurs de ces valeurs. Par exemple 26404 (colonne H) est divisible par 7 et par 23, mais pas par 11, 19 ou 31. Inversement 12515 (colonne G) ou 40295 (colonne I) ne sont pas divisibles par un premier de la liste.

Ainsi le crible permet de connaître très rapidement quels sont les X² – N qui peuvent être divisés par les premiers retenus (comme 26404) et inversement ceux pour lesquels il est inutile de perdre son temps à vouloir les factoriser (comme 12515 et 40295). Le tout sans avoir à calculer les X² – N et sans faire la moindre division !

Attention à cette nuance, car rien ne dit que les X² – N ainsi identifiés sont divisibles uniquement par les premiers retenus et forment donc une relation. Il faudra faire une factorisation classique pour s'en assurer (sans oublier les nombres 2 et 5 non repris dans le crible car, vous le comprenez maintenant, ils ne sont pas représentatifs), mais cela permet d'évincer des calculs un grand nombre de cas infructueux.

Et autre point fort de cette technique, j'ai indiqué en ligne 3 la « valeur de pondération » du crible, soit l'addition des logarithmes des diviseurs retenus. Par exemple l'élément 1 (cellule H1) qui vaut 8 est l'addition de log(7) et log(23), c'est-à-dire ses diviseurs. Sachant que la valeur maximale possible est l'addition des logarithmes de tous les premiers, soit 20, ici 8 représente 40 % de ce maximum (vous pouvez prendre la pondération de votre choix).

Dans la pratique nous évincerons aussi des calculs tous les X² – N du crible qui n'ont pas une pondération significative, et sont donc non prometteurs, quitte à perdre certaines bonnes relations mais peu importe puisque le crible est si rapide qu'obtenir des cas prometteurs n'est plus un problème.

En effet, si dans notre exemple une taille de crible de 1000 éléments n'était pas suffisante pour trouver assez de relations, le crible peut être relancé on prenant comme origine 6944 + 1000 (notez que le décalage de p ne sera même pas à recalculer car il suffit de récupérer la valeur de dépassement de la boucle précédente, i – taille du crible).

Il est possible aussi d'avoir une avance négative, en prenant comme origine 6944 – 1000. Dans ce cas il faudra ajouter à la matrice un facteur « x -1 » pour identifier ces relations. C'est très utile car plus l'on s'éloigne de la racine carrée de N et plus la probabilité de trouver une bonne relation est faible.

Comme nous l'avons vu, lorsque X² – N est factorisé avec les premiers de la liste, la relation est mémorisée. Quand on obtient suffisamment de relations, la recherche de la solution est lancée, en utilisant le pivot de Gauss pour trouver les relations à additionner, après avoir simplifié la matrice…

La matrice en modulo 2 de la factorisation de 48206621 est la suivante :

Les lignes à additionner pour trouver une solution ne sautent pas aux yeux.

Une étape intermédiaire au pivot de Gauss est de simplifier la matrice. Concrètement, cela consiste :

• à supprimer les colonnes où toutes les lignes valent zéro (ici les colonnes des diviseurs -1 et 2) ;
• à supprimer les colonnes où une seule ligne est à 1, en effet la solution doit provenir d'une addition modulo 2 donnant zéro, il est donc impossible qu'une colonne seule puisse être utilisée pour la solution. Et donc supprimer la colonne permet de supprimer la ligne qui l'utilise ;
• à regrouper sur une unique ligne deux lignes dont la colonne vaut 2 (ici la colonne du diviseur 7 vaut 2, donc soit les deux lignes sont incluses dans la solution, soit elles sont toutes les deux exclues car il est impossible de prendre l'une sans prendre l'autre).

Ce qui donne la matrice filtrée suivante :

Les lignes 13364 et 25299 ont fusionné pour donner la ligne 338095836, c'est-à-dire multiplication des « X », concaténation des « diviseurs X² – N », addition modulo 2 des matrices (ou un traitement XOR), ce qui permet de supprimer la colonne du diviseur 7.

Grâce à cette simplification on passe de la matrice d'origine constituée de 7 lignes x 8 colonnes, à une matrice de 6 lignes x 5 colonnes, soit environ deux fois moins de cases.

Ce n'est pas terminé, nous allons poursuivre la simplification de la matrice en procédant à son « échelonnement » …

Lors de mes recherches sur Internet (malheureusement j'ai perdu la source et ne peux citer son auteur) j'ai découvert par hasard une autre méthode de simplification de la matrice, dite « échelonnement », dont voici le principe.

Pour chacune des colonnes, on boucle sur les lignes de la colonne à la recherche du premier élément valant 1, qui devient la ligne de pivot. On boucle alors sur les autres colonnes dont la ligne du pivot vaut 1, et on additionne modulo 2 la matrice de la colonne de référence dans ces colonnes.

Une astuce permet de diminuer la « complexité » de cette méthode (qui est au pire le nombre de colonnes de la matrice au carré fois le nombre de lignes au carré) : pour ne pas avoir à additionner l'intégralité des lignes d'un diviseur avec les lignes de la colonne de référence, il suffit de mémoriser préalablement les lignes de la colonne de référence où la matrice vaut 1 et ne traiter que ces lignes :

Pour chaque colonne :

• Mémoriser les lignes où la matrice vaut 1, la première de ces lignes est le pivot.

• Si le pivot existe alors :

  • Pour chacune des autres colonnes :

    • Si la ligne pivot vaut 1 alors :

      • Pour les lignes mémorisées, faire l'addition 1 modulo 2.

Si l'on reprend la matrice filtrée de la page précédente, la 2^e ligne de la 1^re colonne est le pivot, et ses lignes 5 et 6 valent aussi 1. La 2^e ligne de la 2^e colonne vaut aussi 1, donc les lignes 2, 5 et 6 de la 2^e colonne sont additionnées 1 modulo 2, les trois autres colonnes n'ont pas de 1 en ligne 2 et ne sont pas concernées (on ne touche pas aux valeurs de « X » et aux « diviseurs de X² – N »), ce qui donne :

À la fin du traitement vous obtenez cette matrice qui n'a plus que 7 valeurs à 1 au lieu des 17 d'origine, et libère la colonne du diviseur 31 pour le pivot de Gauss :

Le temps passé à réaliser cet échelonnement sera en partie compensé par celui gagné à générer un pivot de Gauss sur une matrice ainsi condensée. Je ne peux pas le démontrer, mais j'ai l'intime conviction que cette simplification permet au pivot de Gauss de trouver plus facilement une bonne relation sur les grandes matrices.

Et la méthode étant aussi originale qu'efficace elle mérite, à mon avis, d'être mentionnée.

C'est pour ces raisons que je l'ai retenue.

Le fameux pivot de Gauss fonctionne sur un principe qui ressemble à l'échelonnement, avec cette fois une vision horizontale et non plus verticale.

Repartons de la matrice simplifiée après son échelonnement :

On part de la 1^re colonne de la matrice à la recherche de la 1^re ligne dont la valeur est 1, et qui devient la ligne de pivot. En principe la 1^re ligne est le pivot de la 1^re colonne. Ici ce n'est pas le cas puisque le pivot est en 2^e ligne, on procède alors à un échange entre ces deux lignes pour rétablir la situation.

Puis pour toutes les lignes sous le pivot, si la valeur de la colonne de référence vaut 1 alors on additionne la ligne du pivot à la ligne testée, c'est-à-dire multiplication des « X », concaténation des « diviseurs X² – N », addition modulo 2 des matrices (comme pour l'échelonnement, pour améliorer la rapidité des traitements, au lieu d'additionner toutes les colonnes il est préférable de mémoriser préalablement les colonnes de la ligne pivot qui valent 1 et de ne traiter que ces colonnes).

Ce qui donne ceci après le 1^er pivot :

Ce processus est répété sur toutes les autres colonnes de la matrice.

Le 2^e pivot devrait se trouver à la 2^e ligne de la 2^e colonne, ce qui nécessite encore ici un échange.

Puis l'addition de la ligne 2 dans les lignes 4 et 6 (car ces lignes valent 1 dans la 2^e colonne).

Finalement vous obtenez cette matrice après le traitement de toutes les colonnes :

Cette fois ça saute aux yeux, les deux dernières lignes sont des solutions possibles, car l'addition de leurs colonnes donne zéro.

Il reste juste à les tester.

La ligne 5 donne X = (338095836 x 24289) Modulo 48206621 = 11873254 ;
et Y = Racine(31 x 23 x 19…) modulo 48206621 = 11873254 ;
soit PGCD(X – Y, 48206621) = 48206621 ;
et PGCD(X + Y, 48206621) = 1 ;
ce n'est donc pas la solution.

La ligne 6 donne X = (60261 x 24289 x 53131) Modulo 48206621 = 23571483 ;
et Y = Racine(31 x 23 x 23…) modulo 48206621 = 47746801 ;
soit PGCD(X – Y, 48206621) = 9601 ;
et PGCD(X + Y, 48206621) = 5021 ;
et comme 9601 x 5021 = 48206621, c'est la solution.

De fait, si la méthode de factorisation du crible quadratique est très efficace, elle comporte juste un défaut : elle fait manipuler des nombres d'une grandeur impressionnante.

Démonstration

Dans notre exemple où il faut trouver les deux facteurs premiers de 48 206 621, qui est un tout petit nombre avec une matrice de factorisations qui ne fait que sept lignes, vous remarquerez que la racine carrée du produit des diviseurs de X² – N est déjà un nombre de 14 chiffres.

Mais la factorisation d'une clé de 56 chiffres génère une matrice de plus de 10 000 lignes, des relations complexes avec des centaines de diviseurs, dont le produit donne un nombre de plus de 100 000 chiffres !

Or certains langages de programmation ne permettent pas de manipuler de tels nombres avec les fonctions natives, comme le VBA qui est limité à 28 chiffres.

Et pourtant en VBA c'est possible : en programmant les opérations élémentaires (addition, soustraction, multiplication, puissances, etc.) et en travaillant avec des chaînes de caractères.

Vous trouverez dans le code source du fichier joint les principales fonctions nécessaires aux calculs sur les grands nombres en VBA dans le module « Gamma_Calculs_GN_Classiques », voir aussi http://fordom.free.fr/ pour d'autres calculs.

Ces fonctions ont déjà été utilisées dans mon tutoriel sur la cryptologie. Mais il manquait la division des grands entiers…

N'ayant rien trouvé de convaincant sur Internet, j'ai développé une fonction basée sur l'utilisation des logarithmes en m'inspirant en partie de cette discussion (je vous recommande la lecture de l'article passionnant de Wikipédia qui relate l'origine de l'invention des logarithmes, https://fr.wikipedia.org/wiki/Logarithme).

Vous trouverez le code source en VBA dans le module « Gamma_Calculs_GN_Autres ».

Si vous connaissez une méthode plus efficace je suis preneur.

L'épineux problème des calculs sur les grands nombres en VBA est ainsi réglé, mais en contrepartie les temps de traitement sont forcement bien plus longs qu'avec les langages qui les intègrent.

Surtout que le VBA est un langage « interprété », l'exécution est donc moins rapide qu'avec un langage « compilé ».

Alors si la factorisation en VBA d'une clé de 56 chiffres demande plus d'une heure, elle ne devrait prendre que quelques minutes dans un langage adapté.

Mais même si vous programmez en C, ne criez pas victoire trop vite et gardez à l'esprit que retrouver les deux facteurs premiers d'un nombre de 512 bits, soit 155 chiffres, demande des ressources énormes et nécessite l'usage de plusieurs centaines d'ordinateurs.

À ce jour aucune clé RSA de 1024 bits n'a officiellement été cassée.

Par prudence (excessive ?) certains recommandent l'usage de clés RSA de 2048 bits (617 chiffres).

L'autre contrainte rencontrée en VBA est que la taille maximale de la mémoire qu'il est possible d'allouer pour la matrice du pivot de Gauss est limitée à environ 15 000 lignes x 15 000 colonnes. Ce qui n'est pas suffisant pour la factorisation d'un nombre de 60 chiffres. Il faudra alors plafonner la base des facteurs et donc réduire l'efficacité du crible. Cette solution de contournement a ses limites et n'est plus efficace sur des nombres plus grands.

J'ai déjà signalé que plus le crible quadratique s'éloigne de la racine carrée du nombre à factoriser, plus la probabilité de trouver une bonne relation est faible. Dit autrement, plus l'on s'approche du but et plus il devient difficile à atteindre. L'idée ici est de donner un petit coup de pouce au crible quadratique dans la dernière ligne droite.

En factorisant un nombre de 50 chiffres, je me suis aperçu qu'environ 5 % des nombres de la base des facteurs n'avaient jamais été utilisés pour générer les bonnes relations retenues lorsque la collecte avait atteint 90 %. Inversement, ils sont parfois présents dans de mauvaises relations qui ont nécessité des divisions pour rien.

À croire que ces facteurs sont maudits.

Il semble donc peu probable qu'ils servent à trouver de bonnes relations dans les 10 % restants à collecter.

L'idée est alors de supprimer ces nombres « tabous » de la base des facteurs, donc du crible quadratique.

C'est facile, il suffit de leur attribuer une fausse valeur fortement négative à leur logarithme et ainsi forcer une « valeur de pondération » du crible qui sera toujours non significative lorsqu'un nombre tabou est présent.

De plus, la taille de la collecte peut en être diminuée d'autant.

En effet, la règle du pivot de Gauss qui veut qu'il y ait plus de lignes que de colonnes est respectée si l'on part du principe que les colonnes des nombres tabous seront vides, donc supprimées de la matrice.

Je me suis dit que :

• si la base des facteurs est plus petite, n'ont été supprimés que des nombres premiers jamais utilisés et non prometteurs, donc cela ne devrait pas ralentir la méthode ;
• voire au contraire cela pourrait l'accélérer, car certaines divisions inutiles seront épargnées ;
• et puisque la quantité de données à collecter diminue, c'est donc logiquement plus rapide ;
• enfin, la simplification de la matrice sera là aussi en principe plus rapide puisqu'il y aura un peu moins de lignes.

Finalement j'ai obtenu un gain de temps dans la collecte des données du nombre « 15 023 780 778 022 201 893 758 785 029 013 218 879 770 000 000 957 » d'environ 6 % (10min 17s contre 10min 57s, en VBA), de 33 % sur la résolution de la matrice (6min 40s contre 9min 58s), de 8 % sur la recherche de la solution.

On atteint au global près de 17 % de gains.

Voici en synthèse les étapes à suivre.

Génération de la base des facteurs, qui représente la liste des nombres premiers à retenir pour le traitement. Pour cela il faut : calculer le seuil de friabilité selon la formule officielle (éventuellement le doubler pour les petits nombres), trouver des nombres premiers inférieurs à ce seuil et vérifier par le test de Legendre si ces nombres premiers sont à inclure ou non dans la base des facteurs.

En partant de la racine carrée de N, boucle tant que la solution n'est pas trouvée (normalement elle est trouvée quand le nombre de relations collectées correspond à la taille de la base des facteurs).

Génération du crible quadratique (avance positive et négative).

Ne retenir du crible que les éléments significatifs (éventuellement, si le crible quadratique ne donne pas assez d'éléments significatifs, il est possible de baisser le taux de référence).

Pour chaque élément retenu :

• Vérifier que X² - N se divise uniquement avec les nombres premiers de la base des facteurs.
• Si c'est le cas alors mémoriser la relation obtenue (éventuellement si l'on a atteint 90 % de la collecte on peut désactiver de la base des facteurs les nombres jamais utilisés dans les bonnes relations, ce qui permet aussi de diminuer la quantité de données à collecter).

• Si le nombre de relations est au moins égal à la taille de la base des facteurs alors rechercher une solution :

  • simplification de la matrice ;
  • échelonnement de la matrice ;
  • pivot de Gauss sur la matrice.

• Boucle sur les lignes de la matrice dont l'addition des colonnes vaut zéro tant que la solution n'est pas trouvée :

  • calcul de X' = (produit des X) modulo N ;
  • calcul de Y' = Racine(produit des diviseurs) modulo N ;
  • calcul du PGCD(X'-Y', N) et du PGCD(X'+Y', N) pour savoir s'ils divisent N, c'est-à-dire s'ils sont la solution.

Comme indiqué au chapitre VII, avec les ressources matérielles dont je dispose et en utilisant le VBA pour programmer cet algorithme, la taille maximale de la matrice qui mémorise les relations des nombres premiers est d’environ 15 000 lignes x 15 000 colonnes. Ce qui n'est pas suffisant pour la factorisation d'un nombre de plus de 60 chiffres.
Début 2021, deux ans après la publication de cet article, je me suis replongé dans le code source avec pour objectif de trouver une astuce pour augmenter cette limite afin de pouvoir factoriser un nombre de 70 chiffres.

Et j’en ai trouvé trois :

• une pour augmenter le nombre de relations qui peuvent être mémorisées ;
• une pour accélérer le pivot de Gauss ;
• une pour calculer plus rapidement la solution finale.

- Augmenter le nombre de relations qui peuvent être mémorisées :
Reprenez le chapitre « II - Le principe de base » pour bien comprendre...
Dans la version d’origine, une variable d’un octet, soit 8 bits, est utilisée pour mémoriser le nombre de fois qu’un premier est utilisé pour factoriser « X²-N ».
Puis ce nombre est transformé en modulo 2 pour trouver les bonnes relations, donc un bit suffit pour représenter 0 ou 1.
Désormais, une variable de type long (32 bits) est utilisée, où chaque bit représente directement modulo 2 le nombre de fois qu’un premier est utilisé (en pratique, seuls 31 bits sont utilisés, car la variable est signée, il n’y a pas d’équivalence de la déclaration Unsigned long du C en VBA).
En reprenant l’exemple du chapitre, le Bit 0 de la variable correspond au « x2 », le bit 1 au « x3 », le bit 2 au « x5 » et ainsi de suite. Pour lire l’état du bit ou le changer, en VBA, il faut utiliser les opérateurs logiques AND, XOR ou OR.

Exemples avec ce tableau :

L’opération (45 AND 8) / 8 renvoie la valeur du bit 3, qui ici vaut 1, et (45 AND 2) / 2 renvoie 0.
L’opérateur XOR inverse la valeur d’un bit, soit l’équivalent d’une opération modulo 2.
45 XOR 8 donne 37 (le bit 3 passe à 0), et 37 XOR 8 donne 45 (le bit 3 passe à 1).
Pour forcer un bit à 1 il faut utiliser l’opérateur OR. 45 OR 16 = 61 (si le bit est déjà à 1 rien ne change, 45 OR 8 = 45).
Inversement, pour forcer un bit à 0, il faut dans un premier temps le forcer à 1 avec OR, puis l’inverser avec XOR.
Evidemment, si votre langage de programmation permet d’utiliser des variables sur un bit, ces opérateurs logiques ne sont pas utiles.
Travailler sur les bits de variables de type long au lieu d’utiliser des variables d’un octet permet de gagner assez de mémoire pour générer une matrice de 100 000 lignes x 100 000 colonnes au lieu de 15 000 x 15 000.


- Accélérer le pivot de Gauss :
Voir les chapitres « V - Échelonnement de la matrice des relations » et « VI - Le pivot de Gauss ». Lorsqu’en 2015 je présentais l’échelonnement (un traitement de simplification de la matrice qui précède le pivot de Gauss), je concluais en ces termes : « Je ne peux pas le démontrer, mais j'ai l'intime conviction que cette simplification permet au pivot de Gauss de trouver plus facilement une bonne relation sur les grandes matrices ». Mon intuition était bonne !
J’ai découvert qu’en faisant un échelonnement non pas sur la totalité de la matrice, mais uniquement sur les deux derniers tiers, le traitement est beaucoup plus rapide tout en restant aussi efficace.
Cela est dû au fait que les petits nombres premiers, ceux du premier tiers, sont très utilisés et rendent l’échelonnement très long. Sauter cette partie permet de minimiser les additions modulo 2 à faire et donc de diminuer considérablement la durée de ce traitement, sans en perdre son efficacité.
Le pivot de Gauss est alors simplifié et tourne beaucoup plus vite.
Finalement la résolution de la matrice, c’est-à-dire l’échelonnement et le calcul du pivot de Gauss, est dix fois plus rapide.
J'ai également remarqué qu'il est encore plus efficace de partir de la dernière colonne et de décrémenter jusqu'à 1/3 des colonnes.
Cette amélioration, non négligeable, est la plus importante que j’ai découverte. Et elle est valable pour tous les langages de programmation.


- Calculer plus rapidement la solution finale :
C’était pourtant évident, ça sautait aux yeux, et je n’ai rien vu !
Reprenez le chapitre « II - Le principe de base » une dernière fois...
La formule d’origine consiste à calculer le produit des « X », puis son modulo N et enfin son PGCD avec N ; calculer le produit des diviseurs de « X²-N » (noté Y) un sur deux (pour éviter la racine carrée) puis son modulo N, et enfin son PGCD avec N.
Cela donne des produits très grands, de plusieurs centaines de chiffres, et donc des temps de traitements très longs.
Il est bien plus simple de calculer directement le PGCD de chaque nombre et de repartir de ce résultat dans la suite du traitement.
Exemple avec de petits nombres pour comprendre le principe : soit les nombres suivants 200, 15 et 36, pour trouver le PGCD de ce produit avec 234, au lieu de calculer PGDC(108000, 234)=18, il est plus rapide de faire :
PGDC(200, 234)=2 puis PGDC(2*15, 234)=6 puis PGDC(6*36, 234)=18.
Ainsi, le calcul de la solution (sur de grands nombres) est cent fois plus rapide.


En conclusion :
Ces optimisations permettent, en théorie, de factoriser un nombre semi-premier de 70 caractères.
L’objectif est donc atteint, mais j’avoue ne pas avoir tenté l’expérience car d’après mes estimations, en VBA, cela me prendrait deux jours de traitement. Et c’est sans intérêt sachant qu’une simple calculatrice de lycéen donne le résultat en moins de 10 minutes. Si vous n’en avez pas une, vous pouvez télécharger gratuitement sur votre téléphone l’application « Calculatrice scientifique 115 es plus 991 ».
Le but ici est uniquement de comprendre comment marche le crible quadratique.
Une transcription en C serait plus causante en ce qui concerne la durée du traitement.

Il paraît que ceux qui étudient les méthodes de factorisation qui permettent de casser les clés RSA (qui sont la base du cryptage des informations confidentielles que nous échangeons tous les jours) seraient surveillés par les services de renseignements, mais rassurez-vous, vous ne risquez rien en lisant cette documentation, car comme indiqué dans l'introduction, l'algorithme présenté ici est une version personnelle confectionnée à l'aide d'informations glanées sur Internet qui ne contiennent rien de révolutionnaire, les professionnels utilisent des programmes plus complexes, font des traitements en parallèle sur plusieurs ordinateurs, et utilisent des langages de programmation bien plus rapides.

Et je le répète, le but ici est uniquement de comprendre comment marche le crible quadratique.

J'espère que cet objectif a été rempli.

Vous avez dû remarquer que cette documentation ne comporte aucune démonstration mathématique : c'est que je me voyais mal faire un copier/coller de notions que je ne comprends pas. Les forts en maths trouveront leur bonheur sur Internet avec la recherche « crible quadratique ».

Dans la pratique le crible quadratique est réservé à la factorisation de nombres d'au moins 18 chiffres, pour des nombres plus petits vous utiliserez l'algorithme « RhoPollard » donné dans le module VBA du fichier joint.

En cas de besoin d'explications sur le VBA vous pouvez vous référer à mon tutoriel : https://laurent-ott.developpez.com/tutoriels/programmation-excel-vba-tome-1/.

N'hésitez pas à partager vos codes si vous transposez cet algorithme du crible quadratique dans un autre langage de programmation, ou à nous indiquer les temps de traitement que vous obtenez.

Laurent OTT. 2019

CQGN_VBA.xlsm : (à ne plus utiliser) version VBA 32 bits pour EXCEL 2010 compatible 2016.

La feuille « Menu » permet de lancer la recherche de factorisation en deux nombres premiers d'une clé d'un maximum de 60 chiffres ; ou de reprendre le dernier traitement en cours (les données sont automatiquement sauvegardées tous les 250 cribles).

La combinaison des touches [MAJ][Fin] permet d'interrompre un traitement.

Les différents modules :

• Algorithme_CQGN : contient la fonction « CQGN » qui exécute la factorisation d'un nombre par l'algorithme du crible quadratique.
• Algorithme_Pollard : contient la fonction « RhoPollard » qui exécute la factorisation d'un nombre par l'algorithme « Rho de Pollard ».
• Algorithme_Tonelli_Shanks : contient la fonction « Tonelli_Shanks_GN » pour calculer les racines de N modulo un facteur premier.
• Formules_MillerRabin : contient l'algorithme du test de primalité de « Miller-Rabin ».
• Formules_MillerRabin_GN : permet de tester la primalité de grands nombres.
• Gamma_Calculs_GN_Autres : contient principalement les fonctions de calculs pour la division de grands nombres entiers.
• Gamma_Calculs_GN_Classiques : contient diverses fonctions pour les calculs sur des grands nombres, issues de http://fordom.free.fr/.
• Gamma_Calculs_PN : contient diverses fonctions pour les calculs sur des petits nombres.
• NP_Crible : contient un algorithme de type crible d'Ératosthène qui renvoie un million de nombres premiers.

CQGN_2021.xlsm : nouvelle version suite aux optimisations - ajout de juillet 2021. Cette version contient aussi l'algorithme ECM Lenstra pour factoriser un nombre.

Je remercie LittleWhite, Pierre Fauconnier, Chrtophe, pour l’intérêt qu'ils ont porté à cette documentation, pour leur relecture et leurs conseils.

Ainsi que Claude Leloup et escartefigue pour la correction orthographique, Winjerome pour la mise au gabarit et Malick pour la publication.

Sans oublier celles et ceux qui voudront apporter un commentaire...