Montgomery curve group operations

def madd(P,Q,R):
    """Montgomery addition: returns (x_{m+n},z_{m+n}) given P=(x_m,z_m), Q=(x_n,z_n), and R=(x_{m-n},z_{m-n})"""
    a = P[0]-P[1]; b=P[0]+P[1]; c=Q[0]-Q[1]; d=Q[0]+Q[1]
    e=a*d; f=b*c;
    return (R[1]*(e+f)^2,R[0]*(e-f)^2)

def mdbl(P,C):
    """Montgomery doubling: returns (x_{2n},z_{2n}) given P=(x_n,z_n) and C=(A+2)/4 where By^2=x^3+Ax^2+x is the Montgomery curve equation"""
    a=P[0]+P[1]; b=P[0]-P[1]
    c=a^2; d=b^2; e=c-d
    return (c*d,e*(d+C*e))
    
def mmul(P,n,C):
    """Montgomery multiplication: returns (x_n,z_n) given P=(x_1,z_1) and C=(A+2)/4 where By^2=x^3+Ax^2+x is the Montgomery curve equation"""
    Q = [P,mdbl(P,C)]
    b=n.digits(2)
    for i in range(len(b)-2,-1,-1):
        Q[1-b[i]] = madd(Q[1],Q[0],P)
        Q[b[i]] = mdbl(Q[b[i]],C)
    return Q[0]
F=GF(random_prime(2^256,2^255))
A=F.random_element(); B=F.random_element()
C=(A+2)/4
while true:
    x = F.random_element()
    if ((x^3+A*x^2+x)/B).is_square(): break;
P=(x,1);
t=cputime()
for p in primes(0,10000): P=mmul(P,p,C)
print "Montgomery time:", cputime()-t
a4=1/B^2-A^2/(3*B^2); a6=-A^3/(27*B^3)-a4*A/(3*B)
E=EllipticCurve([a4,a6])
P=E.random_element()
t=cputime()
for p in primes(0,10000): P=p*P
print "Weierstrass time:", cputime()-t
Montgomery time: 0.165339 Weierstrass time: 0.688554

Single stage ECM using Montgomery curves

def L(a,c,N):
    z = ceil(exp(c*log(N)^a*log(log(N))^(1-a)))
    return ceil(z)

def ECM_1curve(N,p,M,B1):
    R=Integers(N)
    # generate Montgomery curve using parameterization that guarantees a point of order 3
    c=R(randint(1,10^10))
    a=6*c/(c^2+6)
    b=R(randint(1,10^10))
    A=(-3*a^4-6*a^2+1)/(4*a^3)
    B=(a^2-1)^2/(4*a*b^2)
    C=(A+2)/4
    P=(3*a/4,1)
    a4=1/B^2-A^2/(3*B^2); a6=-A^3/(27*B^3)-a4*A/(3*B)
    F=GF(p)
    E=EllipticCurve([F(a4.lift()),F(a6.lift())])
    print factor(E.cardinality())
    m=ceil(M+2*sqrt(M)+1)
    for p in primes(B1):
        q = p; r=p*q
        while r <= m: q=r; r=p*q
        P=mmul(P,q,C)
        d=gcd(N,P[1].lift())
        if d == N: print "N fail"; return 0
        if d > 1: return d
    return 0
k=50; n=200
q=random_prime(2^(n-k),2^(n-k-1))
p=random_prime(2^k,2^(k-1))
B=2.5*L(1/2,1/sqrt(2),2^k) # includes an empirical fudge factor
print "B =",B
print "Expected number of iterations is about", L(1/2,1/(2*sqrt(2)),2^k)
t=cputime()
i=1
while true:
    print i,
    d = ECM_1curve(p*q,p,2^50,B)
    if d: break
    i += 1
print d
print "ECM time:", cputime()-t
print p

B = 6340.00000000000 Expected number of iterations is about 51 1 2^3 * 3 * 5 * 7 * 1567 * 2531 * 21557 2 2^3 * 3 * 47^2 * 5479 * 247241 3 2^2 * 3^3 * 53 * 12546688591 4 2^4 * 3 * 13 * 115091709323 5 2^5 * 3^4 * 11^2 * 19 * 47 * 256423 6 2^2 * 3^3 * 7 * 94996330883 7 2^2 * 3 * 5984769840877 8 2^3 * 3^5 * 36943028737 9 2^3 * 3^2 * 59 * 16906132519 10 2^4 * 3 * 5 * 17 * 1693 * 10397087 11 2^4 * 3^2 * 13 * 19 * 2019153197 12 2^5 * 3 * 37 * 263 * 751 * 102367 13 2^2 * 3 * 211 * 28363835737 14 2^7 * 3^2 * 131 * 191 * 1259 * 1979 71817240331703 ECM time: 5.743919 71817240331703