## Modern Elliptic Curve Factorization, Part 1

### April 23, 2010

We discussed Hendrik Lenstra’s algorithm for factoring integers using elliptic curves in three previous exercises. After Lenstra published his algorithm in 1987, mathematicians studied the algorithm extensively and made many improvements. In this exercise, and the next, we will study a two-stage version of elliptic curve factorization that features improved elliptic arithmetic and is much faster than Lenstra’s original algorithm. Today’s exercise looks at the improved elliptic arithmetic; we will look at the two-stage algorithm in the next exercise.

If you look at Lenstra’s original algorithm, much of the time is spent calculating modular inverses. Peter Montgomery worked out a way to eliminate the modular inverses by changing the elliptic arithmetic; instead of using the curve y2 = x3 + ax + b, Montgomery uses the curve y2 = x3 + ax2 + x. Additionally, he does a co-ordinate transform to eliminate the y variable, so we don’t need to figure out which of the quadratic roots to follow. Finally, Montgomery stores the numerator and denominator of the a parameter separately, eliminating the modular inverse. We won’t give more of the math; for those who are interested, a good reference is Prime Numbers: A Computational Perspective by Richard Crandall and Carl Pomerance.

Here are the elliptic addition, doubling and multiplication rules for Montgomery’s elliptic curves:

A curve is given by the formula y2 = x3 + (An / Ad – 2)x2 + x (mod n).

A point is stored as a two-tuple containing its x and z co-ordinates.

To add two points P1(x1,z1) and P2(x2,z2), given the point P0(x0,z0) = P1P2 (notice that the curve is not used in the calculation):

x = ((x1z1) × (x2+z2) + (x1+z1) × (x2z2))2 × z0 mod n
z = ((x1z1) × (x2+z2) – (x1+z1) × (x2z2))2 × x0 mod n

To double the point P:

x = (x+z)2 × (xz)2 × 4 × Ad mod n
z = (4 × Ad × (xz)2 + t × An) × t mod n where t = (x+z)2 – (xz)2

Multiplication is done by an adding/doubling ladder, similar to the method used in the Lucas chain in the exercise on the Baillie-Wagstaff primality checker; the Lucas chain is needed so that we always have the point P1P2. For instance, the product 13P is done by the following ladder:

2P = double(P)
3P = add(2P, P, P)
4P = double(2P)
6P = double(3P)
7P = add(4P, 3P, P)
13P = add(7P, 6P, P)

Richard Brent gave a method to find a suitable elliptic curve and a point on the curve given an initial seed s on the half-open range [6, n), with all operations done modulo n:

u = s2 – 5
v = 4s

An = (vu)3 (3u+v)

P = (u3, v3)

Your task today is to write the three functions that add, double, and multiply points on an elliptic curve and the function that finds a curve and a point on the curve; the payoff will be in the next exercise that presents elliptic curve factorization. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.

Pages: 1 2

### 5 Responses to “Modern Elliptic Curve Factorization, Part 1”

1. […] today’s Programming Praxis exercise we have to implement an algorithm for curve factorization. The Scheme […]

2. Remco Niemeijer said

My Haskell solution (see http://bonsaicode.wordpress.com/2010/04/23/programming-praxis-modern-elliptic-curve-factorization-part-1/ for a version with comments):

```add :: Integral a => (a, a) -> (a, a) -> (a, a) -> a -> (a, a)
add (x1, z1) (x2, z2) (x0, z0) n = (f (+) z0, f (-) x0) where
f g m = (g ((x1-z1)*(x2+z2)) ((x1+z1)*(x2-z2)))^2 * m `mod` n

double :: Integral a => (a, a) -> a -> a -> a -> (a, a)
double (x,z) an ad n = (mod x' n, mod z' n) where
x' = 4 * ad * (x-z)^2 * (x+z)^2
z' = (4 * ad * (x-z)^2 + t * an) * t
t = (x+z)^2 - (x-z)^2

multiply :: Integral a => Int -> (a, a) -> a -> a -> a -> (a, a)
multiply 0 _ _  _  _ = (0,0)
multiply 1 p _  _  _ = p
multiply k p an ad n = let half = multiply (div k 2) p an ad n in
if odd k then add half (multiply (div k 2 + 1) p an ad n) p n
else double half an ad n

curve12 :: Integral a => a -> a -> (a, a, (a, a))
curve12 n s = ((v-u)^3 * (3 * u + v) `mod` n,
4 * u^3 * v `mod` n,
(u^3 `mod` n, v^3 `mod` n)) where
u = s * s - 5 `mod` n
v = 4 * s `mod` n
```
3. Mike said

Python version.

```def init( n, seed ):
u = ( seed*seed - 5 ) % n
v = 4*seed % n

num = pow(v-u, 3, n) * (3*u + v) % n
den = 4 * pow(u,3,n) * v % n

A = num, den
P = (pow(u, 3, n), pow( v, 3, n) )
return A, P

def add( p1, p2, p0, n ):
x1,z1 = p1
x2,z2 = p2
x0,z0 = p0

t1 = (x1-z1)*(x2+z2)
t2 = (x1+z1)*(x2-z2)

newx = z0 * pow( t1 + t2, 2, n ) % n
newz = x0 * pow( t1 - t2, 2, n ) % n

return newx, newz

def double( p, A, n ):
x,z = p
An, Ad = A

t1 = pow(x + z, 2, n)
t2 = pow(x - z, 2, n)
t = t1 - t2

newx = t1 * t2 * 4 * Ad % n
newz = (4 * Ad * t2 + t * An) * t % n

return newx, newz

def multiply( m, p, A, n ):
if m == 0:
return 0,0
elif m == 1:
return p
else:
q = double( p, A, n )

if m == 2:
return q

b = 1
while b < m:
b <<= 1
b >>= 2

r = p
while b:
if m&b:
q, r = double(q, A, n), add(q, r, p, n)
else:
q, r = add(r, q, p, n), double(r, A, n)
b >>= 1
return r

n,s = 5569379,4921134
A, P1 = init(n, s)
P2 = double(P1,A,n)
P4 = double(P2, A, n)
P6 = double(P3, A, n)
P13a = multiply(13,P1,A,n)

print P13, P13a, P13 == P13a
```
4. Hi Programmingpraxis
Could you tell me if my implementation is correct or not specially for multiPoint function as i was trying to reduce one more call of mulitpOint function.

```--y^2 = x^3+ax^2+b
data Elliptic = Conelliptic Integer Integer Integer deriving (Eq,Show)
data Point = Conpoint Integer Integer deriving (Eq,Show)

addPoints (Conpoint x_1 z_1) (Conpoint x_2 z_2) (Conpoint 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 (Conpoint x_3 z_3)

doublePoint::Point->Elliptic->Point
doublePoint (Conpoint x z) (Conelliptic an ad n) = (Conpoint (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

multiPoint::Integer->Integer->Point->Elliptic->Point
multiPoint _ 0 _ _ = (Conpoint 0 0)
multiPoint _ 1 pOint _ = pOint
multiPoint n k pOint eCurve =let
p=multiPoint n (div k 2) pOint eCurve
in case even k of
True -> doublePoint p eCurve
_ ->  addPoints p (addPoints p pOint pOint n) pOint  n  -- for 7p i am using 6p+pOint

montgomeryCurve::Integer->Integer->[Integer]
montgomeryCurve s n = let
u=mod (s*s-5) n
v=mod (4*s) n
an=mod (((v-u)^3)*(3*u+v)) n