// Example code from paper 3:
// _Some ternary Diophantine equations of signature (n, n, 2)_ by Nils Bruin
//
print "Examples from: Some ternary Diophantine equations of signature (n, n, 2)";
ei := GetEchoInput();
SetEchoInput(true);
// Section 2: Proof of Proposition 1.3
_<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) ];


// 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;


// 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);

for delta in [ s : s in S2E | Order(s) ne 1 ] do
  print "Looking at delta =", delta @@ toS2E;
  psi := tor(delta @@ toS2E);
  XX := Domain(psi);
  C, CtoXX := Quartic(XX);
  V := RationalPoints(C : Bound := 20);
  print "Points found on quartic:", V;
  W := { psi(CtoXX(v)) : v in V };
  print "points on the curve:", W;
  assert forall{ P : P in W | toS2E(mu(P)) eq delta };
end for;


// 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 };
SetEchoInput(ei);