#This program compute biquadratic forms for hyperelliptic curves of #genus 3 defined by a polynomial of degree 7. These computations are #very long and need a lot of memory (I can do them with magma with a lot #of memory (see biquadratic2). So that it is much simpler to consider an #exemple). I choose this one because torsion points are easy to compute #butthe user can change the function without any problem. Moreover the #matrix which gives addition by a divisor of order 2 is given by the #previous program and is denoted W in the following. Be carefull, this #program is long to execute even with any given f. It took about 45 mn #on a pentium 433. f:=x->x*(x-1)*(x-2)*(x-3)*(x-4)*(x-5)*(x-6):i:='i': for i from 0 to 7 do f.i:=coeff(f(x),x,i); od: #sym3 transform a symetric polynomial in three variables in a polynomial #in the elementary symmetric functions s=x1+x2+x3, b=x1*x2+x1*x3+x2*x3 #and p=x1*x2*x3. s:='s':b:='b':p:='p': sym3:=proc(P) local Q,Q1,Q2,Q3,Q4,R,R1,R2,R3,u,v,w,x,y,t1,t2,s12,b12; Q:=P: t1:=simplify(Q-subs(x1=s2,x2=x1,s2=x2,Q)); t2:=simplify(Q-subs(x1=s2,x3=x1,x2=x3,s2=x2,Q)); if t1=0 and t2=0 then else print("polynomial is not symmetric"):break:fi: while Q <> subs(x1=u,x2=v,x3=w,Q) do #test if Q depends of x1, x2 and x3 R:=subs(x3=0,Q): while R <> subs(x1=x,x2=y,R) do #test if Q depends of x1, x2 R1:=subs(x2=0,R); R2:=R-subs(x1=x1+x2,R1); R3:=simplify(R2/x1/x2); R:=subs(x1=s12,R1)+b12*R3 od: Q2:=subs(s12=x1+x2+x3,b12=x1*x2+x1*x3+x2*x3,R); Q3:=Q-Q2; Q4:=simplify(Q3/(x1*x2*x3)); Q:=subs(s12=s,b12=b,R)+p*Q4 od: Q end: #We first transform the matrix W (constructed with the previous program) #which gives addition by a divisor of order 2. This matrix is given in #terms of the coefficients of g and h, where g is cubic (defining the #divisor of order 2) and h is of degree 4. We first transform to the #coefficients of G and those of F=GH, eliminating those of H. We shall #also put g3=1, g2=-s, g1=b and g0=-p. To eliminate the h's, it is #convenient to define the procedure subsh. subsh:=proc(P) local Q; Q:=P; Q:=subs(h0=f3+s*h1-b*h2+p*h3,Q); Q:=subs(h1=f4+s*h2-b*h3+p*h4,Q); Q:=subs(h2=f5+s*h3-b*h4,Q); Q:=subs(h3=f6+s*h4,Q); Q:=subs(h4=f7,Q); Q:=simplify(Q); Q end: if length(evalm(W))=1 then print("You have probably forgotten to compute the matrix W, you must do it before with the program matrixW or by reading the file matrixWr") fi: #Let us now transform the matrix W for i from 1 to 8 do for j from 1 to 8 do W[i,j]:=simplify(subs(g2=-s,g1=b,g0=-p,g3=1,subsh(W[i,j]))) od od: #We now want g=0 to be a generic divisor of order 2. Since we no longer #need the original uses of x1, x2, x3, we use them as a triplet of #distinct zeros of f. We need to express that they are distinct. B123:=simplify((x1*f(x2)-x2*f(x1)-x3*f(x2)+x3*f(x1)-x1*f(x3)+x2*f(x3)) /((x1-x2)*(x2-x3)*(x3-x1))): B23:=simplify((f(x2)-f(x3))/(x2-x3)): #x1 is a root of f of degree 7. x2 and x3 are other roots of f, so #satisfies equation B123 of degree 5 and B23 of degree 6. #By using B123 and B23 we can reduce a symmetric poly in x1, x2, x3 to #degree at most 4 in x3. By symmetry it is of degree at most 4 in x1 #and x2 and so of the shape s^i*b^j*p^k, with i+j+k at most 4. #The following program takes a symmetric function of x1, x2, x3 and #makes the reduction. simp:=proc(P) local Q; Q:=P; Q:=rem(Q,B123,x1); Q:=rem(Q,B23,x2); Q:=rem(Q,f(x3),x3); Q:=simplify(Q); Q end: #Any polynomial in s, b, p is equal to a linear form in the elements #s^i*b^j*p^k with i+j+k at most 4. We need a routine to do this. simpsp:=proc(P) local Q; Q:=P; Q:=subs(s=x1+x2+x3,b=x1*x2+x1*x3+x2*x3,p=x1*x2*x3,Q); Q:=simp(Q); Q:=sym3(Q); Q end: #Let us now compute the coordinate of {(x1,0),(x2,0),(x3,0)} a0:=1: a1:=x1+x2+x3: a2:=x1*x2+x1*x3+x2*x3: a3:=x1*x2*x3: c0:=-f7*s^3+f7*b*s-f6*s^2+3*f7*p+2*f6*b: c1:=-f7*s^4+3*f7*b*s^2-f6*s^3-f7*b^2-f7*p*s+2*f6*b*s-f5*s^2+2*f5*b: c2:=f7*b*s^3-2*f7*b^2*s+f6*b*s^2+f7*b*p-f6*b^2+f5*b*s-3*f5*p: c3:=f7*b^2*s^2-f7*p*s^3+f7*p*b*s-f7*b^3+f6*b^2*s-f6*p*s^2+f5*b^2-f5*p*s: #We want now to compute the product of all these coordinates k:=[1,s,b,p,c0,c1,c2,c3]: remp:=(i,j)->simpsp(k[i]*k[j]): K2:=matrix(8,8,remp): #Let us now write these products in the vector V. V:=[K2[1,1],K2[1,2],K2[1,3],K2[1,4],K2[1,5],K2[1,6],K2[1,7], K2[2,2],K2[2,3],K2[2,4],K2[2,5],K2[2,6],K2[2,7],K2[2,8], K2[3,3],K2[3,4],K2[3,5],K2[3,6],K2[3,7],K2[3,8], K2[4,4],K2[4,5],K2[4,6],K2[4,7],K2[4,8], K2[5,5],K2[5,6],K2[5,7],K2[5,8], K2[6,6],K2[6,7],K2[6,8], K2[7,7],K2[7,8], K2[8,8]]: #the vector V corresponds to the following products Z:=[A0*A0,A0*A1,A0*A2,A0*A3,A0*C0,A0*C1,A0*C2, A1*A1,A1*A2,A1*A3,A1*C0,A1*C1,A1*C2,A1*C3, A2*A2,A2*A3,A2*C0,A2*C1,A2*C2,A2*C3, A3*A3,A3*C0,A3*C1,A3*C2,A3*C3, C0*C0,C0*C1,C0*C2,C0*C3, C1*C1,C1*C2,C1*C3, C2*C2,C2*C3, C3*C3]: #Hence V[i] is Z[i] expressed with s,b and p. Note that in these vectors #there is not the product a0c3. That is because of the relation #a0c3-a1c2-a2c1-a3c0-2f5*a1*a3+f5*a2^2+2*f6*a2*a3+3*f7*a3^2=0 #In fact opposite to the genus 2 case, all the product are not linearly #independant. Nevertheless, if we replace a0c3 in all the expressions by #a1c2+a2c1+..., the remaining products are independant. Since this #relation always occur (not only for divisors of order 2) we can still #deduce the general case from the 2-torsion case. #We now need a program (called tryon) to convert a linear form in the #s^i*b^j*p^k into a quadratic form in a0,a1,a2,a3,c0,c1,c2,c3 (without #any term in a0c3). For example c3^2 is the only product where p^4 lies. #Hence the term in c3^2 is given by the term of p^4. We can then look at #this problem as a system. Let us construct the matrix M of 35 #coefficients in some s^i*b^j*p^k of each of the 35 products in V. coef:=(cs,cb,cp,es)->(coeff(coeff(coeff(es,s,cs),b,cb),p,cp)): M:=matrix(35,35,0): t:=0: for i from 0 to 4 do for j from 0 to 4-i do for l from 0 to 4-i-j do t:=t+1: for m from 1 to 35 do M[t,m]:=coef(i,j,l,V[m]) od: od od od: Mi:=evalm(M^(-1)): tryon:=proc(P) local Q,R,S,t,i,j,l; Q:=P: R:=matrix(35,1): t:=0: for i from 0 to 4 do for j from 0 to 4-i do for l from 0 to 4-i-j do t:=t+1: R[t,1]:=coef(i,j,l,Q) od od od: #R is the vector of coefficients of s^i*b^j*p^l, so that Mi*R is the #vector of the coefficients of the desired quadratic form in the 35 #products in V. S:=evalm(Z&*evalm(Mi&*R)): Q:=simplify(S[1]): Q end: #We now construct the biquadratic forms. Let us first reduce the #coefficients of the matrix W using simpsp for i from 1 to 8 do for j from 1 to 8 do W[i,j]:=simpsp(W[i,j]) od od: #Let A be a generic divisor [y1,y2,y3,y4,y5,y6,y7,y8]. Unfortunately #we cannot make directly the computations of the products Y[i]*Y[j] #by using Y[i]:=sum('y.i*W[i,j]',j=1..8) because it requires too many #RAM. We can only compute them for little i and j. #So, we need to compute the coefficients of these products one by one. #In other words, we need to compute all the product of the elements of #the matrix W. We write this in a matrix 8*8 (called WW) which #coefficients are 8*8 matrices such that WW[i,j][k,l]=W[i,j]*W[k,l]. #The procedure compare is used to do computations of each W[i,j]*W[k,l] #only one time. compare:=proc(u,v) local a; a:=0: if u[1] ij WW[i,j][k,l]:=tryon(simpsp(expand(W[i,j]*W[k,l]))) fi od od od: print(i,"time=",time()-t) od: print("total time=",time()-T); #Let us now complete the matrix WW by symmetry for i from 1 to 8 do for j from 1 to 8 do for k from 1 to 8 do for l from 1 to 8 do if compare([k,l],[i,j])=1 then WW[i,j][k,l]:=WW[k,l][i,j] fi od od od od: #We can now write the products Y[i]*Y[j] and then deduce the matrix #Phi of the biquadratic forms. Note that, as we choose to replace the #term in a0c3, we have to replace the term in y1y8 #if we want the Y[i]*Y[j] to be symmetric k:='k':i:='i':j:='j':l:='l': Phi:=matrix(8,8,0): for i from 1 to 8 do for j from 1 to 8 do phi:=simplify(sum('sum('WW[i,k][j,l]*y.k',k=1..8)*y.l',l=1..8)): toreplace:=coeff(coeff(phi,y1,1),y8,1): Phi[i,j]:=simplify(phi-toreplace*y1*y8+toreplace*(y2*y7+y3*y6+ y4*y5+2*f5*y2*y4-f5*y3^2-2*f6*y3*y4-3*f7*y4^2)) od od: #We can now check that these Phi[i,j] actually are symmetric in the y's #and the ai,ci although they are treated very differently. test:=matrix(8,8,1): for i from 1 to 8 do for j from 1 to 8 do temp:=Phi[i,j]: temp:=subs(A0=t1,A1=t2,A2=t3,A3=t4,C0=t5,C1=t6,C2=t7,C3=t8, y1=s1,y2=s2,y3=s3,y4=s4,y5=s5,y6=s6,y7=s7,y8=s8,temp): temp:=subs(s1=A0,s2=A1,s3=A2,s4=A3,s5=C0,s6=C1,s7=C2,s8=C3, t1=y1,t2=y2,t3=y3,t4=y4,t5=y5,t6=y6,t7=y7,t8=y8,temp): test[i,j]:=simplify(Phi[i,j]-temp); if test[i,j]=0 then else print("mistake",i,j) fi; od: od: #test must be 0