SetAutoColumns(false);
SetColumns(100);
SetEchoInput(true);
Attach("routines.m");
_<x>:=PolynomialRing(Rationals());
Q:=NumberField(x-1:DoLinearExtension) where x:=PolynomialRing(Rationals()).1;
OQ:=IntegerRing(Q);
Qx:=PolynomialRing(Q);

//Lemma 3.1
f:=x*(x^3+8);
lcf:=LeadingCoefficient(f);
fmonic:=f/lcf;

AQ<theta>:=quo<Qx|f>;
A<theta>:=quo<Parent(x)|f>;

slm,mp:=pSelmerGroup(AQ,3,Support(2*3*OQ));
slmQ,mpQ:=pSelmerGroup(3,Support(2*3*OQ));
Nslmmp:=map<slm->slmQ|a:->mpQ(Norm(a@@mp))>;
H:=Kernel(hom<slm->slmQ|[Nslmmp(a):a in OrderedGenerators(slm)]>);

//Now H is the subgroup of A(3,{2,3}) of cube norm. It is generated by:
[a@@mp:a in OrderedGenerators(H)];

//This code constructs the curve C_delta:
_<y0,y1,y2,y3>:=ProjectiveSpace(Rationals(),3);
YA<y0A,y1A,y2A,y3A>:=PolynomialRing(A,4);
AY<thetaY>,swp:=SwapExtension(YA);
_<y0,y1,y2,y3>:=BaseRing(AY);
P3:=Proj(Parent(y0));
function C(delta)
  _,_,eq1,eq2:=Explode(Eltseq(
      swp(delta*(&+[YA.i*theta^(i-1):i in [1..4]])^3)));
  return Scheme(P3,[eq1,eq2]);
end function;

//select those such that C(delta)(Q3) is non-empty.
deltas:=[delta: delta in [A!Eltseq(a@@mp): a in H] |
           IsLocallySolvable(C(delta),3:AssumeNonsingular,AssumeIrreducible)];
P2<X,Y,Z>:=ProjectiveSpace(Rationals(),2);
m1:=map<A->Rationals()|delta:->Evaluate(delta,0)>;
assert Evaluate(f,m1(theta)) eq 0;
m2:=map<A->Rationals()|delta:->Evaluate(delta,-2)>;
assert Evaluate(f,m2(theta)) eq 0;
//These maps give the elliptic curves E1_delta and E2_delta
E1:=map<A->PowerStructure(CrvEll)|
  delta:->WeierstrassModel(SimplifiedModel(EllipticCurve(C,C![-2,0,1]))) where C:=
                         Curve(P2,d*Y^3-(X^3+8*Z^3)) where
                         d:=PowerFreePart(Norm(delta)/m1(delta),3)>;
E2:=map<A->PowerStructure(CrvEll)|
  delta:->WeierstrassModel(SimplifiedModel(EllipticCurve(C,C![0,0,1]))) where C:=
                         Curve(P2,d*Y^3-X*(X^2-2*X*Z+4*Z^2)) where
                         d:=PowerFreePart(Norm(delta)/m2(delta),3)>;
rkE1:=[r where _,r:=RankBounds(E1(delta)): delta in deltas];
rkE2:=[r where _,r:=RankBounds(E2(delta)): delta in deltas];

//Table 1
[<deltas[i],
     PowerFreePart(Norm(deltas[i])/m1(deltas[i]),3),
     PowerFreePart(Norm(deltas[i])/m2(deltas[i]),3),
     rkE1[i],rkE2[i]>: i in [1..#deltas]];
     
///////////////////////////////////////////////////////////////////////////////

//Lemma 3.3

//First the data in Table 2:
K<alpha>:=NumberField(x^4-2*x^3-2*x+1);
E:=[EllipticCurve([0,0,0,0,2*alpha^3 - 4*alpha^2 - 2*alpha - 6]),
 EllipticCurve([0,0,0,0,-2*alpha^3 + 4*alpha^2 + 2*alpha - 2]),
 EllipticCurve([0,0,0,0,2*alpha^3 - 6*alpha^2 + 2*alpha]),
 EllipticCurve([0,0,0,0,2*alpha^3 - 2*alpha^2 - 6*alpha]),
 EllipticCurve([0,0,0,0,-30*alpha^3 - 22*alpha^2 - 30*alpha]),
 EllipticCurve([0,0,0,0,14*alpha^3 - 6*alpha^2 - 74*alpha + 32])];
G:=[*[E[1]![2,alpha^3-2*alpha^2-alpha],E[1]![2*alpha^2+2*alpha+2,-3*alpha^3-4*alpha^2-5*alpha-2]],
[E[2]![2*alpha,alpha^3 - alpha],E[2]![1,alpha^2 - 2*alpha]],
[E[3]![2,alpha^2 - 3]],
[E[4]![2, 2*alpha^3 - 3*alpha^2 - 2*alpha - 1]],
[E[5]![2*alpha^2 + 2*alpha + 2,6*alpha^3 - alpha^2 + 2*alpha - 5]],
[E[6]![-2*alpha^3+6*alpha^2-4*alpha+2,6*alpha^3-13*alpha^2-4*alpha+5]]*];

MWmaps:=[map<AbelianGroup([0:j in [1..#G[i]]])->E[i]|
  c:->&+[Eltseq(c)[j]*G[i][j]:j in [1..#G[i]]]>:i in [1..#G]];
Es:=E;

delete E,G;

//beware, this might take a while (15 minutes or so?)
//it verifies (via a 2-descent) that the given generators indeed generate a
//group of odd finite index in the full Mordell-Weil group.

for MWmap in MWmaps do
  E:=Codomain(MWmap);
  two:=MultiplicationByMMap(E,2);
  mu:=IsogenyMu(two);
  slm,mp:=SelmerGroup(two);
  assert slm eq
     sub<slm|[mp(mu(MWmap(g))):g in OrderedGenerators(Domain(MWmap))]>;
end for;

//and also check that the groups are of index prime to 3.
assert forall{MWmap:MWmap in MWmaps|IsPSaturated(MWmap,3)};

//////////////////////////////////////////////////////////////////////////////
//Setting up the data structures to store the data involving the application
//of Lemma 2.2:

OK:=IntegerRing(K);

//The group K(3,S) of which the deltas are a member:
Adelt,Adeltmap:=pSelmerGroup(3,Support(6*OK));
//The covered P1
P1<s,t>:=ProjectiveSpace(Rationals(),1);
Qst:=CoordinateRing(P1);
Kst<alphaKst>:=quo<R|Evaluate(DefiningPolynomial(K),R.1)> where
                                               R:=PolynomialRing(Qst);
Pr:=ProjectiveSpace(Rationals(),Degree(K)-1);
PrK<alphaPr>:=quo<R|Evaluate(DefiningPolynomial(K),R.1)> where
                                  R:=PolynomialRing(CoordinateRing(Pr));
AssignNames(~Pr,["y" cat Sprint(i): i in [0..Degree(K)-1]]);
assert Basis(K) eq [alpha^i:i in [0..Degree(K)-1]];
P2K<S,Y,T>:=ProjectiveSpace(K,2);

Ecovs:=[PowerStructure(MapSch)|];
//////////////////////////////////////////////////////////////////////////////
//we'll be gathering the covers E_delta -> P^1 in the list Ecovs

//first, define f and theta and make sure (s-theta*t) divides f
f:=s^4 + 6*s^2*t^2 - 3*t^4;
theta:=alpha^2-2*alpha;
assert Evaluate(f,[theta,1]) eq 0;

g:=(s-Kst!Eltseq(theta)*t);
assert TotalDegree(f) eq TotalDegree(Norm(g));
D:=Rationals()!(f/Norm(g));
deltaset:=[K|a@@Adeltmap:a in Adelt|IsPower(Norm(a@@Adeltmap)*D,3)];

//this function constructs the cover (y0:...:y3):->(s:t)
//between C_delta -> P^1
function Ccov(delta)
  Gvec:=Eltseq(PrK!Eltseq(delta)*(&+[Pr.i*alphaPr^(i-1):i in [1..Degree(K)]])^3);
  stK:=PolynomialRing(K,2);
  gstK:=stK.1-theta*stK.2;
  bs:=[Vector(Eltseq(MonomialCoefficient(gstK,stK.i))):i in [1,2]];
  M:=Matrix(ExtendBasis(bs,Universe(bs)))^(-1);
  I:=Eltseq(Vector(Gvec)*Matrix(Parent(Pr.1),RowSequence(M)));
  C:=Scheme(Pr,[I[3],I[4]]);
  return map<C->P1|[I[1],I[2]]>;
end function;

//the values of delta in Table 3:
deltas:=[1,alpha^3-2*alpha^2-alpha-3,
 -3*alpha^3+5*alpha^2+2*alpha+7,7*alpha^3-11*alpha^2-5*alpha-16];

//Check that these deltas represent the classes from deltaset that are locally
//solvable at 3.
assert {Adeltmap(delta):delta in deltas} eq
   {Adeltmap(delta):delta in deltaset|
        IsLocallySolvable(Domain(Ccov(delta)),3:
                 AssumeIrreducible,AssumeNonsingular)};

//The values s/t(p0) in Table 3:
STp0:=[P1|[1,0],[0,1],[1,1],[-1,1]];
//check that these points actually lift to rational points on C_delta
assert forall{i:i in [1..4]|#RationalPoints(STp0[i]@@Ccov(deltas[i])) gt 0};

//The indices i in Table 3:
idxs:=[1,2,3,4];

//given delta and E_i, compute the cover E_i -> E_delta -> (s:t)
//it is an error if E_delta and E_i are not isomorphic
function Ecov(delta,E)
  cov:=Ccov(delta);
  Edelta:=Curve(P2K,(Y^3-delta*Parent(S)!(Evaluate(f,[S,T])/(S-theta*T))));
  STmap:=map<Edelta->Codomain(cov)|[S,T]>;
  P:=Edelta!Rep(RationalPoints(Flexes(Edelta)));
  Et,EdeltaToEt:=EllipticCurve(Edelta,P);
  bl,EtToE:=IsIsomorphic(Et,E);
  assert bl;
  EdeltaToE:=EdeltaToEt*EtToE;
  return Inverse(EdeltaToE)*STmap;
end function;

//Compute the covers. Note that this line also checks that the entries i in
//Table 3 are correct
Ecovs:=Ecovs cat [Expand(Ecov(deltas[i],Es[idxs[i]])):i in [1..4]];

//////////////////////////////////////////////////////////////////////////////
//Second case. Exactly the same as above.

f:=-s^4 + 6*s^2*t^2 + 3*t^4;
theta:=alpha^3-alpha^2-alpha-2;
assert Evaluate(f,[theta,1]) eq 0;
g:=(s-Kst!Eltseq(theta)*t);
assert TotalDegree(f) eq TotalDegree(Norm(g));
D:=Rationals()!(f/Norm(g));
deltaset:=[K|a@@Adeltmap:a in Adelt|IsPower(Norm(a@@Adeltmap)*D,3)];
function Ccov(delta)
  Gvec:=Eltseq(PrK!Eltseq(delta)*(&+[Pr.i*alphaPr^(i-1):i in [1..Degree(K)]])^3);
  stK:=PolynomialRing(K,2);
  gstK:=stK.1-theta*stK.2;

  bs:=[Vector(Eltseq(MonomialCoefficient(gstK,stK.i))):i in [1,2]];
  M:=Matrix(ExtendBasis(bs,Universe(bs)))^(-1);
  I:=Eltseq(Vector(Gvec)*Matrix(Parent(Pr.1),RowSequence(M)));
  C:=Scheme(Pr,[I[3],I[4]]);
  return map<C->P1|[I[1],I[2]]>;
end function;

deltas:=[1,-4*alpha^3+8*alpha^2+4*alpha+11,alpha^2+alpha+1,alpha^3+alpha^2+alpha];
assert {Adeltmap(delta):delta in deltas} eq
   {Adeltmap(delta):delta in deltaset|
        IsLocallySolvable(Domain(Ccov(delta)),3:
                 AssumeIrreducible,AssumeNonsingular)};
idxs:=[2,1,6,5];
STp0:=[P1|[1,0],[0,1],[1,1],[-1,1]];
assert forall{i:i in [1..4]|#RationalPoints(STp0[i]@@Ccov(deltas[i])) gt 0};

function Ecov(delta,E)
  cov:=Ccov(delta);
  Edelta:=Curve(P2K,(Y^3-delta*Parent(S)!(Evaluate(f,[S,T])/(S-theta*T))));
  STmap:=map<Edelta->Codomain(cov)|[S,T]>;
  P:=Edelta!Rep(RationalPoints(Flexes(Edelta)));
  Et,EdeltaToEt:=EllipticCurve(Edelta,P);
  bl,EtToE:=IsIsomorphic(Et,E);
  assert bl;
  EdeltaToE:=EdeltaToEt*EtToE;
  return Inverse(EdeltaToE)*STmap;
end function;

Ecovs:=Ecovs cat [Expand(Ecov(deltas[i],Es[idxs[i]])):i in [1..4]];

//////////////////////////////////////////////////////////////////////////
// Now we have gathered all covers E_i -> (s:t) to which Chabauty must be
// applied.

//this list gives which prime to use for the corresponding entry Ecov[i]
primes:=[31,11,31,31,11,11,31,31];

//We'll gather the points (s:t) in P^1(Q) covered by the E_delta(K) in this
//set:
Points:={};

for i in [1..8] do
  Ecov:=Ecovs[i];
  // select the appropriate MW-group for Ecov
  bl:=exists(mwmap){mwmap: mwmap in MWmaps| Codomain(mwmap) eq Domain(Ecov)};
  assert bl;
  
  //apply Chabauty
  V,R:=Chabauty(mwmap,Ecov,primes[i]);
  
  //The value V returned is a set of elements of Domain(mwmap) that correspond
  //to the points of E(K) that map to P^1(Q), provided the index of the group
  //Domain(mwmap) in the Mordell-Weil group E(K) is prime to R. Since we've
  //checked 2- and 3- saturation, we verify that that is sufficient:
  assert Seqset(PrimeDivisors(R)) subset {2,3};
  
  //and we map these points to (s:t)-coordinates and put them in Points.
  Points join:={Ecov(mwmap(P)):P in V};
end for;

//so these are the values covered, exactly as claimed in Lemma 3.3.
Points;
