My Weblog

Blog about programming and math

Improved Elliptic Curve prime factorisation

Currently 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 . This code factorise the product of prime numbers in range [1000..2000] withing less than second on my slow computer.Algorithm explained is previous post. Modified Haskell code.

--y^2= x^3+ax+b mod n
import Random
import Control.Monad
import Data.List
import Data.Bits

data Elliptic = Conelliptic Integer Integer Integer deriving ( Eq , Show )
data Point = Conpoint Integer Integer | Identity deriving ( Eq , Show )


powM :: Integer -> Integer -> Integer -> IO Integer
powM a d n = return $ 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 -> IO ( Integer , Integer ) 
calSd n = return ( s , d ) where 
		s = until ( \x -> testBit ( n - 1) ( fromIntegral x ) )  ( +1 ) 0
                d = div ( n - 1 ) (  shiftL 1 ( fromIntegral s )  )


isProbable::Integer->IO Bool
isProbable n 
   | n <= 1 = return False
   | n == 2 = return True
   | even n = return False
   | otherwise	= calSd n >>= (\( s , d ) -> rabinMiller 0 n s d )


rabinMiller::Integer->Integer->Integer->Integer->IO Bool
rabinMiller cnt n s d 
   | cnt>=5 = return True
   | otherwise = 
       randomRIO ( 2 , n - 2 ) >>= 
         (\a ->  powM a d n >>=
	    (\x -> case x == 1 of  
		   True  -> rabinMiller ( cnt + 1 ) n s d
		   _   -> if any ( == pred n ) . take ( fromIntegral s ) . iterate (\e -> mod ( e^2 ) n ) $ x  
		           then rabinMiller ( cnt + 1) n s d 
			    else return False  ) )


--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 u v= helpGcd u v 0 [ [ 1 , 0 ] , [ 0 , 1 ] ] where
	helpGcd u v n m @ ( [ [ a , b ] , [ c , d ] ] ) 
	  | v == 0 = if u<0 then [ - ( ( -1 ) ^ n ) * ( m !! 1 !! 1 ) , - ( ( -1 ) ^ ( n + 1 ) ) * ( m !! 0 !! 1 ) , -u ]
           	      else [ ( ( -1 ) ^ n ) * ( m !! 1 !! 1 ) , ( ( -1 ) ^ ( n + 1 ) ) * ( m !! 0 !! 1 ) , u ]
	  | otherwise = helpGcd u' v' ( n + 1 ) m' where 
		( q , v' ) = quotRem u v 
		t = [ [q , 1 ] , [ 1 , 0 ] ]
		m' = [ [ q * a + b , a ] , [ q * c + d , c ] ] --mult m t 
		u' = v 
  



multiEC :: Elliptic -> Point -> Integer -> IO ( Either Point Integer )
multiEC _ _ 0 = return $ Left Identity
multiEC ecurve point k  = return $ helpEC Identity point k
	where
	  helpEC p _ 0 = Left p
	  helpEC p q n = 
             case (p',q') of
		  ( Left p'' , Left q'' ) -> helpEC p'' q'' (div n 2)
		  ( Right p'' , _ ) -> Right p''
		  ( _ , Right q'' ) -> Right q''
	          where
		   p' = if odd n then addPoints ecurve p q else Left p
	           q' =  addPoints ecurve q q
			    

dscrmntEC a b = 4 * a * a * a + 27 * b * b


randomCurve :: Integer -> IO [Integer]
randomCurve n = 
   randomRIO ( 1 , n ) >>= 
      ( \a -> randomRIO ( 1 , n ) >>= 
         ( \u -> randomRIO (1 , n) >>= 
             (\v -> return [a , mod ( v*v - u*u*u - a*u ) n , n , u , v ] ) ) ) 


factor :: Integer -> Integer -> IO [Integer]
factor 1 _   = return []
factor n m   = 
    isProbable n >>= (\x -> if x then  return [n] 
			     else 
			       randomCurve n >>= (\[ a , b , n , u , v ] -> 
                             	 multiEC ( Conelliptic a b n ) ( Conpoint u v ) m >>= 
                                   ( \p -> case p of 
                                          Left p' -> factor n m 
                                          Right p'->  factor ( div n p' ) m  >>= ( \x -> factor p' m  >>= (\y  -> return $ x ++ y   ) ) ) ) )




solve :: Integer -> IO [ Integer ]
solve n = factor n ( foldl lcm 1 [1..10000] )   

main = liftM read getLine >>=( \n -> solve n  ) >>= print
[user@localhost Haskell]$ time ./Elliptic_Len 
14771291995722167329218042370113255785101138530361517229308539087571387624569184860194674818433642463913959873342852255543111135151457209398370030102228858128758106619882182645353368949658870563576166323307466834284670054906895828663228915811739807554988946409285249055824035034910383070200873326731464685909980630202533115332294811841112677175465781631555282945350555352712075708686158031835188673039526232167090966030801645637
[1759,1997,1291,1637,1669,1889,1907,1823,1231,1579,1597,1327,1427,1499,1831,1951,1213,1753,1993,1487,1429,1601,1319,1447,1721,1847,1697,1873,1699,1663,1061,1879,1087,1093,1117,1423,1933,1493,1303,1021,1109,1973,1361,1811,1123,1543,1249,1783,1867,1301,1789,1901,1657,1693,1801,1091,1609,1913,1277,1787,1483,1283,1511,1381,1877,1297,1181,1571,1103,1229,1039,1471,1289,1373,1861,1069,1987,1481,1279,1559,1193,1129,1747,1949,1733,1489,1063,1201,1307,1187,1871,1031,1741,1433,1321,1013,1567,1627,1549,1399,1553,1523,1051,1439,1723,1259,1667,1621,1777,1151,1367,1019,1097,1451,1163,1223,1999,1979,1931,1583,1619,1217,1153,1009,1409,1613,1459,1171,1531,1453,1237,1607,1049,1709,1033]

real	0m16.424s
user	0m1.905s
sys	0m0.009s
Advertisements

July 18, 2011 - Posted by | Programming | , , ,

No comments yet.

Leave a Reply

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

WordPress.com Logo

You are commenting using your WordPress.com 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 )

Google+ photo

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

Connecting to %s

%d bloggers like this: