Checking associativity for elliptic curves

grouplaw_associativity.sage

S.<a1, a2, a3, a4,a6, x1, x2, x3,y1,y2,y3> = PolynomialRing(QQ,11)
I = S.ideal([y1^2 + a1*x1*y1 + a3*y1 -(x1^3 + a2*x1^2 + a4*x1 + a6), \
             y2^2 + a1*x2*y2 + a3*y2 - (x2^3 + a2*x2^2 + a4*x2 + a6), \
             y3^2 + a1*x3*y3 + a3*y3 -(x3^3 + a2*x3^2 + a4*x3 + a6)])
P = (x1,y1)
Q = (x2,y2)
R = (x3,y3)

def add(P,Q):
    s1 = (Q[1]-P[1]) / (Q[0]-P[0])
    xPQ = s1^2 + a1 *s1 - a2 - P[0] - Q[0]
    yPQ = -a1 * xPQ -a3 - P[1] + s1*(P[0] - xPQ)
    return (xPQ, yPQ)

def double(P):
    s1 = (3*P[0]**2 + 2*a2 *P[0] -a1*P[1]+ a4) / (2 * P[1] + a1*P[0] + a3)
    xPQ = s1^2 + s1*a1 - a2 - 2*P[0]
    yPQ = -a1 * xPQ -a3 - P[1] + s1*(P[0] - xPQ)
    return (xPQ, yPQ)

A1 = add( add(P, Q), R )
B1 = add(P, add(Q, R) )
print( (A1[0] - B1[0]).numerator() in I and (A1[1] - B1[1]).numerator() in I )
# True

A2 = add(double(P),Q)
B2 = add(P,add(P,Q))
print( (A2[0] - B2[0]).numerator() in I and (A2[1] - B2[1]).numerator() in I )
# True

A3 = double(add(P,Q))
B3 = add(P,add(Q,add(P,Q)))
print( (A3[0] - B3[0]).numerator() in I and (A3[1] - B3[1]).numerator() in I )
# True

# Now get the actual polynomials. Here is the worst of them:
f, g, h =(A1[1] - B1[1]).numerator().lift(I)
print(len(f.coefficients()), len(g.coefficients()), len(h.coefficients()))
# Reports (75871, 126838, 132975)