3
\$\begingroup\$

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 ]
     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 paiB₁;
       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())(Z(R)+Z())−α+βδ)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.

\$\endgroup\$
3
  • 1
    \$\begingroup\$ Please add some example code to your question that calls the interesting code. You should also describe "it crashes" more specifically, otherwise you can only expect "then your code must be wrong" as a comment, and that would make your question off-topic on this site. \$\endgroup\$ Jan 25, 2020 at 5:25
  • 1
    \$\begingroup\$ Basically you have to convince us that "the code works as intended according to the best of the author's knowledge". \$\endgroup\$ Jan 25, 2020 at 5:26
  • \$\begingroup\$ I'm pretty sure this code has a bug. Division in python produces floating point numbers, which I'm 95% sure isn't what you intended. \$\endgroup\$ Dec 20, 2021 at 6:57

1 Answer 1

1
\$\begingroup\$

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.