My Weblog

Blog about programming and math

Elliptic curve prime factorisation using Montgomery curve

This is almost similar to previous example except using Montgomery curve rather than Lenstra curve. Programmingpraxis used Montgomery curve to find the prime factorisation. I personally feel that Lenstra algorithm worked fine when factors were small how ever Montgomery return the number as prime. Run this program against a small input (digits 15) and chances are that you will get small digit as itself prime and large input as prime factor. This code uses some bitwise operations , upper bound of number base 2 and Montgomery Ladder algorithm.

--y^2 = x^3+ax^2+x mod n
import Random
import System.Random
import Data.Bits
data Elliptic = Cone Integer Integer Integer deriving (Eq,Show) 
data Point = Conp Integer Integer  deriving (Eq,Show)

addPoints (Conp 0 0) p2 _ _ = p2
addPoints p1 (Conp 0 0) _ _ = p1
addPoints (Conp x_1 z_1) (Conp x_2 z_2) (Conp x_0 z_0) n = 
		x_3= mod ((((x_1-z_1)*(x_2+z_2)+(x_1+z_1)*(x_2-z_2))^2)*z_0) n 
		z_3= mod ((((x_1-z_1)*(x_2+z_2)-(x_1+z_1)*(x_2-z_2))^2)*x_0) n
	in (Conp x_3 z_3)

doublePoint (Conp x z) (Cone an ad n) = Conp (mod x' n) (mod y' n)
					   x'= mod (4*ad*(x+z)^2*(x-z)^2) n
					   y'= mod ((4*ad*(x-z)^2+t*an)*t) n 

ilognew b n = f (div ubound 2)  ubound where
    ubound = until (\e -> b ^ e > n) (* 2) 1 
    f lo hi | lo>=hi = div (lo+hi) 2 
            | 2 ^ mid <= n = f (mid+1) hi
            | otherwise = f lo mid
            where mid = div (lo + hi) 2

montMul _ _ _ 0= (Conp 0 0)
montMul n eCurve pOint k = tmpFun  (Conp 0 0)  pOint (ilognew 2 k  -1) where 
	tmpFun r0 r1  0 = if (.&.) k ( shiftL (1::Integer) 0 ) /=0 then  addPoints r0 r1 pOint n  
			  else  doublePoint r0 eCurve
	tmpFun r0 r1  cnt |(.&.) k ( shiftL (1::Integer) cnt ) /=0  = 
			    tmpFun (addPoints r0 r1 pOint n) (doublePoint r1 eCurve)  (cnt-1)
			  |otherwise =  
			    tmpFun (doublePoint r0 eCurve) (addPoints r0 r1 pOint n)  (cnt-1)

montgomeryCurve::Integer-> IO [Integer]
montgomeryCurve  n = do  
			s<-getStdRandom (randomR (6,n-1))
			let u=mod (s*s-5) n
			let v=mod (4*s) n
			let an=mod (((v-u)^3)*(3*u+v)) n
			let ad=mod (4*u^3*v) n
		        return [an,ad,mod (u*u*u) n,mod (v*v*v) n]

prime = 2:3:5:filter isPrime [6*k+r | k<-[1..],r<-[1,5]]
isPrime n = all (notDiv n ) $ takeWhile (\p->p*p<=n) prime where 
		notDiv n p = mod n p /=0

allp p_1 p_2 = takeWhile (<=p_2) $ dropWhile (<p_1) prime

ellipticFac  b_1 b_2 (Cone an ad n)  (Conp u v) = 
		c= foldl lcm 1 [1..b_1]
		(Conp x z)= montMul n (Cone an ad n) (Conp u v) c
		g=gcd n z
	in case (and [1<g, g<n]) of 
		True -> g
		_ ->   
		   (Conp x' z')= foldl ( montMul n (Cone an ad n) ) (Conp x z) (allp b_1 b_2)
		   g'=gcd (mod (g*z') n) n
	         in case (and [1<g',g'<n]) of 
		     True -> g'
		     _ -> n

factor::Integer->Integer->IO [Integer]
factor 1 _ = return []
factor n cnt | cnt>=5 = return [n]
	     | otherwise = do 
				[an,ad,u,v] <-montgomeryCurve n
				let p=ellipticFac 10000 20000 (Cone an ad n) (Conp u v)
				case (and [1<p,p<n]) of 
					False -> factor n (cnt+1)
					_ -> do 
						ps<-factor (div n p) (0)
						ls<-factor p (0)
						return (ls++ps)
main = do
		let n=read line ::Integer
		lst<-factor n 0
		print lst
		return ()

January 21, 2011 - Posted by | Programming

No comments yet.

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s

%d bloggers like this: