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

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

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

(* Changes a number in Q(sqrt(-7)) to the form a + b*rr, where rr stands for sqrt{-7} *)
ToRR[x_]:=Ret=Rationalize[Chop[N[Re[x],1900],10^(-1800)]]+rr Rationalize[Chop[N[Im[x]/(Sqrt[7]),1900],10^(-1800)]];

Part 1: Find the 6 global sections of H^0(FPP, 5H + D) 
The idea is that we can obtain 5H + D as 12H = (5H + D) + 4H + (3H + D), so we can find quadratic equations vanishing on 4H and 3H + D.

(* We first create points on r_3 = 0 and U0 = 0 *)
NewS3 =((((1-27*rr)*U1*U3)/29+(4*(101-rr)*U4^2)/29+(8*(15+rr)*U0*U5)/29+(8*(1+2*rr)*U2*U6)/29+8*U4*U7-4*U0*U8+((101-rr)*U0*U8)/29+((101-rr)*U2*U9)/58+(I/58)*Sqrt[7]*(101-rr)*U2*U9)/8)/.rr->(Sqrt[-7]);
NewD =( ((35-17*rr)*U0^2+2*(35-17*rr)*(U1*U4+U2*U5+U3*U6)+(21+13*rr)*(U1*U7+U2*U8+U3*U9))/812)/.rr->(Sqrt[-7]);
(* We compute 40 points on s3 = d = 0 and 20 points on U0=0. While we don't need this many points of s3 = d3 = 0 now, they will be useful later. *)
iter = 4;
PointsOn4H= {};
Do[
EqsFPPPlusCut  = Join[EqsFPP, {NewS3, NewD, U2 - 1,RandomSample[Join[-Range[100],Range[100]], 10].Us}];
current = NSolve[EqsFPPPlusCut ==0, 700];
PointsOn4H= Join[PointsOn4H, current[[1;;10]]];
Print[Length[PointsOn4H]];
, {i, 1, iter}];

Solutions=PointsOn4H;
Do[
EqsFPPPlusCut = Join[EqsFPP,{U0, U2-1, RandomSample[Join[-Range[100],Range[100]], 10].Us}];
current = NSolve[EqsFPPPlusCut==0, 700];
Solutions = Join[Solutions, current[[1;;10]]];
Print[Length[Solutions]];
, {i, 1, iter}];
(* Because accuracy can drop quite a bit here, we are making sure that it is still decent. *)
Print[Accuracy[Solutions]];

(* Generate all exponents of monomials of degree 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];
AllEquations = {};

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

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

(*Find the monomials that respect the C7 weighting*)
Monomials=Table[Product[Us[[n]]^(Powers[[n]]/.CorrectPowers[[k]]),{n,1,10}],{k,1,Length[CorrectPowers]}];
(* Find the space of monomials vanishing on random points of the cut  *)
TableOfValues=Table[Monomials[[k]]/.Solutions[[l]],{l,1,Length[Solutions]},{k,1,Length[Monomials]}];
(*Print[Dimensions[TableOfValues]];
Print[MatrixRank[TableOfValues]];*)

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

If[Length[result]>0,
rowResult = RowReduce[Chop[result, 10^(-270)]];
Print[{Length[rowResult],Accuracy[rowResult]}];

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

AppendTo[AllEquations, ResultEquations];
,
Print["No solutions."];
];
];
(* The weight 3 eigenspace has no quadratics, this is to be expected *)
AllEquations = Flatten[AllEquations];
(* Adjust weights - the presence of 4H here gives a -3 on the weights of every quadratic *)
SectionsOf5HD = Join[AllEquations[[4;;6]], AllEquations[[1;;3]]];
(* This ordering makes the weights in the order 1 2 3 4 5 6 *)
(* Print and write to the file. *)
Do [
Print[eq];
, {eq, SectionsOf5HD}];
fname=OpenWrite["Equations/5H+D_Sections.txt"];
Write[fname,Flatten[SectionsOf5HD]];
Close[fname];

Part 2: Produce a map of the FPP into CP^5 and compute equations of the image.

(* we will compute 70 points on the FPP and *)
Upoints={};
numpoints=70;
prec=800;
For[k=1,k≤Round[numpoints/10],k++,
Print[k];
EqsFPPplus= Join[EqsFPP, {U0-1, Sum[Us[[l]] RandomInteger[50],{l,1,10}], Sum[Us[[l]] RandomInteger[50],{l,1,10}]}];
sols=NSolve[EqsFPPplus==0,prec];
Upoints=Join[Upoints,Table[sols[[l]],{l,1,10}]];
];

(* Compute the images of these points in P^5. *)
Zs={Z1,Z2,Z3,Z4,Z5,Z6};
Zpoints=Factor[Table[Zs[[j]]->SectionsOf5HD[[j]]/.rr->(Sqrt[-7])/.Upoints[[i]],{i,1,Length[Upoints]},{j,1,Length[Zs]}]];
{Dimensions[Zpoints],Accuracy[Zpoints]}

(*We will now compute the degree six polynomials in Z which vanish at these points. Wesplit the calculation up by the C_7 weights. *)
deg=6;
Powers={b1,b2,b3,b4,b5,b6};
Zweights={1,2,3,4,5,6};
AllPowers=Solve[Join[Table[Powers[[k]]≥0,{k,1,6}],{Sum[Powers[[k]],{k,1,6}]==deg}],Integers];

C7Monomials = {};
EqsFPPInP5rr = {};

For[mod7=0,mod7≤6,mod7++,
Print["Weight ", mod7];
CorrectPowers={};
Do[
If[Mod[Sum[Zweights[[k]] Powers[[k]],{k,1,Length[Zweights]}]/.AllPowers[[k]],7]==mod7,
AppendTo[CorrectPowers,AllPowers[[k]]];
];
,{k,1,Length[AllPowers]}];

(*Find the monomials that respect the C7 weighting, we also store it in C7Monomials for future use below. *)
Monomials=Table[Product[Zs[[n]]^(Powers[[n]]/.CorrectPowers[[k]]),{n,1,Length[Powers]}],{k,1,Length[CorrectPowers]}];
AppendTo[C7Monomials, Monomials];

(* Find the space of monomials vanishing on random points of the cut  *)
TableOfValues=Table[Monomials[[k]]/.Zpoints[[l]],{l,1,Length[Zpoints]},{k,1,Length[Monomials]}];

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

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

AppendTo[EqsFPPInP5rr, ResultEquations];
	
];
ByteCount[EqsFPPInP5rr]

Now we perform some simplifications of the equations

(* This is a module to convert the equations of the FPP into nice linear combinations (in terms of the coefficient matrix) *)

NiceLinComd[CoeffMatrix_,root_,prec_]:=Module[{BackToCxReduced,CoeffMatrixReIm,i,j,k,u,r,v,SmallLinComb,rk,NewRow1,NewRow2,BackToCxMatrix},(*Building the rational matrix of double the size by looking at real and imaginary parts of rows and (I root) times rows.*)

CoeffMatrixReIm={};
For[i=1,i<=Length[CoeffMatrix],i++,NewRow1={};
NewRow2={};
For[j=1,j<=Length[CoeffMatrix[[i]]],j++,AppendTo[NewRow1,Re[CoeffMatrix[[i]][[j]]]];
AppendTo[NewRow1,Im[CoeffMatrix[[i]][[j]]]/root];
AppendTo[NewRow2,Im[CoeffMatrix[[i]][[j]]] (-root)];
AppendTo[NewRow2,Re[CoeffMatrix[[i]][[j]]]];];
AppendTo[CoeffMatrixReIm,NewRow1];
AppendTo[CoeffMatrixReIm,NewRow2];];

(*Getting rid of denominators.*)
CoeffMatrixReIm=Table[(Denominator[CoeffMatrixReIm[[k]]]/.
List->LCM) CoeffMatrixReIm[[k]],{k,1,Length[CoeffMatrixReIm]}];

(*Saturating the lattice using Smith decomposition.*)
Print["starting Smith decomposition"];
{u,r,v}=SmithDecomposition[CoeffMatrixReIm];
Print["Smith decomposition is done"];
SmallLinComb=Table[LinearSolve[Transpose[v],Table[KroneckerDelta[k,i],{k,1,Length[CoeffMatrixReIm[[1]]]}]],{i,1,MatrixRank[CoeffMatrixReIm]}];

(*Key step:using the lattice reduction algorithm.*)
Print["Starting lattice reduction."];
SmallLinComb=LatticeReduce[SmallLinComb];
Print["Lattice reduction is done."];

(*Getting back to the complex matrix.*)BackToCxMatrix=Table[Table[SmallLinComb[[i]][[2 k-1]]+I root SmallLinComb[[i]][[2 k]],{k,1,Length[SmallLinComb[[1]]]/2}],{i,1,Length[SmallLinComb]}];

(*Picking the first rows that give the correct rank.*)
BackToCxReduced={};
rk=0;
For[i=1,i<=Length[BackToCxMatrix],i++,If[MatrixRank[N[Append[BackToCxReduced,BackToCxMatrix[[i]]],prec]]>rk,rk++;
BackToCxReduced=Append[BackToCxReduced,BackToCxMatrix[[i]]];
(*Print[i," ",rk];*)
If[rk==Length[CoeffMatrix],i=Length[BackToCxMatrix]+1];
];
];

BackToCxReduced];

EQSFPP = EqsFPPInP5rr;
Zs={Z1,Z2,Z3,Z4,Z5,Z6};
SimpleEqsFPP = {};

For[mod7=0,mod7≤6,mod7++,
Print["Weight ", mod7];
EQS=RootReduce[EQSFPP[[mod7+1]]/.rr-> (-7)^(1/2)];
GradedMonoms = C7Monomials[[mod7 + 1]];
CoeffMatrix=Table[Coefficient[EQS[[k]],GradedMonoms[[i]]],{k,1,Length[EQS]},{i,1,Length[GradedMonoms]}];

(* Producing the coefficient matrix of the nice linear combination *)
BackToCxReduced = NiceLinComd[CoeffMatrix, Sqrt[7], 400];
(* Getting the equations in rr *)
coeffInrr = Table[ToRR[BackToCxReduced[[i]]],{i,1,Length[BackToCxReduced]}];
EQSrr=Table[Sum[coeffInrr[[i]][[k]] GradedMonoms[[k]],{k,1,Length[GradedMonoms]}],{i,1,Length[coeffInrr]}];
AppendTo[SimpleEqsFPP, EQSrr];
];

ByteCount[SimpleEqsFPP]
fname=OpenWrite["Equations/5H+D_Embedding.txt"];
Write[fname,Flatten[SimpleEqsFPP]];
Close[fname];

Part 3: Find the 6 global sections of H^0(FPP, 5H) 
Now we will repeat this procedure for 5H. It will turn out that the resulting map here is actually singular.

The plan for 5H is as follows, we can realize 5H as

30H  = 5H + (4H + D) + \sum_{i = 1}^7 (3H + D_i) 

Hence, we can solve for degree 5 monomials vanishing on 4H + D and all of 3H + D_i.

We can reduce some of the degree 5 monomials out by the relations on the FPP. We call the monomials left “useful monomials”. There are exactly 58 of them in each C7 weight.

(* To save time, we have precomputed the useful monomials already. 
*)
UsefulMonomials = {{U7^3*U8^2,U6^2*U8^3,U6^3*U7^2,U5*U7^3*U8,U5*U6^2*U8^2,U5^2*U7^3,U5^2*U6^2*U8,U5^3*U6^2,U4*U7^2*U8^2,U4*U6^3*U7,U4*U5*U7^2*U8,U4*U5^2*U7^2,U4^2*U7*U8^2,U4^2*U6^3,U4^2*U5*U7*U8,U4^2*U5^2*U7,U4^3*U8^2,U4^3*U5*U8,U4^3*U5^2,U3*U7^4,U3*U6^2*U7*U8,U3*U5*U6^2*U7,U3*U4*U7^3,U3*U4*U6^2*U8,U3*U4*U5*U6^2,U3*U4^2*U7^2,U3*U4^3*U7,U3*U4^4,U3^2*U6*U8^2,U3^2*U5*U6*U8,U3^2*U5^2*U6,U3^3*U6*U7,U3^3*U4*U6,U3^4*U8,U3^4*U5,U2*U6*U7*U8^2,U2*U6^4,U2*U5*U6*U7*U8,U2*U5^2*U6*U7,U2*U4*U6*U8^2,U2*U4*U5*U6*U8,U2*U4*U5^2*U6,U2*U3*U8^3,U2*U3*U6*U7^2,U2*U3*U5*U8^2,U2*U3*U5^2*U8,U2*U3*U5^3,U2*U3*U4*U6*U7,U2*U3*U4^2*U6,U2*U3^2*U7*U8,U2*U3^2*U5*U7,U2*U3^2*U4*U8,U2*U3^2*U4*U5,U2^2*U7^2*U8,U2^2*U4^2*U8,U1*U8^4,U8^3*U9^2,U4^2*U9^3},

{U7^2*U8^3,U6*U7^4,U6^3*U7*U8,U5*U7^2*U8^2,U5*U6^3*U7,U5^2*U7^2*U8,U5^3*U7^2,U4*U7*U8^3,U4*U6*U7^3,U4*U6^3*U8,U4*U5*U7*U8^2,U4*U5*U6^3,U4*U5^2*U7*U8,U4*U5^3*U7,U4^2*U8^3,U4^2*U6*U7^2,U4^2*U5*U8^2,U4^2*U5^2*U8,U4^2*U5^3,U4^3*U6*U7,U4^4*U6,U3*U7^3*U8,U3*U6^2*U8^2,U3*U5*U7^3,U3*U5*U6^2*U8,U3*U5^2*U6^2,U3*U4*U7^2*U8,U3*U4*U5*U7^2,U3*U4^2*U7*U8,U3*U4^2*U5*U7,U3*U4^3*U8,U3*U4^3*U5,U3^2*U6^2*U7,U3^2*U4*U6^2,U3^3*U6*U8,U3^3*U5*U6,U3^5,U2*U6*U8^3,U2*U6^2*U7^2,U2*U5*U6*U8^2,U2*U5^2*U6*U8,U2*U5^3*U6,U2*U4*U6^2*U7,U2*U4^2*U6^2,U2*U3*U6*U7*U8,U2*U3*U5*U6*U7,U2*U3*U4*U6*U8,U2*U3*U4*U5*U6,U2*U3^2*U8^2,U2*U3^2*U5*U8,U2*U3^2*U5^2,U2*U3^3*U7,U2*U3^3*U4,U2^2*U7*U8^2,U2^2*U4*U8^2,U4^4*U9,U2*U8^3*U9,U4*U8*U9^3},

{U7*U8^4,U6*U7^3*U8,U6^3*U8^2,U5*U7*U8^3,U5*U6*U7^3,U5*U6^3*U8,U5^2*U7*U8^2,U5^2*U6^3,U5^3*U7*U8,U5^4*U7,U4*U8^4,U4*U6*U7^2*U8,U4*U5*U8^3,U4*U5*U6*U7^2,U4*U5^2*U8^2,U4*U5^3*U8,U4*U5^4,U4^2*U6*U7*U8,U4^2*U5*U6*U7,U4^3*U6*U8,U4^3*U5*U6,U3*U7^2*U8^2,U3*U6^3*U7,U3*U5*U7^2*U8,U3*U5^2*U7^2,U3*U4*U7*U8^2,U3*U4*U6^3,U3*U4*U5*U7*U8,U3*U4*U5^2*U7,U3*U4^2*U8^2,U3*U4^2*U5*U8,U3*U4^2*U5^2,U3^2*U7^3,U3^2*U6^2*U8,U3^2*U5*U6^2,U3^2*U4*U7^2,U3^2*U4^2*U7,U3^2*U4^3,U3^4*U6,U2*U7^4,U2*U6^2*U7*U8,U2*U5*U6^2*U7,U2*U4*U7^3,U2*U4*U6^2*U8,U2*U4*U5*U6^2,U2*U4^2*U7^2,U2*U4^3*U7,U2*U4^4,U2*U3*U6*U8^2,U2*U3*U5*U6*U8,U2*U3*U5^2*U6,U2*U3^2*U6*U7,U2*U3^2*U4*U6,U2*U3^3*U8,U2*U3^3*U5,U2^2*U8^3,U4^3*U8*U9,U8^2*U9^3},

{U5^2*U8^3,U5^2*U6*U7^2,U5^3*U8^2,U5^4*U8,U5^5,U8^5,U6*U7^2*U8^2,U6^4*U7,U5*U8^4,U5*U6*U7^2*U8,U4*U6*U7*U8^2,U4*U6^4,U4*U5*U6*U7*U8,U4*U5^2*U6*U7,U4^2*U6*U8^2,U4^2*U5*U6*U8,U4^2*U5^2*U6,U3*U7*U8^3,U3*U6*U7^3,U3*U6^3*U8,U3*U5*U7*U8^2,U3*U5*U6^3,U3*U5^2*U7*U8,U3*U5^3*U7,U3*U4*U8^3,U3*U4*U6*U7^2,U3*U4*U5*U8^2,U3*U4*U5^2*U8,U3*U4*U5^3,U3*U4^2*U6*U7,U3*U4^3*U6,U3^2*U7^2*U8,U3^2*U5*U7^2,U3^2*U4*U7*U8,U3^2*U4*U5*U7,U3^2*U4^2*U8,U3^2*U4^2*U5,U3^3*U6^2,U2*U7^3*U8,U2*U6^2*U8^2,U2*U5*U7^3,U2*U5*U6^2*U8,U2*U5^2*U6^2,U2*U4*U7^2*U8,U2*U4*U5*U7^2,U2*U4^2*U7*U8,U2*U4^2*U5*U7,U2*U4^3*U8,U2*U4^3*U5,U2*U3*U6^2*U7,U2*U3*U4*U6^2,U2*U3^2*U6*U8,U2*U3^2*U5*U6,U2*U3^4,U2^2*U6*U7*U8,U2^3*U4^2,U4^2*U8^2*U9,U4*U9^4},

{U6*U7*U8^3,U6^2*U7^3,U6^4*U8,U5*U6*U7*U8^2,U5*U6^4,U5^2*U6*U7*U8,U5^3*U6*U7,U4*U6*U8^3,U4*U6^2*U7^2,U4*U5*U6*U8^2,U4*U5^2*U6*U8,U4*U5^3*U6,U4^2*U6^2*U7,U4^3*U6^2,U3*U8^4,U3*U6*U7^2*U8,U3*U5*U8^3,U3*U5*U6*U7^2,U3*U5^2*U8^2,U3*U5^3*U8,U3*U5^4,U3*U4*U6*U7*U8,U3*U4*U5*U6*U7,U3*U4^2*U6*U8,U3*U4^2*U5*U6,U3^2*U7*U8^2,U3^2*U6^3,U3^2*U5*U7*U8,U3^2*U5^2*U7,U3^2*U4*U8^2,U3^2*U4*U5*U8,U3^2*U4*U5^2,U3^3*U7^2,U3^3*U4*U7,U3^3*U4^2,U2*U7^2*U8^2,U2*U6^3*U7,U2*U5*U7^2*U8,U2*U5^2*U7^2,U2*U4*U7*U8^2,U2*U4*U6^3,U2*U4*U5*U7*U8,U2*U4*U5^2*U7,U2*U4^2*U8^2,U2*U4^2*U5*U8,U2*U4^2*U5^2,U2*U3*U7^3,U2*U3*U6^2*U8,U2*U3*U5*U6^2,U2*U3*U4*U7^2,U2*U3*U4^2*U7,U2*U3*U4^3,U2*U3^3*U6,U2^2*U6*U8^2,U2^3*U4*U8,U4*U8^3*U9,U4^3*U9^2,U8*U9^4},

{U7^5,U6*U8^4,U6^2*U7^2*U8,U5*U6*U8^3,U5*U6^2*U7^2,U5^2*U6*U8^2,U5^3*U6*U8,U5^4*U6,U4*U7^4,U4*U6^2*U7*U8,U4*U5*U6^2*U7,U4^2*U7^3,U4^2*U6^2*U8,U4^2*U5*U6^2,U4^3*U7^2,U4^4*U7,U4^5,U3*U6*U7*U8^2,U3*U6^4,U3*U5*U6*U7*U8,U3*U5^2*U6*U7,U3*U4*U6*U8^2,U3*U4*U5*U6*U8,U3*U4*U5^2*U6,U3^2*U8^3,U3^2*U6*U7^2,U3^2*U5*U8^2,U3^2*U5^2*U8,U3^2*U5^3,U3^2*U4*U6*U7,U3^2*U4^2*U6,U3^3*U7*U8,U3^3*U5*U7,U3^3*U4*U8,U3^3*U4*U5,U2*U7*U8^3,U2*U6*U7^3,U2*U6^3*U8,U2*U5*U7*U8^2,U2*U5*U6^3,U2*U5^2*U7*U8,U2*U5^3*U7,U2*U4*U8^3,U2*U4*U6*U7^2,U2*U4*U5*U8^2,U2*U4*U5^2*U8,U2*U4*U5^3,U2*U4^2*U6*U7,U2*U4^3*U6,U2*U3*U7^2*U8,U2*U3*U5*U7^2,U2*U3*U4*U7*U8,U2*U3*U4*U5*U7,U2*U3*U4^2*U8,U2*U3*U4^2*U5,U2^3*U8^2,U8^4*U9,U4^2*U8*U9^2},

{U7^4*U8,U6^2*U7*U8^2,U6^5,U5*U7^4,U5*U6^2*U7*U8,U5^2*U6^2*U7,U4*U7^3*U8,U4*U6^2*U8^2,U4*U5*U7^3,U4*U5*U6^2*U8,U4*U5^2*U6^2,U4^2*U7^2*U8,U4^2*U5*U7^2,U4^3*U7*U8,U4^3*U5*U7,U4^4*U8,U4^4*U5,U3*U6*U8^3,U3*U6^2*U7^2,U3*U5*U6*U8^2,U3*U5^2*U6*U8,U3*U5^3*U6,U3*U4*U6^2*U7,U3*U4^2*U6^2,U3^2*U6*U7*U8,U3^2*U5*U6*U7,U3^2*U4*U6*U8,U3^2*U4*U5*U6,U3^3*U8^2,U3^3*U5*U8,U3^3*U5^2,U3^4*U7,U3^4*U4,U2*U8^4,U2*U6*U7^2*U8,U2*U5*U8^3,U2*U5*U6*U7^2,U2*U5^2*U8^2,U2*U5^3*U8,U2*U5^4,U2*U4*U6*U7*U8,U2*U4*U5*U6*U7,U2*U4^2*U6*U8,U2*U4^2*U5*U6,U2*U3*U7*U8^2,U2*U3*U6^3,U2*U3*U5*U7*U8,U2*U3*U5^2*U7,U2*U3*U4*U8^2,U2*U3*U4*U5*U8,U2*U3*U4*U5^2,U2*U3^2*U7^2,U2*U3^2*U4*U7,U2*U3^2*U4^2,U2^2*U4^3,U1*U4*U8^3,U4*U8^2*U9^2,U9^5}};

(* Sanity check on the Monomials here *)
(* Rank should be 58 for all weights. *)
Table[MatrixRank[UsefulMonomials[[i]]/.Upoints],{i,1,7}]
Table[DeleteDuplicates[Table[Mod[Sum[Exponent[UsefulMonomials[[i]][[k]],Us[[l]]]Mod7Weights[[l]],{l,1,10}],7],{k,1,58}]],{i,1,7}]

(* We also need random points on 3H + D_1. Since we will be looking at one eigenspace at a time, we will not need to include their C_7 translates. *)
(* We compute these points and then points on a section r_3 of 4H *)
iter = 5;
PointsOnDi= {};
Do [
EqsFPPPlusCut  = Join[EqsFPP, { U0+(1/2) (1+I (Sqrt[7])) (U1+U2+U3), U2 - 1,RandomSample[Join[-Range[100],Range[100]], 10].Us}];
current = DeleteDuplicates[NSolve[EqsFPPPlusCut ==0, 500]];
PointsOnDi= Join[PointsOnDi, current[[1;;10]]];
Print[Length[PointsOnDi]];
, {i, 1, iter}];
Solutions = Join[PointsOn4H, PointsOnDi];
{Dimensions[Solutions],Accuracy[Solutions]}

(* Now we will compute the Z_i *)

ans= {};
Do[
Monomials = UsefulMonomials[[mod7+1]];
TableOfValues=Table[Monomials[[k]]/.Solutions[[l]],{l,1,Length[Solutions]},{k,1,Length[Monomials]}];
(*Print[Dimensions[TableOfValues]];
Print[MatrixRank[TableOfValues]];*)

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

If[Length[result]>0,
rowResult = RowReduce[Chop[result, 10^(-200)]];
Print[{Length[rowResult],Accuracy[rowResult]}];

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

AppendTo[ans, ResultEquations];
,
Print["No solutions."];
];
,{mod7,0,6}];

(* The weight 3 eigenspace has no quadratics, this is to be expected *)
(* Adjust weights - the presence of 4H here gives a -3 on the weights of every quadratic *)
SectionsOf5H ={ans[[4]][[1]],ans[[5]][[1]],ans[[6]][[1]],ans[[1]][[1]],ans[[2]][[1]],ans[[3]][[1]]};
(* This ordering makes the weights in the order 1 2 3 4 5 6 *)


(* The following block of code makes an appropriate scaling of the Zs. The idea is that it turns out to be possible to adjust Zs so that ratios of their values at a C3-fixed point are among {(1+(Sqrt[-7]))^(n/3)} for some n. The precise n were found by a little trial and error. The values of "close" must be defined over our quadratic field. 
 *)
SectionsOf5Hclean=SectionsOf5H;
SectionsOf5Hclean[[1]]=SectionsOf5H[[1]]/.rr->(-7)^(1/2);
fixedpoint =NSolve[Join[EqsFPP,{U0-1,U1-U2,U2-U3,U4-U5,U5-U6,U7-U8,U8-U9}]==0,500][[1]];
close=RootReduce[RootApproximant[(SectionsOf5H[[1]]/SectionsOf5H[[2]])/.rr->(-7)^(1/2)/.fixedpoint]/(1 + (-7)^(1/2))^(1/3)]
SectionsOf5Hclean[[2]] =RootReduce[Collect[ close (SectionsOf5H[[2]]/.rr->(-7)^(1/2)),Us]];

close=RootReduce[RootApproximant[(SectionsOf5H[[1]]/SectionsOf5H[[4]])/.rr->(-7)^(1/2)/.fixedpoint]/(1 + (-7)^(1/2))^(2/3)]
SectionsOf5Hclean[[4]] =RootReduce[Collect[ close SectionsOf5H[[4]]/.rr->(-7)^(1/2),Us]];

close=RootReduce[RootApproximant[(SectionsOf5H[[1]]/SectionsOf5H[[3]])/.rr->(-7)^(1/2)/.fixedpoint]/(1 + (-7)^(1/2))^(2/3)]
SectionsOf5Hclean[[3]] =RootReduce[Collect[ close SectionsOf5H[[3]]/.rr->(-7)^(1/2),Us]];

close=RootReduce[RootApproximant[(SectionsOf5H[[1]]/SectionsOf5H[[5]])/.rr->(-7)^(1/2)/.fixedpoint]/(1 + (-7)^(1/2))^(1/3)]
SectionsOf5Hclean[[5]] =RootReduce[Collect[ close SectionsOf5H[[5]]/.rr->(-7)^(1/2),Us]];

close=RootReduce[RootApproximant[(SectionsOf5H[[1]]/SectionsOf5H[[6]])/.rr->(-7)^(1/2)/.fixedpoint]]
SectionsOf5Hclean[[6]] =RootReduce[Collect[ close SectionsOf5H[[6]]/.rr->(-7)^(1/2),Us]];

(* These are ugly, so there is no particular reason to save them. *)

Part 4: Produce a map of the FPP to CP^5

Zs={Z1,Z2,Z3,Z4,Z5,Z6};
Zpoints=Factor[Table[Zs[[j]]->SectionsOf5Hclean[[j]]/.rr->(Sqrt[-7])/.Upoints[[i]],{i,1,Length[Upoints]},{j,1,Length[Zs]}]];
Dimensions[Zpoints]

deg =6;
AllPowers=Solve[Join[Table[Powers[[k]]≥0,{k,1,6}],{Sum[Powers[[k]],{k,1,6}]==deg}],Integers];

EqsFPPInP5rr = {};

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

Monomials={};
For[i=1,i≤ Length[AllPowers],i++,
If[ Mod[Sum[(Powers[[k]]/.AllPowers[[i]])Zweights[[k]],{k,1,6}],7]==mod7,
AppendTo[Monomials,Product[Zs[[k]]^(Powers[[k]]/.AllPowers[[i]]),{k,1,6}]];
];
];
(* Print[Length[Monomials]]; *)(* always 66 *)

(* Find the space of monomials vanishing on random points of the cut  *)
TableOfValues=Table[Monomials[[k]]/.Zpoints[[l]],{l,1,Length[Zpoints]},{k,1,Length[Monomials]}];

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

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

AppendTo[EqsFPPInP5rr, ResultEquations];
	
];

EqsFPPInP5rr=Flatten[EqsFPPInP5rr];

fname=OpenWrite["Equations/5H_Equations.txt"];
Write[fname,EqsFPPInP5rr];
Close[fname];

This surface is actually singular. We suspect that the only 3 singular points here are exactly (1:0:0:0:0:0), (0:1:0:0:0:0), (0:0:0:1:0:0). Here we will show that these are indeed singular points.

q1=Table[Zs[[i]]->KroneckerDelta[i-1],{i,1,6}];
q2=Table[Zs[[i]]->KroneckerDelta[i-2],{i,1,6}];
q3=Table[Zs[[i]]->KroneckerDelta[i-4],{i,1,6}];
(* The following code, together with computation that the ideal defined by 59 sextic equations is prime,  verifies that these three points lie on the surface and that their tangent spaces are  We see that their ranks are 1 - which indicates that these are singular points *)
JacobianFPPSingular=Table[D[EqsFPPInP5rr[[i]],Zs[[k]]],{i,1,Length[EqsFPPInP5rr]},{k,1,6}]/.rr->(Sqrt[-7]);
{DeleteDuplicates[EqsFPPInP5rr/.q1],MatrixRank[JacobianFPPSingular/.q1]}
{DeleteDuplicates[EqsFPPInP5rr/.q2],MatrixRank[JacobianFPPSingular/.q2]}
{DeleteDuplicates[EqsFPPInP5rr/.q2],MatrixRank[JacobianFPPSingular/.q2]}

Part 5: Finding Quadratics Vanishing on Sections of 5H
We will verify here that the 6 sections produced have no common zeroes, hence this is a well-defined map.
To do this, we want to find quadratics vanishing on each of the Zi's, hence we need random points on each sections of 5H.
To obtain the random points, we will generate random cuts on the degree 5 equations for each Zi. The points that do not evaluate to 0 on each of 3H + Di and 4H will be a point on each sections of 5H + D.

NewS3 =(1/464)*((2-54*rr)*U1*U3-8*(-101+rr)*U4^2+240*U0*U5+16*rr*U0*U5+16*U2*U6+32*rr*U2*U6+464*U4*U7-30*U0*U8-2*rr*U0*U8+101*U2*U9+101*I*Sqrt[7]*U2*U9-rr*U2*U9-I*Sqrt[7]*rr*U2*U9)/.rr->(Sqrt[-7]);
(* FirstCut=U0+(1/2) (1+I (Sqrt[7])) (U1+U2+U3); *)
ζ=Exp[2*π*I/7];
Firstcuts={};
For[i=0, i≤6, i++, AppendTo[Firstcuts,U0+(1/2) (1+I (Sqrt[7])) (ζ^(6i)*U1+ζ^(5i)*U2+ζ^(3i)U3)]];
EightCurves = Join[Firstcuts, {NewS3}];

(* Each Zi takes about 10 minutes, so about 20 minutes total. *)
(* To save time, we have pre-computed the points on Z1 and Z5 already. If you want, you could uncomment the code below to run it yourself. *)

Z1points = ToExpression[Import["Dependency\\Zi_for_5H\\z1_points.txt"]];
Accuracy[Z1points]
Z5points = ToExpression[Import["Dependency\\Zi_for_5H\\z5_points.txt"]];
Accuracy[Z5points]

(*
prec=500;
(* one iteration is enough *)
current = NSolve[ Join[EqsFPP, {SectionsOf5Hclean[[1]], U2 - 1,RandomSample[Join[-Range[100],Range[100]], 10].Us}]==0,Us,prec];currentsorted=Sort[Table[{Min[Table[Abs[EightCurves[[m]]/.current[[k]]],{m,1,8}]],current[[k]]},{k,1,Length[current]}]];
Z1points =Table[currentsorted[[k]][[2]],{k,Length[currentsorted]-9,Length[currentsorted]}];
Print["Z1 done"];
Accuracy[Z1points]

(* one iteration is enough *)
current = NSolve[ Join[EqsFPP, {SectionsOf5Hclean[[5]], U2 - 1,RandomSample[Join[-Range[100],Range[100]], 10].Us}]==0,Us,prec];currentsorted=Sort[Table[{Min[Table[Abs[EightCurves[[m]]/.current[[k]]],{m,1,8}]],current[[k]]},{k,1,Length[current]}]];
Z5points =Table[currentsorted[[k]][[2]],{k,Length[currentsorted]-9,Length[currentsorted]}];
Print["Z5 done"];
Accuracy[Z5points]
*)

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];

(* Finding quadratics vanishing on Z1 *)
EquationsZ1 = {};
Do [
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]}];
TableOfValues=Table[Monomials[[k]]/.Z1points[[l]],{l,1,Length[Z1points]},{k,1,Length[Monomials]}];

result = NullSpace[RowReduce[TableOfValues, Tolerance->10^(-Accuracy[TableOfValues]+30)]];
Print[{Dimensions[result],Accuracy[result]}];

rowResult = RowReduce[Chop[result,10^(-Accuracy[result]+30)],Tolerance->10^(-Accuracy[result]+30)];
exactResult = RootApproximant[Chop[rowResult, 10^(-Accuracy[rowResult]+30)],2];
exactResultrr = Table[ToRR[exactResult[[i]]],{i,1,Length[exactResult]}];
newEquations = exactResult.Monomials;
newEquationsrr = exactResultrr.Monomials;
(*
Print[newEquations];
Print[newEquationsrr];
*)
EquationsZ1 = Join[EquationsZ1, newEquationsrr];

, {mod7, 0, 6}];

EquationsZ1

(* Finding quadratics vanishing on Z1 *)
EquationsZ5 = {};
Do [
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]}];
TableOfValues=Table[Monomials[[k]]/.Z5points[[l]],{l,1,Length[Z1points]},{k,1,Length[Monomials]}];

result = NullSpace[RowReduce[TableOfValues, Tolerance->10^(-Accuracy[TableOfValues]+30)]];
Print[{Dimensions[result],Accuracy[result]}];

If[mod7==5, Debugresult=result];

rowResult = RowReduce[Chop[result,10^(-Accuracy[result]+30)],Tolerance->10^(-Accuracy[result]+30)];
exactResult = RootApproximant[Chop[rowResult, 10^(-Accuracy[rowResult]+30)],2];
exactResultrr = Table[ToRR[exactResult[[i]]],{i,1,Length[exactResult]}];
newEquations = exactResult.Monomials;
newEquationsrr = exactResultrr.Monomials;
(*
Print[newEquations];
Print[newEquationsrr];
*)
EquationsZ5 = Join[EquationsZ5, newEquationsrr];

, {mod7, 0, 6}];

EquationsZ5

(* We will now use this to compute the equations for the other {Z_i =0}. We use Magma to confirm that these are indeed curves of correct degree and that {Z_i = 0} have no common points. Another verification is below. It shows that all 90 quadrics together span the space of all quadrics. *)
ZQuadratics = {1, 2, 3, 4, 5, 6};
Order3Aut={U1->U2,U2->U3,U3->U1,U4->U5,U5->U6,U6->U4,U7->U8,U8->U9,U9->U7};

ZQuadratics[[1]] = EquationsZ1;
ZQuadratics[[5]] =EquationsZ5;

ZQuadratics[[2]] = ZQuadratics[[1]]/.Order3Aut;
ZQuadratics[[4]] = ZQuadratics[[2]]/.Order3Aut;

ZQuadratics[[3]] = ZQuadratics[[5]]/.Order3Aut;
ZQuadratics[[6]] = ZQuadratics[[3]]/.Order3Aut;

fname=OpenWrite["Equations/Quadratics_on_Zi.txt"];
Write[fname,ZQuadratics];
Close[fname];

(* Numerical verification. *)
Deg2Monomials=DeleteDuplicates[Flatten[Table[Us[[i]]Us[[j]],{i,1,10},{j,1,10}]]];
ZQuadraticsflat=Flatten[ZQuadratics]/.rr->(-7)^(1/2);
If[MatrixRank[Table[Coefficient[ZQuadraticsflat,Deg2Monomials[[k]]],{k,1,Length[Deg2Monomials]}]]== 55,
Print["No base points."]];