Algorithme ECM de factorisation par les courbes elliptiques

Algorithme ECM – Une invention de Hendrik Lenstra

Laurent Ott

Dans cet article, vous allez apprendre à programmer l’algorithme ECM de Hendrik Lenstra qui utilise les courbes elliptiques pour factoriser un nombre entier.
Cette méthode est bien adaptée à la recherche de « petits facteurs », car sa complexité ne dépend pas de la taille du nombre à factoriser, mais de la taille du plus petit de ses facteurs.

Vous pouvez déposer vos commentaires dans cette discussion.

- Titre : Algorithme ECM de factorisation par les courbes elliptiques

- Auteur : Laurent Ott

- Parution : 24 juin 2021

- Mise à jour : 24 juin 2021

- Licence : Le contenu de cet article est rédigé par laurent_ott et est mis à disposition selon les termes de la Licence Creative Commons Attribution - Pas d'Utilisation Commerciale 3.0 non transposé. Les logos Developpez.com, en-tête, pied de page, css, et look & feel de l'article sont Copyright ® 2013 Developpez.com.

Voir également :

- Tome 1 - Des bases de la programmation à l'algorithme de classement rapide QuickRanking

- Tome 2 - Des bases de la programmation en mode graphique à la programmation d'un jeu d'arcade en VBA et Microsoft Excel

- Tome 3 - Problème du Voyageur de commerce - Une méthode d'approximation basée sur des principes simples

- Tome 4 - Un algorithme de chiffrement/déchiffrement des cellules pour une meilleure confidentialité

- Tome 5 - Sentinelle - Une application qui veille sur vos classeurs sensibles - Exemples d'utilisations des tableaux de données et des requêtes SQL en VBA

- Tome 6 - Apprendre la programmation de nouvelles fonctions Excel pour l'utilisateur

- Comprendre la méthode de factorisation du Crible Quadratique

- Programmation en VBA de menus personnalisés pour Excel

- Manipuler les données des bases Access depuis Excel

- Transfert des fichiers volumineux sous Outlook

- Algorithme ECM de factorisation par les courbes elliptiques

- Un classeur Excel multi-utilisateur

- Compresser/décompresser des fichiers au format ZIP

- Fonctions en VBA pour gérer les Tableaux Structurés d’Excel

- Des fonctions en VBA pour générer des courriels par les balises HTML

- MsgBox_Perso : une boîte de message personnalisée

- Utiliser une collection ArrayList à plusieurs dimensions

- Fonctions en VBA pour gérer de gros volumes de données

- Les cryptosystèmes RSA et AES-256

Sommaire

- I. Avant-propos

- II. Introduction

- III. Présentation du principe de la méthode de factorisation par les courbes elliptiques

- III-A. Les règles de l’addition ECM

- III-B. Pseudocode de l’addition ECM

- III-C. Les règles de la multiplication ECM

- III-D. Pseudocode de la multiplication ECM

- IV. Présentation de l’algorithme ECM de Lenstra

- V. Conclusion

- VI. Des liens qui m’ont été utiles

- VII. Remerciements

- VIII. Annexe – Algorithmes des fonctions utilisées

- IX. Annexe – Algorithme ECM Lenstra en BASIC

- X. Annexe – Algorithme de factorisation utilisant ECM Lenstra

I - Avant-propos

Factoriser un nombre entier N consiste à retrouver les nombres premiers dont le produit donne N.
Par exemple 40 est composé de quatre facteurs premiers : 2 x 2 x 2 x 5. Qui s’écrit aussi : 2^3 x 5.

Au quotidien la décomposition en facteurs premiers permet, par exemple, de réduire une fraction pour la simplifier :

Ou de simplifier une racine :

Pour réaliser ces factorisations, de tête, j’ai utilisé une méthode qui consiste à tester la division de N par la liste des nombres premiers que je connais.
Méthode qui évidemment atteint vite ses limites…
Les mathématiciens ont donc recherché des procédés plus efficaces, entre autres celle dite ECM (sigle anglais de Elliptic Curve factorization Method), une invention de Hendrik Lenstra en 1985.

Notez que lorsqu’un entier ne se décompose qu’en deux facteurs, on dit qu’il est semi-premier. Comme 15 qui se factorise en 3 x 5.
Cette propriété est utilisée en cryptographie moderne, car s’il est facile de calculer le produit de deux grands nombres premiers, l’opération inverse pour retrouver les semi-premiers d’un grand nombre est beaucoup plus difficile.
Ce n’est pas le sujet de cette documentation, car l’algorithme ECM est destiné à des nombres de taille plus petite que ceux utilisés en cryptologie.

II - Introduction

Nous allons étudier dans cette documentation l’algorithme ECM de Lenstra, mais pas les mathématiques qui expliquent comment et pourquoi cette méthode de factorisation par les courbes elliptiques fonctionne, d’une part parce que cela nécessite des connaissances que je n’ai pas, d’autre part parce que vous trouverez sur Internet tout ce dont vous avez besoin.

A contrario lors de mes recherches, je n’ai rien trouvé qui décrive cet algorithme de façon facilement compréhensible pour une personne ayant de faibles notions en mathématiques.
Alors j’espère combler ici ce manque et permettre, à ceux qui le souhaitent, de programmer dans leur langage préféré cette méthode de factorisation très intéressante, car elle est, théoriquement, bien adaptée à la recherche de « petits facteurs » puisque sa complexité ne dépend pas de la taille du nombre à factoriser, mais de la taille du plus petit de ses facteurs.

III - Présentation du principe de la méthode de factorisation par les courbes elliptiques

« En mathématiques, une courbe elliptique est un cas particulier de courbe algébrique, munie entre autres propriétés d'une addition géométrique sur ses points ».

Par exemple cette courbe elliptique, en rouge, est la représentation de l’équation « y² = x^3 + 1 ».

Cette équation admet les points aux coordonnées suivantes : (2, 3), (0, 1), (–1, 0), (0, –1), (2, –3).

Vérifions :
2^3 + 1 = 3^2 ou (-3)2.
0^3 + 1 = 1^2 ou (-1)2.
-1^3 + 1 = 0^2.

Maintenant, prenez le point P1 sur cette courbe elliptique. Pour trouver P2, il suffit d’additionner P1 à lui-même. Soit P1 + P1 = P2.

Nous utiliserons la formule ci-dessous pour réaliser cette « addition », où (Px,Py) sont les coordonnées du premier point, (Qx,Qy) celles du second, et (Rx,Ry) les coordonnées du nouveau point.
P1 + P1 = P2 donne (Px,Py) + (Qx,Qy) = (Rx,Ry)

Soit ici, P1 + P1, avec (Px,Py) = (2,-3) et (Qx,Qy) = (2,-3) donne (0,-1) qui est bien le point P2.
Pour un troisième point, noté P3, nous additionnons P1 (2,-3) à P2 (0,-1), soit P1 + P2 = P3, et obtenons les coordonnées du point P3 du graphique (-1,0).
Pour le point P4, nous pouvons faire au choix, soit P1 + P3 = P4 soit P2 + P2 = P4, comme dans une addition classique : 1 + 3 = 4 ou 2 + 2 = 4. Et obtenons (0,1).
Et donc pour un cinquième point, P5, nous pouvons faire soit P1 + P4 = P5 soit P2 + P3 = P5.
Vous remarquerez que P6 n’existe pas (ce qui est normal) en effet P1 + P5 = (0,-1) soit les coordonnées du point P3.

Mais à quoi ça sert ?
Dans notre cas, avec l’équation « b = y² - x^3 – ax », cela va nous servir pour factoriser un nombre N, car une règle mathématique dit que dans une courbe elliptique de cette forme, il existe un point qui permet de trouver un diviseur de N.

À force d’additionner, vous parcourez les points de la courbe elliptique et finissez donc par tomber sur ce point particulier, c’est-à-dire trouver une solution à la factorisation de N.
Les règles d’addition de l’équation « b = y² - x^3 – ax » sont légèrement différentes de celles présentées ici.

III-A - Les règles de l’addition ECM

Pour additionner deux points « P » et « Q » sur une courbe elliptique de la forme b = y² - x^3 – ax, et déterminer ainsi un troisième point « R » (de coordonnées Rx et Ry), nous avons besoin des variables suivantes :

- les coordonnées du point P, notées ici Px et Py ;

- les coordonnées du point Q, notées ici Qx et Qy ;

- une variable « a », que l’on verra dans le chapitre suivant ;

- le nombre N que l’on cherche à factoriser.

La première étape consiste à calculer la « pente », que l’on notera « λ ».
Deux formules différentes sont à l’œuvre suivant que Px est égal ou non à Qx :

L’inverse modulo est le nombre qui donne le modulo 1 à N. Voir en annexe si votre langage de programmation n’a pas cette fonction.

Avant de calculer l’inverse modulo N, nous pouvons évaluer si le dividende ou le diviseur de la pente donnent un diviseur de N, c’est-à-dire une solution à notre problème :
soit X le dividende ou le diviseur ;
en calculant le plus grand commun dénominateur de X et N. X’ = PGCD(X, N) ;
si X’ est différent de 1, c’est qu’il divise N. Nous avons donc trouvé deux facteurs de N (X’ et N/X’) que nous retournons et le traitement est terminé.

La pente nous permet de calculer la deuxième étape, la valeur de Rx :
Rx = (λ² - Px – Qx) modulo N

Et cette valeur est elle-même utilisée dans la troisième et dernière étape pour calculer Ry :
Ry = (λ(Px – Rx) – Py) modulo N

Attention au piège de la gestion des modulos négatifs (et des inverses modulos négatifs) : si le résultat est négatif, il faut l’ôter de N.
Par exemple -25 modulo 20 = -5. Puis -5 + 20 = 15. Donc -25 modulo 20 = 15.
L’addition des points P et Q a soit constitué le point R de coordonnées Rx et Ry, soit trouvé une solution.

III-B - Pseudocode de l’addition ECM

Algorithme Add_ECM pour additionner deux points.

Entrées : Px, Py (les coordonnées du premier point), Qx, Qy (les coordonnées du second point), a, N (le nombre à factoriser).
Sorties : renseigne les variables Rx, Ry (nouveau point), retourne Vrai si une solution est trouvée.

Si Px = Qx alors
Diviseur = (3 * Px * Px + a)
Dividende = 2 * Py
Si non
Diviseur = Qy - Py
Dividende = Qx - Px
Fin Si

Test si le diviseur ou le dividende apportent une solution :
X = PGCD(Diviseur, N)
Si X <> 1 alors Rx = X et Ry = N / X et sortir de la fonction en retournant Vrai

X = PGCD(Dividende, N)
Si X <> 1 alors Rx = X et Ry = N / X et sortir de la fonction en retournant Vrai

Calcul de la pente :
S = Abs(Diviseur * InverseModulo(Dividende, N))
S = S Modulo N

Si le diviseur ou le dividende est négatif, alors il faut gérer un inverse modulo négatif :
Si (Diviseur < 0 et Dividende > 0) ou (Diviseur > 0 et Dividende < 0) alors S = N – S

Calcul des coordonnées du nouveau point :
Rx = (S * S - Px – Qx) modulo N
Si Rx < 0 alors Rx = Rx + N
Ry = (S * (Px - Rx) – Py) Modulo N
Si Ry < 0 alors Ry = Ry + N

Fin de la fonction

III-C - Les règles de la multiplication ECM

La multiplication d’un point P pour retourner le point R va s’appuyer sur une succession d’additions.

Par exemple pour multiplier le point « P » 59 fois :
nous additionnons P1 à lui-même pour avoir P2 ;
puis : P2 + P2 = P4, P4 + P4 = P8, P8 + P8 = P16, P16 + P16 = P32 ;
inutile d’additionner P32 + P32, car cela donne P64 qui est trop grand.
Notez que tous ces résultats sont mémorisés dans une table pour la suite du traitement qui va s’en servir.
P32 + P16 = P48
P48 + P8 = P56
P56 + P2 = P58
P58 + P1 = P59

Voilà, inutile de faire 58 additions, ce qui accélère les traitements.

Mais à quoi ça sert de faire une multiplication ?
Imaginez une formule magique qui détermine le point d’une courbe qui a une grande probabilité de retourner un diviseur : la multiplication permet d’accéder rapidement à ce point en évitant une cascade d’additions inutiles. Et donc de gagner du temps.

III-D - Pseudocode de la multiplication ECM

Algorithme Mult_ECM pour multiplier un point.

Entrées : Px, Py (les coordonnées du premier point), Qx, Qy (les coordonnées du second point), a, N (le nombre à factoriser), Nb (le nombre de fois qu’il faut multiplier).
Sorties : renseigne les variables Rx, Ry (nouveau point), retourne Vrai si une solution est trouvée.

Addition du premier point pour trouver le deuxième :
Si Add_ECM(Px, Py, Px, Py, a, N, Rx, Ry) = Vrai alors renseigner Rx et Ry et retourner Vrai.

Mémorisation des deux premiers points :
Mx(1) = Px: My(1) = Py: Mz(1) = 1
Mx(2) = Rx: My(2) = Ry: Mz(2) = 2
i = 2
k = 1

Boucle tant que l’on peut additionner les carrés :
Tant que Mz(i) * 2 < Nb
Si Add_ECM(Mx(i), My(i), Mx(i), My(i), a, N, Rx, Ry) = Vrai alors renseigner Rx et Ry et retourner Vrai
i = i + 1
Mx(i) = Rx: My(i) = Ry: Mz(i) = Mz(i - 1) * 2
k = Mz(i)
Fin de Tant que

Affine en additionnant avec la plus grande valeur possible :
Faire tant que k < Nb
i = i - 1
Si i < 1 alors i = 1
Si Mz(i) + k <= Nb alors
Qx = Rx: Qy = Ry
Si Add_ECM(Mx(i), My(i), Qx, Qy, a, N, Rx, Ry) = Vrai alors renseigner Rx et Ry et retourner Vrai
k = k + Mz(i)
Fin Si
Fin de Faire

Fin de la fonction

IV - Présentation de l’algorithme ECM de Lenstra

Pour vous présenter l’algorithme ECM de Lenstra, j’ai repris ci-dessous le pseudocode issu de ce lien :

Que je vais exprimer en d’autres termes.

Soit N le nombre à factoriser que vous passez en argument de la fonction ECM.

(1) Tirez au hasard trois nombres entiers compris entre 1 et N – 1 qui alimenteront les variables x, y et a (cette dernière étant la variable utilisée dans l’addition et la multiplication que nous avons étudiées dans les pages précédentes).

(2) Calculez la variable b avec b = y² - x^3 – ax.

(3) Calculez g qui est le plus grand commun dénominateur entre (4a^3 + 27b²) et N.

(4) Si g est égal à 1 alors le traitement continue, si g est égal à N alors retournez en (1), et donc si g est supérieur à 1 et inférieur à N, c’est que vous avez trouvé (par hasard) un diviseur de N et que vous pouvez terminer le traitement en retournant ce diviseur.

(9) Les variables x et y vont servir pour les calculs à venir, car elles représentent les coordonnées d’un point « A » sur une courbe elliptique.

(10-a) Calculez B un « seuil de friabilité », par exemple j’ai pris celui-ci arbitrairement :

B = Entier(Exp((1/2) x Racine(Log(N) x Log(Log(N)))))

Ou exprimé autrement, quand N vaut 1234567890 :

B = Entier(2,718281828459 Puissance(0,5 x Racine(Log(1234567890) x Log(Log(1234567890))))) = 54.

Je n’ai rien lu sur un éventuel calcul qui donnerait la taille optimale de B. Il semble qu’il est plus efficace d’avoir un petit B et relancer plusieurs fois la génération d’une courbe elliptique, que l’inverse.

(10-b) Bouclez sur les nombres premiers compris entre 2 et B inclus. On note p le nombre premier testé.

(11) Calculez e maximal tel que pe est inférieur ou égal à B.

Par exemple pour p = 2 et B = 54 alors p^5 = 32 est la valeur maximale, car p^6 dépasse 54 et donc e vaut 5.

(13) Multipliez le point « A », pe fois. Comme vu au chapitre précédent, cette multiplication permettra peut-être de trouver un facteur de N, ou dans le cas contraire, retourne un nouveau point.

C’est ce nouveau point qui devient le point de départ « A ».

(16) Lorsque tous les nombres premiers ont été testés, l’algorithme retourne qu’aucune solution n’a été trouvée.

La solution se trouve rarement au premier passage et vous repartirez en (1) plusieurs fois, éventuellement en augmentant B, tant qu’une solution n’est pas trouvée et que le temps que vous souhaitez y consacrer n’est pas écoulé.

V - Conclusion

L’algorithme ECM de Hendrik Lenstra (1985) qui utilise les courbes elliptiques pour factoriser un nombre entier est bien adapté à la recherche de « petits facteurs », car sa complexité ne dépend pas de la taille du nombre à factoriser, mais de la taille du plus petit de ses facteurs.

Il vient donc en complément d’autres algorithmes (par exemple RhoPollard P-1) à utiliser pour la factorisation d’un entier, pour tenter d’éliminer rapidement ses petits facteurs.

C’est pourquoi il me semblait utile de fournir de plus amples informations sur cet algorithme à ceux que le sujet intéresse et qui ne sont pas mathématiciens.

Surtout que cet algorithme est bien plus simple à programmer que l’algorithme du Crible Quadratique, qui est destiné à la factorisation de nombres semi-premiers plus grands.

Vous trouverez en annexe des informations complémentaires sur les algorithmes pour réaliser certains calculs que vos langages de programmation n’intègrent peut-être pas.
Ainsi qu’un exemple complet en BASIC de l’algorithme ECM suivi d’un exemple d’algorithme de factorisation permettant de retrouver tous les facteurs d’un nombre en utilisant ECM.

Bonne programmation.
Laurent Ott - 2021.

VI - Des liens qui m’ont été utiles

Sur Internet vous trouverez de nombreux articles sur la factorisation par les courbes elliptiques et l’ECM, j’ai retenu ceux-ci qui m’ont été utiles pour comprendre comment fonctionne l’algorithme.

Le pseudocode repris dans cette documentation est issu de ce lien.

Un article de vulgarisation : https://www.apprendre-en-ligne.net/crypto/rsa/251_088_096.pdf.

Ou cet article : http://images.math.cnrs.fr/Le-rang-des-courbes-elliptiques.html.

Un article pour les mathématiciens.

Une synthèse de plusieurs algorithmes.

Articles Wikipédia : Décomposition en produit de facteurs premiers, et Courbe elliptique.

VII - Remerciements

Je remercie dourouc05 et gaby277 pour leurs apports, et Claude Leloup pour la correction orthographique.

Ainsi que toute l’équipe de Developpez.com qui participe à la maintenance du site.

VIII - Annexe – Algorithmes des fonctions utilisées

Vous trouverez ici des algorithmes pour réaliser certains calculs que votre langage de programmation n’intègre peut-être pas.

Modulo
Le modulo retourne le reste d’une division :

'------------------------------------------------------------------------------------------------
Function Modulo(D, N)
'------------------------------------------------------------------------------------------------
Modulo = D - N * Int(D / N)
End Function

Inverse Modulo
L’inverse modulo (D, N) retrouve la valeur de x dans : Dx modulo N = 1

'------------------------------------------------------------------------------------------------
Function InversMod(D, N)
'------------------------------------------------------------------------------------------------
A = N: B = d: X = 1: Y = 0

While (B <> 0)
T = B
Q = Int(A / T)
B = A - Q * T
A = T
T = X
X = Y - Q * T
Y = T
Wend

InversMod = Y
If Y < 0 Then InversMod = Y + N
End Function

PGCD
Cette fonction retourne le plus grand commun dénominateur, ou 1 s’il n’y en a pas :

'------------------------------------------------------------------------------------------------
Function PGCD(NbA, NbB)
'------------------------------------------------------------------------------------------------
If NbA = 0 And NbB = 0 Then PGCD = 0: Exit Function 'Pas de solution
A = NbA: B = NbB
If NbA < NbB Then R = B: B = A: A = R 'inverse les valeurs

Do While Abs(B) >= 1
R = A - Int(A / B + 0.5) * B
A = B
B = R
Loop

PGCD = Abs(A)
End Function

IX - Annexe – Algorithme ECM Lenstra en BASIC

Voici une fonction en BASIC basée sur l’algorithme ECM Lenstra, que vous adapterez à votre langage.

Ici sont passés en arguments : N le nombre à factoriser, RetA et RetB deux variables qui contiendront les facteurs trouvés, et Nb_Passages le nombre de passages souhaité.

'---------------------------------------------------------------------------------------
Function ECM_Lenstra(N, RetA, RetB, Nb_Passages)
'---------------------------------------------------------------------------------------
Dim X, Y, A, B, d, g, i, Prime, p, Np

Prime = Array(2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271)

B = 271

Do
Do
' Prendre a, x, y au hasard:
Randomize Timer
i = Sqr(Sqr(N))
A = Int(Rnd() * i)
X = Int(Rnd() * i)
Y = Int(Rnd() * i)

d = (Y ^ 2) - (X ^ 3) - (A * X)
d = Modulo(d, N)

' Tester si c'est une solution:
g = PGCD(4 * A ^ 3 + 27 * d ^ 2, N)

If g > 1 And g < N Then
RetA = g: RetB = N / g
If RetA > RetB Then i = RetA: RetA = RetB: RetB = i
ECM_Lenstra = True ' Debug.Print RetA & " * " & RetB & " = "; RetA * RetB
Exit Function
End If

Loop While g = N

' Boucler sur une liste de nombres premiers:
For p = LBound(Prime) To UBound(Prime)
Np = Prime(p)
i = 1
While Np ^ (i + 1) < B
i = i + 1
Wend
i = Np ^ i
If Mult_ECM(X, Y, A, N, i, RetA, RetB) = True Then
If RetA <> 1 And RetB <> 1 Then
If RetA > RetB Then i = RetA: RetA = RetB: RetB = i
ECM_Lenstra = True ' Debug.Print RetA & " * " & RetB & " = " & RetA * RetB
Exit Do
End If
End If
X = RetA
Y = RetB
Next p

' S'il reste des passages à faire alors relancer la boucle:
Nb_Passages = Nb_Passages - 1

Loop While Nb_Passages > 0

End Function

'---------------------------------------------------------------------------------------
Function Add_ECM(Px, Py, Qx, Qy, A, N, Rx, Ry)
'---------------------------------------------------------------------------------------
Dim Diviseur, Dividende, X, S

If Px = Qx Then
Diviseur = (3 * Px * Px + A)
Dividende = 2 * Py
Else
Diviseur = Qy - Py
Dividende = Qx - Px
End If

X = PGCD(Diviseur, N)
If X <> 1 And X <> N Then
Rx = X: Ry = N / X: Add_ECM = True: Exit Function
End If

X = PGCD(Dividende, N)
If X <> 1 And X <> N Then
Rx = X: Ry = N / X: Add_ECM = True: Exit Function
End If

S = Abs(Diviseur * InversMod(Dividende, N))
S = Modulo(S, N)
If (Diviseur < 0 And Dividende > 0) Or (Diviseur > 0 And Dividende < 0) Then S = N - S

Rx = Modulo(S * S - Px - Qx, N)
If Rx < 0 Then Rx = Rx + N

Ry = Modulo(S * (Px - Rx) - Py, N)
If Ry < 0 Then Ry = Ry + N

End Function

'---------------------------------------------------------------------------------------
Function Mult_ECM(Px, Py, A, N, Nb, Rx, Ry)
'---------------------------------------------------------------------------------------
Dim Qx, Qy, X, S
Dim i, k
Dim Mx(1 To 9), My(1 To 9), Mz(1 To 9)

If Add_ECM(Px, Py, Px, Py, A, N, Rx, Ry) = True Then Mult_ECM = True: Exit Function
Mx(1) = Px: My(1) = Py: Mz(1) = 1
Mx(2) = Rx: My(2) = Ry: Mz(2) = 2
i = 2
k = 1

While Mz(i) * 2 < Nb
If Add_ECM(Mx(i), My(i), Mx(i), My(i), A, N, Rx, Ry) = True Then Mult_ECM = True: Exit Function
i = i + 1
Mx(i) = Rx: My(i) = Ry: Mz(i) = Mz(i - 1) * 2
k = Mz(i)
Wend

Do While k < Nb
i = i - 1
If i < 1 Then i = 1
If Mz(i) + k <= Nb Then
Qx = Rx: Qy = Ry
If Add_ECM(Mx(i), My(i), Qx, Qy, A, N, Rx, Ry) = True Then Mult_ECM = True: Exit Do
k = k + Mz(i)
End If
Loop

End Function

X - Annexe – Algorithme de factorisation utilisant ECM Lenstra

Toujours en BASIC, cette fonction utilise l’algorithme ECM Lenstra pour factoriser intégralement un nombre. La fonction est récursive pour décomposer les facteurs retournés tant qu’ils ne sont pas premiers. Il utilise à cet effet le test de primalité Miller Rabin. Vous noterez que le code fourni considère que 2 n’est pas un nombre premier (ce qui peut se comprendre puisqu’il est pair) alors qu’il est, dans notre cas, le plus petit facteur à retourner. D’où la correction opérée.

Passez en argument le nombre à factoriser.

La fonction affiche les facteurs premiers trouvés.

'---------------------------------------------------------------------------------------
Function Factorisation(N)
'---------------------------------------------------------------------------------------
Dim RetA, RetB, TestA, TestB

If N > 1 Then

If ECM_Lenstra(N, RetA, RetB, 100) = True Then

If RetA = 2 Then TestA = True Else TestA = TestMillerRabin(RetA)
If RetB = 2 Then TestB = True Else TestB = TestMillerRabin(RetB)

If TestA = True Then Debug.Print ">>> " & RetA
If TestB = True Then Debug.Print ">>> " & RetB

If TestA = False Then Call Factorisation(RetA)
If TestB = False Then Call Factorisation(RetB)

Else

If TestMillerRabin(N) = True Or N = 2 Then Debug.Print ">>> " & N

End If

End If

End Function

'---------------------------------------------------------------------------------------
Function TestMillerRabin(N)
'---------------------------------------------------------------------------------------
' Test de primalité MILLER RABIN. Retourne Vrai si n est un nombre premier
'---------------------------------------------------------------------------------------
Dim A, BaseMaxi
Dim Base

' Optimisation nb base à tester:
Base = Array(2, 3, 5, 7, 11, 13, 17, 19, 23)
BaseMaxi = Int(4 + Log(N) / 8 - Abs(4 - Log(N) / 8))

TestMillerRabin = False
If IsMultiple(N, 2) = True Then Exit Function

' Algo:
A = 0
Do
If IsMultiple(N, Base(A)) = False Then
If MillerRabin(N, Base(A)) = False Then Exit Function
End If
A = A + 1
Loop Until A > BaseMaxi

TestMillerRabin = True
End Function

'---------------------------------------------------------------------------------------
Function MillerRabin(N, A)
'---------------------------------------------------------------------------------------
Dim h As Long, d As Variant, i As Long

If Modulo(N, 2) = 0 Or N < 3 Then Exit Function

' Calcul de n-1 = 2^h * d avec d impair:
Do
h = h + 1
d = (N - 1) / 2 ^ h
Loop Until d = Int(d) And Modulo(d, 2) = 1

' Test a^d=1 mod n:
MillerRabin = (ExpoMod(A, d, N) = 1)
If MillerRabin Then Exit Function

' Teste s'il existe un 0 <= i <= h-1 tel que a^(h*2^i)=-1 mod n:
For i = 0 To h - 1
MillerRabin = (ExpoMod(A, d * 2 ^ i, N) = N - 1)
If MillerRabin Then Exit For
Next i

End Function

'---------------------------------------------------------------------------------------
Function ExpoMod(N, Exp, NbModulo)
'---------------------------------------------------------------------------------------
' EXPONENTIATION MODULAIRE RAPIDE : Nb^Expo MOD Modulo.
'---------------------------------------------------------------------------------------
Dim Nb, Expo, NModulo

Nb = N: Expo = Exp: NModulo = NbModulo
ExpoMod = 1
Do
If Modulo(Expo, 2) = 1 Then
ExpoMod = MODProd(Nb, ExpoMod, NModulo)
Expo = (Expo - 1) / 2
Nb = MODProd(Nb, Nb, NModulo)
End If

If Modulo(Expo, 2) = 0 Then
ExpoMod = MODProd(ExpoMod, 1, NModulo)
Expo = Expo / 2
Nb = MODProd(Nb, Nb, NModulo)
End If

Loop Until Expo = 0
End Function

'---------------------------------------------------------------------------------------
Function IsMultiple(Nb1, Nb2)
'---------------------------------------------------------------------------------------
' Teste si Nb1 est multiple de Nb2.
'---------------------------------------------------------------------------------------
If Nb2 = 0 Then
IsMultiple = True
Else
IsMultiple = ((Int(Nb1 / Nb2) = Nb1 / Nb2) And Nb1 <> 0)
End If
End Function

'---------------------------------------------------------------------------------------
Function MODProd(Nb1, Nb2, NbModulo)
'---------------------------------------------------------------------------------------
' Renvoie le modulo du produit "nb1*nb2 MOD Modulo"
'---------------------------------------------------------------------------------------
MODProd = Modulo(Nb1 * Nb2, NbModulo)
End Function