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::Point->Point->Point->Integer->Point 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 = let 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::Point->Elliptic->Point doublePoint (Conp x z) (Cone an ad n) = Conp (mod x' n) (mod y' n) where x'= mod (4*ad*(x+z)^2*(x-z)^2) n y'= mod ((4*ad*(x-z)^2+t*an)*t) n where t=(x+z)^2-(x-z)^2 ilognew::Integer->Integer->Int 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::Integer->Elliptic->Point->Integer->Point 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::Integer->Integer->[Integer] allp p_1 p_2 = takeWhile (<=p_2) $ dropWhile (<p_1) prime --http://programmingpraxis.com/2010/04/27/modern-elliptic-curve-factorization-part-2/ ellipticFac::Integer->Integer->Elliptic->Point->Integer ellipticFac b_1 b_2 (Cone an ad n) (Conp u v) = let 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 _ -> let (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 line<-getLine let n=read line ::Integer lst<-factor n 0 print lst return ()
No comments yet.
Leave a comment