////////////////////////////////////////////////////////////////////////////////////
// This file handles the 4 genus 3 hyperelliptic curves
////////////////////////////////////////////////////////////////////////////////////

P<x> := PolynomialRing(Rationals());

eqns := [
x^7 + 4*x^6 - 7*x^5 - 8*x^4 + 7*x^3 + 4*x^2 - x,
x^8 - 4*x^7 - 12*x^6 + 28*x^5 + 38*x^4 - 28*x^3 - 12*x^2 + 4*x + 1
];
  
// X556. Each twist maps to a genus 2 curve whose Jacobian has rank 2

for d in {1,-1,2,-2} do
    Hd := QuadraticTwist(HyperellipticCurve((x - 1)*(x + 1)*(x^4 + 4*x^3 - 6*x^2 - 4*x + 1)),d);
    J := Jacobian(Hd);
    if RankBound(J) eq 1 then
      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 ]);       
      Chabauty(bas[1]);
    else
     Chabauty0(Jacobian(SimplifiedModel(Hd)));
    end if;
end for;    

// returns that the only rational points satisfy xyz eq 0
    A1<t> := AffineSpace(Rationals(),1);
    X556 := HyperellipticCurve(eqns[1]);

    _,y := IsSquare(Evaluate(eqns[1],0)); {@X556![0,y],X556![0,-y]@};    
    [X556![a[1],0] :  a in RationalPoints(Scheme(A1,Evaluate(eqns[1],t)))];
    PointsAtInfinity(X556); // one point
    // 4 points total, all known points



// X558. 

for d in {1,-1,2,-2} do
    Hd := QuadraticTwist(HyperellipticCurve((x^2 + 2*x - 1)*(x^4 - 4*x^3 - 6*x^2 + 4*x + 1)),d);
    IsLocallySolvable(Hd,2);
    J := Jacobian(Hd);
    if RankBound(J) eq 1 then
      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 ]);       
      Chabauty(bas[1]);
    else
     Chabauty0(Jacobian(SimplifiedModel(Hd)));
    end if;
end for;    
// returns that the only rational points satisfy xz eq 0 or x eq +-1
  
  
    A1<t> := AffineSpace(Rationals(),1);    
    X558 := HyperellipticCurve(eqns[2]);

    bool,y := IsSquare(Evaluate(eqns[2],0));   if bool then {@X558![0,y],X558![0,-y]@}; end if;
    bool,y := IsSquare(Evaluate(eqns[2],-1)); if bool then {@X558![-1,y],X558![-1,-y]@}; end if;
    bool,y := IsSquare(Evaluate(eqns[2],1));   if bool then {@X558![1,y],X558![1,-y]@}; end if;
    [X558![a[1],0] :  a in RationalPoints(Scheme(A1,Evaluate(eqns[2],t)))];
    PointsAtInfinity(X558); // one point
    // 8 points total, all known points
