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.
[…] i am trying to understand monad. I wrote elliptic curve prime factorisation some times ago and now using monads, i wrote a bit improved elliptic curve prime factorisation . […]
Pingback by Improved Elliptic Curve prime factorisation « My Weblog | July 18, 2011 |
Was testing your code with 103695158367006753023 = 10183081967 x 10183081969. It doesn’t factor this number. ??
Comment by Kevin | April 1, 2013 |
@Kevin I think you are getting correct factorization result.
Comment by tiwari_mukesh | April 1, 2013 |
The above python code some times produce correct result and some times not. I will try to see what is going wrong in python. Try improved Elliptic curve factorization code.
Mukeshs-MacBook-Pro:Programming mukeshtiwari$ python Elliptic.py
103695158367006753023
10183081969
10183081967
Mukeshs-MacBook-Pro:Programming mukeshtiwari$ python Elliptic.py
103695158367006753023
103695158367006753023
Mukeshs-MacBook-Pro:Programming mukeshtiwari$ python Elliptic.py
103695158367006753023
10183081967
10183081969
Comment by tiwari_mukesh | April 1, 2013 |
Did you find out what may be causing this issue in python? I have followed your python example (can’t make out much out of monad version yet) and created an equivalent code in java. When it comes to larger factors I don’t get the proper factorization very often (like 1 out of 10 times). Is it possible that it has something to do with random numbers generation? Greetings and thanks for your effort.
Comment by apm | April 8, 2013 |
I think it’s because of bound on m ( m=int(math.factorial(2000)) ). If you numbers are bigger then try to increase the m by 4000 or 5000. Also try the Improved Haskell version. I think the only difference between python and Haskell code is m=int(math.factorial(2000)) and foldl lcm 1 [1..10000]. Hope it helps.
Comment by tiwari_mukesh | April 8, 2013 |