How to implement birthday paradox continuation of elliptic curve factorization algorithm

Question

I have already implemented Lenstra's algorithm for factoring integers using elliptic curves; it is shown below, or you can run it at http://ideone.com/QEDmMY. Beware that my code is optimized for simplicity and clarity, not for speed.

I would like to implement the second stage of elliptic curve factorization using the birthday paradox continuation, which I read about in Richard Brent's paper "Some Integer Factorization Algorithms Using Elliptic Curves."

However, I don't understand Brent's algorithm, or perhaps just his notation. In Section 6, what is $Q_j^2$ in the definition of $Q_{j+1}$? Likewise, what is $Q_j^2 * Q$? Since $Q$ is just a point on the elliptic curve, how can it be squared or multiplied by another point? And how is $r$ computed? Brent discusses $r$ but never defines it.

Many thanks,

Phil

# lenstra's algorithm

from random import randint
from fractions import gcd

def primes(n):
    b, p, ps = [True] * (n+1), 2, []
    for p in xrange(2, n+1):
        if b[p]:
            ps.append(p)
            for i in xrange(p, n+1, p):
                b[i] = False
    return ps

def bezout(a, b):
    if b == 0: return 1, 0, a
    q, r = divmod(a, b)
    x, y, g = bezout(b, r)
    return y, x-q*y, g

def add(p, q, a, b, m):
    if p[2] == 0: return q
    if q[2] == 0: return p
    if p[0] == q[0]:
        if (p[1] + q[1]) % m == 0:
            return 0, 1, 0 # infinity
        n = (3 * p[0] * p[0] + a) % m
        d = (2 * p[1]) % m
    else:
        n = (q[1] - p[1]) % m
        d = (q[0] - p[0]) % m
    x, y, g = bezout(d, m)
    if g > 1: return 0, 0, d # failure
    z = (n*x*n*x - p[0] - q[0]) % m
    return z, (n * x * (p[0] - z) - p[1]) % m, 1

def mul(k, p, a, b, m):
    r = (0,1,0)
    while k > 0:
        if p[2] > 1:
            return p
        if k % 2 == 1:
            r = add(p, r, a, b, m)
        k = k // 2
        p = add(p, p, a, b, m)
    return r

def lenstra(n, limit):
    g = n
    while g == n:
        q = randint(0, n-1), randint(0, n-1), 1
        a = randint(0, n-1)
        b = (q[1]*q[1] - q[0]*q[0]*q[0] - a*q[0]) % n
        g = gcd(4*a*a*a + 27*b*b, n)
    if g > 1: return g # lucky factor
    for p in primes(limit):
        pp = p
        while pp < limit:
            q = mul(p, q, a, b, n)
            if q[2] > 1:
                return gcd(q[2], n)
            pp = p * pp
    return False

Answer 1

On the first line of page 4 Brent introduces a multiplicative notation for the group operation on the elliptic curve. This is a bit unusual in the context of elliptic curves, but it is difficult to think what else it might mean. Therefore his

1. $Q_j^2$ is the double of the point $Q_j$, many authors use the notation $2Q_j$ for this. IIRC $[2]Q_j$ is also in use, but I may be wrong.
2. Similarly his $Q_j^2*Q$ is more commonly denoted by $2Q_j+Q$.

In other words it sounds something like a step in a double-and-add -algorithm. I didn't read all the paper, so don't know for sure.

So as long as you have implemented the usual point doubling and point addition routines for the elliptic curve, you should be fine.