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) = P1 – P2 (notice that the curve is not used in the calculation):
x = ((x1–z1) × (x2+z2) + (x1+z1) × (x2–z2))2 × z0 mod n
z = ((x1–z1) × (x2+z2) – (x1+z1) × (x2–z2))2 × x0 mod nTo double the point P:
x = (x+z)2 × (x–z)2 × 4 × Ad mod n
z = (4 × Ad × (x–z)2 + t × An) × t mod n where t = (x+z)2 – (x–z)2Multiplication 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 P1–P2. 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 = 4sAn = (v–u)3 (3u+v)
Ad = 4u3vP = (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.
Programming Praxis 
[…] today’s Programming Praxis exercise we have to implement an algorithm for curve factorization. The Scheme […]
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` nPython 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) P3 = add(P2,P1,P1,n) P4 = double(P2, A, n) P6 = double(P3, A, n) P7 = add(P4,P3,P1,n) P13 = add(P6,P7,P1,n) P13a = multiply(13,P1,A,n) print P13, P13a, P13 == P13aHi 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::Point->Point->Point->Integer->Point 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 ad=mod (4*u^3*v) n in [an,ad,n,mod (u*u*u) n,mod (v*v*v) n][…] 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 […]