(* Code takes several hours to run *)

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

filename = "Dependency/EqsFPP.txt";
EqsFPPwithrr=ToExpression[Import[filename]];

EqsFPP = EqsFPPwithrr/.rr-> (-7)^(1/2);

We elucidate the group relations in the Picard group of 3H + Torsion explicitly. First, we lift the two cuts in 6H that produce a nonreduced scheme over F_11 to the number field Q[sqrt(-7)]. These correspond to two elements in the said group. Then we use the C_7 action to produce the the cuts corresponding to the elements in their orbits. Finally, we obtain the group relations explicitly by determining the dimension of the space of quadratics (i.e. the dimension of the subspace of H^0(12H)) that vanishes on triples of these cuts. 
This following part calculates explicit representatives of 3H + Torsion as non-reduced cuts in 6H with coefficients in Q[sqrt(-7)].

(*The setup. cut1 and cut2 are the two cuts that give unreduced schemes, computed via magma*)

Us={U0,U1,U2,U3,U4,U5,U6,U7,U8,U9};
prime= 11;
rootvalue=rr->2;
cut1=U0+7(U1+U2+U3)+0(U4+U5+U6)+0(U7+U8+U9);
cut2=U0+8(U1+U2+U3)+7(U4+U5+U6)+7(U7+U8+U9);
cuts={cut1,cut2};

PointsMod11={
{{U0+3*U9,U1+3*U9,U2+5*U9,U3+U9,U4+9*U9,U5+8*U9,U6+10*U9,U7+3*U9,U8+5*U9},{U0+8*U9,U1+5*U9,U2+6*U9,U3+2*U9,U4+3*U9,U5+5*U9,U6,U7+7*U9,U8+U9},{U0+3*U9,U1+9*U9,U2+6*U9,U3+5*U9,U4,U5+8*U9,U6+6*U9,U7+U9,U8+4*U9},{U0+6*U9,U1+2*U9,U2+6*U9,U3+10*U9,U4+9*U9,U5+7*U9,U6+5*U9,U7+9*U9,U8+6*U9},{U0+3*U8,U1+2*U8,U2+10*U8,U3+8*U8,U4+5*U8,U5+10*U8,U6+U8,U7+3*U8,U9},{U0+10*U9,U1+2*U9,U2+7*U9,U3+10*U9,U4+U9,U5+4*U9,U6+8*U9,U7+2*U9,U8+4*U9},{U0+2*U9,U1+7*U9,U2+6*U9,U3+4*U9,U4+4*U9,U5,U6+9*U9,U7+3*U9,U8+8*U9},{U0+3*U9,U1+8*U9,U2+2*U9,U3+10*U9,U4+U9,U5+5*U9,U6+10*U9,U7,U8+3*U9}},
{{U0+8*U9,U1,U2,U3+4*U9,U4,U5+9*U9,U6+6*U9,U7+7*U9,U8},
{U0+2*U9,U1,U2+4*U9,U3+U9,U4+U9,U5+4*U9,U6+7*U9,U7,U8+5*U9},{U0+10*U9,U1+10*U9,U2+9*U9,U3+4*U9,U4+6*U9,U5+U9,U6+10*U9,U7+2*U9,U8+3*U9},{U0+2*U8,U1+4*U8,U2+U8,U3,U4+4*U8,U5+7*U8,U6+U8,U7+5*U8,U9},{U0+6*U9,U1+U9,U2+9*U9,U3+6*U9,U4+5*U9,U5+6*U9,U6+8*U9,U7+4*U9,U8+6*U9},{U0+4*U9,U1+6*U9,U2+4*U9,U3+8*U9,U4+4*U9,U5+9*U9,U6+7*U9,U7+4*U9,U8+3*U9}}};

(*Lifting the points and cuts to a larger ring.*)

Uns={Un0,Un1,Un2,Un3,Un4,Un5,Un6,Un7,Un8,Un9};
Coeffs=Table[ToExpression[StringJoin["c",ToString[m]]],{m,1,2000}];
maxdepth=20;
slsroot=Solve[x^2+7==0,Modulus->(prime^(maxdepth+1))];
Mod[x/.slsroot,11]

GoodCuts={};
For[numcut=1,numcut≤2,numcut++,

Print["cut number ",numcut];

(* building pairs of tangent vectors at each point *)
pointscut=PointsMod11[[numcut]];
ptssub={};
tanvecs={};
For[k=1,k≤Length[pointscut],k++,
pt=Solve[Join[Table[pointscut[[k]][[l]]==0,{l,1,9}],{U0==1}],Modulus->prime][[1]];
AppendTo[ptssub,pt];

tanvector=Table[Us[[k]]->( Us[[k]]/.pt)+ t Uns[[k]] ,{k,1,10}];
Eqs=Join[Coefficient[((EqsFPPwithrr/.rootvalue)/.U0->1)/.tanvector,t],{Un0}];
(* We set up two orthogonal tangent vectors. *)
sls1=Solve[Join[Table[0==Eqs[[m]],{m,1,Length[Eqs]}],{Un1==0,Un2==1}],Modulus->prime][[1]];
sls2=Solve[Join[Table[0==Eqs[[m]],{m,1,Length[Eqs]}],{Un1==1,Un2==0}],Modulus->prime][[1]];
AppendTo[tanvecs,{tanvector/.sls1,tanvector/.sls2}];
];

currentcut=cuts[[numcut]];
currentptssub=ptssub;
currenttanvecs=tanvecs;

For[depth=1,depth≤maxdepth,depth++,
If[Mod[depth,10]==5,Print[depth]];
rootvalue=rr-> Mod[x/.slsroot[[1]],prime^(depth+1)];
EqsFPPnewdepth=PolynomialMod[(EqsFPPwithrr/.U0->1)/.rootvalue,prime^(depth+1)];

count=1;
(* We set up the cuts, points and tangent vectors for p^(depth+1). *)
newcut=currentcut + prime^depth (a1 (U1+U2+U3)+a2 (U4+U5+U6)+a3(U7+U8+U9));
newptssub={};
newtanvecs={};
For[k=1,k≤Length[currentptssub],k++,
newpt=(Us/.currentptssub[[k]]) + prime^depth  
{0,0,Coeffs[[count]] , Coeffs[[count+1]] ,  Coeffs[[count+2]] ,  Coeffs[[count+3]] ,
 Coeffs[[count+4]] ,  Coeffs[[count+5]] ,  Coeffs[[count+6]] ,  Coeffs[[count+7]]  };
count+=8;
AppendTo[newptssub,Table[Us[[j]]->newpt[[j]],{j,1,10}]];
newvec1=newpt + Coefficient[Us/.currenttanvecs[[k]][[1]],t]t +  prime^depth  t 
{0,0,0,Coeffs[[count]], Coeffs[[count+1]] ,  Coeffs[[count+2]] ,  Coeffs[[count+3]] ,
 Coeffs[[count+4]] ,  Coeffs[[count+5]] ,  Coeffs[[count+6]] };
count+=7;
newvec2=newpt + Coefficient[Us/.currenttanvecs[[k]][[2]],t]t +  prime^depth  t 
{0,0,0,Coeffs[[count]], Coeffs[[count+1]] ,  Coeffs[[count+2]] ,  Coeffs[[count+3]] ,
 Coeffs[[count+4]] ,  Coeffs[[count+5]] ,  Coeffs[[count+6]] };
count+=7;
AppendTo[newtanvecs,{Table[Us[[j]]->newvec1[[j]],{j,1,10}],Table[Us[[j]]->newvec2[[j]],{j,1,10}]}];
];

(* Solve for the cuts, points and tangent vectors for p^(depth+1).*)
Eqs={};

For[l=1,l≤Length[tanvecs],l++,
For[k=1,k≤84,k++,
tbz=(EqsFPPnewdepth[[k]]/.newtanvecs[[l]][[1]])/prime^depth;
AppendTo[Eqs,PolynomialMod[Factor[tbz/.t->0],prime]];
AppendTo[Eqs,PolynomialMod[Factor[Coefficient[tbz,t]],prime]];
tbz=(EqsFPPnewdepth[[k]]/.newtanvecs[[l]][[2]])/prime^depth;
AppendTo[Eqs,PolynomialMod[Factor[Coefficient[tbz,t]],prime]];
];
tbz=(newcut/.newtanvecs[[l]][[1]])/prime^depth;
AppendTo[Eqs,PolynomialMod[Factor[tbz/.t->0],prime]];
AppendTo[Eqs,PolynomialMod[Factor[Coefficient[tbz,t]],prime]];
];

sls=Solve[Table[0==Eqs[[k]],{k,1,Length[Eqs]}],Modulus->prime][[1]];
(*repeat the steps for a higher power of p.*)
currentcut=newcut/.sls;
currentptssub=newptssub/.sls;
currenttanvecs=newtanvecs/.sls;

];
AppendTo[GoodCuts,currentcut];
];

(* We now find the minimal polynomial of the coefficients. *)

(*The first cut*)

Print["cut 1"];
r7=Mod[rr/.rootvalue,11^21];
(*Assuming that the minimal polynomial has degree 2, it has the form u_1 + u_2*t + u_3*t^2 + Sqrt[-7]*(u_4 + u_5*t + u_6*t^2) = 0. 
We further know that the coefficient in GoodCuts satisfies this equation mod 11^21. So we compute a nullspace basis of 
u_1 + u_2*t + u_3*t^2 + Sqrt[-7]*(u_4 + u_5*t + u_6*t^2) +u_7 * 11^21 = 0 where we let t be the coefficient in GoodCuts.*)
whatis=Coefficient[GoodCuts[[1]],U1];
thebasis={1,r7,whatis,11^21};
smalllincomb=LatticeReduce[NullSpace[{thebasis}]][[1]];
smallpoly1={1,(-7)^(1/2),t}.{smalllincomb[[1]],smalllincomb[[2]],smalllincomb[[3]]};

(* Solving for the coefficients of the first cut. *)
FirstCut=RootReduce[U0+(U1+U2+U3)(t/.Solve[smallpoly1==0][[1]])];
Print[FirstCut];

(*The second cut*)
Print["cut 2"];
whatis=Coefficient[GoodCuts[[2]],U1];
thebasis={1,r7,whatis,11^21};
smalllincomb=LatticeReduce[NullSpace[{thebasis}]][[1]];
smallpoly21={1,(-7)^(1/2),t}.{smalllincomb[[1]],smalllincomb[[2]],smalllincomb[[3]]};


whatis=Coefficient[GoodCuts[[2]],U4];
thebasis={1,r7,whatis,11^21};
smalllincomb=LatticeReduce[NullSpace[{thebasis}]][[1]];
smallpoly22={1,(-7)^(1/2),t}.{smalllincomb[[1]],smalllincomb[[2]],smalllincomb[[3]]};


whatis=Coefficient[GoodCuts[[2]],U7];
thebasis={1,r7,whatis,11^21};
smalllincomb=LatticeReduce[NullSpace[{thebasis}]][[1]];
smallpoly23={1,(-7)^(1/2),t}.{smalllincomb[[1]],smalllincomb[[2]],smalllincomb[[3]]};

(* Solving for the coefficients of the second cut. *)
SecondCut=RootReduce[U0+(U1+U2+U3)(t/.Solve[smallpoly21==0][[1]])
      +(U4+U5+U6)(t/.Solve[smallpoly22==0][[1]]) 
+(U7+U8+U9)(t/.Solve[smallpoly23==0][[1]]) ];
Print[SecondCut];

(* Check that they produce nonreduced curves. *)
sols=NSolve[Join[Table[0==EqsFPP[[k]],{k,1,84}],{0==FirstCut,U0==1,U1==RandomInteger[10]+1}],50];
Length[DeleteDuplicates[sols]]
sols=NSolve[Join[Table[0==EqsFPP[[k]],{k,1,84}],{0==SecondCut,U0==1,U1==RandomInteger[10]+1}],50];
Length[DeleteDuplicates[sols]]

In this part we investigate the group relations of the torsion divisors. First we compute points on elements of |3H+T| where T is a nontrivial torsion divisor.

(* Obtain the orbit under the action of C_7. The weights are (0, 6, 5, 3, 1, 2, 4, 1, 2, 4).*)
FirstCut=U0+(1/2) (1+I (Sqrt[7])) (U1+U2+U3);
Firstcuts=Table[ToExpression[StringJoin["Firstcut",ToString[i]]],{i,1,7}];
Firstsolns=Table[ToExpression[StringJoin["Firstsoln",ToString[i]]],{i,1,7}];
ζ=Exp[2*π*I/7];
For[i=1, i<=7, i++, Firstcuts[[i]]=U0+(1/2) (1+I (Sqrt[7])) (ζ^(6i)*U1+ζ^(5i)*U2+ζ^(3i)U3)];
SecondCut = U0+(-5+I (Sqrt[7])) (U1+U2+U3)+(4-4 I (Sqrt[7])) (U4+U5+U6)-4 (U7+U8+U9);
Secondcuts=Table[ToExpression[StringJoin["Secondcut",ToString[i]]],{i,1,7}];
Secondsolns=Table[ToExpression[StringJoin["Secondsoln",ToString[i]]],{i,1,7}];
For[i=1, i<=7, i++, Secondcuts[[i]]=U0+(-5+I (Sqrt[7])) (ζ^(6i)*U1+ζ^(5i)*U2+ζ^(3i)*U3)+(4-4 I (Sqrt[7])) (ζ^(i)*U4+ζ^(2i)*U5+ζ^(4i)*(U6))-4 (ζ^(i)*U7+ζ^(2i)*U8+ζ^(4i)*U9)];

(*calculating some points on this cut*)
For[i=1,i<=7,i++,
Firstsolns[[i]]=Table[DeleteDuplicates[NSolve[Join[Table[0==EqsFPP[[k]],{k,1,84}],{0==Firstcuts[[i]],U0==1,RandomInteger[10] U0+RandomInteger[10] U1+RandomInteger[10] U2+RandomInteger[10] U3+RandomInteger[10] U4+RandomInteger[10] U5+RandomInteger[10] U6+RandomInteger[10] U7+RandomInteger[10] U8+RandomInteger[10] U9==1}],50]], {j, 1, 4}]];
For[i=1, i<=7, i++, Firstsolns[[i]]=Join[Firstsolns[[i]][[1]],Firstsolns[[i]][[2]],Firstsolns[[i]][[3]],Firstsolns[[i]][[4]]]];
Print["C7 translates of the first cut are done."];
(*Each time it produces at most 18 distinct solutions after removing duplicates so at most 72 points for each orbit which contains 7 cuts.*)

(*Do the same for the second cut*)
For[i=1,i<=7,i++,
Secondsolns[[i]]=Table[DeleteDuplicates[NSolve[Join[Table[0==EqsFPP[[k]],{k,1,84}],{0==Secondcuts[[i]],U0==1,RandomInteger[10] U0+RandomInteger[10] U1+RandomInteger[10] U2+RandomInteger[10] U3+RandomInteger[10] U4+RandomInteger[10] U5+RandomInteger[10] U6+RandomInteger[10] U7+RandomInteger[10] U8+RandomInteger[10] U9==1}],50]], {j, 1, 4}]];
For[i=1, i<=7, i++, Secondsolns[[i]]=Join[Secondsolns[[i]][[1]],Secondsolns[[i]][[2]],Secondsolns[[i]][[3]],Secondsolns[[i]][[4]]]];
Print["C7 translates of the second cut are done."];

(*Then we obtain some points for the U0=0 cut.*)
SolnU0=Table[DeleteDuplicates[NSolve[Join[Table[0==EqsFPP[[k]],{k,1,84}],{0==U0,U1==1,RandomInteger[10] U0+RandomInteger[10] U1+RandomInteger[10] U2+RandomInteger[10] U3+RandomInteger[10] U4+RandomInteger[10] U5+RandomInteger[10] U6+RandomInteger[10] U7+RandomInteger[10] U8+RandomInteger[10] U9==1}],50]],{j,1,4}];
SolnU0 =Join[SolnU0[[1]],SolnU0[[2]],SolnU0[[3]],SolnU0[[4]]];

(*Put all cuts together.*)
Solns=Append[Join[Firstsolns, Secondsolns], SolnU0];

Here we calculate the dimension of the space of H^0(12H) that vanishes on three different cuts. The zero of a section in this space is a quadratic with unknown coefficient, which if evaluated at the points on the cuts should be 0. So we substitute in the numerically calculated points to get matrices and compute the their rank. We record the triple when the rank is 55 since it means that the triple sums up to the trivial bundle. This produces the following explicit addition table. We denote by D the cut U0=0, D1 through D7 the orbit of the first cut and D8 through D14 that of the second cut.
In the above calculations 7 represents the first cut and 14 represents the second cut, but we change it here with 7 going to D1, 1 going to D2, etc and 14 going to D8,  8 going to D9, etc.

Triples={};
MatRows=DeleteDuplicates[Flatten[Table[Us[[i]]Us[[j]],{i,1,10},{j,1,10}]]];
For[i=1, i<=Length[Solns]-2, i++,
Print[i];
For[j=i+1, j<= Length[Solns]-1, j++,
For[k = j+1, k<= Length[Solns], k++,
CMat=Join[Table[MatRows/.Solns[[i]][[l]],{l,1,Length[Solns[[i]]]}], Table[MatRows/.Solns[[j]][[l]],{l,1,Length[Solns[[j]]]}],  Table[MatRows/.Solns[[k]][[l]],{l,1,Length[Solns[[k]]]}]];
RMat=MatrixRank[CMat];
If[RMat==  55,AppendTo[Triples,{i,j,k}]];
]]];

(* Now we compute the table. *)
For[i=1, i<= Length[Triples],i++,
For[j=1, j<= Length[Triples[[i]]], j++,
If[Triples[[i]][[j]]!=15 ,
If[Triples[[i]][[j]]==7,
Triples[[i]][[j]]=1,
If[Triples[[i]][[j]]==14,
Triples[[i]][[j]]=8,
Triples[[i]][[j]]=Triples[[i]][[j]]+1;]
]
]
]
]
Mult=Table["+", 16, 16];
For[i=1, i<=Length[Solns],i++,
For[j =1, j<=Length[Solns],j++,
For[k=1, k<=Length[Triples], k++,
If[i ==j, 
Mult[[i+1]][[j+1]]=0,
If[MemberQ[Triples[[k]],i]&& MemberQ[Triples[[k]],j],
For[m=1, m<=3,m++,
If[Triples[[k]][[m]]!=i && Triples[[k]][[m]]!=j,
Mult[[i+1]][[j+1]]=If[Triples[[k]][[m]] !=15, 
StringJoin["D",ToString[Triples[[k]][[m]]]], "D"]]]]
]]]]
For[i=1, i<=14, i++,
Mult[[1]][[i+1]]=StringJoin["D",ToString[i]]; 
Mult[[i+1]][[1]]=StringJoin["D",ToString[i]];]
Mult[[1]][[16]]="D";
Mult[[16]][[1]]="D";
Grid[Mult]

We henceforth denote D_8 through D_14 as D_1+D through D_7+D.

In this next part, we will compute for explicit quadratics vanishing on each unique section of 3H + Torsion.

Us = {U0, U1, U2, U3, U4, U5, U6, U7, U8, U9};
(* C7 weights of U0, ..., 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};

(* These are the cuts obtained from the file 3H_Torsion_Cuts.nb *)

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

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

Computing 28 Quadratics on 3H + D
This has a C7 eigenspace splitting

Cut = {InvariantCut, U2 - 1};
Solutions = {};
iter = 3;
Do[
Print["Iteration - ", n];
EqsFPPPlusCut = Union[EqsFPP, Cut,{RandomSample[Join[-Range[100],Range[100]], 10].Us}];
current = NSolve[EqsFPPPlusCut==0, {U0, U1, U2, U3, U4, U5, U6, U7, U8, U9},500];
(*Print[current];*)
Solutions = Join[Solutions, current];
, {n, 1, iter}];

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

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]}];
(* Print[Monomials];*)

(* 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]]; *)

rowResult = RowReduce[Chop[result, 10^(-300)]];
(* Print[Accuracy[rowResult]]; *)

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

AppendTo[AllEquations, ResultEquations];
(* Print["-----------------------------------------"];*)
	
];

(* Here we are adjusting the quadratics to be generated by the C3 action *)
AllEquationsC3 = Range[1, 7];
AllEquationsC3[[0+1]] = AllEquations[[0+1]];

(* Equations in weight 2, 4 can be obtained as C3 translates of weight 1 *)
AllEquationsC3[[1+1]] = AllEquations[[1+1]];
AllEquationsC3[[2+1]] = AllEquationsC3[[1+1]]/.Order3Aut;
AllEquationsC3[[4+1]] = AllEquationsC3[[2+1]]/.Order3Aut;

(* Equations in weight 6, 5 can be obtained as C3 translates of weight 3 *)
AllEquationsC3[[3+1]] = AllEquations[[3+1]];
AllEquationsC3[[6+1]] = AllEquationsC3[[3+1]]/.Order3Aut;
AllEquationsC3[[5+1]] = AllEquationsC3[[6+1]]/.Order3Aut;

Do [
Print[eq];
, {eq, AllEquationsC3}];

fname=OpenWrite["Equations/3H_D.txt"];
Write[fname, Flatten[AllEquationsC3]];
Close[fname];

Length[Flatten[AllEquationsC3]]

Computing 28 Quadratics on 3H + D1

Cut = {FirstCut, U2 - 1};
Solutions = {};
iter = 3;
Do[
Print["Iteration - ", n];
EqsFPPPlusCut = Union[EqsFPP, Cut,{RandomSample[Join[-Range[100],Range[100]], 10].Us}];
current = NSolve[EqsFPPPlusCut==0, {U0, U1, U2, U3, U4, U5, U6, U7, U8, U9},500];
(*Print[current];*)
Solutions = Join[Solutions, current];
, {n, 1, iter}];
(* *)

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

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

rowResult = RowReduce[Chop[result, 10^(-300)]];
Print[Accuracy[rowResult]];

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

fname=OpenWrite["Equations/3H_D1.txt"];
Write[fname, ResultEquations];
Close[fname];

Length[ResultEquations]

Computing 28 Quadratics on 3H + D + D1

Cut = {SecondCut, U2 - 1};
Solutions = {};
iter = 3;
Do[
Print["Iteration - ", n];
EqsFPPPlusCut = Union[EqsFPP, Cut,{RandomSample[Join[-Range[100],Range[100]], 10].Us}];
current = NSolve[EqsFPPPlusCut==0, {U0, U1, U2, U3, U4, U5, U6, U7, U8, U9},500];
(*Print[current];*)
Solutions = Join[Solutions, current];
, {n, 1, iter}];

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

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

rowResult = RowReduce[Chop[result, 10^(-300)]];
Print[Accuracy[rowResult]];

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

fname=OpenWrite["Equations/3H_D+D1.txt"];
Write[fname, ResultEquations];
Close[fname];

Length[ResultEquations]