My Weblog

Blog about programming and math

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.

January 17, 2011 - Posted by | Programming

6 Comments »

  1. […] 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 | Reply

  2. Was testing your code with 103695158367006753023 = 10183081967 x 10183081969. It doesn’t factor this number. ??

    Comment by Kevin | April 1, 2013 | Reply

    • @Kevin I think you are getting correct factorization result.

      import Data.Bits
      import Control.Parallel.Strategies
       
      powM :: Integer -> Integer -> Integer ->  Integer
      powM a d n = powM' a d n where
        powM' a d n
              | d == 0 = 1
              | d == 1 = mod a n
              | otherwise = mod q n  where
                 p = powM'  ( mod ( a^2 ) n ) ( shiftR d 1 ) n
                 q = if (.&.) d 1 == 1 then mod ( a * p ) n else p
       
      calSd :: Integer -> ( Integer , Integer )
      calSd n = ( s , d ) where
            s = until ( \x -> testBit ( n - 1) ( fromIntegral x ) )  ( +1 ) 0
            d = div ( n - 1 ) (  shiftL 1 ( fromIntegral s )  )
       
       
      rabinMiller::Integer->Integer->Integer->Integer-> Bool
      rabinMiller  n s d a
         | n == a = True
         | otherwise =  case x == 1 of
                True -> True
                _ ->   any ( == pred n ) . take ( fromIntegral s )
                            . iterate (\e -> mod ( e^2 ) n ) $ x 
              where
                    x = powM a d n
          
       
      isPrime::Integer-> Bool
      isPrime n
         | n <= 1 = False
         | n == 2 = True
         | even n = False
         | otherwise  = all ( == True ) . map ( rabinMiller n s d ) $  [ 2 , 3 , 5 , 7 , 11 , 13 , 17 ] where
                      ( s , d ) = calSd n 
      
      *Main> isPrime 10183081967
      Loading package containers-0.5.0.0 ... linking ... done.
      Loading package parallel-3.2.0.3 ... linking ... done.
      True
      *Main> isPrime 10183081969
      True
      *Main> 
      
      

      Comment by tiwari_mukesh | April 1, 2013 | Reply

  3. 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 | Reply

  4. 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 | Reply

    • 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 | Reply


Leave a reply to Kevin Cancel reply