Factorisation par courbe elliptique

Il y a quelques années, lorsque nous penchions sur le Green Code, nous avions trouvé une méthode de factorisation intéressante appelée communément appelé ECM (pour Elliptic Curve Method). Etant curieux, nous avons voulu percer les zones ombres de cette méthode. Notre article précédent nous a permis de définir les bases sur les courbes elliptiques. Cela nous sera utile pour cet article.

Hendrik Lenstra a, en 1985, découvert/inventé un algorithme de factorisation qui utilise des courbes elliptiques. Cet algorithme a la particularité de permettre de trouver  les petits facteurs premiers d’un nombre aussi grand qu’il soit.  Et là où la plus part des algorithmes de factorisation dépendent de la taille du nombre à factoriser celui-ce dépend principalement de la taille du facteur à trouver.

L’idée de la factorisation par courbe elliptique consiste à choisir une courbe elliptique modulo N (le nombre à factoriser) et de tenter de multiplier un point P plusieurs fois sur cette courbe.

Comme nous l’avons vu, pour le calcul de la somme de deux points sur une courbe elliptique il faut calculer la pente m. Cette pente contient une division qui en arithmétique modulaire n’est rien d’autre que le calcul d’un inverse. Or cet inverse n’existe pas toujours dans le cas où N est composé. Ainsi si notre addition de deux points échoue, nous avons probablement trouvé un diviseur de N.

Petit exemple avec 91

Plus concrêtement, essayons de factoriser le nombre  N=91 avec la courbe E(1,1) :

On choisit un point de la courbe : par exemple  P (37,37) est bien sur la courbe car :

37^2 \mod 91 = 1369 \mod 91 = 4
37^3+37+1 \mod 91 = 50653 +37 +1 \mod 91 = 50691 \mod 91 = 4

Puis on l’additionne à lui-même plusieurs fois : Au début tout se passe bien : P = (37,37) ;  2P = (56,15) ; 3P = (49,20) ; 4P = (86,40)
Mais pour 5P = 4P + P = (86,40)+(37,37), on doit calculer m =  \frac{40 -37}{86-37} \mod 91 = \frac{3}{49} \mod 91 Or 49 n’a pas d’inverse modulo 91.

Comme le calcul échoue, on sait que \gcd(91,49) > 1 on a donc trouvé un diviseur de N.

Notez que bien que le point 5P n’existe pas, le point 6P existe bel et bien sur la courbe :(23,58).

Les calculs échouent pour k = 5 ( et tous ces multiples) car 7 est un diviseur de 91 et que 5 est le nombre d’éléments de E(1,1) modulo 7.

Projetons-nous

Si l’on veut comprendre d’où vient ce lien, il faut comprendre les liens des courbes E(a,b) modulo 91 et modulo 7.

Chaque point de la courbe modulo 91 peut être « projeter » sur la courbe modulo 7 simplement en réduisant ces coordonnées modulo 7.

Tout d’abord remarquons que sur la courbe E(1,1) modulo 7, il n’y a que 5 éléments.

sage: Z = Zmod(7)
sage: E = EllipticCurve(Z, (1,1))
sage: E.points()
[(0 : 1 : 0), (0 : 1 : 1), (0 : 6 : 1), (2 : 2 : 1), (2 : 5 : 1)]
sage: E.cardinality()
5
sage:

Si l’on reprend le point P_{91} (37,37) sur la courbe modulo 91 et qu’on le réduit modulo 7 alors on tombe sur P_7 (2,2), mais cela fonctionne aussi avec les multiples de P : 2*P_{91}  = (56,15) qui se projete Q_{7}  = (0,1) = 2P_7.

Mais c’est Absurde

On peut continuer de vérifier les additions pour les multiples suivants. Vient alors le cas 5*P_{91} et raisonnons par l’absurde.

Si  nous supposons qu’un tel point existe, alors on pourrait calculer m donc ces coordonnées seraient finis. Et nécessairement la réduction modulo 7 de ces coordonnées le seraient aussi (ie celles de 5*P_{7} ).

Or on sait que l’ordre de cette courbe est 5 et donc d’après Lagrange, 5 fois n’importe quel point est égal à l’élément neutre. Dans le cas des courbes elliptiques, l’élément neutre a des coordonnées infinies.

D’où l’absurdité.

Quels paramètres choisir ?

Revenons au cas général. Si nous souhaitons utiliser cette technique pour factoriser un entier quelconque, plusieurs questions se posent :

Y-a-t-il un point P à choisir en particulier ?

La réponse à cette question est assez simple : pour l’instant nous n’avons pas de moyen de déterminer à l’avance si un point sera meilleurs qu’un autre. Le point est donc souvent choisi au hazard.

Quelle est la valeur optimale pour k (le nombre de fois qu’on multiplie le point P) ?

Pour répondre à cette question, posons-nous en une autre quelle sont les propriétés que nous souhaitons avoir pour k ? De plus, nous avons vu plus haut dans notre exemple que 5*P n’existait pas (ni ses multiples) mais que 6*P lui existait. Autrement dit, nous cherchons un k qui soit au moins égal à l’ordre de la courbe modulo le facteur cherché. Bien entendu, nous n’avons aucune idée de l’ordre en question.

Soyons naïfs :

Du coup une approche naïve pourrait être de prendre k = r!  avec r l’ordre maximum possible pour notre facteur :

Si nous cherchons les facteurs plus petit que 1000000 on sait que l’ordre maximum sera 1000000+2\sqrt(1000000)+1 = 1002001.

Mais ça fait quand même beaucoup…

Soyons attentifs !

Nous pourrions donc utiliser non plus la factorielle mais la primorielle. Cela réduit significitivement notre nombre. Cependant cela reste incomplet et nous pouvons « probablement » mieux faire.

Cela reste incomplet car si l’ordre pour le facteur chercher est par exemple 9. La multiplication par k = 11# (primorielle 11 = 2*3*5*7*11), ne donnera pas le résultat escompté car nous n’avons pas de puissance des nombres premiers.

Idéalement, il faut choisir une borne B et calculer le produit de toutes les puissances de nombres premiers inférieurs à B. Pour 27, nous aurions  k = 2^4*3^3*5*2*7*11*13*17*19*23. Ce nombre est aussi le plus petit commun multiple (ppcm ou lcm chez Shakespeare)  de tous les nombres entre 2 et 23 (le dernier nombre premier précédent 27)

Mais le nombre k est toujours très (trop) grand pour une utilisation en pratique

Par ailleurs, les lecteurs assidus de ce blog auront peut-être reconnu un problème se rapproche d’une notion déjà abordée ici. En effet, les plus attentifs auront probablement fait le lien avec la friabilité introduite lorsque nous parlions de l’algorithme p-1 de Pollard et l’optimisation de k.

Les valeurs de a et b varient chacune de 0 à N, Il y a un peu moins de N^2 (il faut enlever les cas où \Delta = 4*a^3 + 27*b^2 mod N = 0 ) courbes possibles. Mais leurs ordres se répartissent sur 4*\sqrt{N}. Nécessairement, il y a plusieurs courbes qui possèdent le même ordre.

Cependant on doit aussi s’intéresser aux ordres super-friables. Grâce à un petit script, nous obtenons les résultats suivants qui donnent une  tendance.

P Ordres min et max Nombre d’ordres différents Valeurs B-smooth avec B = Sqrt(P) Ratio
11 (5, 18) 13 0 0
101 (81, 122) 41 4 0.097
1009 (946, 1073) 127 22 0.173
10009 (9809, 10210) 401 83 0.206
100003 (99371, 100636) 1265 314 0.248
1000003 (998003, 1002004) 4001 1036 0.258
10000019 (9993695, 10006344) 12649 3371 0.266
100000007 (99980007, 100020008) 40001 10961 0.274
1000000007 (999936762, 1000063253) 126491 35247 0.278

En effet, il semble que plus P augmente, plus le nombre de nombre B-superfriable (avec B = sqrt(B)) augmente. Nos calculs montrent que si nous avons un P proche de 10^6 alors nous avons 25% des ordres sont 1000-superfriables. Cela donne une précision intéressante.

Pour la suite, soyons probabilistes !

En utilisant la notion de friabilité et le fait que nous pouvons tester plusieurs courbes, il est peut-être possible de drastiquement réduire le nombre k. En effet, avec un peu de chance, nous pouvons tomber sur une courbe d’ordre B-super-friable, avec un B de taille raisonnable.

Puisque le nombre total de courbes est très grand, il n’est pas possible d’effectuer des recherches sur toutes les courbes, nous devons nous limiter à effectuer des calculs sur un sous-ensemble des courbes mais combien précisement ?  Il y a deux aspects à prendre en compte : il nous faut d’un coté avoir suffisament de courbes pour avoir une chance (probabilité) raisonnable de trouver un facteur de la taille choisie et d’un autre coté, il nous faut avoir un minimum de courbe pour que les temps de calculs soient compétitifs par rapport aux autres méthodes de recherche de facteur.

Nous avons donc réalisé une implémentation de l’algorithme ECM en utilisant l’outils SageMath. Le script donc le code est ci-dessous nous a permis d’étudier les valeurs optimales du nombre de courbe et de B en fonction de la taille du facteur à trouver. Le résultat est donné dans le tableau ci-aprés.

Taille du        B    Nb Moyen           Nb de courbe pour les 10 tests
facteur                        Nb Médian    réalisés pour chaque B     

1.00E+05	20	1322.8	1086	[849, 3276, 1932, 73, 1323, 2471, 599, 414, 187, 2104]
1.00E+05	50	41.5	24	[72, 25, 3, 23, 13, 18, 23, 140, 65, 33]
1.00E+05	100	18.3	15.5	[15, 11, 35, 5, 16, 35, 9, 14, 20, 23]
1.00E+05	200	5.8	3.5	[13, 3, 1, 12, 14, 3, 1, 5, 4, 2]
1.00E+05	500	3.1	2.5	[3, 2, 9, 1, 2, 3, 5, 4, 1, 1]
1.00E+05	1000	2	1	[5, 1, 1, 1, 2, 1, 2, 5, 1, 1]
1.00E+05	1500	2.2	2	[2, 1, 1, 2, 4, 5, 2, 2, 2, 1]
1.00E+06	20	5263.1	5603.5	[416, 26, 10000, 10000, 2149, 10000, 7829, 5929, 1004, 5278]
1.00E+06	50	156.4	128.5	[376, 31, 97, 160, 69, 168, 97, 358, 28, 180]
1.00E+06	100	21.6	12.5	[14, 10, 6, 72, 11, 1, 1, 25, 36, 40]
1.00E+06	200	12.8	7	[29, 21, 1, 4, 15, 1, 8, 4, 6, 39]
1.00E+06	500	4	3.5	[4, 3, 4, 3, 2, 7, 2, 4, 2, 9]
1.00E+06	1000	3.5	3.5	[4, 3, 2, 5, 4, 6, 6, 2, 1, 2]
1.00E+06	1500	3.1	2.5	[5, 3, 1, 1, 2, 5, 7, 1, 2, 4]
1.00E+07	20	10000	10000	[10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000]
1.00E+07	50	1204.3	802	[183, 794, 930, 2831, 372, 810, 55, 1747, 3974, 347]
1.00E+07	100	210.1	204	[67, 315, 71, 373, 10, 8, 555, 275, 133, 294]
1.00E+07	200	32.7	14.5	[10, 18, 3, 123, 1, 35, 83, 32, 11, 11]
1.00E+07	500	15.8	12	[14, 14, 6, 2, 7, 40, 5, 43, 17, 10]
1.00E+07	1000	3.6	3	[1, 1, 3, 6, 6, 4, 3, 3, 3, 6]
1.00E+07	1500	6.7	5	[5, 1, 8, 16, 5, 13, 3, 5, 4, 7]
1.00E+08	20	10000	10000	[10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000]
1.00E+08	50	2408.1	1555	[2199, 1213, 2309, 1724, 3939, 1386, 10000, 506, 385, 420]
1.00E+08	100	460.3	521	[599, 603, 127, 203, 868, 549, 1020, 90, 493, 51]
1.00E+08	200	53.2	50	[9, 51, 14, 49, 147, 67, 91, 3, 84, 17]
1.00E+08	500	23.1	16.5	[8, 25, 6, 18, 83, 7, 16, 13, 38, 17]
1.00E+08	1000	15.6	15	[4, 38, 17, 15, 12, 3, 22, 27, 3, 15]
1.00E+08	1500	9.4	3.5	[6, 1, 1, 2, 1, 4, 54, 13, 9, 3]
1.00E+09	20	10000	10000	[10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000]
1.00E+09	50	10000	10000	[10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000]
1.00E+09	100	3262.9	2229	[2496, 1305, 3913, 272, 282, 9391, 2991, 1962, 10000, 17]
1.00E+09	200	284.2	186.5	[92, 631, 123, 799, 169, 399, 168, 204, 233, 24]
1.00E+09	500	43	43.5	[41, 46, 49, 19, 4, 56, 88, 5, 119, 3]
1.00E+09	1000	25.7	18	[85, 15, 3, 21, 9, 36, 35, 13, 1, 39]
1.00E+09	1500	18.3	11.5	[62, 8, 15, 5, 19, 5, 5, 30, 26, 8]
1.00E+11	200	4203.5	2883.5	[479, 6265, 10000, 10000, 491, 6260, 3211, 2556, 1825, 948]
1.00E+11	500	188.9	203	[263, 238, 93, 362, 22, 229, 342, 57, 177, 106]
1.00E+11	1000	168.2	137.5	[134, 144, 223, 230, 567, 84, 71, 141, 80, 8]
1.00E+11	1500	109	58.5	[8, 40, 319, 35, 252, 57, 160, 60, 37, 122]
1.00E+12	200	7652.1	10000	[10000, 10000, 1126, 1425, 10000, 10000, 6572, 10000, 10000, 7398]
1.00E+12	500	1156.9	872	[403, 1414, 3639, 954, 1297, 303, 487, 663, 1619, 790]
1.00E+12	1000	288.3	186	[165, 207, 72, 728, 140, 324, 398, 628, 125, 96]
1.00E+12	1500	82.4	81	[74, 12, 12, 253, 91, 93, 1, 35, 88, 165]
1.00E+13	200	9562.1	10000	[10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 5621]
1.00E+13	500	3742.1	3837.5	[1038, 8005, 5419, 4126, 4256, 1254, 3449, 5933, 3549, 392]
1.00E+13	1000	746.6	645	[577, 1545, 389, 713, 1412, 1110, 85, 522, 921, 192]

Ici notons que nous avons une approche très naïve. Nous utilisons beaucoup de fonction aléatoire pour générer nos points et nos courbes. Cela donne que des résultats bien loin de standard moderne dont le raisonnement est hors de notre portée. Mais les mathématiciens sont des gens plein de ressources et ont trouvé des moyens pour encore rafiner et  abaisser la borne B et le nombre de courbes en remarquant que certains courbes notamment ont des propriétés intéressantes.

Y-a-t-il une courbe (des coefficients a et b) à choisir en particulier ?

Cette dernière question est en effet importante. Les coefficients a et b permettent d’influer sur l’ordre de la courbe et comme le rappelle le théorème de Haas, l’ordre d’une courbe modulo p varie entre p+1 -2\sqrt{p} et p+1+2\sqrt{p}. L’expérience montre que l’ordre de la courbe est assez régulièrement réparti pour l’ensemble des valeurs du couple (a,b) mais que certains ordres ont plus de succes que d’autre. Ci-dessous nous avons fait une recherche exhaustive du nombre de courbe par ordre pour différentes valeurs de nombres premiers

17: [4, 8, 24, 8, 16, 24, 28, 8, 32, 8, 28, 24, 16, 8, 24, 8, 4] 17 17
53: [39, 26, 104, 52, 104, 130, 52, 26, 260, 52, 117, 104, 156, 78, 156, 78, 156, 104, 117, 52, 260, 26, 52, 130, 104, 52, 104, 26, 39] 29 29
79: [52, 156, 78, 156, 91, 156, 156, 312, 78, 390, 78, 156, 156, 364, 117, 156, 234, 390, 234, 156, 117, 364, 156, 156, 78, 390, 78, 312, 156, 156, 91, 156, 78, 156, 52] 35 35
101: [25, 50, 300, 100, 100, 250, 300, 100, 400, 150, 500, 200, 200, 200, 600, 150, 200, 400, 375, 100, 700, 100, 375, 400, 200, 150, 600, 200, 200, 200, 500, 150, 400, 100, 300, 250, 100, 100, 300, 50, 25] 41 41
127: [126, 63, 336, 147, 252, 252, 630, 189, 252, 378, 504, 315, 252, 126, 1008, 504, 378, 252, 504, 189, 756, 273, 630, 273, 756, 189, 504, 252, 378, 504, 1008, 126, 252, 315, 504, 378, 252, 189, 630, 252, 252, 147, 336, 63, 126] 45 45
173 [129, 86, 516, 86, 516, 602, 344, 258, 1032, 172, 516, 602, 1032, 430, 688, 430, 516, 860, 516, 258, 2064, 344, 559, 430, 860, 430, 1204, 430, 860, 430, 559, 344, 2064, 258, 516, 860, 516, 430, 688, 430, 1032, 602, 516, 172, 1032, 258, 344, 602, 516, 86, 516, 86, 129] 53 53
257: [64, 128, 896, 256, 768, 1024, 768, 256, 1024, 384, 1536, 896, 768, 512, 2816, 640, 512, 1280, 1792, 896, 2048, 384, 768, 640, 1536, 1024, 3072, 512, 512, 1664, 1984, 512, 2048, 512, 1984, 1664, 512, 512, 3072, 1024, 1536, 640, 768, 384, 2048, 896, 1792, 1280, 512, 640, 2816, 512, 768, 896, 1536, 384, 1024, 256, 768, 1024, 768, 256, 896, 128, 64] 65 65

Même si aucune formule actuellement ne permet de déterminer l’ordre d’une courbe en fonction de a,b et p, En autre des mathématiciens en utilisant des notions que nous détaillerons pas ici (comprendre qu’on n’a pas vraiment compris) ont trouvé des sous-ensembles de courbes dont l’ordre est au moins toujours divisibles par 12 (Paramétrage de Sumaya). Cela donne un certain avantage et permet de réduire le nombre de courbes à tester.

Mais l’algorithme a encore et toujours été raffiné : tant sur l’implémentation (pour gagner en vitesse de calcul) que sur les choix des courbes et que celui de k.

Ainsi les dernieres valeurs recommandées lors de l’utilisation de GMP-ECM (la référence pour la factorisation par courbes elliptiques) sont les suivantes :

Digits optimal B1 expected curves
(default parameters for GMP-ECM 6)
expected curves
(default parameters for GMP-ECM 7)
20 11,000 86 107
25 50,000 214 261
30 250,000 430 513
35 1,000,000 910 1,071
40 3,000,000 2,351 2,753
45 11,000,000 4,482 5,208
50 43,000,000 7,557 8,704
55 110,000,000 17,884 20,479
60 260,000,000 42,057 47,888
65 850,000,000 69,471 78,923
70 2,900,000,000 102,212 115,153
75 7,600,000,000 188,056 211,681
80 25,000,000,000 265,557 296,479

Voici quelques lectures intéressantes pour ce qui voudraient aller plus loin :

En espérant que vous avez réussi à suivre, sinon n’hésitez pas à laisser un commentaire .

 

Implémentation naïve d’ECM avec Sage:


#!/usr/bin/env sage

from numpy import median
from itertools import groupby
from sage.all import *

'''
Naive implemenation of ECM method
for test and education purposes
'''
def factorECM(k,m,verbose):
	condition1 = True
	condition2 = True
	Z = Zmod(m)
	curvecount = 0
	x = 0
	y = 1
	while curvecount  0):
			curvecount += 1
			#The curve is defined modulo M
			E = EllipticCurve(Z, (a,b))
			if(verbose):
				print "Curve", a ,b
			'''For some reasons, Sage does not allow us to
			use E.random_point() when M is composite
			so we have to compute a point on the curve
			We simply try various value of x and hope to find one
			'''
			while condition2 :
				x = ZZ.random_element(0,m)
				y2 = x**3 + a*x + b % m
				try:
					y = mod(y2,m).sqrt(extend=False) 		

				except:
					continue
				break	

			P = E(x,y)
			#We compute k*P	on catch the error if any, if not we'll try another curve
			try:
				Q = k*P
			except ZeroDivisionError as err:
				d = Integer(err.args[0].split()[2])
				if(verbose):
					print "After testing ", curvecount, " curve(s), we have a winner : "
					print "\t Point P :"  , P
					print "\t on E(",a,",",b,")"
					print "\t In K*P computation, d=", d ,"has no inverse modulo", m
					print "\t So we found one factor gcd(m,d) =", gcd(m,d)
				print "m =", m , " = ", gcd(m,d),"x", m/gcd(m,d)
				return curvecount

	if(verbose):
		print "No factor found with the given parameter : try increase B"
	return curvecount

def main(argv):
	if len(sys.argv) != 4:
	    print("Usage: %s     " % sys.argv[0])
	    print("Naive implementation of ECM for factorization")
	    print("v : Try to factor  with ")
	    print("q : Do some benchmark with random composite and predefine B")

	    sys.exit(1)

	#The number we want to factor
	m = sage_eval(sys.argv[2])

	#super-smooth border B
	B = sage_eval(sys.argv[3])

	if sys.argv[1] == 'v':
		#We compute the number k
		k = LCM(2..previous_prime(B))
		factorECM(k,m,True)
		return 0
	if sys.argv[1] == 'q':
		Blist = (200,500,1000,1500)
		for i in range (11,16):
			rg = 10**i
			for B in Blist:
				nbcurve = list()
				k = LCM(2..previous_prime(B))
				for l in range(10):
					p = next_prime(ZZ.random_element(rg,rg*10))
					q = next_prime(ZZ.random_element(rg*321,rg*14311))
					nbcurve.append(factorECM(k,p*q,False))
				print rg, B,
				print float(sum(nbcurve) / len(nbcurve)), float(median(nbcurve)),
				print nbcurve

if __name__ == "__main__":
   main(sys.argv[1:])
<span id="mce_SELREST_start" style="overflow:hidden;line-height:0;">&#65279;</span>

Publicités

A propos JoMendes

Amateur de mathématiques et d'hexadécimal. Je m'intéresse de près ou de loin suivant mon niveau à tous les sujets de sécurité de l'information.
Cet article, publié dans Cryptographie Asymétrique, Outils, RSA, est tagué , , , , , , , . Ajoutez ce permalien à vos favoris.

Laisser un commentaire

Entrez vos coordonnées ci-dessous ou cliquez sur une icône pour vous connecter:

Logo WordPress.com

Vous commentez à l'aide de votre compte WordPress.com. Déconnexion /  Changer )

Photo Google+

Vous commentez à l'aide de votre compte Google+. Déconnexion /  Changer )

Image Twitter

Vous commentez à l'aide de votre compte Twitter. Déconnexion /  Changer )

Photo Facebook

Vous commentez à l'aide de votre compte Facebook. Déconnexion /  Changer )

Connexion à %s