/* --------------------------------------------------------------------------------------
Please direct any comments to jmorrow4692@gmail.com
 -------------------------------------------------------------------------------------- */

//Prelim info

A2<u,v> := AffineSpace(Rationals(),2);
C := Curve(A2,  u^3*v^2 - u^3*v - v^3 + 6*v^2 - 3*v - 1);
PointSearch(C,1000);
assert (Genus(C) eq 4);
PC<U,V,Z> := ProjectiveClosure(C);
ptsPC:= PointSearch(PC,1000);
assert (IsHyperelliptic(PC) eq false);
aut := AutomorphismGroup(PC);
Quo := CurveQuotient(aut);
assert (Genus(Quo) eq 0);
assert (HasRationalPoint(Conic(Quo)) eq true);

//Looks like Jac_PC is simple over Q

bad := [3];

p:= 1;
while p le 20 do
  p := NextPrime(p); 
  p := p in bad select NextPrime(p) else p;
  "*******************************";
  p;  
  Cp<x,y> := Curve(Reduction(C,p));
  num<A> := Numerator(ZetaFunction(Cp));
  Factorization(num);
end while;

/*Looks like Jac_PC splits into two 2-dimensional abelian varieties over a quadratic extension (Time = 42.140) */

p:=1;
while p le 5 do 
  p := NextPrime(p); 
  p := p in bad select NextPrime(p) else p;
  "*******************************";
  p;  
  Cp<x,y> := Curve(Reduction(C,p^2));
  num<A> := Numerator(ZetaFunction(Cp));
  Factorization(num);
end while;


/* We find a number field that is unramified outside of 3, and study our curve over said number field */

P<x> := PolynomialRing(Rationals());
K := NumberField(x^2 - x + 1);
O:= MaximalOrder(K);
A1<UU,VV,ZZ> := ProjectiveSpace(K,2);
CC := Curve(A1 , Evaluate(DefiningPolynomials(PC)[1],[UU,VV,ZZ]));

/* Check to see CC maps to any lower genus curves (which it should) */

autK := AutomorphismGroup(CC);
QK0 := CurveQuotient(autK);
QK1 := CurveQuotient(AutomorphismGroup(CC,[autK.1]));
QK2 := CurveQuotient(AutomorphismGroup(CC,[autK.1*autK.1])); 
QK3 := CurveQuotient(AutomorphismGroup(CC,[autK.2])); 
QK4 := CurveQuotient(AutomorphismGroup(CC,[autK.2*autK.2])); 
QK5,phi5 := CurveQuotient(AutomorphismGroup(CC,[autK.1*autK.2])); 
QK6,phi6 := CurveQuotient(AutomorphismGroup(CC,[autK.1*autK.1*autK.2])); 
QK7 := CurveQuotient(AutomorphismGroup(CC,[autK.1*autK.2*autK.2])); 

assert (HasRationalPoint(ChangeRing(Conic(QK0),Rationals())) eq true);
assert (HasRationalPoint(ChangeRing(Conic(QK1),Rationals())) eq true);
assert (HasRationalPoint(ChangeRing(Conic(QK2),Rationals())) eq true);
assert (HasRationalPoint(ChangeRing(Conic(QK3),Rationals())) eq true);
assert (Genus(QK5) eq 2);
assert (Genus(QK6) eq 2) and (Genus(QK7) eq 2) and (QK6 eq QK7 eq true);

/* Great so we have found factors of the Jac_PC(K)! Let's see how these hyper elliptic curves QK5,QK6 behave */

J5 := Jacobian(QK5);
assert (RankBound(J5) eq 0);
J6 := Jacobian(QK6);
assert (RankBound(J6) eq 0);

//Even better! Now we want to compute what the torsion looks like

//Determining the probably order of Jac_QK5

torsOrders := {@@};
bool := false; p := 3;
while  p le 10 do
  p :=  NextPrime(p);     
  g:= ChangeRing(QK5,O);  	
  R,phi:=ResidueClassField(Decomposition(O,p)[1][1]);
  G<X,Y,Z> := BaseChange(g,phi); 
  boo, G<X,Y,Z>,_ := IsHyperelliptic(G);
  numH := Numerator(ZetaFunction(G));  
  torsOrders :=   torsOrders join  {@ Evaluate(numH,1) @};
  [p,GCD(torsOrders)];
  J := Jacobian(G);
  pts := RationalPoints(G);
 [Order(J!(Divisor(pts[i]) - Divisor(pts[1])))  : i in [2..#pts] ];
"****************************";
end while;

//Determining the probably order of Jac_QK6

torsOrders := {@@};
bool := false; p := 3;
while  p le 10 do
  p :=  NextPrime(p);     
  g:= ChangeRing(QK6,O);  	
  R,phi:=ResidueClassField(Decomposition(O,p)[1][1]);
  G<X,Y,Z> := BaseChange(g,phi); 
  boo, G<X,Y,Z>,_ := IsHyperelliptic(G);
  numH := Numerator(ZetaFunction(G));  
  torsOrders :=   torsOrders join  {@ Evaluate(numH,1) @};
  [p,GCD(torsOrders)];
  J := Jacobian(G);
  pts := RationalPoints(G);
 [Order(J!(Divisor(pts[i]) - Divisor(pts[1])))  : i in [2..#pts] ];
"****************************";
end while;


/* By looking mod 7, we deduce that Jac ~ Z/3 for each of the hyper elliptic curves QK5,QK6. We want to try and find a divisor class that has 3-torsion */

pts5 := PointsAtInfinity(QK5);
pts6 := PointsAtInfinity(QK6);

//Checking that difference of cusps is 3-torsion on QK6

assert (IsPrincipal(Divisor(pts6[1]) - Divisor(pts6[2])) eq false);
assert (IsPrincipal( 2*(Divisor(pts6[1]) - Divisor(pts6[2]))) eq false);
assert (IsPrincipal( 3*(Divisor(pts6[1]) - Divisor(pts6[2]))) eq true);


//Checking that the difference of cusps is 3-torsion on QK5

assert (IsPrincipal( (Divisor(pts5[1]) - Divisor(pts5[2]))) eq false);
assert (IsPrincipal( 2*(Divisor(pts5[1]) - Divisor(pts5[2]))) eq false);
assert (IsPrincipal( 3*(Divisor(pts5[1]) - Divisor(pts5[2]))) eq true);

/* Hence difference of the cusps for each hyper elliptic curve generates the torsion of Jac_PC(K).This tells us that the only points on C come from pullbacks of these two points.
Below, we check that the two points at infinity pull back to the known 2 singular points on CC and hence the 2 singular points on PC */

for P in pts5 do
  Points(P@@phi5);
end for;
for P in pts6 do
  Points(P@@phi6);
end for;

