Modern Elliptic Curve Factorization, Part 2
April 27, 2010
The elliptic curve factorization algorithm is remarkably similar to Pollard’s p-1 algorithm:
define p-1(n, B1, B2)
|
1 |
define ecf(n, B1, B2, C)
|
| q = 2 | 2 | Qx,z ∈ C // random point |
| for p ∈ π(2 … B1) | 3 | for p ∈ π(2 … B1) |
a = ⌊logp B1⌋
|
4 |
a = ⌊logp B1⌋
|
q = qpa (mod n)
|
5 | Q = pa ⊗C Q |
g = gcd(q, n)
|
6 |
g = gcd(Qz, n)
|
| if (1 < g < n) return g | 7 | if (1 < g < n) return g |
| for p ∈ π(B1 … B2) | 8 | for p ∈ π(B1 … B2) |
q = qp (mod n)
|
9 | Q = p ⊗C Q |
| 10 |
g = g × Qz (mod n)
|
|
g = gcd(q, n)
|
11 |
g = gcd(g, n)
|
| if (1 < g < n) return g | 12 | if (1 < g < n) return g |
| return failure | 13 | return failure |
Lines 3 and 8 loop over the same lists of primes. Line 4 computes the same maximum exponent. Lines 6 and 7, and lines 11 and 12, compute the same greatest common divisor. And line 13 returns the same failure. The structures of the two algorithms are identical. Also, note that both algorithms use variables that retain their meaning from the first stage (lines 2 to 7) to the second stage (lines 8 to 12).
Lines 2, 5, 9 and 10 differ because p-1 is based on modular arithmetic and ecf is based on elliptic arithmetic (note that the ⊗C operator, otimes-sub-C, indicates elliptic multiplication on curve C). The most important difference is on the first line, where ecf has an extra argument C indicating the choice of elliptic curve. If the p-1 method fails, there is nothing to do, but if the ecf method fails, you can try again with a different curve.
Use of Montgomery’s elliptic curve parameterization from the previous exercise is a huge benefit. Lenstra’s original algorithm, which is the first stage of the modern algorithm, required that we check the greatest common divisor at each step, looking for a failure of the elliptic arithmetic. With the modern algorithm, we compute the greatest common divisor only once, at the end, because Montgomery’s parameterization has the property that a failure of the elliptic arithmetic propagates through the calculation. And without inversions, Montgomery’s parameterization is also much faster to calculate. Brent’s parameterization of the randomly-selected curve is also important, as it guarantees that the order of the curve (the number of points) will be a multiple of 12, which has various good mathematical properties that we won’t enumerate here.
There is a simple coding trick that can speed up the second stage. Instead of multiplying each prime times Q, we iterate over i from B1 + 1 to B2, adding 2Q at each step; when i is prime, the current Q can be accumulated into the running solution. Again, we defer the calculation of the greatest common divisor until the end of the iteration.
Your task is to write a function that performs 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 
Two questions:
in line 6 of ecf, what is ‘q’
in line 11 of ecf, is that supposed to be g = gcd(*q*,n) like in the p-1 function?
In line 6, q should read Q-sub-z (it’s the z parameter of the point Q). I’ll fix that.
In line 11, g is the g that is initialized on line 6 and accumulated on line 11.
Since no one posted the source code for this so here is mine.How ever running this algorithm and Lenstra implementation, Lenstra seems to be good.
--y^2 = x^3+ax^2+x mod n import Random import System.Random import Data.Bits data Elliptic = Cone Integer Integer Integer deriving (Eq,Show) data Point = Conp Integer Integer deriving (Eq,Show) addPoints::Point->Point->Point->Integer->Point addPoints (Conp 0 0) p2 _ _ = p2 addPoints p1 (Conp 0 0) _ _ = p1 addPoints (Conp x_1 z_1) (Conp x_2 z_2) (Conp 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 (Conp x_3 z_3) doublePoint::Point->Elliptic->Point doublePoint (Conp x z) (Cone an ad n) = Conp (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 --logupperbound inspired by programming praxis --https://programmingpraxis.com/2010/05/07/integer-logarithms/ ilognew::Integer->Integer->Int ilognew b n = f (div ubound 2) ubound where ubound = until (\e -> b ^ e > n) (* 2) 1 f lo hi | lo>=hi = div (lo+hi) 2 | 2 ^ mid <= n = f (mid+1) hi | otherwise = f lo mid where mid = div (lo + hi) 2 --http://en.wikipedia.org/wiki/Elliptic_curve_point_multiplication --montogomery ladder addition but do a binary search to find the 2^p>k so all bits from p-1 will be set montMul::Integer->Elliptic->Point->Integer->Point montMul _ _ _ 0= (Conp 0 0) montMul n eCurve pOint k = tmpFun (Conp 0 0) pOint (ilognew 2 k -1) where tmpFun r0 r1 0 = if (.&.) k ( shiftL (1::Integer) 0 ) /=0 then addPoints r0 r1 pOint n else doublePoint r0 eCurve tmpFun r0 r1 cnt |(.&.) k ( shiftL (1::Integer) cnt ) /=0 = tmpFun (addPoints r0 r1 pOint n) (doublePoint r1 eCurve) (cnt-1) |otherwise = tmpFun (doublePoint r0 eCurve) (addPoints r0 r1 pOint n) (cnt-1) montgomeryCurve::Integer-> IO [Integer] montgomeryCurve n = do s<-getStdRandom (randomR (6,n-1)) let u=mod (s*s-5) n let v=mod (4*s) n let an=mod (((v-u)^3)*(3*u+v)) n let ad=mod (4*u^3*v) n return [an,ad,mod (u*u*u) n,mod (v*v*v) n] prime = 2:filter rabinMiller [3,5..] fun2mk :: Integer->(Integer,Integer) fun2mk n = let (s,d)= f 0 (n-1) where f s d | mod d 2 ==0 = f (s+1) (div d 2) | otherwise = (s,d) in (s,d) rabinMiller :: Integer->Bool rabinMiller n | n<2 = False | n==2 = True | even n = False | b0==1 || b0== n' = True | otherwise = iter (tail b) where n'= n-1 (k,m)=fun2mk n a= head ( take 1 $ randomRs (2,n-2) (mkStdGen 3) :: [Integer]) b0 = mod (a^m) n b= take (fromIntegral k) $ iterate (\x-> mod (x*x) n) b0 iter []= False iter (x:xs) | x==1 = False | x==n'= True | otherwise = iter xs allp::Integer->Integer->[Integer] allp p_1 p_2 = takeWhile (<=p_2) $ dropWhile (<p_1) prime --https://programmingpraxis.com/2010/04/27/modern-elliptic-curve-factorization-part-2/ ellipticFac::Integer->Integer->Elliptic->Point->Integer ellipticFac b_1 b_2 (Cone an ad n) (Conp u v) = let c= foldl lcm 1 [1..b_1] (Conp x z)= montMul n (Cone an ad n) (Conp u v) c g=gcd n z in case (and [1<g, g<n]) of True -> g _ -> let (Conp x' z')= foldl ( montMul n (Cone an ad n) ) (Conp x z) (allp b_1 b_2) g'=gcd (mod (g*z') n) n in case (and [1<g',g'<n]) of True -> g' _ -> n factor::Integer->Integer->IO [Integer] factor 1 _ = return [] factor n cnt | cnt>=5 = return [n] | otherwise = do [an,ad,u,v] <-montgomeryCurve n let p=ellipticFac 10000 20000 (Cone an ad n) (Conp u v) case (and [1<p,p<n]) of False -> factor n (cnt+1) _ -> do ps<-factor (div n p) (cnt+1) ls<-factor p (cnt+1) return (ls++ps) main = do line<-getLine let n=read line ::Integer lst<-factor n 0 print lst return ()Lenstra Implementation
[…] is almost similar to 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 […]
A solution in Python, borrowing the functions from https://programmingpraxis.com/2010/04/23/modern-elliptic-curve-factorization-part-1/#comment-1192:
from random import randrange # Recursive sieve of Eratosthenes def primegen(): yield 2; yield 3; yield 5; yield 7; yield 11; yield 13 ps = primegen() # yay recursion p = ps.next() and ps.next() q, sieve, n = p**2, {}, 13 while True: if n not in sieve: if n < q: yield n else: next, step = q + 2*p, 2*p while next in sieve: next += step sieve[next] = step p = ps.next() q = p**2 else: step = sieve.pop(n) next = n + step while next in sieve: next += step sieve[next] = step n += 2 def add(p1, p2, p0, n): # Add two points p1 and p2 given point P0 = P1-P2 modulo n x1,z1 = p1; x2,z2 = p2; x0,z0 = p0 t1, t2 = (x1-z1)*(x2+z2), (x1+z1)*(x2-z2) return (z0*pow(t1+t2,2,n) % n, x0*pow(t1-t2,2,n) % n) def double(p, A, n): # double point p on A modulo n x, z = p; An, Ad = A t1, t2 = pow(x+z,2,n), pow(x-z,2,n) t = t1 - t2 return (t1*t2*4*Ad % n, (4*Ad*t2 + t*An)*t % n) def multiply(m, p, A, n): # multiply point p by m on curve A modulo 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 *= 2 b /= 4 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 /= 2 return r def ilog(x, b): # greatest integer l such that b**l <= x. # TODO: is there a faster way to compute this? Binsearch, perhaps? l = 0 while x >= b: x /= b l += 1 return l def ecm(n, B1=10, B2=20): # TODO: determine the best defaults for B1 and B2 and the best way to increment them and iters iters = 1 while True: iters *= 2 for _ in xrange(iters): seed = randrange(6, n) u, v = (seed**2 - 5) % n, 4*seed % n p = pow(u, 3, n) Q, C = (pow(v-u,3,n)*(3*u+v) % n, 4*p*v % n), (p, pow(v,3,n)) pg = primegen() p = pg.next() while p <= B1: Q, p = multiply(p**ilog(B1, p), Q, C, n), pg.next() g = gcd(Q[1], n) if 1 < g < n: return g while p <= B2: # "There is a simple coding trick that can speed up the second stage. Instead of multiplying each prime times Q, # we iterate over i from B1 + 1 to B2, adding 2Q at each step; when i is prime, the current Q can be accumulated # into the running solution. Again, we defer the calculation of the greatest common divisor until the end of the # iteration. TODO: implement this trick and compare performance." Q = multiply(p, Q, C, n) g *= Q[1] g %= n p = pg.next() g = gcd(g, n) if 1 < g < n: return g # This seed failed. Try again with a new one. # These bounds failed. Increase and try again. B1 *= 2 B2 *= 2 print ecm(20639383 * 117876683047)I cannot wrap my head around the concept of Elliptic Curves let alone the ECM but I think I found a P195 factor. Is that a record as your article implies the current record stands at P73?
Here is the C227 composite and here are its factors using Yet Another Factorization Utility 1.34 from http://yafu.sourceforge.net
72862001761201652067512910912311143529912811111098135112781186477227938869603430735445831828875855453895937353829184560496255789629222413141111181212305252442828205640315040464229193625221719263020152111881958811118395473635456
=
2 x
2 x
2 x
2 x
2 x
2 x
2 x
7 x
5527 x
137885520290732817351584663 x
106705047172471266468660928563317852379901416010404723328566172957737395755840658022032569277618563231217065120661299797289022426791599936582159383554005399760617203071477939525944722428007054711
What am I missing ???
Also I currently I am factoring this number:
65453637459381111885198811211520302619172225361929424640503140562028284452523012121811111413242229967855624960451829383537598953548575881828345547330346069889322777641187811213598110111128995243111123109129752061651201762002867
Again a C227
Thank you for any feedback,
Email: pureprimes@yahoo.com
Web: http://www.heliwave.com
It seems the numbers have been truncated so I will break them into multiple 50-digit lines, sorry for the double post.
72862001761201652067512910912311143529912811111098
13511278118647722793886960343073544583182887585545
38959373538291845604962557896292224131411111812123
05252442828205640315040464229193625221719263020152
111881958811118395473635456
=
2 x
2 x
2 x
2 x
2 x
2 x
2 x
7 x
5527 x
137885520290732817351584663 x
10670504717247126646866092856331785237990141601040
47233285661729577373957558406580220325692776185632
31217065120661299797289022426791599936582159383554
005399760617203071477939525944722428007054711
Is the above P195 a record for ECMs?
And here is the number being currently factored:
65453637459381111885198811211520302619172225361929
42464050314056202828445252301212181111141324222996
78556249604518293835375989535485758818283455473303
46069889322777641187811213598110111128995243111123
109129752061651201762002867
Thank you for any feedback
You found seven factors of 2, one factor of 7, one factor of 5527, and one 27-digit factor of 137885520290732817351584663. You did not find the remaining P195 co-factor, it is just “left over” after all the factors that you did find. Sorry, not a record.