Inversionless Lenstra ECM running way slower than it should

Question

So I was coding Inversionless ECM using Prime Numbers: A computational perspective and when I finished and tried running my program it ran way, way slower than it is supposed to. I have no idea what I am doing wrong. The algorithm I'm coding is here and is algorithm 7.4.4(Ctrl-f it)

1. [Choose criteria]
    B₁ = 10000;  // Stage-one limit (must be even).
    B₂ = 100 B₁; // Stage-two limit (must be even).
    D = 100;    // Total memory is about 3 D size-n integers.
2. [Choose random curve Eσ]
    Choose random σ∈[6, n−1];  // Via Theorem 7.4.3.
    u = (σ²−5) mod n;
    v = 4σ mod n;
    C = ((v−u)³(3 u+v)/(4 u³v)−2) mod n;
      // Note: C determines curve y²=x³+Cx²+x,
      // yet, C can be kept in the form num/den.
    Q = [u³mod n : v³mod n];   // Initial point is represented [X : Z].
3. [Perform stage one]
    for(1≤i≤π(B₁)) { // Loop over primes pi.
      Find largest integer a such that pai≤B₁;
      Q = [pai]Q; // Via Algorithm 7.2.7, and perhaps use FFT enhancements (see text following).
    }
    g = gcd(Z(Q), n); // Point has form Q = [X(Q) : Z(Q)].
    if(1<g<n) return g;  // Return a nontrivial factor of n.
4. [Enter stage two]  // Inversion-free stage two.
    S₁ = doubleh(Q);
    S₂ = doubleh(S₁);
    for(d∈[1,D]) {  // This loop computes Sd=[2 d]Q.
      if(d>2) Sd = addh(Sd−1,S₁,Sd−2);
      βd = X(Sd)Z(Sd)mod n;  // Store the XZ products also.
    }
    g = 1;
    B = B₁−1;  // B is odd.
    T = [B−2 D]Q;  // Via Algorithm 7.2.7.
    R = [B]Q;  // Via Algorithm 7.2.7.
    for(r=B; r<B₂; r=r+2 D) {
      α = X(R)Z(R)mod n;
      for(prime q∈[r+2, r+2 D]) {  //Loop over primes.
        δ = (q−r)/2;  // Distance to next prime.
            // Note the next step admits of transform enhancement.
        g = g((X(R)−X(Sδ))(Z(R)+Z(Sδ))−α+βδ)mod n;
      }
      (R, T) = (addh(R, SD, T), R);
    }
    g = gcd(g, n);
    if(1<g<n) return g;  // Return a nontrivial factor of n.
5. [Failure]
   goto [Choose random curve...]; // Or increase B₁, B₂limits, etc.

Here is my code:

from random import randrange
from primesieve import primes
from math import floor, log

def gcd(a, b):
    if a == b: return a
    while b > 0: 
        a, b = b, a % b
    return a  



def addh(xp, zp, xq, zq, x0, z0, n):
    t = (xp - zp) * (xq + zq)
    v = (xp + zp) * (xq - zq)
    addx, addz = (t + v), (t - v)
    addx, addz = addx * addx, addz * addz
    addx, addz = addx * z0, addz * x0
    if addx >= n:
        addx = addx % n
    if addz >= n:
        addz = addz % n
    return (addx, addz)

def doubleh(xp, zp, a24, n):
    t, v = (xp + zp), (xp - zp)
    t, v = t * t, v * v
    u = t - v
    addx = t * v
    addz = u * (v + a24 * u)
    if addx >= n:
        addx = addx % n
    if addz >= n:
        addz = addz % n
    return (addx, addz)


def montladder(n, X, Z, a24, r):
    if n == 1:
        return (X, Z)
    if n == 2:
        return doubleh(X, Z, a24, r)
    U, V = X, Z
    T, W = doubleh(X, Z, a24, r)
    bk = bin(n)
    for nj in bk[2:]:
        if nj == 1:
            U, V = addh(T, W, U, V, X, Z, r)
            T, W = doubleh(T, W, a24, r)
        else:
            T, W = addh(U, V, T, W, X, Z, r)
            U, V = doubleh(U, V, a24, r)
    if bk[-1] == 1:
        return addh(U, V, T, W, X, Z, r)
    return doubleh(U, V, a24, r)


def fastECM(n, B1 = 10000, B2 = 100, D = 100):
    # Criteria
    B2 = B1 * B2
    S = [0] * (D * 2 + 1)
    be = [0] * (D+1)
    # Choose random Curve Eo
    g = 1
    while g == n or g == 1:
        o = randrange(6, n-1)
        u = (o**2 - 5) % n
        v = 4 * o % n
        t1 = (v - u)**3
        t2 = (3 * u + v)
        t3 = 4 * (u**3) * v
        C = ((t1 * t2 / t3) - 2) % n
        a24 = (C+2) / 4
        Q = (pow(u, 3, n), pow(v, 3, n))

        # Perform Stage one

        for pi in primes(B1):
            a = int(log(B1, pi))
            Q = montladder(pi**a, Q[0], Q[1], a24, n)
        g = gcd(Q[1], n)
        if g > 1 and g < n: return g

        # Perform Stage two
        S[1] = doubleh(Q[0], Q[1], a24, n)
        S[2] = doubleh(S[1][0], S[1][1], a24, n)
        be[1] = S[1][0] * S[1][1] % n
        be[2] = S[2][0] * S[2][1] % n
        for d in range(3, D+1):
            S[d] = addh(S[d-1][0], S[d-1][1], S[1][0],S[1][1], S[d-2][0], S[d-2][1], n)
            be[d] = (S[d][0] * S[d][1]) % n
        g = 1
        B = B1 - 1
        T = montladder(B - 2 * D, Q[0], Q[1], a24, n)
        R = montladder(B, Q[0], Q[1], a24, n)
        r = B
        for r in range(B1, B2, 2 * D):
            alph = R[0] * R[1] % n
            for q in primes(r + 2, r + 2 * D):
                spec = (q - r) // 2
                t1 = (R[0] - S[spec][0])
                t2 = R[1] + S[spec][1]
                g = (g * (t1 * t2 - alph + be[spec])) % n
            R, T = addh(R[0], R[1], S[D][0], S[D][1], R[0], R[1], n), R
        g = gcd(g, n)
    return g

Basically this is taking minutes for 6 digit numbers and anything above my terminal just crashes. I'm not sure if this is the right place to ask this so if it isn't please tell me where I can.

Answer 1

Just some general thoughts since I am not familiar with the ECM algorithm you are using:

The standard implementation of Python is interpreted. The Python executive is switched to for every line of your source. Thus idiomatic Python will use Numpy or other libraries to express the algorithm as calls to matrix routines or as operations on complex values or quaternions. Thus the interpreter is called less often and code spends more of its time running in compiled routines. In contrast your code looks more like Fortran.

Other options are using Cython, which compiles Python to machine code, and writing compiled functions for routines that are called often inside loops.

The standard module profile can tell you how often a line is executed and identify where your code is spending its time.

You don't say which version of Python you are running with, but if you are using a 32 bit version of Python 2 you can get an immediate increase in speed and stability by upgrading to Python 3.