//////////////////////////////////////////////////////////////////
// This code determines the rational points on X650
//////////////////////////////////////////////////////////////////

// X650 is y^4 = (x^2 - 2xz - z^2)(x^2 + z^2)


 gramMatrix := function(Q)
   a := Evaluate(Q,[1,0,0]);    b := Evaluate(Q,[0,1,0]);  c := Evaluate(Q,[0,0,1]); 
   d := (Evaluate(Q,[1,1,0]) - a - b)/2;   e := (Evaluate(Q,[1,0,1]) - a - c)/2;   f := (Evaluate(Q,[0,1,1]) - b - c)/2;
  return Matrix(3,[a,d,e,d,b,f,e,f,c]);
 end function;


   P2<x,y,z> := ProjectiveSpace(Rationals(),2);    
   Q := [(x^2 - 2*x*z - z^2), (y^2), (x^2 + z^2)];
   C := Curve(P2,Q[1]*Q[3]-Q[2]^2);
    _<t> := PolynomialRing(Rationals());
   H := [QuadraticTwist(HyperellipticCurve(f),d)
     where f is -Determinant(M[1] + 2*t*M[2] + t^2*M[3])
     where M is [gramMatrix(Q[i]) : i in [1..3]]
         : d in [1,-1,2,-2]];


   P4<u,v,w,r,s> := ProjectiveSpace(Rationals(),4);
   D := [Curve(P4,[Evaluate(Q[1],[u,v,w]) -d*r^2, Evaluate(Q[2],[u,v,w]) 
               -d*r*s, Evaluate(Q[3],[u,v,w]) -d*s^2]) : d in [1,-1,2,-2]];
   
// Check for 2-adic points 
 for d in D do  IsEmpty(d(pAdicField(2,20))); end for; 
    // only d = 1 has 2 adic points
    
         
         
// maps to Hyperelliptics
  D := D[1];         
  aut := AutomorphismGroup(D); 
  H1<X,Y,Z>, prDH1 := CurveQuotient(AutomorphismGroup(D,[aut.1*aut.2]));
 
 H2 := HyperellipticCurve(-t^5 + 2*t);

 prH12 :=  map<H1->H2 | [X^2,X*Y*Z,Z^2]>;
 // Do Chabauty on H2
    J := Jacobian(H2);
      RankBound(J); // returns 1
    heightconst := HeightConstant(J : Effort:=2, Factor);
      pts := Points(J : Bound := 1000);
      min := Minimum({Height(P) : P in pts | Order(P) eq 0}); 
      bool := exists(PJ){pt : pt in pts | Height(pt) eq min};
      LogarithmicBound := Height(PJ) + heightconst;  // Bound on h(Q)
      AbsoluteBound := Ceiling(Exp(LogarithmicBound));
      PtsUpToAbsBound := Points(J : Bound:=AbsoluteBound );
      bas, ht := ReducedBasis([ pt : pt in PtsUpToAbsBound ]);       
    ptsH2 :=    Chabauty(bas[1]);

// Compute preimages of rational points on H2
  ptsH1 :=  {@qt  : qt in RationalPoints(pt @@ prH12), pt in ptsH2@};
  ptsD :=   {@qt  : qt in RationalPoints(pt @@ prDH1), pt in ptsH1@}; 
  ptsD; // returns only the known rational points.   
  
  
  
