# define a polynomial ring for the variables
QQx.<a,b,l,x_1,y_1,x_2,y_2> = QQ[]
QQE.<X,Y> = QQx[]
L = Y - (l*(X-x_1) + y_1)
# the line equation above is L(X,Y) = 0, but I need to isolate on the left-hand side Y = something,
# the following does that in a blind way: divide the terms that are not multiple of Y by the coefficient of Y
Ly = -L.coefficient({Y:0})/L.coefficient({Y:1})
E = Y^2 - X^3 - a*X - b # short Weierstrass form

# evaluate E at Ly = the value of Y in the line
Eq = E([X, Ly]); print("E(X, value of Y in the line L) = {}".format(Eq))
#-X^3 + l^2*X^2 + (-2*l^2*x_1 + 2*l*y_1 - a)*X + l^2*x_1^2 - 2*l*y_1*x_1 + y_1^2 - b
print("L intersects E has solutions Eq, and\nEq mod (X-x_1) = {} = {} mod E([x_1,y_1])".format(Eq % (X-x_1), (Eq % (X-x_1)) % E([x_1,y_1])))
assert (Eq % (X-x_1)) % E([x_1,y_1]) == 0
print("simplify Eq by X-x1")
#-x_1^3 + y_1^2 - a*x_1 - b
Eq2 = Eq // (X-x_1) ; print("Eq2 = Eq // (X-x_1) = {}".format(Eq2))
#-X^2 + (l^2 - x_1)*X - l^2*x_1 + 2*l*y_1 - x_1^2 - a
Eq3 = Eq2 % (X-x_2); print("Eq3 = Eq2 % (X-x_2) = {}".format(Eq3))
#-l^2*x_1 + l^2*x_2 + 2*l*y_1 - x_1^2 - x_1*x_2 - x_2^2 - a
# replace l by its value
Eq3 = QQx(Eq3) # coerce in the coefficient ring
lambda_add = (y_2-y_1)/(x_2-x_1)
Eq4 = Eq3(l=lambda_add); print("Eq4 = Eq3(l = {})".format(lambda_add))

N = Eq4.numerator(); print("N = Eq4.numerator() = {} = {} mod (E(x_1,y_1), E(x_2,y_2))".format(N, N % (E(x_1,y_1)) % E(x_2,y_2)))
# x_1^3 - x_2^3 + a*x_1 - y_1^2 - a*x_2 + y_2^2
D = Eq4.denominator();  print("N = Eq4.denominator() = {}".format(D))
# -x_1 + x_2
assert (N % (E(x_1,y_1))) % E(x_2,y_2) == 0
print("Simplify Eq by (X-x2) because Eq % (X-x2) = 0 mod the curve equation, when replacing l by its value")
Eq5 = Eq2 // (X-x_2); print("Eq5 = Eq2 // (X-x_2) = {}".format(Eq5))
print("The equation is now linear in X, that is a solution x_3 for X can be deduced easily")
x_3 = -Eq5.coefficient({X:0})/Eq5.coefficient({X:1}); print("x_3 = {}".format(x_3))
y0_3 = QQx(Ly([x_3,Y])); print("knowing x_3, gets y0_3 with L(x_3,y0_3) = 0\ny0_3 = {}".format(y0_3))
# in all generallity, y_3 is not -y0_3
y_3 = -y0_3; print("finally, y_3 = -y0_3 = {}".format(y_3))


# now the doubling
print("\ndoubling formula")
L = E.derivative(X)(x_1,y_1)*(X-x_1) + E.derivative(Y)(x_1,y_1)*(Y-y_1)
print("partial f / partial x (x_1,y_1) = {}".format(E.derivative(X)(x_1,y_1)))
print("partial f / partial y (x_1,y_1) = {}".format(E.derivative(Y)(x_1,y_1)))
lambda_dbl = -L.coefficient({X:1})/L.coefficient({Y:1}); print("lambda_dbl = {}".format(lambda_dbl))

Lyl = l*(X-x_1) + y_1
Ly = lambda_dbl*(X-x_1) + y_1
Eq = E([X, Lyl]); print("E(X, value of Y in the langent L) = {}".format(Eq))
# -X^3 + l^2*X^2 + (-2*l^2*x_1 + 2*l*y_1 - a)*X + l^2*x_1^2 - 2*l*x_1*y_1 + y_1^2 - b
assert QQx(Eq([x_1,y_1])) % QQx(E([x_1,y_1])) == 0
assert QQx(Eq % (X - x_1)) % E([x_1,y_1]) == 0
print("L intersects E has solutions Eq, and\nEq mod (X-x_1) = {} = {} mod E([x_1,y_1])".format(Eq % (X-x_1), (Eq % (X-x_1)) % E([x_1,y_1])))
assert (Eq % (X-x_1)) % E([x_1,y_1]) == 0
print("simplify Eq by X-x1")
Eq2 = Eq // (X-x_1) ; print("Eq2 = Eq // (X-x_1) = {}".format(Eq2))
Eq3 = Eq2 % (X-x_1); print("Eq3 = Eq2 mod (X-x_1) = {}".format(Eq3))
if Eq3 != 0:
    Eq4 = QQx(Eq3)(l=lambda_dbl)
    print("Eq3(l=lambda_dbl) = {}".format(Eq4))
    if Eq4 != 0:
        print("Eq4(l=lambda_dbl) % E([x_1,y_1]) = {}".format(Eq4 % QQx(E([x_1,y_1]))))
# -3*x_1^2 + 2*l*y_1 - a
assert QQx(Eq3)(l=lambda_dbl) == 0
Eq5 = Eq2 // (X-x_1) ; print("Eq5 = Eq2 // (X-x_1) = {}".format(Eq5))
# -X + l^2 - 2*x_1
x_4 = -Eq5.coefficient({X:0})/Eq5.coefficient({X:1}); print("x_4 = {}".format(x_4))
# l^2 - 2*x_1
y_4 = QQx(Lyl([x_4,Y])); print("knowing x_4, gets y0_4 with L(x_4,y0_4) = 0\ny0_4 = {}".format(y_4))
y_4 = -y_4; print("finally, y_4 = -y0_4 = {}".format(y_4))