#Let X={[x1,y1],[x2,y2],[x3,y3]} be a generic divisor and
#A={[O1,0],[O2,0],[O3,0]} an element of order 2. We write
#g(x)=(x-O1)(x-O2)(x-O3) and f(x)=g(x)*h(x). This maple session computes
#the matrix which gives addition by the divisor A
h:=x->h4*x^4+h3*x^3+h2*x^2+h1*x+h0:
g:=x->g3*x^3+g2*x^2+g1*x+g0:
f:=g*h:
#As in Stubbs's thesis, we define the coordinates in the Kummer of X.
a0:=1:
a1:=x1+x2+x3:
a2:=x1*x2+x2*x3+x1*x3:
a3:=x1*x2*x3:
#If we define
delta:=(x1-x2)*(x2-x3)*(x3-x1):
H0:=-f7*a1^3+f7*a2*a1-f6*a1^2+3*f7*a3+2*f6*a2:
H1:=-f7*a1^4+3*f7*a2*a1^2-f6*a1^3-f7*a2^2-f7*a3*a1+2*f6*a2*a1-f5*a1^2+
2*f5*a2:
H2:=f7*a2*a1^3-2*f7*a2^2*a1+f6*a2*a1^2+f7*a3*a2-f6*a2^2+f5*a2*a1-3*f5*a3:
H3:=f7*a2^2*a1^2-f7*a1^3*a3+f7*a1*a2*a3-f7*a2^3+f6*a2^2*a1-f6*a3*a1^2+
f5*a2^2-f5*a3*a1:
b0:=(x1*y2-x2*y1-x3*y2+x3*y1-x1*y3+x2*y3)/delta:
b1:=(x3^2*y2-x3^2*y1+x2^2*y1+y3*x1^2-y2*x1^2-y3*x2^2)/delta:
#In our case, we have
f7:=coeff(f(x),x,7):
f6:=coeff(f(x),x,6):
f5:=coeff(f(x),x,5):
#We have
c0:=a0*b0^2+H0:
c1:=a1*b0^2+2*a0*b0*b1+H1:
c2:=a0*b1^2-a2*b0^2+H2:
c3:=a1*b1^2+2*a2*b0*b1+a3*b0^2+H3:
#We first want to find the quartic M and the element gamma such that the
#equation (X-gamma)Y=M(X) passes trought the points [x1,y1], [x2,y2],
#[x3,y3], [O1,0], [O2,0] and [O3,0]. The 3 remaining points in the
#intersection of the curve and (X-gamma)Y=M(X) will give the opposite of
#the sum of X and A (which is the same in the kummer than the sum of X
#and A). For this, we have to solve the system given by
#(xi-gamma)yi=g(xi)(alpha xi+beta) for i=1,2,3.
#In other words, if we define
System1:=[[x1*G(x1),G(x1),y1],[x2*G(x2),G(x2),y2],[x3*G(x3),G(x3),y3]]:
System2:=[x1*y1,x2*y2,x3*y3]:
#the solution (alpha,beta,gamma) is given by System1^(-1)*System2. If d
#is the determinant of System1, this solution can be written
#(alpha=a/d,beta=b/d,gamma=A/d)
with(linalg):
d:=-det(System1):
sol:=simplify(evalm(d*System1^(-1)&*System2)):
a:=sol[1]:
b:=sol[2]:
A:=sol[3]:
#We have now to compute the squares of these quantities in order to
#compute the residual intersection between (X-gamma)Y=M(X) and the curve.
#But we have to eliminate the term in y, as for example the term in y1y2.
#The term in yi^2 can of course be replaced by f(xi). Let us now express
#the terms in yiyj with the ai and the ci using b0^2, b0b1 and b1^2 (where
#y1y2, y1y3 and y2y3 appear). For this, we first compute delta^2*b0^2.
#Note that with our notations, we have
#a0*b0^2=c0-H0, (1)
#a1*b0^2+2*a0*b0*b1=c1-H1, (2)
#a0*b1^2-a2*b0^2=c2-H2. (3)
#So that the following must be zero
simplify(a0*b0^2-(c0-H0));
simplify(a1*b0^2+2*a0*b0*b1-(c1-H1));
simplify(a0*b1^2-a2*b0^2-(c2-H2));
#In the following, we prefer keeping the ci in generic form. In the same
#way, we use the function F,G,H instead of f,g,h and we will replace them
#later.
c0:='c0':c1:='c1':c2:='c2':c3:='c3':
b0b0:=simplify(expand(delta^2*b0^2)):
#We then replace yi^2 by F(xi)=G(xi)*H(xi). This can be done using the
#following procedure
replace:=proc(f)
local f2,f3,f4;
f2:=simplify(f+coeff(f,y1,2)*(G(x1)*H(x1)-y1^2)):
f3:=simplify(f2+coeff(f2,y2,2)*(G(x2)*H(x2)-y2^2)):
f4:=simplify(f3+coeff(f3,y3,2)*(G(x3)*H(x3)-y3^2)):
f4:
end:
b0b0:=replace(b0b0):
#We extract from this the term without any yi
b0b0_y:=subs(y1=0,y2=0,y3=0,b0b0):
#Let us do the same thing for b0b1 and b1^2
b0b1:=simplify(expand(delta^2*b0*b1)):
b0b1:=replace(b0b1):
b0b1_y:=subs(y1=0,y2=0,y3=0,b0b1):
b1b1:=simplify(expand(delta^2*b1^2)):
b1b1:=replace(b1b1):
b1b1_y:=subs(y1=0,y2=0,y3=0,b1b1):
#Let us now consider equations (1), (2) and (3) as a system whose
#variables are 2*(x3-x2)*(x1-x3)*y1*y2, 2*(x3-x2)*(x2-x1)*y1*y3 and
#2*(x1-x3)*(x2-x1)*y2*y3
System3:=simplify([
[coeff(coeff(a0*b0b0,y1,1),y2,1)/(2*(x3-x2)*(x1-x3)),
coeff(coeff(a0*b0b0,y1,1),y3,1)/(2*(x3-x2)*(x2-x1)),
coeff(coeff(a0*b0b0,y2,1),y3,1)/(2*(x1-x3)*(x2-x1))],
[coeff(coeff(a1*b0b0+2*a0*b0b1,y1,1),y2,1)/(2*(x3-x2)*(x1-x3)),
coeff(coeff(a1*b0b0+2*a0*b0b1,y1,1),y3,1)/(2*(x3-x2)*(x2-x1)),
coeff(coeff(a1*b0b0+2*a0*b0b1,y2,1),y3,1)/(2*(x1-x3)*(x2-x1))],
[coeff(coeff(a0*b1b1-a2*b0b0,y1,1),y2,1)/(2*(x3-x2)*(x1-x3)),
coeff(coeff(a0*b1b1-a2*b0b0,y1,1),y3,1)/(2*(x3-x2)*(x2-x1)),
coeff(coeff(a0*b1b1-a2*b0b0,y2,1),y3,1)/(2*(x1-x3)*(x2-x1))]]):
System4:=simplify([
(c0-H0)*delta^2-a0*b0b0_y,
(c1-H1)*delta^2-(a1*b0b0_y+2*a0*b0b1_y),
(c2-H2)*delta^2-(a0*b1b1_y-a2*b0b0_y)]):
#System3 is a Van der Mond matrix with determinant equals to delta, so
#that we can invert it and then obtain the expression of y1y2, y1y3 and
#y2y3 in terms of the ai and the ci. Note that, at first, we prefer write
#System4 in a generic form
sol3:=simplify(evalm(System3^(-1)&*[e0,e1,e2])):
#The intersection of the curve and (dX-A)Y=G(X)(aX+b) is then given by
#(dX-A)^2F(X)=G(X)^2(aX+b)^2 which is of degree 9 in X. G(X) is in the
#two terms (it corresponds to the points [O1,0], [O2,0] and [O3,0]). So
#that we define L which is a polynomial in X of degree 6.
L:=(d*X-A)^2*H(X)-G(X)*(a*X+b)^2:
#We replace yi^2 by F(xi)=G(xi)*H(xi) in this expression
L:=replace(L):
#We now replace yiyj thanks to system 3 and 4
Ly1y2:=simplify(coeff(coeff(L,y1,1),y2,1)/(x1-x3)/(x3-x2)/2):
Ly1y3:=simplify(coeff(coeff(L,y1,1),y3,1)/(x2-x1)/(x3-x2)/2):
Ly2y3:=simplify(coeff(coeff(L,y2,1),y3,1)/(x1-x3)/(x2-x1)/2):
L_y:=subs(y1=0,y2=0,y3=0,L):
#So that the following must be 0
simplify(L-Ly1y2*2*(x3-x2)*(x1-x3)*y1*y2-Ly1y3*2*(x3-x2)*(x2-x1)*y1*y3-
Ly2y3*2*(x1-x3)*(x2-x1)*y2*y3-L_y);
L:=simplify(L_y+Ly1y2*sol3[1]+Ly1y3*sol3[2]+Ly2y3*sol3[3]):
#Note that L is divisible by G(x1)G(x2)G(x3). Since we are looking for
#roots of L, we can divide L by G(x1)G(x2)G(x3).
L:=simplify(L/G(x1)/G(x2)/G(x3)):
#We can then write L=P+P0*e0+P1*e1+P2*e2 with
P:=simplify(subs(e0=0,e1=0,e2=0,L)):
P0:=simplify(coeff(L,e0,1)):
P1:=simplify(coeff(L,e1,1)):
P2:=simplify(coeff(L,e2,1)):
#Our goal is now to express that x1, x2 and x3 are in the intersection.
#In other words, L is divisible by (X-x1)(X-x2)(X-x3). For this we have
#to write polynomials G and H with their coefficients
P:=simplify(subs(G=g,H=h,P)/((X-x1)*(X-x2)*(X-x3))):
P0:=simplify(subs(G=g,H=h,P0)/((X-x1)*(X-x2)*(X-x3))):
P1:=simplify(subs(G=g,H=h,P1)/((X-x1)*(X-x2)*(X-x3))):
P2:=simplify(subs(G=g,H=h,P2)/((X-x1)*(X-x2)*(X-x3))):
#P0, P1 and P2 are very simple. These polynomial are all of degree 3 in
#X, as required. Remember we want to obtain polynomials linear in the ai
#and the ci. As c0 lies only in e0, it is natural to define
p0:=simplify(P0-a3*h4*g(X)):
p1:=simplify(P1-a2*h4*g(X)):
p2:=simplify(P2-a1*h4*g(X)):
#So that p0, p1 and p2 do not depend of the ai. Hence p0 is a good
#candidate to be the coefficient of c0. We have with these notations
N:=subs(G=g,H=h,subs(e0=System4[1],e1=System4[2],e2=System4[3],
P+(p0+a3*h4*g(X))*e0+(p1+a2*h4*g(X))*e1+(p2+a1*h4*g(X))*e2)):
#We can then divide this expression by delta^2 in order to obtain p0 as
#the coefficient of c0.
N:=simplify(N/delta^2):
#But we want this polynomial to be linear in the ai and the ci, and there
#is in this expression a term in a3*c0. To eliminate it, we use the
#relation on the Kummer
#a3*(c0-H0)+a2*(c1-H1)+a1*(c2-H2)=a0*c3*(c3-H3) where
H3:=f7*a2^2*a1^2-f7*a1^3*a3+f7*a1*a2*a3-f7*a2^3+f6*a2^2*a1-f6*a3*a1^2+
f5*a2^2-f5*a3*a1:
coeff_a3c0:=simplify((coeff(N,c0,1)-subs(x1=0,coeff(N,c0,1)))/a3):
#this coefficient is h4*g(X) and is the same for a2*c1 and a1*c2
N:=simplify(N-coeff_a3c0*(a3*(c0-H0)+a2*(c1-H1)+a1*(c2-H2)-a0*(c3-H3))):
#By this way, N contains only linear terms in c0, c1, c2 and c3. The
#remaining terms are linear in the ai
M:=subs(c0=0,c1=0,c2=0,c3=0,N):
#We compute the terms in a0, a1, a2 and a3
ta0:=simplify(subs(x1=0,x2=0,x3=0,M)):
ta1:=simplify(subs(x2=0,x3=0,M-ta0)/x1):
ta2:=simplify(subs(x3=0,M-ta0-ta1*a1)/x1/x2):
ta3:=simplify((M-ta0-ta1*a1-ta2*a2)/a3):
#So that the following must be 0
simplify(N-(ta0*a0+ta1*a1+ta2*a2+ta3*a3+p0*c0+p1*c1+p2*c2+h4*g(X)*c3));
#ta0, ta1, ta2, ta3 are polynomials in X of degree 3 with coefficient in
#Z[g0,..,g3,h0,..,h4]. Finally N is a polynomial in X of degree 3 and
#is a linear combination of the ai and the ci. The roots of this
#polynomial define the sum of the element X and the divisor A of order 2.
#These roots are of course awful to compute but we are only interested in
#the coordinates in the Kummer of this sum. The first four coordinates
#are projectively equal to the coefficients of the polynomial N.
#If xi_i(X) denotes the (i+1)-th coordinate in the Kummer of X, we can
#write xi_i(X+A)=sum_{j=0}^7 wij*xi_j(X) for i=0..3.
#The coefficient wij are given by
w01:=coeff(ta0,X,3): w11:=-coeff(ta0,X,2):
w02:=coeff(ta1,X,3): w12:=-coeff(ta1,X,2):
w03:=coeff(ta2,X,3): w13:=-coeff(ta2,X,2):
w04:=coeff(ta3,X,3): w14:=-coeff(ta3,X,2):
w05:=coeff(p0,X,3): w15:=-coeff(p0,X,2):
w06:=coeff(p1,X,3): w16:=-coeff(p1,X,2):
w07:=coeff(p2,X,3): w17:=-coeff(p2,X,2):
w08:=coeff(h4*g(X),X,3): w18:=-coeff(h4*g(X),X,2):
w21:=coeff(ta0,X,1): w31:=-coeff(ta0,X,0):
w22:=coeff(ta1,X,1): w32:=-coeff(ta1,X,0):
w23:=coeff(ta2,X,1): w33:=-coeff(ta2,X,0):
w24:=coeff(ta3,X,1): w34:=-coeff(ta3,X,0):
w25:=coeff(p0,X,1): w35:=-coeff(p0,X,0):
w26:=coeff(p1,X,1): w36:=-coeff(p1,X,0):
w27:=coeff(p2,X,1): w37:=-coeff(p2,X,0):
w28:=coeff(h4*g(X),X,1): w38:=-coeff(h4*g(X),X,0):
#We then compute the 4 last lines of the matrix W using that the square
#of W is a scalar matrix. We put
A:=[[w01,w02,w03,w04],
[w11,w12,w13,w14],
[w21,w22,w23,w24],
[w31,w32,w33,w34]]:
B:=[[w05,w06,w07,w08],
[w15,w16,w17,w18],
[w25,w26,w27,w28],
[w35,w36,w37,w38]]:
AA:=simplify(evalm(A^2)):
iB:=evalm(B^(-1)):
DD:=simplify(evalm(-iB&*A&*B)):
c:=cte*resultant(g(x),h(x),x):
Id:=[[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]]:
cIAA:=c*Id-AA:
C:=simplify(evalm(iB&*cIAA)):
#Finally, we have the matrix W
W:=matrix(8,8):
for i from 1 to 4 do
for j from 1 to 4 do
W[i,j]:=A[i,j];
W[i,j+4]:=B[i,j];
W[i+4,j]:=C[i,j];
W[i+4,j+4]:=DD[i,j]
od
od:
#The following must be the identity matrix
simplify(evalm(W^2/c));
#We have now to compute the constant c. We think it's not far from the
#resultant of g and h. That's the reason why we have chosen c to be
#cte*resultant(g,h). Hence we have to compute cte.
#For this, we use the following remark:
#if Wi is the matrix of the addition of the 2-torsion element Ai, and
#Wi^2=ci*I (the ci are not equal a priori because the constant c depends
#of g and h), there is a constant C such that W1(A2)=C*W2(A1). In the
#following, Ai is defined by the polynomial Gi, and Ki represents its
#coordinates in the Kummer.
#kumcoord computes the Kummer coordinates of an element of 2-torsion
#given by the polynomial G1
kumcoord:=proc(G1,H1)
local a0,a1,a2,a3,b0,b1,c0,c1,c2,c3,f5,f6,f7;
a0:=1:a1:=-coeff(G1(x),x,2)/coeff(G1(x),x,3):
a2:=coeff(G1(x),x,1)/coeff(G1(x),x,3):
a3:=-coeff(G1(x),x,0)/coeff(G1(x),x,3):
b0:=0:b1:=0:
f5:=coeff(G1(x)*H1(x),x,5):
f6:=coeff(G1(x)*H1(x),x,6):
f7:=coeff(G1(x)*H1(x),x,7):
c0:=a0*b0^2-f7*a1^3+f7*a1*a2-f6*a1^2+3*f7*a3+2*f6*a2:
c1:=a1*b0^2+2*a0*b0*b1-f7*a1^4+3*f7*a2*a1^2-f6*a1^3-f7*a2^2-f7*a3*a1+
2*f6*a2*a1-f5*a1^2+2*f5*a2:
c2:=a0*b1^2-a2*b0^2+f7*a2*a1^3-2*f7*a2^2*a1+f6*a2*a1^2+f7*a3*a2-f6*a2^2+
f5*a1*a2-3*f5*a3:
c3:=a1*b1^2+2*a2*b0*b1+a3*b0^2+f7*a2^2*a1^2-f7*a3*a1^3+f7*a3*a2*a1-
f7*a2^3+f6*a2^2*a1-f6*a3*a1^2+f5*a2^2-f5*a3*a1:
[a0,a1,a2,a3,c0,c1,c2,c3]:
end:
G1:=x->(gg3*x^2+gg2*x+gg1)*(x-e1):
H1:=x->(x-e2)*(hh4*x^3+hh3*x^2+hh2*x+hh1):
K1:=kumcoord(G1,H1):
G2:=x->(gg3*x^2+gg2*x+gg1)*(x-e2):
H2:=x->(x-e1)*(hh4*x^3+hh3*x^2+hh2*x+hh1):
K2:=kumcoord(G2,H2):
#We have
W1:=simplify(subs(
g0=coeff(G1(x),x,0),g1=coeff(G1(x),x,1),g2=coeff(G1(x),x,2),
g3=coeff(G1(x),x,3),h0=coeff(H1(x),x,0),h1=coeff(H1(x),x,1),
h2=coeff(H1(x),x,2),h3=coeff(H1(x),x,3),h4=coeff(H1(x),x,4),
cte=cte1,evalm(W))):
W2:=simplify(subs(g0=coeff(G2(x),x,0),g1=coeff(G2(x),x,1),
g2=coeff(G2(x),x,2),g3=coeff(G2(x),x,3),h0=coeff(H2(x),x,0),
h1=coeff(H2(x),x,1),h2=coeff(H2(x),x,2),h3=coeff(H2(x),x,3),
h4=coeff(H2(x),x,4),cte=cte2,evalm(W))):
W1K2:=simplify(evalm(W1&*K2)):
W2K1:=simplify(evalm(W2&*K1)):
#The constant C (such that W1*K2=C*W2*K1) can be obtained by
C:=evalm(W2K1[2]/W1K2[2]):
#The fifth and sixth line of W1*K2-C*W2*K1 then gives two equations in
#cte1 and cte2
eq1:=simplify(W1K2[5]-C*W2K1[5]):
eq2:=simplify(W1K2[6]-C*W2K1[6]):
System5:=[[coeff(eq1,cte1,1),coeff(eq1,cte2,1)],
[coeff(eq2,cte1,1),coeff(eq2,cte2,1)]]:
System6:=-[subs(cte1=0,cte2=0,eq1),subs(cte1=0,cte2=0,eq2)]:
sol5:=simplify(evalm(System5^(-1)&*System6)):
#we then obtain that cte is equal to h4. We can now replace the constant
#cte in the matrix W.
W:=simplify(subs(cte=h4,evalm(W))):
#By this way the matrix W is entirely determined.