# My Weblog

## 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 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  ) )

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
```