#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.