#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