//////////////////////////////////////////////////////////////////
// This code determines the rational points on X520
//////////////////////////////////////////////////////////////////
// We need to handle the d = 1 twist of the genus 3 2-cover of X520

// X520 is y^2 =  (x^2 + 2)*(x^4 + 8*x^3 + 4*x^2 - 16*x + 4)
// C_d is  dy1^2 = x^2 +2 , dy2^2 =  (x^4 + 8*x^3 + 4*x^2 - 16*x + 4)
// D_{d1,d2} is  d1y1^2 = x^2 + 2, d1d2y2^2 = x^2 + (2*K.1 + 4)*x - 2, d2y3^2 = x^2 + (-2*K.1 + 4)*x - 2
//
// D maps to the genus 0 curves 
//   H1_{d1,d2}: d1^2d2 y1^2 = (x^2 + 2)*(x^2 + (2*K.1 + 4)*x - 2)
//   H2_{d1,d2}:    d1 *d2 y1^2 = (x^2 + 2)*(x^2 + (-2*K.1 + 4)*x - 2)


_<t> := PolynomialRing(Rationals());
  K :=SplittingField(t^2 - 2);
  Q<x> := PolynomialRing(K);
  f := x^4 + 8*x^3 + 4*x^2 - 16*x + 4; // Discriminant is 2^21
  OK := Integers(K);
  pp := Factorization(2*OK)[1][1];
 
sel2, msel2 := pSelmerGroup(2,{pp});
  msel2 := Inverse(msel2);

P1 := ProjectiveSpace(Rationals(),1);

possibleCoords:= {@@};

for dd in sel2 do 
  print "******************************";  
  d := msel2(dd);
  print "Working on the twist of the genus 5 cover by", d;

// D:
// y1^2 = (x^2 +  2)
// d*y2^2 = (x^2 + (2*K.1 + 4)*x - 2) 
// d*y3^2 = (x^2 + (-2*K.1 + 4)*x - 2)

  H := QuadraticTwist(HyperellipticCurve((x^2 + 2)*(x^2 + (-2*K.1 + 4)*x - 2))  , d);
    f := HyperellipticPolynomials(H);
  pts := RationalPoints(H : Bound := 100);
  if #pts eq 0 then
    print "No visible K-points.";
    print "IsLocallySolvable?";
    IsLocallySolvable(H,pp);
  else
  print "There are visible K-points. Lets try ellitpic Chabauty";
   e,m :=   EllipticCurve(H,pts[1]);
   minv := Extend(Inverse(m));
   success, G, GtoEK := PseudoMordellWeilGroup(e);

// elliptic chabauty map
  eqns := DefiningEquations(minv);
  xmap := Extend(map<e -> P1 | [eqns[1],eqns[3]]>);   

// Elliptic Chabauty

//  visiblePts := [<i,j> :   i in [0..1], j in [-10..10] |     xmap(GtoEK(i*G.1 + j*G.2))[1] in Rationals()];
  print "Does it work?";

  exists(p){ p : p in SequenceToSet(PrimesUpTo(30)) diff {2}  |  not p eq 2 and N eq #V   where   N, V, R, C is   Chabauty(GtoEK, xmap, p)};
  N, V, R, C :=   Chabauty(GtoEK, xmap, p); 
  posPts := {minv((GtoEK(a))) : a in V};
  print "The points on H with rational x cooordinate are",  posPts;

  possibleCoords:= possibleCoords join {@[Integers()!pt[1],Integers()!pt[3]] : pt in posPts@};
  
// works in every case that has local points.

  end if;

end for;
  
print "The possible pairs [x:z] are", possibleCoords; 
zeroCoords := {@pt : pt in possibleCoords | pt[2] eq 0@};
print "Of these, ", zeroCoords, "and", {pt : pt in possibleCoords diff zeroCoords | 
                    IsSquare( (x^2 + 2)*(x^4 + 8*x^3 + 4*x^2 - 16*x + 4)) where x is pt[1]/pt[2]},
                    "correspond to points on X520";

                    
//  The above code computes the lifts of all other points too. 
//  There are only the known points                    