Us = {U0, U1, U2, U3, U4, U5, U6, U7, U8, U9};
Mod7Weights={0,6,5,3,1 ,2,4,1,2,4}; 

(*Order 3 automorphism*)
Order3Aut={U1->U2,U2->U3,U3->U1,U4->U5,U5->U6,U6->U4,U7->U8,U8->U9,U9->U7};

(*Import file of equations for FPP*)
SetDirectory[NotebookDirectory[]];
FileNames[]

filename = "Dependency/EqsFPP.txt";
EqsFPP=ToExpression[Import[filename]]/.rr->(Sqrt[-7]);

ToRR[x_]:=Ret=Rationalize[Chop[N[Re[x],1900],10^(-300)]]+rr Rationalize[Chop[N[Im[x]/(Sqrt[7]),1900],10^(-300)]];

Computing Divisors in |4H + Torsion| and Quadratics vanishing on 4H + Torsion, for nontrivial torsion elements.
This parts needs results from the 5H calculations. Here’s the idea. 

For a torsion divisor T, we can realize 4H + T as follows:
12H = (4H + T) + (3H + T) + 5H
We have computed in 5H-Torsion.nb what the quadratics vanishing on sections of 5H are. This allows us to find (again) multiple points on a divisor in  |5H|.
To find sections on 4H + T, we can compute a bunch of random points on 3H + T and 5H, and then solve for quadratic monomials vanishing on it. 
To find quadratics vanishing on 4H + T, we can use the sections computed above and find points that don’t vanish on 3H + T and the quadratics of 5H.

We can use C_7 action to reduce the question to computing divisors in |4H+D| , |4H+D_1|, |4H+D+D_1|. In the former case we have C7:C3 action, but in the latter cases we only have C3
action.

Only one of the Zi’s is really needed. For our case, we will choose Z2.

Z2 = {(-5/32-rr/32)*U0^2+U1*U4+(U3*U6)/2+(1/8+rr/8)*U1*U7+(-3/16+rr/16)*U2*U8,(-1/8+(3*rr)/8)*U2*U5+(5/4+rr/4)*U3*U6+(-5/8-rr/8)*U1*U7+(-1/16+(3*rr)/16)*U2*U8+U3*U9,U4^2+(5/16+rr/16)*U0*U5+(1/4+rr/4)*U4*U7+(-3/32+rr/32)*U7^2+(3/32-rr/32)*U0*U8,U1*U3+(-1/8-rr/8)*U7^2,(5/16+rr/16)*U0*U5+U2*U6+(-3/4+rr/4)*U4*U7-U7^2/4+(5/32+rr/32)*U0*U8,(-3/8-rr/8)*U1*U2+2*U5^2+U0*U6+(1/16+(3*rr)/16)*U3*U7+U5*U8,-(U1*U2)/4+U3*U4+(3/8-rr/8)*U0*U6+(3/16+rr/16)*U3*U7,(-5/4+rr/4)*U3^2+(3/4+rr/4)*U2*U7+2*U5*U9+U8*U9,(1/8+rr/8)*U0*U1+(5/16-rr/16)*U3^2+U5*U6+(-1/16-(3*rr)/16)*U2*U7,-(U2*U3)+(-1/2-rr/2)*U1*U5+(5/16+rr/16)*U0*U7+U6*U9,(3/8-rr/8)*U2*U3+(1/4+rr/4)*U1*U5+U6^2-(U0*U7)/4+(1/8+rr/8)*U1*U8,(-3/8+rr/8)*U2^2+(5/32+rr/32)*U0*U3-U1*U6+(7/16+(3*rr)/16)*U5*U7+U4*U8+(U7*U8)/2+(1/8+rr/8)*U1*U9,(3/16-rr/16)*U2^2+(-1/64+(3*rr)/64)*U0*U3+U4*U5+(5/32+rr/32)*U5*U7+(-1/16-rr/16)*U1*U9,(-1/4-rr/4)*U3*U5+U4*U6+(1/4+rr/4)*U6*U7-(U3*U8)/4+(U7*U9)/8,U1^2+U3*U5+(-3/2-rr/2)*U6*U7+(1/4-rr/4)*U3*U8+(-3/8-rr/8)*U7*U9}/.rr->(Sqrt[-7]);

FirstCut=U0+(1/2) (1+I (Sqrt[7])) (U1+U2+U3);
SecondCut = U0+(-5+I (Sqrt[7])) (U1+U2+U3)+(4-4 I (Sqrt[7])) (U4+U5+U6)-4 (U7+U8+U9);

(* Computing points on D, D1, and D + D1 *)

PointsOn3HD ={};
Do[
EqsFPPPlusCut = Union[EqsFPP, {U2-1,U0,RandomSample[Join[-Range[100],Range[100]], 10].Us}];
current = DeleteDuplicates[NSolve[EqsFPPPlusCut==0, {U0, U1, U2, U3, U4, U5, U6, U7, U8, U9},500]];
PointsOn3HD= Join[PointsOn3HD, current];
, {n, 1, 3}];
{"Length and accuracy of points on |3H+D| ",Length[PointsOn3HD],Accuracy[PointsOn3HD]}


PointsOn3HD1 ={};
Do[
EqsFPPPlusCut = Union[EqsFPP, {U2-1,FirstCut,RandomSample[Join[-Range[100],Range[100]], 10].Us}];
current = DeleteDuplicates[NSolve[EqsFPPPlusCut==0, {U0, U1, U2, U3, U4, U5, U6, U7, U8, U9},500]];
PointsOn3HD1 = Join[PointsOn3HD1, current];
, {n, 1, 3}];
{"Length and accuracy of points on |3H+D_1| ",Length[PointsOn3HD1],Accuracy[PointsOn3HD1]}


PointsOn3HDD1 ={};
Do[
EqsFPPPlusCut = Union[EqsFPP, {U2-1,SecondCut,RandomSample[Join[-Range[100],Range[100]], 10].Us}];
current = DeleteDuplicates[NSolve[EqsFPPPlusCut==0, {U0, U1, U2, U3, U4, U5, U6, U7, U8, U9},500]];
PointsOn3HDD1 = Join[PointsOn3HDD1, current];
, {n, 1, 3}];
{"Length and accuracy of points on |3H+D+D_1| ",Length[PointsOn3HDD1],Accuracy[PointsOn3HDD1]}

Z2Points ={};
Do[
EqsFPPPlusCut = Union[EqsFPP,Z2, {U2-1,RandomSample[Join[-Range[100],Range[100]], 10].Us}];
current = DeleteDuplicates[NSolve[EqsFPPPlusCut==0, {U0, U1, U2, U3, U4, U5, U6, U7, U8, U9},500]];
Z2Points = Join[Z2Points, current];
, {n, 1, 3}];
{"Length and accuracy of points on |Z2| ",Length[Z2Points],Accuracy[Z2Points]}

The case for 4H + D. This does have a C7 splitting.

deg=2;

Powers={a1,a2,a3,a4,a5,a6,a7,a8,a9,a10};
AllPowers=Solve[Join[Table[Powers[[k]]≥0,{k,1,10}],{Sum[Powers[[k]],{k,1,10}]==deg}],Integers];

Solutions = Join[Z2Points, PointsOn3HD];

Equations = {};
Equationsrr = {};

For[mod7=0,mod7≤6,mod7++,
Print["mod7=",mod7];

CorrectPowers={};
Do[
If[Mod[Sum[Mod7Weights[[k]] Powers[[k]],{k,1,10}]/.AllPowers[[k]],7]==mod7,
AppendTo[CorrectPowers,AllPowers[[k]]];
];
,{k,1,Length[AllPowers]}];


Monomials=Table[Product[Us[[n]]^(Powers[[n]]/.CorrectPowers[[k]]),{n,1,10}],{k,1,Length[CorrectPowers]}];
(* Print[Monomials]; *)

TableOfValues=Table[Monomials[[k]]/.Solutions[[l]],{l,1,Length[Solutions]},{k,1,Length[Monomials]}];

(* 
Print["matrix done"];
Print[Dimensions[TableOfValues]];
Print[MatrixRank[TableOfValues]];
Print[Accuracy[TableOfValues]]; 
*)

result =NullSpace[Chop[RowReduce[TableOfValues,Tolerance->10^(-270)],10^(-270)]];
(* Print["nullspace done"];
Print[Dimensions[result]];
Print[Accuracy[result]]; *)

If[ Length[result]==0,Print["Nothing in this weight"]];
If[ Length[result]>0,
rowResult = RowReduce[result, Tolerance->10^(-250)];
(* Print["Rowreduce done"];
Print[Accuracy[rowResult]];  *)

exactResult = RootApproximant[Chop[rowResult, 10^(-270)]];

rowResultrr =Table[ToRR[exactResult[[i]]],{i,1,Length[exactResult]}];
newEquationsrr = rowResultrr.Monomials;
newEquations = exactResult.Monomials;
Print[newEquations];

AppendTo[Equations, newEquations];
AppendTo[Equationsrr, newEquationsrr];
];
];
sections = Flatten[{Equations[[1]], Equations[[2]], Equations[[3]]}];
sectionsrr = Flatten[{Equationsrr[[1]], Equationsrr[[2]], Equationsrr[[3]]}];

fname=OpenWrite["Equations/4HD.txt"];
Write[fname,sectionsrr];
Close[fname];

(* We want to find points on the sections above that do not vanish on Z2 and D. We only need to do this with the first section since the other 2 sets can be obtained by a C3 action.

To save time, we have precomputed it already. The code to compute the points are commented below. *)


PointsOn4HTor=ToExpression[Import["Dependency/points_on_4HD.txt"]];

(*
iter = 5;
(* The other curves Z2 and D *)
OtherCurves = Join[Z2, {U0}];
PointsOn4HTor= {};

Do [
Print["Iteration: ", i];

RandomCoeff = RandomSample[Join[-Range[100],Range[100]], 10];
RandomCut = RandomCoeff.Us;
Print[RandomCut];
	
InputEquations = Join[EqsFPP, {sections[[1]], U2 - 1,RandomCut}];
Print[Dimensions[InputEquations]];

current = NSolve[InputEquations==0, {U0, U1, U2, U3, U4, U5, U6, U7, U8, U9},500];
Print["Solve done"];


sortedPoints = Sort[current,Min[Map[Norm][OtherCurves/.#1]] > Min[Map[Norm][OtherCurves/.#2]]&];
PointsOn4HTor = Join[PointsOn4HTor, sortedPoints[[1;;10]]];

Print[Length[PointsOn4HTor]];

, {i, 1, iter}];
*)

Equations = {};
Equationsrr = {};

For[mod7=0,mod7≤6,mod7++,
Print["mod7=",mod7];

CorrectPowers={};
Do[
If[Mod[Sum[Mod7Weights[[k]] Powers[[k]],{k,1,10}]/.AllPowers[[k]],7]==mod7,
AppendTo[CorrectPowers,AllPowers[[k]]];
];
,{k,1,Length[AllPowers]}];

(* Print[Length[CorrectPowers]]; *)

Monomials=Table[Product[Us[[n]]^(Powers[[n]]/.CorrectPowers[[k]]),{n,1,10}],{k,1,Length[CorrectPowers]}];
TableOfValues=Table[Monomials[[k]]/.PointsOn4HTor[[l]],{l,1,Length[PointsOn4HTor]},{k,1,Length[Monomials]}];

(* Print["matrix done"];
Print[Dimensions[TableOfValues]];
Print[MatrixRank[TableOfValues]];
Print[Accuracy[TableOfValues]]; *)

result =NullSpace[Chop[RowReduce[TableOfValues,Tolerance->10^(-270)],10^(-270)]];
(* Print["nullspace done"];
Print[Dimensions[result]];
Print[Accuracy[result]];*)


rowResult = RowReduce[result, Tolerance-> 10^(-250)];
(* Print["Rowreduce done"];
Print[Accuracy[rowResult]]; *)

exactResult = RootApproximant[Chop[rowResult, 10^(-250)]];

rowResultrr =Table[ToRR[exactResult[[i]]],{i,1,Length[exactResult]}];
newEquationsrr = rowResultrr.Monomials;
newEquations = exactResult.Monomials;
Print[newEquations];

Equations=Join[Equations, newEquations];
Equationsrr=Join[Equationsrr, newEquationsrr];

];

fname=OpenWrite["Equations/4HD_one_section.txt"];
Write[fname,Equationsrr];
Close[fname];

(* We will also check if there base points. *)
Equations2=Equations/.Order3Aut;
Equations3=Equations2/.Order3Aut;
AllEquations = Union[Equations,Equations2,Equations3];
AllMonoms=DeleteDuplicates[Flatten[Table[Us[[i]]Us[[j]],{i,1,10},{j,1,10}]]];
mr=MatrixRank[Table[Coefficient[AllEquations,AllMonoms[[k]]],{k,1,Length[AllMonoms]}]]
(* answer is 54 *)
(* The only way to have a nontrivial intersection of codimension one space of quadrics is if there is only one intersection point. 
However, this contradicts the fact that the base locus is C7:C3 invariant and there are no points with this stabilizer. *)

The case of 4H + D1.

(* Generate all powers whose degree sum up to 2 *)
deg=2;

Powers={a1,a2,a3,a4,a5,a6,a7,a8,a9,a10};
AllPowers=Solve[Join[Table[Powers[[k]]≥0,{k,1,10}],{Sum[Powers[[k]],{k,1,10}]==deg}],Integers];

(* There is no C7 eigenspace splitting in this case *)
Monomials=Table[Product[Us[[n]]^(Powers[[n]]/.AllPowers[[k]]),{n,1,10}],{k,1,Length[AllPowers]}];
(* Print[Monomials]; *)

Solutions = Join[Z2Points, PointsOn3HD1];
TableOfValues=Table[Monomials[[k]]/.Solutions[[l]],{l,1,Length[Solutions]},{k,1,Length[Monomials]}];

(* Print["matrix done"];
Print[Dimensions[TableOfValues]];
Print[MatrixRank[TableOfValues]];
Print[Accuracy[TableOfValues]]; *)

result =NullSpace[Chop[RowReduce[TableOfValues,Tolerance->10^(-270)],10^(-270)]];
(* Print["nullspace done"];
Print[Dimensions[result]];
Print[Accuracy[result]]; *)

rowResult = RowReduce[Chop[result, 10^(-270)],Tolerance->  10^(-250)];
(* Print["Rowreduce done"];
Print[Accuracy[rowResult]]; *)

exactResult = RootApproximant[Chop[rowResult, 10^(-250)]];

rowResultrr =Table[ToRR[exactResult[[i]]],{i,1,Length[exactResult]}];
Equationsrr = rowResultrr.Monomials;
Equations = exactResult.Monomials;

Print[Equationsrr];

fname=OpenWrite["Equations/4HD1.txt"];
Write[fname,Equationsrr];
Close[fname];

(* We want to find points on the sections above that do not vanish on Z2 and D1. We will do this with the first section and hopefully get the other two with the group action. To save time, we have precomputed it already. The code to compute the points are commented below. *)

PointsOn4HTor=ToExpression[Import["Dependency/points_on_4HD1.txt"]];

(*

sections = Equations;

iter = 5;
OtherCurves = Join[Z2, {FirstCut}];
PointsOn4HTor= {};

Do [
Print["Iteration: ", i];

RandomCoeff = RandomSample[Join[-Range[100],Range[100]], 10];
RandomCut = RandomCoeff.Us;
Print[RandomCut];
	
InputEquations = Join[EqsFPP, {sections[[1]], U2 - 1,RandomCut}];
Print[Dimensions[InputEquations]];

current = NSolve[InputEquations==0, {U0, U1, U2, U3, U4, U5, U6, U7, U8, U9},500];
Print["Solve done"];


sortedPoints = Sort[current,Min[Map[Norm][OtherCurves/.#1]] > Min[Map[Norm][OtherCurves/.#2]]&];
PointsOn4HTor = Join[PointsOn4HTor, sortedPoints[[1;;10]]];



Print[Length[PointsOn4HTor]];

, {i, 1, iter}];

*)

Solutions = PointsOn4HTor;
TableOfValues=Table[Monomials[[k]]/.Solutions[[l]],{l,1,Length[Solutions]},{k,1,Length[Monomials]}];

(* Print[Dimensions[TableOfValues]];
Print[MatrixRank[TableOfValues]];
Print[Accuracy[TableOfValues]]; *)

result =NullSpace[Chop[RowReduce[TableOfValues,Tolerance->10^(-270)],10^(-270)]];
(* Print[Dimensions[result]];
Print[Accuracy[result]];*)

rowResult = RowReduce[result, Tolerance-> 10^(-300)];
(* Print["Rowreduce done"];
Print[Accuracy[rowResult]]; *)

exactResult = RootApproximant[Chop[rowResult, 10^(-300)]];
rowResultrr =Table[ToRR[exactResult[[i]]],{i,1,Length[exactResult]}];
ResultEquations = exactResult.Monomials;
ResultEquationsrr = rowResultrr.Monomials;

fname=OpenWrite["Equations/4HD1_one_section.txt"];
Write[fname,ResultEquationsrr];
Close[fname];

(* quick check of base point freeness gave 54, which was not enough to make a conclusion in this case *)

The case of 4H+D+D1

(* Generate all powers whose degree sum up to 2 *)
deg=2;

Powers={a1,a2,a3,a4,a5,a6,a7,a8,a9,a10};
AllPowers=Solve[Join[Table[Powers[[k]]≥0,{k,1,10}],{Sum[Powers[[k]],{k,1,10}]==deg}],Integers];

(* There is no C7 eigenspace splitting in this case *)
Monomials=Table[Product[Us[[n]]^(Powers[[n]]/.AllPowers[[k]]),{n,1,10}],{k,1,Length[AllPowers]}];
(* Print[Monomials]; *)

Solutions = Join[PointsOn3HDD1,Z2Points];
TableOfValues=Table[Monomials[[k]]/.Solutions[[l]],{l,1,Length[Solutions]},{k,1,Length[Monomials]}];


result =NullSpace[Chop[RowReduce[TableOfValues,Tolerance->10^(-270)],10^(-270)]];
(* Print["nullspace done"];
Print[Dimensions[result]];
Print[Accuracy[result]]; *)

rowResult = RowReduce[result,Tolerance->  10^(-250)];
(* Print["Rowreduce done"];
Print[Accuracy[rowResult]]; *)

exactResult = RootApproximant[Chop[rowResult, 10^(-250)]];

rowResultrr =Table[ToRR[exactResult[[i]]],{i,1,Length[exactResult]}];
Equationsrr = rowResultrr.Monomials;
Equations = exactResult.Monomials;

Print[Equationsrr];

fname=OpenWrite["Equations/4HDD1.txt"];
Write[fname,Equationsrr];
Close[fname];

(* We want to find points on the sections above that do not vanish on Z2 and D + D1. We only need to do this with the first section since the other 2 sets can be obtained by a C3 action.To save time, we have precomputed it already. The code to compute the points are commented below. *)
PointsOn4HTor=ToExpression[Import["Dependency/points_on_4HDD1.txt"]];
(*

sections = Equations;

iter = 5;
OtherCurves = Join[Z2, {SecondCut}];
PointsOn4HTor= {};

Do [
Print["Iteration: ", i];

RandomCoeff = RandomSample[Join[-Range[100],Range[100]], 10];
RandomCut = RandomCoeff.Us;
Print[RandomCut];
	
InputEquations = Join[EqsFPP, {sections[[1]], U2 - 1,RandomCut}];
Print[Dimensions[InputEquations]];

current = NSolve[InputEquations==0, {U0, U1, U2, U3, U4, U5, U6, U7, U8, U9},500];
Print["Solve done"];


sortedPoints = Sort[current,Min[Map[Norm][OtherCurves/.#1]] > Min[Map[Norm][OtherCurves/.#2]]&];
PointsOn4HTor = Join[PointsOn4HTor, sortedPoints[[1;;10]]];

Print[Length[PointsOn4HTor]];

, {i, 1, iter}];

*)

Solutions = PointsOn4HTor;
TableOfValues=Table[Monomials[[k]]/.Solutions[[l]],{l,1,Length[Solutions]},{k,1,Length[Monomials]}];

(* Print[Dimensions[TableOfValues]];
Print[MatrixRank[TableOfValues]];
Print[Accuracy[TableOfValues]]; *)

result =NullSpace[Chop[RowReduce[TableOfValues,Tolerance->10^(-270)],10^(-270)]];
(* Print[Dimensions[result]];
Print[Accuracy[result]];*)

rowResult = RowReduce[result, Tolerance-> 10^(-300)];
(* Print["Rowreduce done"];
Print[Accuracy[rowResult]]; *)

exactResult = RootApproximant[Chop[rowResult, 10^(-300)]];
rowResultrr =Table[ToRR[exactResult[[i]]],{i,1,Length[exactResult]}];
ResultEquations = exactResult.Monomials;
ResultEquationsrr = rowResultrr.Monomials;

fname=OpenWrite["Equations/4HDD1_one_section.txt"];
Write[fname,ResultEquationsrr];
Close[fname];

Checking that H^0(2H+torsion) = 0.

(* Import points on one section of 4H + D1, and one of 4H + D + D1 respectively. *)

PointsOn4HD1=ToExpression[Import["Dependency/points_on_4HD1.txt"]];
PointsOn4HDD1=ToExpression[Import["Dependency/points_on_4HDD1.txt"]];

N[Norm[Det[Us/.PointsOn4HD1[[16;;25]]]]]
N[Norm[Det[Us/.PointsOn4HDD1[[16;;25]]]]]

(* The fact these are nonzero implies that these points span P^9 and thus there are no sections of 6H which contain these sections of 4H+D1,4H+D+D1. *)