///////////////////////////////////
//NOTE:
//
//This script has been tested to run correctly on MAGMA 2.11-12. Note that there
//is a dependency on a package file "bugfix.m", available from the same location
//as this file.

Attach("bugfix.m");
ClearVerbose( );
SetAutoColumns(false);
SetColumns(100);
SetEchoInput(true);

////////////////////////////////////////////////////////////////////////////
//Case: x^4+y^4=D*z^2
//
//Section 2, Proof of Proposition 1.

_<x>:=PolynomialRing(Rationals());
Dset:={2,3,5,6,10,11,13,17};
{D:D in Dset| IsLocallySolvable(HyperellipticCurve(D*(x^4+1)),2)};

////////////////

C2:=HyperellipticCurve(2*(x^4+1));
p0:=C2![1,2];
E,C2toE:=EllipticCurve(C2,p0);
RankBound(E);
#TorsionSubgroup(E);

////////////////

C17:=HyperellipticCurve(17*(x^4+1));
p0:=C17![2,17];
E,C17toE:=EllipticCurve(C17,p0);
RankBound(E);

////////////////

G,GtoE:=MordellWeilGroup(E);
G;
[Inverse(C17toE)(GtoE(g)):g in OrderedGenerators(G)];

////////////////////////////////////////////////////////////////////////////
//Case: x^5+y^5=D*z^2
//
//Section 4: The equation x^5+y^5=Dz^2

_<x>:=PolynomialRing(Rationals());
f:=x^5+1;
A<theta>:=quo<Parent(x)|f>;

////////////////

PA<Z0A,Z1A,Z2A,Z3A,Z4A>:=PolynomialRing(A,5);
AP<thetaP>,swap:=SwapExtension(PA);
_<Z0,Z1,Z2,Z3,Z4>:=BaseRing(AP);
P4:=Proj(BaseRing(AP));
P1:=ProjectiveSpace(Rationals(),1);

////////////////

function Cdelta(delta)
  g:=delta*(&+[PA.i*theta^(i-1): i in [1..5]])^2;
  Q:=Coefficients(swap(g));
  Crv:=Scheme(P4,[Q[3],Q[4],Q[5]]);
  phi:=map<Crv->P1|[Q[1],-Q[2]]>;
  return Crv,phi;
end function;

////////////////

Q:=NumberField(x-1:DoLinearExtension);
OQ:=IntegerRing(Q);
Qx<xQ>:=PolynomialRing(Q);
AQ:=quo<Qx|Polynomial(Q,f)>;
AQtoA:=hom<AQ->A|[theta]>;
K<zeta>:=NumberField(x^4 - x^3 + x^2 - x + 1);
OK:=IntegerRing(K);
Aa,toAa:=AbsoluteAlgebra(AQ:Fields:={Q,K});

////////////////

function DeltaForD(D)
  S:=Support(D*Discriminant(f)*OQ);
  slmA,slmAmap:=pSelmerGroup(AQ,2,S);
  slmQ,slmQmap:=pSelmerGroup(2,S);
  slmNorm:=map<slmA->slmQ|a:->slmQmap(Norm(a@@slmAmap))>;
  slmSquareNorm:=Kernel(hom<slmA->slmQ|
          [slmNorm(a):a in OrderedGenerators(slmA)]>);
  classD:=slmAmap(D);
  return {AQtoA( (d-classD)@@slmAmap):d in slmSquareNorm};
end function;

Dset:={2,3,5,6,10,11,13,17};
BigDelta:=&join{DeltaForD(D):D in Dset};
Delta:=BigDelta;
Delta:={delta:delta in Delta| IsLocallySolvable(Cdelta(delta),2:
                     AssumeIrreducible,AssumeNonsingular)};
Delta:={delta:delta in Delta| IsLocallySolvable(Cdelta(delta),5:
                     AssumeIrreducible,AssumeNonsingular)};
Delta:={delta:delta in Delta| IsLocallySolvable(Cdelta(delta),11:
                     AssumeIrreducible,AssumeNonsingular)};

////////////////

m1:=hom<A->Rationals()|-1>;
{PowerFreePart(Norm(delta)/m1(delta),2):delta in Delta};

////////////////

E5:=HyperellipticCurve(5*(x^4-x^3+x^2-x+1));
p0:=E5![-1,5];
ell:=EllipticCurve(E5,p0);
RankBound(ell);
#TorsionSubgroup(ell);

////////////////

Delta:={delta:delta in Delta|PowerFreePart(Norm(delta)/m1(delta),2) ne 5};
#Delta;

////////////////

m2:=hom<A->K|zeta>;
slmK,slmKmap:=pSelmerGroup(2,Support(2*3*5*11*13*17*OK));
dset:={K|(slmKmap(Norm(delta)/m2(delta)))@@slmKmap:delta in Delta};
dset;

////////////////

Kd<d>:=RationalFunctionField(K);
KdX<X>:=PolynomialRing(Kd);
FEd1:=(X^4+zeta*X^3+zeta^2*X^2+zeta^3*X+zeta^4)/d;
Ed1:=HyperellipticCurve(FEd1);
Ed2,toEd2:=EllipticCurve(Ed1,Ed1![-1,0]);
umap:=map<Ed1->P1|[Ed1.1,Ed1.3]>;
FEd:=X^3+(-3*zeta^3-zeta+1)*d*X^2+(-zeta^2-zeta-1)*d^2*X;
Ed<xE,yE,zE>:=EllipticCurve(FEd);
bl,toEd:=IsIsomorphic(Ed2,Ed);
u:=Expand(Inverse(toEd2*toEd)*umap);
u:Minimal;

////////////////////////////////////////////////////////////////////////////
//The following code is not literally mentioned in the article.
//The computations here verify the results in the table at the end of Section 4.
//The methods used are explained in Section 7. The examples occurring in that
//Section occur below.

d:=1;
Ed<X,Y,Z>:=EllipticCurve([0,-3*zeta^3-zeta+1*d,0,(-zeta^2-zeta-1)*d^2,0]);
u:=map<Ed->P1|[-X + (zeta^3 - 1)*d*Z,X+(-zeta^3-zeta)*d*Z]>;
two:=MultiplicationByMMap(Ed,2);
mu:=IsogenyMu(two:Fields:={K});
S2Ed,S2map:=SelmerGroup(two);
gs:=[Ed|[zeta^3,0],[zeta^3+zeta-1,zeta^2],[zeta^3+zeta,zeta^2+1]];
assert S2Ed eq sub<S2Ed|[S2map(mu(g)):g in gs]>;
G:=AbelianGroup([2,4,0]);
mwmap:=map<G->Ed|g:->&+[c[i]*gs[i]:i in [1..#gs]] where c:=Eltseq(g)>;
N,V,R:=Chabauty(mwmap,u,31);
assert N eq #V and PrimeDivisors(R) subset {2};
{EvaluateByPowerSeries(u,mwmap(p)):p in V};

////////////////

d:=-2*zeta^2-2;
Ed<X,Y,Z>:=EllipticCurve([0,(-3*zeta^3-zeta+1)*d,0,(-zeta^2-zeta-1)*d^2,0]);
u:=map<Ed->P1|[-X + (zeta^3 - 1)*d*Z,X+(-zeta^3-zeta)*d*Z]>;
two:=MultiplicationByMMap(Ed,2);
mu:=IsogenyMu(two:Fields:={K});
S2Ed,S2map:=SelmerGroup(two);
gs:=[Ed![0,0],Ed![2-2*zeta^3,0],Ed![4-4*zeta^3,4*zeta^2+4],
   Ed![-8*zeta^3 + 4*zeta^2 - 4*zeta + 8,8*zeta^3 + 4*zeta - 8]];
assert S2Ed eq sub<S2Ed|[S2map(mu(g)):g in gs]>;
G:=AbelianGroup([2,2,0,0]);
mwmap:=map<G->Ed|g:->&+[c[i]*gs[i]:i in [1..#gs]] where c:=Eltseq(g)>;
N,V,R:=Chabauty(mwmap,u,31);
assert N eq #V and PrimeDivisors(R) subset {2};
{EvaluateByPowerSeries(u,mwmap(p)):p in V};

////////////////

d:=2*zeta^3-2*zeta^2-2;
Ed<X,Y,Z>:=EllipticCurve([0,(-3*zeta^3-zeta+1)*d,0,(-zeta^2-zeta-1)*d^2,0]);
u:=map<Ed->P1|[-X + (zeta^3 - 1)*d*Z,X+(-zeta^3-zeta)*d*Z]>;
two:=MultiplicationByMMap(Ed,2);
mu:=IsogenyMu(two:Fields:={K});
S2Ed,S2map:=SelmerGroup(two);
gs:=[Ed![0,0],Ed![-2*zeta^3 - 2*zeta + 2,0],Ed![-2*zeta^3,-4*zeta^2],
    Ed![-2*zeta^3 - 4*zeta + 4,-4*zeta^3 + 4*zeta]];
assert S2Ed eq sub<S2Ed|[S2map(mu(g)):g in gs]>;
G:=AbelianGroup([2,2,0,0]);
mwmap:=map<G->Ed|g:->&+[c[i]*gs[i]:i in [1..#gs]] where c:=Eltseq(g)>;
N,V,R:=Chabauty(mwmap,u,31);
assert N eq #V and PrimeDivisors(R) subset {2};
{EvaluateByPowerSeries(u,mwmap(p)):p in V};

////////////////

d:=-3*zeta^3-7*zeta^2-8*zeta-9;
Ed<X,Y,Z>:=EllipticCurve([0,(-3*zeta^3-zeta+1)*d,0,(-zeta^2-zeta-1)*d^2,0]);
u:=map<Ed->P1|[-X + (zeta^3 - 1)*d*Z,X+(-zeta^3-zeta)*d*Z]>;
two:=MultiplicationByMMap(Ed,2);
mu:=IsogenyMu(two:Fields:={K});
S2Ed,S2map:=SelmerGroup(two);
gs:=[Ed![0,0],
     Ed![-17*zeta^3 + 8*zeta^2 - 5*zeta + 15,0],
     Ed![-4*zeta^3 + 14*zeta^2 + 8*zeta + 15,-147*zeta^3 - 89*zeta + 89],
     Ed![-5*zeta^3 + 6*zeta^2 + 9, -69*zeta^3 + 17*zeta^2 - 34*zeta + 60],
     Ed![-6*zeta^3 + 3*zeta^2 - zeta + 6,57*zeta^3 - 16*zeta^2 + 27*zeta - 51]];
assert S2Ed eq sub<S2Ed|[S2map(mu(g)):g in gs]>;
G:=AbelianGroup([2,2,0,0,0]);
mwmap:=map<G->Ed|g:->&+[c[i]*gs[i]:i in [1..#gs]] where c:=Eltseq(g)>;
N,V,R:=Chabauty(mwmap,u,191);
assert N eq #V and PrimeDivisors(R) subset {2};
{EvaluateByPowerSeries(u,mwmap(p)):p in V};

////////////////

d:=15*zeta^3-5*zeta^2+8*zeta-17;
Ed<X,Y,Z>:=EllipticCurve([0,(-3*zeta^3-zeta+1)*d,0,(-zeta^2-zeta-1)*d^2,0]);
u:=map<Ed->P1|[-X + (zeta^3 - 1)*d*Z,X+(-zeta^3-zeta)*d*Z]>;
two:=MultiplicationByMMap(Ed,2);
mu:=IsogenyMu(two:Fields:={K});
S2Ed,S2map:=SelmerGroup(two);
gs:=[Ed![0,0],
     Ed![-9*zeta^3 - 8*zeta^2 - 7*zeta - 3,0],
     Ed![10*zeta^3 - 14*zeta^2 - zeta - 22,-89*zeta^2 - 58*zeta - 89],
     Ed![zeta^3 - 6*zeta^2 - 3*zeta - 6,-26*zeta^3 + 43*zeta^2 + 9*zeta + 60],
     Ed![-3*zeta^3 - 3*zeta^2 - 3*zeta - 2,24*zeta^3 - 35*zeta^2 - 6*zeta - 51]
];
assert S2Ed eq sub<S2Ed|[S2map(mu(g)):g in gs]>;
G:=AbelianGroup([2,2,0,0,0]);
mwmap:=map<G->Ed|g:->&+[c[i]*gs[i]:i in [1..#gs]] where c:=Eltseq(g)>;
N,V,R:=Chabauty(mwmap,u,191);
assert N eq #V and PrimeDivisors(R) subset {2};
{EvaluateByPowerSeries(u,mwmap(p)):p in V};

////////////////////////////////////////////////////////////////////////////
//Section 6: Mordell-Weil groups of Elliptic curves

_<x>:=PolynomialRing(Rationals());
K<zeta>:=NumberField(x^4-x^3+x^2-x+1);
OK:=IntegerRing(K);
d:=2*zeta^3-2*zeta^2-2;
E<X,Y,Z>:=EllipticCurve([0,(-3*zeta^3-zeta+1)*d,0,(-zeta^2-zeta-1)*d^2,0]);
two:=MultiplicationByMMap(E,2);
mu,tor:=IsogenyMu(two);
S2E,toS2E:=SelmerGroup(two);S2E;

////////////////

Etors,EtorsMap:=TorsionSubgroup(E);
sub<S2E|[toS2E(mu(EtorsMap(g))):g in OrderedGenerators(Etors)]>;

////////////////

V:=RationalPoints(E:Bound:=5);
sub<S2E|[toS2E(mu(P)):P in V]> eq S2E;

////////////////

gs:=[E![0,0],
     E![-2*zeta^3 - 2*zeta + 2,0],
     E![-2*zeta^3,-4*zeta^2],
     E![-2*zeta^3 - 4*zeta + 4,-4*zeta^3 + 4*zeta]];
assert S2E eq sub<S2E|[toS2E(mu(g)):g in gs]>;
G:=AbelianGroup([2,2,0,0]);
mwmap:=map<G->E|g:->&+[c[i]*gs[i]:i in [1..#gs]] where c:=Eltseq(g)>;

////////////////

success,G,mwmap:=PseudoMordellWeilGroup(E);
assert success;

////////////////

phi:=TwoIsogeny(E![0,0]);
Sphi,toSphi:=SelmerGroup(phi);
phihat:=DualIsogeny(phi);
Sphihat,toSphihat:=SelmerGroup(phihat);
4*#S2E, #Sphi, #Sphihat, #TwoTorsionSubgroup(E);

////////////////

//We need to fix the notation of delta to get around the randomness in the
//algorithm
delta:=S2E!toS2E([1,1/2*(5*zeta^3-10*zeta^2+5*zeta),1/4*(-4*zeta^3+11*zeta^2-11*zeta+4)]);
psi:=tor(delta@@toS2E);
XX:=Domain(psi);
C,CtoXX:=Quartic(XX);

////////////////

V:=RationalPoints(C:Bound:=20);
assert #V gt 0;
P:=psi(CtoXX(Rep(V)));P;
assert delta eq toS2E(mu(P));

////////////////////////////////////////////////////////////////////////////
//Section 7: Chabauty methods using elliptic curves

P1:=ProjectiveSpace(Rationals(),1);
u:=map<E->P1|[-X + (zeta^3 - 1)*d*Z,X+(-zeta^3-zeta)*d*Z]>;
V3:=RelevantCosets(mwmap,u,Support(3*OK));
Lambda3:=Kernel(V3[1]);
GmodLambda3:=Codomain(V3[1]);
V3;

////////////////

V11:=RelevantCosets(mwmap,u,Support(11*OK));
V3i11:=CosetIntersection(V3,V11:Weak);
V3i11;

////////////////

P0:=mwmap(G.3+G.4);
u(P0);
theta:=ChabautyEquations(P0,u,mwmap,Support(3*OK));
PrintToPrecision(theta[1],1);"";PrintToPrecision(theta[2],1);"";PrintToPrecision(theta[3],1);

////////////////

N,V,R,C:=Chabauty(mwmap,u,3:Aux:={7});
assert N eq #V;
assert #C[2] eq 0;
R;
V;
{EvaluateByPowerSeries(u,mwmap(P)):P in V};

////////////////

d:=-3*zeta^3-7*zeta^2-8*zeta-9;
E<X,Y,Z>:=EllipticCurve([0,(-3*zeta^3-zeta+1)*d,0,(-zeta^2-zeta-1)*d^2,0]);
P1:=ProjectiveSpace(Rationals(),1);
u:=map<E->P1|[-X + (zeta^3 - 1)*d*Z,X+(-zeta^3-zeta)*d*Z]>;
success,G,mwmap:=PseudoMordellWeilGroup(E);
assert success;
[mwmap(P):P in OrderedGenerators(G)];
N,V,R,C:=Chabauty(mwmap,u,3);
C31,R31:=RelevantCosets(mwmap,u,Support(31*OK));
R:=LCM(R,R31);
Cnew:=CosetIntersection(C,C31:Weak);
assert #Cnew[2] eq 0;
R;
V;
{EvaluateByPowerSeries(u,mwmap(P)):P in V};

////////////////////////////////////////////////////////////////////////////
//Section 8: The equations x^n+y^n=Dz^2 for n=6,7,9,11,13,17
//
//The computations in this section do not occur literally in the article,
//but verify the computations sketched in Section 8.

////////////////////////////////////////////////////////////////////////////
//Case: x^6+y^6=D*z^2

Cs:=[HyperellipticCurve(D*(x^6+1)):D in [2,3,5,6,10,11,13,17]];
Cs:=[C:C in Cs|IsLocallySolvable(C,2) and IsLocallySolvable(C,7)];

f:=2*x^2+2;g:=(x^4-x^2+1);
deltas:={1,2,3,6};
deltas:={delta:delta in deltas|
   IsLocallySolvable(HyperellipticCurve(delta*f),2) and
   IsLocallySolvable(HyperellipticCurve(delta*g),2)};
assert deltas eq {1};
C:=HyperellipticCurve(g);
E:=EllipticCurve(C);
assert RankBound(E) eq 0;
assert #TorsionSubgroup(E) eq 8;
V:=RationalPoints(C:Bound:=1);
assert #V eq 8;
V;

////////////////////////////////////////////////////////////////////////////
//Case: x^7+y^7=D*z^2

QQ:=RationalField();

f:=x^7+1;
g:=ExactQuotient(f,x+1);
deltas:=&join[{QQ!(d@@mp):d in s} where s,mp:=pSelmerGroup(2,Support(7*D*OQ)):
               D in [2,3,5,6,10,11,13,17]];
deltas:={d:d in deltas | d gt 0};
deltas:={d:d in deltas|IsLocallySolvable(HyperellipticCurve(d*g),2)};
deltas:={d:d in deltas|IsLocallySolvable(HyperellipticCurve(d*g),7)};
deltas:={d:d in deltas|IsLocallySolvable(HyperellipticCurve(d*g),11)};
deltas;

C:=HyperellipticCurve(7*g);
J:=Jacobian(C);
assert TwoSelmerGroupData(J) eq 1;
gen:=elt<J|x^2 - 3/2*x + 1, 7/4*x - 7/2, 2>;
B:=Ceiling(Exp(CanonicalHeight(gen)/4+HeightConstant(J:Factor,Effort:=10))+1);
assert RationalPoints(J:Bound:=B) subset {n*gen: n in [-1..1]};
Chabauty(gen,3); //only 1 pair of points (there are no rational WS points)
RationalPoints(C:Bound:=100);

C:=HyperellipticCurve(1*g);
K<beta>:=NumberField(x^3 + x^2 - 2*x - 1);
OK:=IntegerRing(K);
KX<X>:=PolynomialRing(K);
Qpol:=X^2+beta*X+1;
Rpol:=ExactQuotient(Polynomial(K,g),Qpol);
slm,mp:=pSelmerGroup(2, Support(2*Resultant(Rpol,Qpol)*OK));
slmQ,mpQ:=pSelmerGroup(2, Support(Norm(2*Resultant(Rpol,Qpol))*OQ));
H:=Kernel(hom<slm->slmQ|[mpQ(Norm(a@@mp)):a in OrderedGenerators(slm)]>);
p:=2*OK;
Kp,toKp:=Completion(K,p);
KpX:=PolynomialRing(Kp);
toKpX:=hom<KX->KpX|toKp*Coercion(Kp,KpX),KpX.1>;
deltas:=[a:a in H | HasRatXpoint([toKpX(d*Rpol),toKpX(d*Qpol)]) where d:=a@@mp];
[K|delta@@mp:delta in deltas];

P1<s,t>:=ProjectiveSpace(Rationals(),1);
D:=HyperellipticCurve(Rpol);
Dcov:=map<D->P1|[D.1,D.3]>;
E,DtoE:=EllipticCurve(D);
bl,G,mwmap:=PseudoMordellWeilGroup(E);
assert bl;
[mwmap(g):g in OrderedGenerators(G)];
Ecov:=Expand(Inverse(DtoE)*Dcov);
N,V,R:=Chabauty(mwmap,Ecov,3);
assert N eq #V and PrimeDivisors(R) subset {2};
{EvaluateByPowerSeries(Ecov,mwmap(P)):P in V};

////////////////////////////////////////////////////////////////////////////
//Case: x^9+y^9=D*z^2

g1:=x^6-x^3+1;
g2:=x^2-x+1;
deltas:=&join[{<Integers()!(d1@@mp),Integers()!(d2@@mp)>:d1,d2 in s} where 
          s,mp:=pSelmerGroup(2,Support(3*D*OQ)):
               D in [2,3,5,6,10,11,13,17]];
deltas:={d:d in deltas| d[1] gt 0 and d[2] gt 0};
deltas:={d:d in deltas| Seqset(PrimeDivisors(GCD(d[1],d[2]))) subset {3}};
Qp:=pAdicField(2);
QpX<X>:=PolynomialRing(Qp);
deltas:={d:d in deltas| HasRatXpoint([QpX!(d[1]*g1),QpX!(d[2]*g2)])};
Qp:=pAdicField(3);
QpX<X>:=PolynomialRing(Qp);
deltas:={d:d in deltas| HasRatXpoint([QpX!(d[1]*g1),QpX!(d[2]*g2)])};
Qp:=pAdicField(5);
QpX<X>:=PolynomialRing(Qp);
deltas:={d:d in deltas| HasRatXpoint([QpX!(d[1]*g1),QpX!(d[2]*g2)])};
{d[1]:d in deltas};

C:=HyperellipticCurve(1*g1);
J:=Jacobian(C);
assert TwoSelmerGroupData(J) eq 1;
gen:=elt<J|x^2 - x, -1, 2>;
B:=Ceiling(Exp(CanonicalHeight(gen)/4+HeightConstant(J:Effort:=10,Factor))+1);
assert CanonicalHeight(ReducedBasis(RationalPoints(J:Bound:=B))[1]) gt
       CanonicalHeight(gen)/2;
Chabauty(gen,11); //gives 6 points
RationalPoints(C:Bound:=1000);

C:=HyperellipticCurve(3*g1);
J:=Jacobian(C);
assert TwoSelmerGroupData(J) eq 1;
gen:=elt<J|x^2 - x + 1, -3, 2>;
assert CanonicalHeight(ReducedBasis(RationalPoints(J:Bound:=B))[1]) gt
       CanonicalHeight(gen)/2;
Chabauty(gen,11); //gives 2 points
RationalPoints(C:Bound:=1000);

////////////////////////////////////////////////////////////////////////////
//Case: x^11+y^11=D*z^2

f:=x^11+1;
g:=ExactQuotient(f,x+1);
deltas:=&join[{QQ!(d@@mp):d in s} where s,mp:=pSelmerGroup(2,Support(11*D*OQ)):
               D in [2,3,5,6,10,11,13,17]];
deltas:={d:d in deltas | d gt 0};

deltas:={d:d in deltas|IsLocallySolvable(HyperellipticCurve(d*g),2)};
deltas:={d:d in deltas|IsLocallySolvable(HyperellipticCurve(d*g),3)};
deltas:={d:d in deltas|IsLocallySolvable(HyperellipticCurve(d*g),5)};
deltas:={d:d in deltas|IsLocallySolvable(HyperellipticCurve(d*g),7)};
deltas:={d:d in deltas|IsLocallySolvable(HyperellipticCurve(d*g),11)};

deltas;

C:=HyperellipticCurve(g);
D:=HyperellipticCurve(x^6 + x^5 - 6*x^4 - 5*x^3 + 9*x^2 + 5*x - 2);
CtoD:=map<C->D|[C.1^2+C.3^2,C.2*(C.1+C.3),C.1*C.3]>;

// delta=1
f:=(x^6 + x^5 - 6*x^4 - 5*x^3 + 9*x^2 + 5*x - 2);
K<beta>:=NumberField(x^5 - x^4 - 4*x^3 + 3*x^2 + 3*x - 1);
OK:=IntegerRing(K);
KX<X>:=PolynomialRing(K);
Qpol:=(X+2)*(X-beta);
Rpol:=ExactQuotient(f,Qpol);
deltas:={d@@mp: d in s} where s,mp:=pSelmerGroup(2,Support(2*11*OK));

Kp,toKp:=Completion(K,2*OK);
KpX:=PolynomialRing(Kp);
toKpX:=hom<KX->KpX|toKp*Coercion(Kp,KpX),KpX.1>;
d1:=[d: d in deltas| IsSquare(toKp(Evaluate(d*Rpol,-2)))];
deltas:={K|d: d in deltas| d in d1 or HasRatXpoint([toKpX(d*Qpol),toKpX(d*Rpol)])};
deltas; //we find {1,beta^2-beta-2}

HC:=HyperellipticCurve(Rpol);
P1<s,t>:=ProjectiveSpace(Rationals(),1);
HCcov:=map<HC->P1|[HC.1,HC.3]>;
E,toE:=EllipticCurve(HC);
bl,G,mwmap:=PseudoMordellWeilGroup(E);
assert bl;
Ecov:=Expand(Inverse(toE)*HCcov);
N,V,R:=Chabauty(mwmap,Ecov,3);
assert N eq #V and PrimeDivisors(R) subset {2};
{EvaluateByPowerSeries(Ecov,mwmap(a)):a in V};

HC:=HyperellipticCurve((beta^2-beta-2)*Rpol);
HCcov:=map<HC->P1|[HC.1,HC.3]>;
E,toE:=EllipticCurve(HC);
bl,G,mwmap:=PseudoMordellWeilGroup(E);
assert bl and IsPSaturated(mwmap,3) and IsPSaturated(mwmap,5);
Ecov:=Expand(Inverse(toE)*HCcov);
N,V,R:=Chabauty(mwmap,Ecov,3);
assert N eq #V and PrimeDivisors(R) subset {2,3,5};
{EvaluateByPowerSeries(Ecov,mwmap(a)):a in V};

//delta=11

f:=11*(x^6 + x^5 - 6*x^4 - 5*x^3 + 9*x^2 + 5*x - 2);
K<beta>:=NumberField(x^5 - x^4 - 4*x^3 + 3*x^2 + 3*x - 1);
OK:=IntegerRing(K);
KX<X>:=PolynomialRing(K);
Qpol:=11*(X+2)*(X-beta);
Rpol:=ExactQuotient(f,Qpol);
deltas:={d@@mp: d in s} where s,mp:=pSelmerGroup(2,Support(2*11*OK));

Kp,toKp:=Completion(K,2*OK);
KpX:=PolynomialRing(Kp);
toKpX:=hom<KX->KpX|toKp*Coercion(Kp,KpX),KpX.1>;
d1:=[d: d in deltas| IsSquare(toKp(Evaluate(d*Rpol,-2)))];
deltas:={K|d: d in deltas| d in d1 or HasRatXpoint([toKpX(d*Qpol),toKpX(d*Rpol)])};

Kp,toKp:=Completion(K,(beta^2+beta-2)*OK);
KpX:=PolynomialRing(Kp);
toKpX:=hom<KX->KpX|toKp*Coercion(Kp,KpX),KpX.1>;
d1:=[d: d in deltas| IsSquare(toKp(Evaluate(d*Rpol,-2)))];
deltas:={K|d: d in deltas| d in d1 or HasRatXpoint([toKpX(d*Qpol),toKpX(d*Rpol)])};
deltas;
//{1,beta^3 - 2*beta^2 - beta + 1,beta^4 - beta^3 - 3*beta^2 + beta + 2}

//we have already done delta=1;

HC:=HyperellipticCurve((beta^3 - 2*beta^2 - beta + 1)*Rpol);
P1<s,t>:=ProjectiveSpace(Rationals(),1);
HCcov:=map<HC->P1|[HC.1,HC.3]>;
E,toE:=EllipticCurve(HC);
T,Tm:=TwoTorsionSubgroup(E);
[RankBound(E:Isogeny:=TwoIsogeny(Tm(t))):t in T | Order(t) ne 1];
Ecov:=Expand(Inverse(toE)*HCcov);
T,Tm:=TorsionSubgroup(E);
{EvaluateByPowerSeries(Ecov,Tm(t)):t in T};

HC:=HyperellipticCurve((beta^4 - beta^3 - 3*beta^2 + beta + 2)*Rpol);
P1<s,t>:=ProjectiveSpace(Rationals(),1);
HCcov:=map<HC->P1|[HC.1,HC.3]>;
E,toE:=EllipticCurve(HC);
bl,G,mwmap:=PseudoMordellWeilGroup(E:Isogeny:=2);
Ecov:=Expand(Inverse(toE)*HCcov);
N,V,R:=Chabauty(mwmap,Ecov,3);
assert N eq #V and PrimeDivisors(R) subset {2};
{EvaluateByPowerSeries(Ecov,mwmap(a)):a in V};

Uvals:={P1|[-2,1],[-1,1],[1,1],[2,1],[1,0]};
XtoU:=map<P1->P1|[s^2+t^2,s*t]>;
Xvals:=&join{RationalPoints(p@@XtoU):p in Uvals};
Xvals;

////////////////////////////////////////////////////////////////////////////
//Case: x^13+y^13=D*z^2

f:=x^13+1;
g:=ExactQuotient(f,x+1);
deltas:=&join[{QQ!(d@@mp):d in s} where s,mp:=pSelmerGroup(2,Support(13*D*OQ)):
               D in [2,3,5,6,10,11,13,17]];
deltas:={d:d in deltas | d gt 0};

deltas:={d:d in deltas|IsLocallySolvable(HyperellipticCurve(d*g),2)};
deltas:={d:d in deltas|IsLocallySolvable(HyperellipticCurve(d*g),3)};
deltas;
//we find delta in {1,13}


//delta=1
K<beta>:=NumberField(x^3 + x^2 - 4*x + 1);
OK:=IntegerRing(K);
KX<X>:=PolynomialRing(K);

Rpol:=X^4+beta*X^3+(beta^2+beta-1)*X^2+beta*X+1;
Qpol:=ExactQuotient(g,Rpol);

slm,mp:=pSelmerGroup(2,Support(2*13*OK));
slmQ,mpQ:=pSelmerGroup(2, Support(Norm(2*Resultant(Rpol,Qpol))*OQ));
H:=Kernel(hom<slm->slmQ|[mpQ(Norm(a@@mp)):a in OrderedGenerators(slm)]>);

deltas:={K|d@@mp: d in H};
Kp,toKp:=Completion(K,2*OK);
KpX:=PolynomialRing(Kp);
toKpX:=hom<KX->KpX|toKp*Coercion(Kp,KpX),KpX.1>;
deltas:={K|d: d in deltas| HasRatXpoint([toKpX(d*Qpol),toKpX(d*Rpol)])};
deltas;
//only delta=1 is left;

HC:=HyperellipticCurve(Rpol);
P1<s,t>:=ProjectiveSpace(Rationals(),1);
HCcov:=map<HC->P1|[HC.1,HC.3]>;
E,toE:=EllipticCurve(HC);
bl,G,mwmap:=PseudoMordellWeilGroup(E);
Ecov:=Expand(Inverse(toE)*HCcov);
N,V,R:=Chabauty(mwmap,Ecov,3);
assert N eq #V and PrimeDivisors(R) subset {2};
{EvaluateByPowerSeries(Ecov,mwmap(a)):a in V};

//delta=13
C:=HyperellipticCurve(13*g);
D:=HyperellipticCurve(13*(x^6 - x^5 - 5*x^4 + 4*x^3 + 6*x^2 - 3*x - 1));
CtoD:=map<C->D|[C.1^2+C.3^2,C.2,C.1*C.3]>;
J:=Jacobian(D);
assert TwoSelmerGroupData(J) eq 1;
gen:=elt<J|x^2 - 1/3*x - 1/3, 13/9*x - 26/9, 2>;
B:=CanonicalHeight(gen)/4+HeightConstant(J:Factor,Effort:=10);
V:=RationalPoints(J:Bound:=Ceiling(Exp(B)));
assert forall{v: v in V| IsZero(v) or CanonicalHeight(v) gt
                   CanonicalHeight(gen)/2};
Chabauty(gen,5); //only 1 U-coordinate
RationalPoints(D:Bound:=100);

Uvals:={P1|[-2,1]};
XtoU:=map<P1->P1|[s^2+t^2,s*t]>;
Xvals:=&join{RationalPoints(p@@XtoU):p in Uvals};
Xvals;

////////////////////////////////////////////////////////////////////////////
//Case: x^17+y^17=D*z^2

f:=x^17+1;
g:=ExactQuotient(f,x+1);
deltas:=&join[{QQ!(d@@mp):d in s} where s,mp:=pSelmerGroup(2,Support(17*D*OQ)):
               D in [2,3,5,6,10,11,13,17]];
deltas:={d:d in deltas | d gt 0};

deltas:={d:d in deltas|IsLocallySolvable(HyperellipticCurve(d*g),2)};
deltas:={d:d in deltas|IsLocallySolvable(HyperellipticCurve(d*g),3)};
deltas:={d:d in deltas|IsLocallySolvable(HyperellipticCurve(d*g),5)};
deltas:={d:d in deltas|IsLocallySolvable(HyperellipticCurve(d*g),11)};
deltas:={d:d in deltas|IsLocallySolvable(HyperellipticCurve(d*g),13)};
deltas;
//we find delta in {1,17}

K<beta>:=NumberField(x^4 + x^3  - 6*x^2 - x + 1);
OK:=IntegerRing(K);
KX<X>:=PolynomialRing(K);

Rpol:=X^4+beta*X^3+1/2*(-beta^3+6*beta+1)*X^2+beta*X+1;
Qpol:=ExactQuotient(g,Rpol);

slm,mp:=pSelmerGroup(2,Support(2*17*OK));
slmQ,mpQ:=pSelmerGroup(2, Support(Norm(2*Resultant(Rpol,Qpol))*OQ));
H:=Kernel(hom<slm->slmQ|[mpQ(Norm(a@@mp)):a in OrderedGenerators(slm)]>);
pi17:=1/2*(-3*beta^3 + 20*beta - 7);
assert Norm(pi17) eq 17;
s17:=mp(pi17);
deltas:={K|d@@mp: d in H} join{K|(d+s17)@@mp:d in H};

Kp,toKp:=Completion(K,(-beta^3 - beta^2 + 6*beta)*OK);
KpX:=PolynomialRing(Kp);
toKpX:=hom<KX->KpX|toKp*Coercion(Kp,KpX),KpX.1>;
deltas:={K|d: d in deltas| HasRatXpoint([toKpX(d*Qpol),toKpX(d*Rpol)])};

Kp,toKp:=Completion(K,1/2*(beta^3 - 6*beta + 5)*OK);
KpX:=PolynomialRing(Kp);
toKpX:=hom<KX->KpX|toKp*Coercion(Kp,KpX),KpX.1>;
deltas:={K|d: d in deltas| HasRatXpoint([toKpX(d*Qpol),toKpX(d*Rpol)])};
deltas;
//deltas = {1,beta^3 + 2*beta^2 - 3*beta + 1}

HC:=HyperellipticCurve(Rpol);
P1<s,t>:=ProjectiveSpace(Rationals(),1);
HCcov:=map<HC->P1|[HC.1,HC.3]>;
E,toE:=EllipticCurve(HC);
bl,G,mwmap:=PseudoMordellWeilGroup(E);
assert bl;
Ecov:=Expand(Inverse(toE)*HCcov);
N,V,R:=Chabauty(mwmap,Ecov,3);
assert N eq #V and PrimeDivisors(R) subset {2};
{EvaluateByPowerSeries(Ecov,mwmap(a)):a in V};

HC:=HyperellipticCurve((beta^3 + 2*beta^2 - 3*beta + 1)*Rpol);
P1<s,t>:=ProjectiveSpace(Rationals(),1);
HCcov:=map<HC->P1|[HC.1,HC.3]>;
E,toE:=EllipticCurve(HC,RationalPoints(HC,-1)[1]);
assert RankBound(E) eq 0;
T,Tm:=TorsionSubgroup(E);
Ecov:=Expand(Inverse(toE)*HCcov);
{EvaluateByPowerSeries(Ecov,Tm(t)):t in T};
