Modern Elliptic Curve Factorization, Part 2

April 27, 2010

Here is our version, which relies on Montgomery’s elliptic arithmetic from the previous exercise:

(define (ec-factor N B1 B2 S)
  (let-values (((An Ad Q) (curve12 N S)))
    (let stage1 ((p 2) (Q Q))
      (if (< p B1)
          (stage1 (next-prime p) (multiply (expt p (ilog p B1)) Q An Ad N))
          (let ((g (gcd (cdr Q) n))) (if (< 1 g n) g
            (let ((QQ (double Q An Ad N))
                  (R (multiply (- B1 1) q An Ad n))
                  (T (multiply (+ B1 1) q An Ad n)))
              (let stage2 ((p (next-prime B1)) (g g) (k (+ B1 1)) (R R) (T T))
                (cond ((< B2 p) (let ((g (gcd g n))) (if (< 1 g n) g #f)))
                      ((= k p) (stage2 (next-prime p) (modulo (* g (cdr T)) N)
                                       (+ k 2) t (add T QQ R N)))
                      (else (stage2 p g (+ k 2) t (add T QQ R N))))))))))))

As an example, consider the factorization of the 99th repunit, the number that consists of 99 consecutive 1s, which is computed as (10^99-1)/9. The factors 3, 3, 37, 67, 199, 397, 21649, 34849, 333667, and 513239 are found by trial division or the rho method, leaving a 70-digit cofactor. The elliptic curve method finds a 19-digit factor 1344628210313298373 using the lucky curve with seed 85398978821858534277806684680913743239224200284672112948006322676318, leaving a remaining 51-digit prime co-factor 362853724342990469324766235474268869786311886053883, thus completing the factorization:

> (ec-factor 250411029487608433927757216420194735307370292135006432870660366092481700801
53000 2120000 85398978821858534277806684680913743239224200284672112948006322676318)

Here’s another example. On September 14, 1998, Conrad Curry found a 53-digit factor of 2677-1 = 1943118631 · 531132717139346021081 · 978146583988637765536217 · P53 · P?? using bounds B1 = 11000000 and B2 = 100B1 with lucky curve S = 8689346476060549; with benefit of hindsight, we are able to reduce those bounds substantially:

> (define (mersenne n) (- (expt 2 n) 1))
> (define x (/ (mersenne 677) 1943118631 531132717139346021081 978146583988637765536217))
> (ec-factor x 9000000 16000000 8689346476060549)

Curry’s accomplishment stood as a record for about fifteen months (that’s a long time in the factoring business). The current record elliptic-curve factor is a 73-digit factor of 21181-1 found by the team of Joppe Bos, Thorsten Kleinjung, Arjen Lenstra, and Peter Montgomery using a program running partly on a cluster of PlayStation 3 computers.

With the addition of the second stage, our elliptic curve factorization program is reasonably powerful; it has been used to find several 25-digit factors, which, notwithstanding the occasional record factorization as mentioned above, is about the limit of elliptic curve factorization. But elliptic curve factorization is generally slow, and our specific implementation is slower still, as we eschewed several algorithmic tricks in the interest of simplicity. There are fast implementations of elliptic curve factorization, but they are quite complicated; for instance, the GNU implementation contains 263 files and weighs in at 4.5MB of code, not counting the required big-integer library.

In addition to the elliptic arithmetic from the previous exercise, ec-factor uses ilog from the Standard Prelude and all of the machinery from the next-prime exercise. You can run the program at

Pages: 1 2

9 Responses to “Modern Elliptic Curve Factorization, Part 2”

  1. Mike said

    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?

  2. programmingpraxis said

    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.

  3. tiwari_mukesh said

    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 (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 = 
    		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 (Conp x z) (Cone an ad n) = Conp (mod x' n) (mod y' n)
    					   x'= mod (4*ad*(x+z)^2*(x-z)^2) n
    					   y'= mod ((4*ad*(x-z)^2+t*an)*t) n 
    --logupperbound inspired by programming praxis
    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
    --montogomery ladder addition but do a binary search to find the 2^p>k so all bits from p-1 will be set 
    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) 
    			  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)
    			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 p_1 p_2 = takeWhile (<=p_2) $ dropWhile (<p_1) prime
    ellipticFac  b_1 b_2 (Cone an ad n)  (Conp u v) = 
    		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
    		_ ->   
    		   (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
    		let n=read line ::Integer
    		lst<-factor n 0
    		print lst
    		return ()
  4. […] 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 […]

  5. Lucas A. Brown said

    A solution in Python, borrowing the functions from

    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 = and
        q, sieve, n = p**2, {}, 13
        while True:
            if n not in sieve:
                if n < q: yield n
                    next, step = q + 2*p, 2*p
                    while next in sieve: next += step
                    sieve[next] = step
                    p =
                    q = p**2
                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
            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 =
                while p <= B1: Q, p = multiply(p**ilog(B1, p), Q, C, n),
                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 =
                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)
  6. Ali Adams said

    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

    2 x
    2 x
    2 x
    2 x
    2 x
    2 x
    2 x
    7 x
    5527 x
    137885520290732817351584663 x

    What am I missing ???

    Also I currently I am factoring this number:

    Again a C227

    Thank you for any feedback,

  7. Ali Adams said

    It seems the numbers have been truncated so I will break them into multiple 50-digit lines, sorry for the double post.

    2 x
    2 x
    2 x
    2 x
    2 x
    2 x
    2 x
    7 x
    5527 x
    137885520290732817351584663 x

    Is the above P195 a record for ECMs?

    And here is the number being currently factored:


    Thank you for any feedback

  8. programmingpraxis said

    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.

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s


Get every new post delivered to your Inbox.

Join 701 other followers

%d bloggers like this: