# My Weblog

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

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 d!=1 : return [Ret,d]
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
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