Elliptic Curve prime factorisation
Lenstra algorithm explains about Lenstra elliptic curve prime factorisation and wiki page about Elliptic curve is good enough to start about elliptic curves.
Here is algorithm which i borrowed from Primality Testing and Integer Factorization in Public-Key Cryptography . Its very good book specially for those who wants to implement the different algorithm about prime numbers. This implementation is highly influenced by Basic Number Theory in Haskell. I borrowed multiEC from this library. Excellent library of number theory in Haskell.
[1] (Choose an Elliptic Curve) Choose a random pair (E, P ),
where E is an elliptic curve y^2 = x^3 +ax+b over Z/nZ, and
P (x, y) ∈ E(Z/nZ) is a point on E. That is, choose a, x, y ∈ Z/nZ at random, and set b ← y^2 − x^3 − ax. If gcd(4a^3 + 27b^2 , n) = 1, then E is not an elliptic curve, start all over and choose another pair
(E, P ).
[2] (Choose an Integer k) Just as in the “p−1” method, select a positive integer k that is divisible by many prime powers, for example, k = lcm(1, 2, . . . , B) or k = B! for a suitable bound B; the larger B is the more likely the method will succeed in producing a factor, but the longer the method will take to work.
[3] (Calculate kP ) Calculate the point kP ∈ E(Z/nZ). We use the following formula to compute
P3 (x3 , y3 ) = P1 (x1 , y1 ) + P2 (x2 , y2 ) modulo n:
(x3 , y3 ) = (s^2 − x1 − x2) mod n, s(x1 − x3 ) − y1 mod n
where
s=3x1^2+a/2y1 mod n if P1=P2
s=y1-y2/x1-x2 mod n
The computation of kP mod n can be done in O(log k) doublings and
additions.
[4] (Computing GCD) If kP ≡ OE (mod n), then compute d = gcd(m2 , n), where m2=2y1 or x1-x2 else go to Step [1] to make a new choice for “a” or even for a new pair (E, P ).
[5] (Factor Found?) If 1 < d < n, then d is a nontrivial factor of n, output d and go to Step [7].
[6] (Start Over?) If d is not a nontrivial factor of n and if you still wish to try more elliptic curves, then go to Step [1] to start all over again, else go to Step [7].
[7] (Exit) Terminate the algorithm.
--y^2= x^3+ax+b mod n
import Random
data Elliptic = Conelliptic Integer Integer Integer deriving (Eq,Show)
data Point = Conpoint Integer Integer | Identity deriving (Eq,Show)
--add points of elliptic curve
addPoints::Elliptic->Point->Point-> Either Point Integer
addPoints _ Identity p_2 = Left p_2
addPoints _ p_1 Identity = Left p_1
addPoints (Conelliptic a b n) (Conpoint x_p y_p) (Conpoint x_q y_q)
| x_p/=x_q = let
[u,v,d]= extended_gcd (x_p-x_q) n
s=mod ((y_p-y_q)*u) n
x_r=mod (s*s-x_p-x_q) n
y_r=mod (-y_p-s*(x_r-x_p)) n
in case ((Conpoint x_r y_r),d) of
(_,1)->Left (Conpoint x_r y_r)
(_,d')->Right d'
| otherwise = if mod (y_p+y_q) n ==0 then Left Identity
else let
[u,v,d]= extended_gcd (2*y_p) n
s=mod ((3*x_p*x_p+a)*u) n
x_r=mod (s*s-2*x_p) n
y_r=mod (-y_p-s*(x_r-x_p)) n
in case ((Conpoint x_r y_r),d) of
(_,1)->Left (Conpoint x_r y_r)
(_,d')->Right d'
extended_gcd::Integer->Integer->[Integer]
extended_gcd a b=
let
[u,v,d]=exgcd a b
where
exgcd a b | mod a b == 0 = [0,1,b]
| otherwise = [y,x-y*(div a b),z]
where
[x,y,z] = exgcd b (mod a b)
in if d<0 then [-u,-v,-d] else [u,v,d]
multiEC::Elliptic->Point->Integer->Either Point Integer
multiEC _ _ 0 = Left Identity
multiEC ecurve point k | k>0 = fun Identity point k
where
fun p _ 0 = Left p
fun p q n = let
p'=if odd n then addPoints ecurve p q else Left p
q'= addPoints ecurve q q
in case (p',q') of
(Left p'',Left q'')-> fun p'' q'' (div n 2)
(Right p'',_)-> Right p''
(_,Right q'')->Right q''
dscrmntEC a b = 4 * a * a * a + 27 * b * b
randomCurve ::Integer-> IO [Integer]
randomCurve n= do
a<-getStdRandom (randomR (1,n))
u<-getStdRandom (randomR (1,n))
v<-getStdRandom (randomR (1,n))
let b=mod (v*v-u*u*u-a*u) n
return [a,b,n,u,v]
factor::Integer->Integer->Integer->IO [Integer]
factor 1 _ _ = return []
factor n m cnt | cnt>=5 = return [n]
| otherwise = do
[a,b,n,u,v]<-randomCurve n
let p=multiEC (Conelliptic a b n) (Conpoint u v ) m
case p of
Left p' -> factor n m (cnt+1)
Right p'-> do
ps<-factor (div n p') m 0 --start again factoring
ls<-factor p' m 0 --start again factoring
return (ls++ps)
main = do
line<-getLine
let n=read line ::Integer
m<-factor n (foldl lcm 1 [1..10000]) 0
print m
return ()
Here is my python implementation which is faster than Haskell. Kind of surprising for me. When i started learning haskell then some where it was given that because of type checking haskell is faster than python.
#!/usr/local/bin/python
# -*- coding: utf-8 -*-
import math
import random
#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=random.randrange(N)
u=random.randrange(N)
v=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]
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
if __name__=="__main__":
n=input()
for p in prime:#preprocessing
while n%p==0:
print p
n/=p
m=int(math.factorial(2000))
while n!=1:
k=ellipticFactor(n,m)
n/=k
print k
It could be possible that these codes may have indentation problem due to formatting and both languages are indentation sensitive.