Lenstra’s Elliptic Curve Factorization Method

Lenstra’s Elliptic Curve Factorization Method

In December 2019, a team led by Paul Zimmermann of INRIA announced the factorization of the largest ever RSA modulus (RSA-240):

RSA-240 = 124620366781718784065835044608106590434820374651678805754818788883289666801188210855036039570272508747509864768438458621054865537970253930571891217684318286362846948405301614416430468066875699415246993185704183030512549594371372159029236099

RSA-240 = 509435952285839914555051023580843714132648382024111473186660296521821206469746700620316443478873837606252372049619334517
        × 244624208838318150567813139024002896653802092578931401452041221336558477095178155258218897735030590669041302045908071447

The factorization involved factorizing a 795-bit integer into its prime number factors. It caused industry experts to define that RSA would only be safe with at least 2,048-bit keys.

The security of several public key methods rely on the difficulty in factorizing a modulus created from the multiplication of large prime numbers. RSA is a good example of this, and where we take two large prime numbers (p and q), and multiply them to create a modulus (N). The difficulty is then to be able to find p and q, if we know N. The two core methods we can use for this factorization are the general number field sieve (GNFS) method and ECM (Elliptic Curve Method).

Hendrik W. Lenstra Jr.

Hendrik W. Lenstra Jr. become a Professor at the University of Amsterdam in 1979, and then, in 1987, he appointed to the University of California, Berkeley. One of his most famous PhD students is Daniel J. Bernstein [here], who is famous for producing standards such as Curve 25119, ChaCha20 and Poly1305. In the year he was appointed to Berkley, he outlined a method to factorize integrations using elliptic curve methods [1]:

ECM

With his method, we define the moving from a point P on an elliptic curve to 2P. For this we find the tangent to the point P and use this to find 2P. This tangent will be defined with a slope (s) and with a gradient in the form of a/b. For this we must find the modular inverse of b. If we cannot find the modular inverse of b, a factor of N is gcd(N,b), and where gcd() is the greatest common denominator.

If our curve is y²=x³+ax+b, the slope (s) will be:

s=(3x²+a)/(2y)

as defined in differentiation.

In detail, we first pick a random point P=(x0,y0) on an elliptic curve of:

y²=x³+Ax+B

We also select a random value of A and then calculate:

B=y0²−x0³−Ax (mod N)

For two points on the curve: P=(Px,Py) and Q=(Qx,Qy) , the slope of the line between them will be:

s=(Py−Qy)(Px−Qx)

With this we have s in the form of a/b(modN) .

Next we define:

R=P+Q=(Rx,−Ry)

and using:

s=(Px−Qx)(Py−Qy)(mod N)

Rx=s2−Px−Qx(mod N)

Ry=Py+s(Rx−Px)(mod N)

If Px=Qx and Py=−Qy, we define as 0 [Identity], else we calculate R=P+P=2P=(Rx,−Ry):

s=(3Px²+A)/(2Py) (mod N)

Rx=s2−2Px(modN)

Ry=P+s(Rx−Px) (mod N)

Here are some test runs:

• p=15 (Factor: 3 and 5 Try!
• p=6,161 (Factors: 61 x 101) Try!
• p=32,128 (Factors: 64 x 502) Try!
• p=53,421 (Factor: 3 x 17,807) Try!
• p=55,440 (Factors: 2 x 2 x 2 x 2 x 3 x 3 x 5 x 7 x 11) Try!
• p=100,001 (Factor: 3 and 17,807) Try!
• p=999,999 (Factor: 3 x 3 x 3 x 7 x 11 x 13 x 37) Try!
• p=455,839 (Factor: 559 and 761) Try!

Coding

An outline of the code used is [code]:

#!/usr/local/bin/python
# -*- coding: utf-8 -*-
import math
import random 
import sys

#y^2=x^3+ax+b mod n 

prime=[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 ]

# ax+by=gcd(a,b). This function returns [gcd(a,b),x,y]. Source Wikipedia
def extended_gcd(a,b):   
	x,y,lastx,lasty=0,1,1,0
	while b!=0:
		q=a/b
		a,b=b,a%b
		x,lastx=(lastx-q*x,x)
		y,lasty=(lasty-q*y,y)
	if a<0: 
		return (-a,-lastx,-lasty)
	else:
		return (a,lastx,lasty)

# pick first a point P=(u,v) with random non-zero coordinates u,v (mod N), then pick a random non-zero A (mod N), 
# then take B = u^2 - v^3 - Ax (mod N).
# http://en.wikipedia.org/wiki/Lenstra_elliptic_curve_factorization

def randomCurve(N):
	A,u,v=random.randrange(N),random.randrange(N),random.randrange(N)
        B=(v*v-u*u*u-A*u)%N
        return [(A,B,N),(u,v)]

	# Given the curve y^2 = x^3 + ax + b over the field K (whose characteristic we assume to be neither 2 nor 3), and points
	# P = (xP, yP) and Q = (xQ, yQ) on the curve, assume first that xP != xQ. Let the slope of the line s = (yP - yQ)/(xP - xQ); since K 	       # is a field, s is well-defined. Then we can define R = P + Q = (xR, - yR) by
	#	s=(xP-xQ)/(yP-yQ) Mod N		
	#	xR=s^2-xP-xQ	Mod N
	#	yR=yP+s(xR-xP)	Mod N 
	# If xP = xQ, then there are two options: if yP = -yQ, including the case where yP = yQ = 0, then the sum is defined as 0[Identity]. 	       # thus, the inverse of each point on the curve is found by reflecting it across the x-axis. If yP = yQ != 0, then R = P + P = 2P = 	      # (xR, -yR) is given by
	#	s=3xP^2+a/(2yP)	Mod N	
	#	xR=s^2-2xP	Mod N
	#	yR=yP+s(xR-xP)	Mod N
	#	http://en.wikipedia.org/wiki/Elliptic_curve#The_group_law''')
	
def addPoint(E,p_1,p_2):
	if p_1=="Identity": return [p_2,1]
	if p_2=="Identity": return [p_1,1]
	a,b,n=E
	(x_1,y_1)=p_1
	(x_2,y_2)=p_2
	x_1%=n
	y_1%=n
	x_2%=n
	y_2%=n
	if x_1 != x_2 :
		d,u,v=extended_gcd(x_1-x_2,n)
		s=((y_1-y_2)*u)%n
		x_3=(s*s-x_1-x_2)%n
		y_3=(-y_1-s*(x_3-x_1))%n
	else:
		if (y_1+y_2)%n==0:return ["Identity",1]
		else:
			d,u,v=extended_gcd(2*y_1,n)
			s=((3*x_1*x_1+a)*u)%n
			x_3=(s*s-2*x_1)%n
			y_3=(-y_1-s*(x_3-x_1))%n

	return [(x_3,y_3),d]

	# http://en.wikipedia.org/wiki/Elliptic_curve_point_multiplication
	#	Q=0 [Identity element]
	#	while m:
	#		if (m is odd) Q+=P
	#		P+=P
	#		m/=2
	#	return Q')

def mulPoint(E,P,m):
	Ret="Identity"
	d=1
	while m!=0:
		if m%2!=0: Ret,d=addPoint(E,Ret,P)
		if d!=1 : return [Ret,d]  # as soon as i got anything otherthan 1 return 
		P,d=addPoint(E,P,P)
		if d!=1 : return [Ret,d]
		m>>=1
	return [Ret,d]

def ellipticFactor(N,m,times=5):
	for i in xrange(times):
		E,P=randomCurve(N);
		Q,d=mulPoint(E,P,m)
		if d!=1 : return d
	return N

first=True
n=455839
if (len(sys.argv)>1):
        n=int(sys.argv[1])

print n,'=',

for p in prime:
	while n%p==0:
		if (first==True):
			print p,
			first=False
		else:
			print 'x',p,
		n/=p

m=int(math.factorial(2000))

while n!=1:
	k=ellipticFactor(n,m)
	n/=k
	if (first==True):
		print k,
		first=False
	else:
		print 'x',k,

Conclusions

The great public key methods of the past, may not scale well into the 2020s, especially as factorization becomes more powerful, and is actually beating Moore’s Law. If you are interested, here are other methods for factorization:

• Difference of squares. Diffsq. This factorizes using the difference of squares method.
• Factors of integers. Factors. Determine factors of an integer.
• Pollard’s ρ method (Factoring integers). Pollard. The Pollard ρ method factorises integers.
• Simplify ap(modN) Go. Simplify the a^p (modN) operation.
• Smooth numbers. Go. Outline of smooth numbers.
• Quadratic residue (mod p). Go. Outline of quadratic residue (mod p).
• Quadratic residue (mod N) and where N=pq. Go. Outline of quadratic residue (mod N) with N made-up of two prime numbers.
• Dixon. Go. Dixon Method.

References

[1] Lenstra Jr, H. W. (1987). Factoring integers with elliptic curves. Annals of mathematics, 649–673.