################################################################################
################################################################################
##
## tow2cov.g
##
## Date: 26 May 2001
## Author: Nils Bruin
##
## This file scripts the main computations referred to in the article
##
## Towers of 2-covers of hyperelliptic curves by Nils Bruin and Victor Flynn
##
## This file needs KASH 2.2 and ALGAE
##
################################################################################
################################################################################

P1:=30*x;P2:=x-1;P3:=x+1;P4:=x^2-5*x+1;
F:=P1*P2*P3*P4;

#the curve we are considering is given by C: y^2=F(x) 

xcoors:=[infinity,0,6,1/2,-1/8,1/3,5,1,-1];

# First we determine which twists of the covering collection are everywhere
# locally solvable.

Norms1:=FilterTwists(P1,P2*P3*P4);
Norms2:=FilterTwists(P2,P1*P3*P4);
Norms3:=FilterTwists(P3,P1*P2*P4);
Norms4:=FilterTwists(P4,P1*P2*P3);
DL1:=Filtered(Cartesian(Norms1,Norms2,Norms3,Norms4),u->IntIsSquare(Product(u)));
BadRatPlaces:=PlaceSupport(PolyDisc(P1*P2*P3*P4));
DL2_0:=Filtered(DL1,d->ForAll(BadRatPlaces,
         plc->HypEllLocRatXTest(plc,[d[1]*P1,d[2]*P2,d[3]*P3,d[4]*P4])));
DL2_1:=Filtered(DL2_0,d->ForAll(BadRatPlaces,
         plc->HypEllLocRatXTest(plc,[d[1]*d[2]*P1*P2,d[3]*P3*d[4]*P4])));
DL2_2:=Filtered(DL2_1,d->ForAll(BadRatPlaces,
         plc->HypEllLocRatXTest(plc,[d[1]*d[3]*P1*P3,d[2]*P2*d[4]*P4])));
DL2_3:=Filtered(DL2_2,d->ForAll(BadRatPlaces,
         plc->HypEllLocRatXTest(plc,[d[2]*d[3]*P2*P3,d[1]*P1*d[4]*P4])));
DL2_4:=Filtered(DL2_3,d->ForAll(BadRatPlaces,
         plc->HypEllLocRatXTest(plc,[d[1]*P1,d[2]*P2*d[3]*P3,d[4]*P4])));
DL2_5:=Filtered(DL2_4,d->ForAll(BadRatPlaces,
         plc->HypEllLocRatXTest(plc,[d[2]*P2,d[1]*P1*d[3]*P3,d[4]*P4])));
DL2_6:=Filtered(DL2_5,d->ForAll(BadRatPlaces,
         plc->HypEllLocRatXTest(plc,[d[3]*P3,d[1]*P1*d[2]*P2,d[4]*P4])));
DL2:=DL2_6;

# Here we've exhausted all we can do over Q. DL2 now
# contains a list of twists we want to consider further

O:=OrderMaximal(x^2-5*x+1);
OrderClassGroup(O);
alpha:=XOrderPrimElt(O);

P41:=x-alpha;
P42:=P4/P41;

BadOPlaces:=PlaceSupport(PolyDisc(P1*P2*P3*P4)+RingZero(O));
Sort(BadOPlaces,function(a,b)return a.primebelow<b.primebelow;end);
tmp:=SquareBasis(BadOPlaces);
tmp2:=Cartesian(List([1..Length(tmp)],i->[0,1]));
K2S:=List(tmp2,e->Product([1..Length(tmp)],i->tmp[i]^e[i]));
tmp:=List(DL2,d->[d[1],d[2],d[3],Filtered(K2S,dk->IntIsSquare(Norm(dk)*d[4]))]);
DL3:=[];
for tmp2 in tmp do
  Append(DL3,List(tmp2[4],d->[tmp2[1],tmp2[2],tmp2[3],d,Norm(d)/d]));
od;

#DL3 essentially contains the same info as DL2, but split over O.
#In fact, we could initialise DL3 with the info obtained from a
#2-descent on C. That is hardly more efficient than this, though.

rlrts:=List([P1,P2,P3,P41,P42],
  u->Re(XEltCon(-PolyToList(u)[2]/PolyToList(u)[1],1)));
perm:=[1,2,3,4,5];SortParallel(rlrts,perm);
idx:=[1,2,3,4,5];SortParallel(perm,idx);
posvals:=List([[true,true,true,true,true],
               [true,true,true,false,false],
               [true,false,false,false,false]],L->List(idx,i->L[i]));

# We use the real places ...

DL4:=Filtered(DL3,d->List(d,e->Re(XEltCon(e,1))>0) in posvals);
DL5_1:=Filtered(DL4,d->ForAll(BadOPlaces,
     plc->HypEllLocRatXTest(plc,
              [d[4]*P41,d[5]*P42*d[1]*P1*d[2]*P2*d[3]*P3])));
DL5_2:=Filtered(DL5_1,d->ForAll(BadOPlaces,
     plc->HypEllLocRatXTest(plc,
              [d[4]*P41*d[1]*P1,d[5]*P42*d[2]*P2*d[3]*P3])));
DL5_3:=Filtered(DL5_2,d->ForAll(BadOPlaces,
     plc->HypEllLocRatXTest(plc,
              [d[4]*P41*d[2]*P2,d[5]*P42*d[1]*P1*d[3]*P3])));
DL5_4:=Filtered(DL5_3,d->ForAll(BadOPlaces,
     plc->HypEllLocRatXTest(plc,
              [d[4]*P41*d[3]*P3,d[5]*P42*d[1]*P1*d[2]*P2])));
DL5_5:=Filtered(DL5_4,d->ForAll(BadOPlaces,
     plc->HypEllLocRatXTest(plc,
              [d[4]*P41,d[1]*P1*d[3]*P3,d[5]*P42*d[2]*P2])));
DL5_6:=Filtered(DL5_5,d->ForAll(BadOPlaces,
     plc->HypEllLocRatXTest(plc,
              [d[4]*P41*d[3]*P3,d[2]*P2,d[5]*P42*d[1]*P1])));

deltas:=DL5_6;

# We're finally there.

List(xcoors{[2..9]},x0->[x0,Filtered(deltas,
  d->ForAll([d[1]*P1,d[2]*P2,d[3]*P3,d[4]*P41,d[5]*P42],
      P->EltIsSquare(Eval(P,x0))))[1]]);

# We match up the twists with x-coordinates (just to check)

# and we construct the elliptic subcovers

E1a:=List(deltas,d->Quar(Eval(P1*P2*P3,x)/d[1]/d[2]/d[3]/30^2));
E1b:=List(deltas,d->Quar(Eval(P1*P2*P3,x-1)/d[1]/d[2]/d[3]/30^2));
E1c:=List(deltas,d->Quar(Eval(P1*P2*P3,x+1)/d[1]/d[2]/d[3]/30^2));

# We do the descents on the elliptic curves to get the rank bounds

E1ranks:=List([1..9],
  i->Minimum([EllRankBound(E1a[i]),
              EllRankBound(E1b[i]),
	      EllRankBound(E1c[i])]));
E2:=List(deltas,d->Quar(Eval(P1*P4,x)/d[1]/Norm(d[4])/900));
E2ranks:=List(E2,EllRankBound);

E3:=List(deltas,d->Quar(Eval(P2*P4,x+1)/d[2]/Norm(d[4])));
E3ranks:=List(E3,EllRankBound);

E4:=List(deltas,d->Quar(Eval(P3*P4,x-1)/d[3]/Norm(d[4])));
E4ranks:=List(E4,EllRankBound);

pl5:=Eval(P1*P2*P4,1/(x+1))*(x+1)^4;
E5:=List(deltas,d->Quar(Eval(pl5,x)/d[1]/d[2]/Norm(d[4])/30^2));
E5ranks:=List(E5,EllRankBound);

pl6:=Eval(P1*P3*P4,1/(x-1))*(x-1)^4;
E6:=List(deltas,d->Quar(Eval(pl6,x)/d[1]/d[3]/Norm(d[4])/30^2));
E6ranks:=List(E6,EllRankBound);

pl7:=Eval(P2*P3*P4,6/(x-3)+1)*(x-3)^4/2^2/3^2;
E7:=List(deltas,d->Quar(Eval(pl7,x)/d[2]/d[3]/Norm(d[4])/3^2));
E7ranks:=List(E7,EllRankBound);

E8a:=List(deltas,d->Quar(Eval(P2*P3*P41,x-1)/d[2]/d[3]/d[4]));
E8b:=List(deltas,d->Quar(Eval(P2*P3*P41,x+1)/d[2]/d[3]/d[4]));
E8c:=List(deltas,d->Quar(Eval(P2*P3*P41,x+PolyRoots(P41)[1])/d[2]/d[3]/d[4]));
E8ar:=List(E8a,EllRankBound);
E8br:=List(E8b,EllRankBound);
E8cr:=List(E8c,EllRankBound);
E8ranks:=List([1..9],
  i->Minimum([EllRankBound(E8a[i]),
              EllRankBound(E8b[i]),
	      EllRankBound(E8c[i])]));
E9a:=List(deltas,d->Quar(Eval(P1*P3*P41,x)/d[2]/d[5]/900));
E9b:=List(deltas,d->Quar(Eval(P1*P3*P41,x-1)/d[2]/d[5]/900));
E9c:=List(deltas,d->Quar(Eval(P1*P3*P41,x+PolyRoots(P41)[1])/d[2]/d[5]/900));
E9ranks:=List([1..9],
  i->Minimum([EllRankBound(E9a[i]),
              EllRankBound(E9b[i]),
	      EllRankBound(E9c[i])]));
E10a:=List(deltas,d->Quar(Eval(P1*P2*P41,x)/d[3]/d[5]/900));
E10b:=List(deltas,d->Quar(Eval(P1*P2*P41,x+1)/d[3]/d[5]/900));
E10c:=List(deltas,d->Quar(Eval(P1*P2*P41,x+PolyRoots(P41)[1])/d[3]/d[5]/900));
E10ranks:=List([1..9],
  i->Minimum([EllRankBound(E10a[i]),
              EllRankBound(E10b[i]),
	      EllRankBound(E10c[i])]));
pl11:=x^4*Eval(P1*P2*P3*P41,1/x)/(1-alpha)^2;
E11a:=List(deltas,d->Quar(Eval(pl11,x-1)/d[5]/10));
E11b:=List(deltas,d->Quar(Eval(pl11,x+1)/d[5]/10));
E11c:=List(deltas,d->Quar(Eval(pl11,x-(alpha-5))/d[5]/10));
E11ranks:=List([1..9],
  i->Minimum([EllRankBound(E11a[i]),
              EllRankBound(E11b[i]),
	      EllRankBound(E11c[i])]));
[E1ranks,E2ranks,E3ranks,E4ranks,E5ranks,E6ranks,E7ranks,E8ranks,E9ranks,
E10ranks,E11ranks];

# The ranks that we get ...
#[ [ 1, 1, 2, 1, 1, 1, 1, 0, 2 ],
#  [ 1, 1, 1, 2, 2, 1, 1, 1, 0 ], 
#  [ 1, 1, 2, 1, 2, 1, 1, 0, 1 ],
#  [ 2, 1, 1, 1, 2, 2, 2, 1, 0 ], 
#  [ 1, 1, 2, 1, 1, 1, 1, 0, 1 ],
#  [ 1, 1, 2, 2, 2, 1, 2, 1, 2 ], 
#  [ 1, 1, 2, 1, 1, 1, 1, 1, 1 ],
#  =============================
#  [ 2, 2, 2, 1, 2, 1, 2, 2, 1 ], 
#  [ 2, 1, 2, 1, 2, 1, 1, 1, 1 ],
#  [ 0, 1, 3, 2, 3, 1, 1, 0, 1 ], 
#  [ 2, 0, 2, 0, 0, 0, 2, 1, 1 ] ]

EllTorBound_Nplaces:=20;
List([E1a,E2,E3,E4,E5,E6,E7,E8a,E9a,E10a,E11a],u->List(u,EllTorBound));
# and, just to check, the torsion bounds ...
#[ [ 4, 4, 4, 4, 4, 4, 4, 4, 4 ],
#  [ 4, 4, 4, 4, 4, 4, 4, 4, 4 ], 
#  [ 2, 2, 2, 2, 2, 2, 2, 2, 2 ],
#  [ 2, 2, 2, 2, 2, 2, 2, 2, 2 ], 
#  [ 2, 2, 2, 2, 2, 2, 2, 2, 2 ],
#  [ 2, 2, 2, 2, 2, 2, 2, 2, 2 ], 
#  [ 2, 2, 2, 2, 2, 2, 2, 2, 2 ],
#  =============================
#  [ 4, 4, 4, 4, 4, 4, 4, 4, 4 ], 
#  [ 4, 4, 4, 4, 4, 4, 4, 4, 4 ],
#  [ 4, 4, 4, 4, 4, 4, 4, 4, 4 ], 
#  [ 4, 4, 4, 4, 4, 4, 4, 4, 4 ] ]

#The one curve where we can do the EC-chabauty and really need it:

E9a:=List(deltas,d->Quar(Eval(P1*P3*P41,x)/d[2]/d[5]/900));
ec:=E9a[7];
EllAddHint(ec,5);
EllAddHint(ec,-1);
EllGensMod2(ec);
EllGenInit([EllXtoPnt(ec,5),EllXtoPnt(ec,0),EllXtoPnt(ec,-1)],1);
delta:=deltas[7];
cov:=QuarCov(P1*P3*P41/delta[2]/delta[5],[0,1,0],ec);
EllCovChab(cov,37);

#The remaining case, where we have to repeat the procedure:

delta:=deltas[3];
d1:=delta[1]*delta[3]*delta[4];
d2:=delta[1]*delta[3]*delta[5];
d0:=Norm(delta[4]);
T:=PolyToList(P4)[2];
Rpol:=P1*P3;
Fpol:=Rpol*(T+4*Norm(d1)*y0);
parm:=[(7-14*z+6*z^6)/(z^2-7),(-z^2+7*z-7)/(z^2-7)];

############################################################################
## Best to restart KASH here
############################################################################

# same exercise as before ...
P1:=42*x;P2:=x^2-14*x+21;P3:=x^2+14*x+21;
F:=P1*P2*P3;
xcoors:=[infinity,0,-7,21];
Norms1:=FilterTwists(P1,P2*P3);
Norms2:=FilterTwists(P2,P1*P3);
Norms3:=FilterTwists(P3,P1*P2);
DL1:=Filtered(Cartesian(Norms1,Norms2,Norms3),u->IntIsSquare(Product(u)));
BadRatPlaces:=PlaceSupport(PolyDisc(P1*P2*P3));
DL2_0:=Filtered(DL1,d->ForAll(BadRatPlaces,
         plc->HypEllLocRatXTest(plc,[d[1]*P1,d[2]*P2,d[3]*P3])));
DL2_1:=Filtered(DL2_0,d->ForAll(BadRatPlaces,
         plc->HypEllLocRatXTest(plc,[d[1]*P1,d[2]*P2*d[3]*P3])));
DL2_2:=Filtered(DL2_1,d->ForAll(BadRatPlaces,
         plc->HypEllLocRatXTest(plc,[d[1]*P1*d[2]*P2,d[3]*P3])));
DL2_3:=Filtered(DL2_2,d->ForAll(BadRatPlaces,
         plc->HypEllLocRatXTest(plc,[d[1]*P1*d[3]*P3,d[2]*P2])));
DL2:=DL2_3;
O:=OrderMaximal(x^2-7);
OrderClassGroup(O);
beta:=Elt(O,[0,1]);
P21:=x-(7+2*beta);
P22:=P2/P21;
P31:=x-(-7+2*beta);
P32:=P3/P31;
BadOPlaces:=PlaceSupport(PolyDisc(P1*P2*P3)+RingZero(O));
Sort(BadOPlaces,function(a,b)return a.primebelow<b.primebelow;end);
tmp:=SquareBasis(BadOPlaces);
tmp2:=Cartesian(List([1..Length(tmp)],i->[0,1]));
K2S:=List(tmp2,e->Product([1..Length(tmp)],i->tmp[i]^e[i]));
tmp:=List(DL2,d->[d[1],Filtered(K2S,dk->IntIsSquare(Norm(dk)*d[2])),
                       Filtered(K2S,dk->IntIsSquare(Norm(dk)*d[3]))]);
DL3_1:=[];
for tmp2 in tmp do
  Append(DL3_1,List(tmp2[2],d->[tmp2[1],d,Norm(d)/d,tmp2[3]]));
od;
DL3_2:=[];
for tmp2 in DL3_1 do
  Append(DL3_2,List(tmp2[4],d->[tmp2[1],tmp2[2],tmp2[3],d,Norm(d)/d]));
od;
#Note that there is more than 2 real embeddings, but they don't give
#more information, since P2 has only positive roots and
#P3 has only negative roots.
DL3:=DL3_2;
rlrts:=List([P1,P21,P22,P31,P32],
  u->Re(XEltCon(-PolyToList(u)[2]/PolyToList(u)[1],1)));
perm:=[1,2,3,4,5];SortParallel(rlrts,perm);
idx:=[1,2,3,4,5];SortParallel(perm,idx);
posvals:=List([[true,true,true,true,true],
               [true,true,true,false,false],
               [true,false,false,false,false]],L->List(idx,i->L[i]));
DL4:=Filtered(DL3,d->List(d,e->Re(XEltCon(e,1))>0) in posvals);
DL5_1:=Filtered(DL4,d->ForAll(BadOPlaces,
     plc->HypEllLocRatXTest(plc,
              [d[2]*P21,d[1]*P1*d[3]*P22*d[4]*P31*d[5]*P32])));
DL5_2:=Filtered(DL5_1,d->ForAll(BadOPlaces,
     plc->HypEllLocRatXTest(plc,
              [d[4]*P31,d[1]*P1*d[2]*P21*d[3]*P22*d[5]*P32])));
DL5_3:=Filtered(DL5_2,d->ForAll(BadOPlaces,
     plc->HypEllLocRatXTest(plc,
              [d[2]*P21*d[4]*P31,d[1]*P1*d[3]*P22*d[5]*P32])));
DL5_4:=Filtered(DL5_3,d->ForAll(BadOPlaces,
     plc->HypEllLocRatXTest(plc,
              [d[2]*P21*d[5]*P32,d[1]*P1*d[3]*P22*d[4]*P31])));
DL5_5:=Filtered(DL5_4,d->ForAll(BadOPlaces,
     plc->HypEllLocRatXTest(plc,
              [d[2]*P21*d[4]*P31*d[5]*P32,d[1]*P1*d[3]*P22])));
deltas:=DL5_5;

# so delta now contains everywhere locally solvable twists

List(xcoors{[2..4]},x0->[x0,Filtered(deltas,
  d->ForAll([d[1]*P1,d[2]*P21,d[3]*P22,d[4]*P31,d[5]*P32],
      P->EltIsSquare(Eval(P,x0))))[1]]);

# and they pair up with the known x-coordinates

ec1:=Quar(P2*P3);
E1:=List(deltas,d->Quar(d[1]*P2*P3,Quar(EllYsqr(ec1,x)/d[1])));
E1ranks:=List(E1,EllRankBound);
#[ 0, 0, 1, 1 ]

E2:=List(deltas,d->Quar(Eval(P1*P2/d[1]/Norm(d[2])*14^2,x/14)));
E2ranks:=List(E2,EllRankBound);
#[ 0, 1, 1, 2 ]

#Note that x -> -x  interchanges E2 and E3 (at least geometrically)

E3:=List(deltas,d->Quar(Eval(P1*P3/d[1]/Norm(d[4])*14^2,x/14)));
E3ranks:=List(E3,EllRankBound);
#[ 1, 1, 1, 1 ]

E4:=List(deltas,d->Quar(Eval(P1*P21*P31/d[1]/d[2]/d[4]/14,x)));
E4ranks:=List(E4,EllRankBound);
#[ 2, 1, 2, 1 ]

E5:=List(deltas,d->Quar(Eval(P1*P21*P32/d[1]/d[2]/d[5]/14,x)));
E5ranks:=List(E5,EllRankBound);
#[ 1, 1, 1, 2 ]

E6:=List(deltas,d->Quar(Eval(P21*P22*P31/d[1]/d[5],x)));
E6a:=List(deltas,d->Quar(Eval(P21*P22*P31/d[1]/d[5],x+PolyRoots(P21)[1])));
E6aranks:=List(E6a,EllRankBound);
#[ 0, 1, 1, 1 ]

ec7a:=Quar((x^3+(22-8*beta)*x^2+(224-88*beta)*x)/deltas[1][5]);
E7:=List(deltas,d->Quar(P1*P2*P31/d[5]));
E7a:=List(deltas,d->Quar(P1*P2*P31/d[5],Quar(EllYsqr(ec7a,x)/d[5])));
E7aranks:=List(E7a,EllRankBound);
E7b:=List(deltas,d->Quar(P1*P2*P31/d[5],Quar(EllYsqr(ec,x-(8-4*beta))/d[5])));
E7branks:=List(E7b,EllRankBound);
E7c:=List(deltas,d->Quar(P1*P2*P31/d[5],Quar(EllYsqr(ec,x-(14-4*beta))/d[5])));
E7cranks:=List(E7c,EllRankBound);
#[ 1, 1, 1, 2 ]

ec8a:=Quar((x^3+(33+12*beta)*x^2+(504+198*beta)*x));
E8:=List(deltas,d->Quar(P1*P21*P3/d[3]));
E8a:=List(deltas,d->Quar(P1*P21*P3/d[3],Quar(EllYsqr(ec8a,x)/d[3])));
#[ 1, 0, 1, 1 ]

E9:=List(deltas,d->Quar(Eval(P21*P31*P32/d[1]/d[3],x)));
E9a:=List(deltas,d->Quar(Eval(P21*P31*P32/d[1]/d[3],x+PolyRoots(P21)[1])));
E9aranks:=List(E9a,EllRankBound);
E9b:=List(deltas,d->Quar(Eval(P21*P31*P32/d[1]/d[3],x+PolyRoots(P31)[1])));
E9branks:=List(E9b,EllRankBound);
E9c:=List(deltas,d->Quar(Eval(P21*P31*P32/d[1]/d[3],x+PolyRoots(P32)[1])));
E9cranks:=List(E9c,EllRankBound);
#[ 0, 1, 1, 2 ]

ec:=E6a[3];
EllAddHint(ec,(14-10*beta));
EllGensMod2(ec);
EllGenInit([EllXtoPnt(ec,(14-10*beta)),
    EllXtoPnt(ec,-4*beta),EllXtoPnt(ec,-14)],1);
cov:=QuarCov(P21*P22*P31/deltas[3][1]/deltas[3][5],21,ec);
cov:=EllCov(ec,cov.f);
EllCovChab(cov,19);

ec:=E6a[4];
EllAddHint(ec,(-14-2*beta));
EllAddHint(ec,-4*beta);
EllAddHint(ec,-14);
EllGenInit([EllXtoPnt(ec,(-14-2*beta)),
    EllXtoPnt(ec,-4*beta),EllXtoPnt(ec,-14)],1);

cov:=QuarCov(P21*P22*P31/deltas[4][1]/deltas[4][5],-7,ec);
cov:=EllCov(ec,cov.f);
EllCovChab(cov,19,5);