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

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

(* We have precomputed 100 points on FPP with high accuracy using NSolve Mathematica command. *)
Upoints=ToExpression[Import["Dependency/PointsOnFPP.txt"]];
Dimensions[Upoints]
Accuracy[Upoints]

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

In this section, we will try to construct equations of W and then Z.
The overall idea going from our base FPP to W is as follows:

1. There is an unique G-invariant section (up to scaling) of 6H + D which contains dim 10 spaces of quadratic relations. 

2. Let’s call U10, ..., U19 for the 10 new quadratics here, then define P0, ..., P19 such that 
P0 = U0, ..., P9 = U9  and  P10 = U10/sqrt{U10}, ..., P19 = U19/sqrt{U10}

3.  We can compute the surface W as quadratics in Pi’s satisfying the following “odd” and “even” relations:

Odd relation - $sum_{ i \in [0,9], j\in [10,19]} c_{ij} P_i P_j = 0$
which is equivalent to $\sum_{ i \in [0,9], j\in [10,19]} c_{ij} U_i U_j = 0$ 

Even relation - $\sum_{ i<=j \in [0,9] ]} c_{ij} P_i P_j + \sum_{ i<=j \in [10,19] ]} c_{ij} P_i P_j  = 0$
which is equivalent to $U * \sum_{ i<=j \in [0,9] ]} c_{ij} U_i U_j + \sum_{ i<=j \in [10,19] ]} c_{ij} U_i U_j  = 0$

Step 1 - Computing the G-invariant Section
The G-invariant section is the unique weight 0 section of 6H + D. Observe that
6H + D + (4H + D) + (4H + D) + (4H + D) = 18H
So we can realize this section as a cubic function with the 3 different sections of 4H + D and then find 10 quadratics vanishing on this section.

Us = {U0, U1, U2, U3, U4, U5, U6, U7, U8, U9};
Order3Aut={U1->U2,U2->U3,U3->U1,U4->U5,U5->U6,U6->U4,U7->U8,U8->U9,U9->U7};

(* Precomputed random points on the three sections of 4H + D *)
fname = "Dependency/points_on_4HD_first.txt";
PointsOn4HD=ToExpression[Import[fname]];

fname1 = "Dependency/points_on_4HD_second.txt";
PointsOn4HD=Join[PointsOn4HD, ToExpression[Import[fname1]]];

fname2 = "Dependency/points_on_4HD_third.txt";
PointsOn4HD=Join[PointsOn4HD, ToExpression[Import[fname2]]];

Dimensions[PointsOn4HD]
Accuracy[PointsOn4HD]

(* TODO *)

deg=3;
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];
Mod7Weights={0,6,5,3,1 ,2,4,1,2,4}; 

(* This is a subset of monomials which span the space of sections of 18H on Y, sorted by mod7 weight *)
UsefulMonomials ={{U6 U7 U8,U5 U6 U7,U4 U6 U8,U4 U5 U6,U3 (U8^2),U3 U5 U8,U3 (U5^2),(U3^2) U7,(U3^2) U4,U2 (U7^2),U2 U4 U7,U2 (U4^2),(U2^2) U6,U1 (U6^2),U1 U2 U3,(U1^2) U8,U7 U8 U9,U4 U8 U9,U4 U5 U9},{U6 (U8^2),U5 U6 U8,(U5^2) U6,U3 U6 U7,U3 U4 U6,(U3^2) U8,(U3^2) U5,U2 U7 U8,U2 U5 U7,U2 U4 U8,U2 U4 U5,(U2^3),U1 (U7^2),U1 U4 U7,U1 (U4^2),U1 U2 U6,(U1^2) U3,(U8^2) U9,U5 U8 U9},{(U6^2) U7,U4 (U6^2),U3 U6 U8,U3 U5 U6,(U3^3),U2 (U8^2),U2 U5 U8,U2 (U5^2),U2 U3 U7,U2 U3 U4,U1 U7 U8,U1 U5 U7,U1 U4 U8,U1 U4 U5,U1 (U2^2),(U1^2) U6,U6 U7 U9,U3 U8 U9,U4 (U9^2)},{(U7^3),(U6^2) U8,U5 (U6^2),U4 (U7^2),(U4^2) U7,(U4^3),(U3^2) U6,U2 U6 U7,U2 U4 U6,U2 U3 U8,U2 U3 U5,U1 (U8^2),U1 U5 U8,U1 (U5^2),U1 U3 U7,U1 U3 U4,(U1^2) U2,U6 U8 U9,U2 U4 U9,U8 (U9^2)},{(U7^2) U8,U5 (U7^2),U4 U7 U8,U4 U5 U7,(U4^2) U8,(U4^2) U5,U3 (U6^2),U2 U6 U8,U2 U5 U6,U2 (U3^2),(U2^2) U7,(U2^2) U4,U1 U6 U7,U1 U4 U6,U1 U3 U8,U1 U3 U5,(U1^3),U3 U6 U9,U2 U8 U9},{U7 (U8^2),(U6^3),U5 U7 U8,(U5^2) U7,U4 (U8^2),U4 U5 U8,U4 (U5^2),U3 (U7^2),U3 U4 U7,U3 (U4^2),U2 U3 U6,(U2^2) U8,(U2^2) U5,U1 U6 U8,U1 U5 U6,U1 (U3^2),U1 U2 U7,U1 U2 U4,(U6^2) U9,(U9^3)},{(U8^3),U6 (U7^2),U5 (U8^2),(U5^2) U8,(U5^3),U4 U6 U7,(U4^2) U6,U3 U7 U8,U3 U5 U7,U3 U4 U8,U3 U4 U5,U2 (U6^2),(U2^2) U3,U1 U3 U6,U1 U2 U8,U1 U2 U5,(U1^2) U7,(U1^2) U4,(U4^2) U9,U2 (U9^2)}};

(* Compute the weight 0 section of 6H + D *)
Equations = {};
Equationsrr = {};

Monomials = UsefulMonomials[[1]];

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

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

result = NullSpace[TableOfValues, Tolerance->10^(-270)];
(*Print["nullspace done"];*)
Print["Dimension ",Length[result]," solution space."];
(*Print[Accuracy[result]];*)


rowResult = RowReduce[result, Tolerance->10^(-270)];
(*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];

weight0section = newEquationsrr[[1]];

(* We calculate points on the weight 0 section of 6H + D. This was already precomputed. Feel free to uncomment the following code block to run it yourself *)

filename = "Dependency/points_on_6HD.txt";
PointsOn6HD=ToExpression[Import[filename]];

Dimensions[PointsOn6HD]
Accuracy[PointsOn6HD]

(*
filename = "Dependency/4HD_one_section.txt";
OpenRead[filename];
OneSection4HD=Read[filename]/.rr->(Sqrt[-7]);
Close[filename];

sectionsOn4HD = Join[OneSection4HD, OneSection4HD/.Order3Aut, OneSection4HD/.Order3Aut/.Order3Aut];

 (* Sanity check that this is C3 invariant *)
Print[Sort[sectionsOn4HD] === Sort[sectionsOn4HD/.Order3Aut]];
Print[Length[sectionsOn4HD]];

iter = 4;
OtherCurves = sectionsOn4HD;
PointsOn6HD= {};

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

RandomCoeff = RandomSample[Join[-Range[100],Range[100]], 10];
RandomCut = RandomCoeff.Us;
Print[RandomCut];
	
InputEquations = Join[EqsFPP, {weight0section, 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]]&];
PointsOn6HD = Join[PointsOn6HD, sortedPoints[[1;;10]]];

Print[Length[PointsOn6HD]];

, {i, 1, iter}];

*)

We are computing the quadratics vanishing on the invariant section of 6H+D.

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];
Mod7Weights={0,6,5,3,1 ,2,4,1,2,4};

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

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

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

result = NullSpace[TableOfValues, Tolerance->10^(-270)];
(* Print[{"Nullspace done",Dimensions[result],Accuracy[result]}]; *)
rowResult = RowReduce[Chop[result, 10^(-270)]];
(* Print[{"RowReduce done",Dimensions[rowResult],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];
]; 

aux=Flatten[Equationsrr];
(* Here we ensure that the 10 sections are generated by the C3 action *)
result = Range[1, 10];
result[[1]] = aux[[1]];
result[[2]] = aux[[2]];
result[[3]] = aux[[3]];
result[[4]] = aux[[2]]/.Order3Aut;
result[[5]] = aux[[3]]/.Order3Aut;
result[[6]] = aux[[6]];
result[[7]] = result[[4]]/.Order3Aut;
result[[8]] = result[[5]]/.Order3Aut;
result[[10]] = result[[6]]/.Order3Aut;
result[[9]] = result[[10]]/.Order3Aut;

(* Here we rearrange the order of the quadratics so that their weights are 0 6 5 3 1 2 4 1 2 4 in that order.*)

ginvariant = Range[1, 10];
ginvariant[[1]] = result[[1]];
ginvariant[[2]] =result[[10]];
ginvariant[[3]] =result[[9]];
ginvariant[[4]] =result[[6]];
ginvariant[[5]] =result[[2]];
ginvariant[[6]] =result[[4]];
ginvariant[[7]] =result[[7]];
ginvariant[[8]] =result[[3]];
ginvariant[[9]] =result[[5]];
ginvariant[[10]] =result[[8]];

(* Here we scale the first section U10 by this constant. This will be important in the construction of B^2/G4 later as it would ensure that the coefficients of equations generated lie in the number field Q(sqrt{-7}). *)
ginvariant[[1]] = (((49+13*rr)/56))*ginvariant[[1]];

Print["The U_{10}:"];
Print[RootReduce[ginvariant[[1]]/.rr->(Sqrt[-7])]];

(* Record the output.*)
fname=OpenWrite["Equations/ginvariant_section_c3.txt"];
Write[fname,Flatten[ginvariant]];
Close[fname];

ginvariant = ginvariant/.rr->(Sqrt[-7]);

Step 2: Now we construct the surface W.

Us = {U0, U1, U2, U3, U4, U5, U6, U7, U8, U9, U10, U11, U12, U13, U14, U15, U16, U17, U18, U19};
Ps = {P0, P1, P2, P3, P4, P5, P6, P7, P8, P9, P10, P11, P12, P13, P14, P15, P16, P17, P18, P19};

Mod7Weights={0,6,5,3,1,2,4,1,2,4,0,6,5,3,1,2,4,1,2,4};

ToMod7Weights={U0->0,U1->6,U2->5,U3->3,U4->1 ,U5->2,U6->4,U7->1,U8->2,U9->4,
U10-> 0,U11-> 6,U12-> 5,U13-> 3,U14->1,U15->2, U16->4,U17->1, U18->2, U19->4}; 

Order3Aut={P1->P2,P2->P3,P3->P1,P4->P5,P5->P6,P6->P4,P7->P8,P8->P9,P9->P7, 
P11->P12, P12->P13, P13->P11, P14->P15, P15->P16, P16->P14, P17->P18, P18->P19, P19->P17};

(* Here we define the rules of how U10 ... U19 translate to quadratics *)
conversion = {};
Do [
AppendTo[conversion, Us[[i+1]]->ginvariant[[i-9]]];
, {i, 10, 19}];

Do [
Print[term];
, {term, conversion}];

ListCoeffs=Table[ToExpression[StringJoin["c",ToString[i], ToString[j]]],{i,0,19}, {j, 0, 19}]; (* Formal symbols of coefficients cij *)
Dimensions[ListCoeffs]

We are now computing odd equations.

AllMonomials = Flatten[Table[Us[[i]]Us[[j]],{i,1,10}, {j, 11, 20}]];
AllMonomialsinPs = Flatten[Table[Ps[[i]] Ps[[j]],{i,1,10}, {j, 11, 20}]];MonomialWeights = Flatten[Table[Mod[Mod7Weights[[i]]+Mod7Weights[[j]],7],{i,1,10}, {j, 11, 20}]];

Equations = {};
Equationsrr = {};
EquationsinPs = {};

For[mod7=0,mod7≤6,mod7++,
Print["mod7=",mod7];
Monomials = {};
MonomialsinPs = {};
Do[
If[MonomialWeights[[i]]== mod7,  
AppendTo[Monomials, AllMonomials[[i]]];
AppendTo[MonomialsinPs, AllMonomialsinPs[[i]]];
];
, {i, 1, Length[AllMonomials]}];
(*
Print[Monomials];
Print[Length[Monomials]];
*)

TableOfValues=Table[(Monomials/.conversion)/.Upoints[[l]],{l,1,Length[Upoints]}];

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

result = NullSpace[TableOfValues];

(*Print["nullspace done"];
Print[Dimensions[result]];
Print[Accuracy[result]];*)

rowResult = RowReduce[result, Tolerance-> 10^(-Round[Accuracy[result] ]+10)];

(*Print["Rowreduce done"];
Print[Accuracy[rowResult]];*)

exactResult = RootApproximant[Chop[rowResult, 10^(-400)]];
rowResultrr =Table[ToRR[exactResult[[i]]],{i,1,Length[exactResult]}];
newEquationsrr = rowResultrr.Monomials;
newEquations = exactResult.Monomials;
newEquationsPs = rowResultrr.MonomialsinPs;
(*Print[newEquations];*)

AppendTo[Equations, newEquations];
AppendTo[Equationsrr, newEquationsrr];
AppendTo[EquationsinPs, newEquationsPs];
];

oddResult= Flatten[EquationsinPs];
Length[Flatten[oddResult]]

We are now computing even equations.

U = U10;
AllMonomials =Join[ Flatten[Table[Us[[i]]Us[[j]],{i,1,10}, {j, 1, i}]],Flatten[Table[U^(-1)Us[[i]]Us[[j]],{i,11,20}, {j, 11, i}]]];
AllMonomialsinPs =Join[ Flatten[Table[Ps[[i]]Ps[[j]],{i,1,10}, {j, 1, i}]],Flatten[Table[Ps[[i]]Ps[[j]],{i,11,20}, {j, 11, i}]]];
MonomialWeights =Join[ Flatten[Table[Mod7Weights[[i]]+Mod7Weights[[j]],{i,1,10}, {j, 1, i}]],Flatten[Table[Mod7Weights[[i]]+Mod7Weights[[j]],{i,11,20}, {j, 11, i}]]];
Equations = {};
Equationsrr = {};
EquationsinPs = {};

For[mod7=0,mod7≤6,mod7++,
Print["mod7=",mod7];
Monomials = {};
MonomialsinPs = {};
Do[
If[Mod[MonomialWeights[[i]],7]== mod7,  
AppendTo[Monomials, AllMonomials[[i]]];
AppendTo[MonomialsinPs, AllMonomialsinPs[[i]]];
];
, {i, 1, Length[AllMonomials]}];
(*
Print[Monomials];
Print[Length[Monomials]];
*)

TableOfValues=Table[(Monomials/.conversion)/.Upoints[[l]],{l,1,Length[Upoints]}];
(*
Print["matrix done"];
Print[Dimensions[TableOfValues]];
Print[MatrixRank[TableOfValues]];
Print[Accuracy[TableOfValues]];
*)

result = NullSpace[TableOfValues];
(*
Print["nullspace done"];
Print[Dimensions[result]];
Print[Accuracy[result]];
*)

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

exactResult = RootApproximant[Chop[rowResult, 10^(-400)]];
rowResultrr =Table[ToRR[exactResult[[i]]],{i,1,Length[exactResult]}];
newEquationsrr = rowResultrr.Monomials;
newEquations = exactResult.Monomials;
newEquationsPs = rowResultrr.MonomialsinPs;
(*Print[newEquations];*)

AppendTo[Equations, newEquations];
AppendTo[Equationsrr, newEquationsrr];
AppendTo[EquationsinPs, newEquationsPs];
];

evenResult= Flatten[EquationsinPs];
Length[Flatten[evenResult]]

result = Join[oddResult, evenResult];
Length[result]
fname=OpenWrite["Equations/W_equations.txt"];
Write[fname,result];
Close[fname];

Now we will try to construct equations of Z.

Order3Ps={P1->P2,P2->P3,P3->P1,P4->P5,P5->P6,P6->P4,P7->P8,P8->P9,P9->P7, 
P11->P12, P12->P13, P13->P11, P14->P15, P15->P16, P16->P14, P17->P18, P18->P19, P19->P17};

PsToUs =Join[Table[Ps[[i]]-> Us[[i]],{i,1,10}],Table[Ps[[i]]-> Us[[i]]/(U10)^(1/2),{i,11,20}]];

(*Import file of equations for FPP*)
filename = "Dependency/EqsFPP.txt";
EqsFPP=ToExpression[Import[filename]]/.rr->(Sqrt[-7]);

(*Import file of equations for W*)
filename = "Equations/W_equations.txt";
EqsW=ToExpression[Import[filename]]/.rr->(Sqrt[-7]);

Step 1: We compute s1, s2, s3, and s4

ListCoeffs=Table[ToExpression[StringJoin["c",ToString[i]]],{i,1,20}]; (* Formal symbols of coefficients ci *)


(* Weight (1, 3) *)
a1=1;
a2 = 3;

(* Here we set up the general forms of s3, s4, s5, and s6 where s5^2 - s6^2 = s1s2, we can recover s1 = s5 + s6 and s2 = s5 - s6 *)
count = 1;

s3Weight = Mod[2*a1, 7];
s3Index = Select[Flatten[Position[Mod7Weights, s3Weight]], # ≤10 &];
s3Var = Function[x, Ps[[x]]]/@s3Index;
s3 = ListCoeffs[[count;;count+Length[s3Var]-1]].s3Var;

count += Length[s3Var];

s4Weight = Mod[2*a2, 7];
s4Index = Select[Flatten[Position[Mod7Weights, s4Weight]], # ≤10 &];
s4Var = Function[x, Ps[[x]]]/@s4Index;
s4 = ListCoeffs[[count;;count+Length[s4Var]-1]].s4Var;

count += Length[s4Var];

s5Weight = Mod[a1+a2, 7];
s5Index = Select[Flatten[Position[Mod7Weights, s5Weight]], # ≤10 &];
s5Var = Function[x, Ps[[x]]]/@s5Index;
s5 = ListCoeffs[[count;;count+Length[s5Var]-1]].s5Var;

count += Length[s5Var];

s6Weight = Mod[a1+a2, 7];
s6Index = Select[Flatten[Position[Mod7Weights, s6Weight]], # > 10 &];
s6Var = Function[x, Ps[[x]]]/@s6Index;
s6 = ListCoeffs[[count;;count+Length[s6Var]-1]].s6Var;

count += Length[s6Var] - 1;
cS = ListCoeffs[[1;;count]];

quadric = (s5^2 - s6^2) - s3*s4;
current =  quadric/.PsToUs
current = current/.conversion;

(* Here we add relations on the s3, s4, s5, and s6 *)
relations = {};
Do [
AppendTo[relations, current/.points];
, {points, Upoints}];

(* From https://mathematica.stackexchange.com/questions/88762/how-to-get-a-list-of-monomials-of-a-polynomial-without-coefficients*)
fn2[x_,var_]:=Collect[List@@Expand@x,var,1&]
Monomials = fn2[relations[[3]], cS];

Clear[a1, a2];

coeffMatrix = Table[Coefficient[relations[[k]],Monomials], {k, 1, Length[Upoints]}];

{Dimensions[coeffMatrix],MatrixRank[coeffMatrix],Accuracy[coeffMatrix]}

reducedCoeff =RowReduce[Chop[coeffMatrix, 10^(-800)]];
(*Accuracy[reducedCoeff]*)
exactCoeff=RootApproximant[Chop[reducedCoeff,10^(-800)], 2];

sol=Solve[exactCoeff.Monomials == 0]
output=quadric/.sol[[1]];

(* Now we can find s1 and s2 *)
(* We have two free variables c1 and c7 and we will set them to 1 *)
s1 = (s5+s6)/.sol[[1]]/.{c1->1, c7->1};
s2 = (s5-s6)/.sol[[1]]/.{c1->1, c7->1};
s3=s3/.sol[[1]]/.{c1->1, c7->1};
s4=s4/.sol[[1]]/.{c1->1, c7->1};

sections = {s1, s2, s3, s4};
RootReduce[Collect[sections,Ps] ]

sectionsUs = sections/.PsToUs/.conversion;

Step 2: Here we construct r_{1, 1}, r_{-1, 1}, and r_{0, 0}
Recall that s1 = r_{1,1} r_{-1,3}; s2 = r_{-1,1} r_{1,3}; s3 = r_{1,1} r_{-1,1}; s4 = r_{1,3} r_{-1,3}

(* Here we compute r_1, 1 *)

(* Random points on r_1, 1. We already precomputed it but feel free to use the code block below *)

fname="Dependency/r11_points.txt";
PointsCommon=ToExpression[Import[fname]];

{Dimensions[PointsCommon],Accuracy[PointsCommon]}

(* 
iter = 30;
PointsCommon={};
For[k=1,k≤iter,k++,
Print["Iteration ", k];
current = NSolve[Join[{0==s1,0==s3,1==P2, RandomInteger[{-300, 300}]==P3},Table[0==EqsW[[i]],{i,1,Length[EqsW]}]],500];

If [Length[current] ≠ 0,
 AppendTo[PointsCommon,current[[1]]];
];
];
*)

(* TBD *)
UsefulMonomials={{P3 P6,P2 P8,P2 P5,P19 P3,P1 P7,P1 P4,P18 P2,P1 P17,P16 P3,P15 P2,P1 P14,P13 P6,P13 P19,P13 P16},{(P6^2),P2 P3,P19 P6,(P19^2),P1 P8,P1 P5,P1 P18,P16 P6,P16 P19,(P16^2),P1 P15,P13 P2,P12 P3,P12 P13,P11 P8,P19 P9},{(P7^2),P4 P7,(P4^2),P2 P6,P19 P2,P1 P3,P17 P7,P17 P4,(P17^2),P16 P2,P14 P7,P14 P4,P1 P13,P12 P6,P12 P19,P12 P16},{P7 P8,P5 P7,P4 P8,P4 P5,(P2^2),P1 P6,P1 P19,P18 P7,P18 P4,P17 P8,P17 P5,P17 P18,P1 P16,P15 P7,P15 P4,P15 P17},{(P8^2),P5 P8,(P5^2),P3 P7,P3 P4,P1 P2,P18 P8,P18 P5,(P18^2),P17 P3,P15 P8,P15 P5,P14 P3,P13 P7,P13 P4,P13 P17},{P6 P7,P4 P6,P3 P8,P3 P5,P19 P7,P19 P4,(P1^2),P18 P3,P17 P6,P17 P19,P16 P7,P16 P4,P16 P17,P15 P3,P14 P6,P14 P19},{P6 P8,P5 P6,(P3^2),P2 P7,P2 P4,P19 P8,P19 P5,P18 P6,P18 P19,P17 P2,P16 P8,P16 P5,P16 P18,P15 P6,P15 P19,P14 P2}};

deg=2;

Powers={a1,a2,a3,a4,a5,a6,a7,a8,a9,a10, a11, a12, a13, a14, a15, a16, a17, a18, a19, a20};
AllPowers=Solve[Join[Table[Powers[[k]]≥0,{k,1,20}],{Sum[Powers[[k]],{k,1,20}]==deg}],Integers];

Equations = {};
Equationsrr = {};

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

Monomials = UsefulMonomials[[mod7 + 1]];
(*Print[Monomials];*)

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

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

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


rowResult = RowReduce[Chop[result, 10^(-270)]];
(*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];


];

fname=OpenWrite["Equations/r11_equations.txt"];
Write[fname,Flatten[Equationsrr]];
Close[fname];
Length[Flatten[Equationsrr]](* 56 equations *)

(* Here we compute r_-1, 1 *)
(* Random points on r_-1, 1. We already precomputed it but feel free to use the code block below *)

fname="Dependency/r-11_points.txt";
PointsCommon=ToExpression[Import[fname]];

{Dimensions[PointsCommon],Accuracy[PointsCommon]}

(* 
iter = 30;
PointsCommon={};
For[k=1,k≤iter,k++,
Print["Iteration ", k];
current = NSolve[Join[{0==s2,0==s3,1==P2, RandomInteger[{-300, 300}]==P3},Table[0==EqsW[[i]],{i,1,Length[EqsW]}]],500];

If [Length[current] ≠ 0,
 AppendTo[PointsCommon,current[[1]]];
];
];
*)

(* TBD *)

Equations = {};
Equationsrr = {};

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

Monomials = UsefulMonomials[[mod7 + 1]];
(*Print[Monomials];*)

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

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

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


rowResult = RowReduce[Chop[result, 10^(-270)]];
(*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];


];

fname=OpenWrite["Equations/r-11_equations.txt"];
Write[fname,Flatten[Equationsrr]];
Close[fname];
Length[Flatten[Equationsrr]](* 56 equations *)

(* Here we compute r_0, 0. This is given by U0 = 0 *)

fname="Dependency/r00_points.txt";
PointsCommon=ToExpression[Import[fname]];

Dimensions[PointsCommon]
Accuracy[PointsCommon]

(* Points on r_0, 0 *)
(* 
iter = 20;
PointsR00={};
PointsR00Us={};

For[k=1,k≤iter,k++,
Print["Iteration ", k];
RandomCoeff = RandomSample[Join[-Range[100],Range[100]], 10];
Print[RandomCoeff];
RandomCut = 0;

Do [
RandomCut = RandomCut + RandomCoeff[[i]]*Us[[i]], {i, 1, 10}];
Print[RandomCut];

current = NSolve[Join[{0==U0, U3==1, RandomCut == 0},Table[0==EqsFPP[[i]],{i,1,Length[EqsFPP]}]],500];
(*Print[current];*)
AppendTo[PointsR00Us,current[[1]]];
];

Do [
currentPoint  = PointsR00Us[[i]];

Do [
AppendTo[currentPoint, Us[[j]]->ginvariant[[j-10]]/.currentPoint]
, {j, 11, 20}];

newPoint = {};
Do [
AppendTo[newPoint, Ps[[j]]->Us[[j]]/.currentPoint];
, {j, 1, 10}];

Do [
AppendTo[newPoint, Ps[[j]]->(Us[[j]]/(U10)^(1/2))/.currentPoint];
, {j, 11, 20}];

AppendTo[PointsR00, newPoint];
, {i, 1, Length[PointsR00Us]}];

*)

Equations = {};
Equationsrr = {};

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

Monomials = UsefulMonomials[[mod7 + 1]];
(*Print[Monomials];*)

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

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

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


rowResult = RowReduce[Chop[result, 10^(-270)]];
(*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];
];

fname=OpenWrite["Equations/r00_equations.txt"];
Write[fname,Flatten[Equationsrr]];
Close[fname];
Length[Flatten[Equationsrr]](* 56 equations *)

Note that z0 = r_0,0, z1 = r_1, 1 and z4 = r_-1, 1.
We can obtain z2 and z3 as C3 translates of z1. And we can obtain z5 and z6 as C3 translates of z4. 

For our computations in Step 3 later, we also need random points on z2, z3, z5, and z6. We can collect points on them by computing random cuts with C3 translates of the equations we computed above.

We already precomputed points on them in the folder Dependency/ZPoints (including z0, z1, z4 we already computed). Here we show an example of how we would compute points on z2.

(*

filename = "Equations/r11_equations.txt";
OpenRead[filename];
Z1Equations=Read[filename]/.rr->(Sqrt[-7]);
Close[filename];

Z2Equations = Z1Equations/.Order3Ps;

iter = 20;
PointsZ2= {};
For[k=1,k≤iter,k++,
Print["Iteration ", k];

RandomCoeff = RandomSample[Join[-Range[100],Range[100]], 20];
Print[RandomCoeff];
RandomCut = 0;

Do [
RandomCut = RandomCut + RandomCoeff[[i]]*Ps[[i]], {i, 1, 20}];
Print[RandomCut];

current = NSolve[Join[Z2Equations, {1==P2, RandomCut==0},Table[0==EqsW[[i]],{i,1,Length[EqsW]}]],500];

If [Length[current] ≠ 0,
 AppendTo[PointsZ2,current[[1]]];
];
];

*)

Step 3: Here we compute the rational functions r_1, ..., r_7

(* Here we set up some preliminary data *)

(* Points on zi *)
filename = "Dependency/ZPoints/z0.txt";
PointsZ0=ToExpression[Import[filename]];
Dimensions[PointsZ0]

filename = "Dependency/ZPoints/z1.txt";
PointsZ1=ToExpression[Import[filename]];
Dimensions[PointsZ1]

filename = "Dependency/ZPoints/z2.txt";
PointsZ2=ToExpression[Import[filename]];
Dimensions[PointsZ2]

filename = "Dependency/ZPoints/z3.txt";
PointsZ3=ToExpression[Import[filename]];
Dimensions[PointsZ3]

filename = "Dependency/ZPoints/z4.txt";
PointsZ4=ToExpression[Import[filename]];
Dimensions[PointsZ4]

filename = "Dependency/ZPoints/z5.txt";
PointsZ5=ToExpression[Import[filename]];
Dimensions[PointsZ5]

filename = "Dependency/ZPoints/z6.txt";
PointsZ6=ToExpression[Import[filename]];
Dimensions[PointsZ6]

(* We also want the Jacobian of the equations of W. This is to add relations for zi^2 *)
JacobianW=Factor[Table[D[EqsW[[i]],Ps[[k]]],{i,1,Length[EqsW]},{k,1,Length[Ps]}]];
Dimensions[JacobianW]

(* C3 Useful Monomials - explanation TBD *)
(* C3 invariant useful monomials - to keep the action nice later *)
UsefulMonomialsC3 ={{P1*P4,P2*P5,P3*P6,P1*P7,P2*P8,P3*P9,P1*P17,P18*P2,P19*P3,P1*P14,P15*P2,P16*P3,P0^2,P0*P10},{P6^2,P2*P3,P19*P6,P19^2,P1*P8,P1*P5,P1*P18,P16*P6,P16*P19,P16^2,P1*P15,P13*P2,P12*P3,P12*P13,P11*P8,P19*P9},{P4^2,P1*P3,P17*P4,P17^2,P2*P9,P2*P6,P19*P2,P14*P4,P14*P17,P14^2,P16*P2,P11*P3,P1*P13,P11*P13,P12*P9,P17*P7},{P5*P7,P4*P5,P2^2,P1*P9,P1*P6,P18*P7,P18*P4,P17*P5,P17*P18,P1*P19,P15*P7,P15*P4,P15*P17,P14*P5,P14*P18,P1*P16},{P5^2,P1*P2,P18*P5,P18^2,P3*P7,P3*P4,P17*P3,P15*P5,P15*P18,P15^2,P14*P3,P1*P12,P11*P2,P11*P12,P13*P7,P18*P8},{P4*P9,P4*P6,P1^2,P3*P8,P3*P5,P17*P9,P17*P6,P19*P4,P17*P19,P18*P3,P14*P9,P14*P6,P14*P19,P16*P4,P16*P17,P15*P3},{P6*P8,P5*P6,P3^2,P2*P7,P2*P4,P19*P8,P19*P5,P18*P6,P18*P19,P17*P2,P16*P8,P16*P5,P16*P18,P15*P6,P15*P19,P14*P2}};

WeightZeroMonomials = UsefulMonomialsC3[[0+1]];
WeightFiveMonomials = UsefulMonomialsC3[[5+1]];
WeightSixMonomials = UsefulMonomialsC3[[6+1]];

(* C3 translates of s3 *)
s3a=s3/.Order3Ps
s3b=s3a/.Order3Ps

3.1 - We can compute r1 = z6 z1^2 z0 / s3’’ P0

(* Now we solve the coefficients from the relations *)
ListCoeffs=Table[ToExpression[StringJoin["c",ToString[i]]],{i,1,100}];
Quadratic = ListCoeffs[[1;;Length[WeightSixMonomials]]].WeightSixMonomials;

(* Singular relations for z1^2 *)
JacobianQuad = Table[D[Quadratic, Ps[[k]]], {k, 1, Length[Ps]}];
SingularRelations = {};
Do [
input =NullSpace[JacobianW/.pt];
TanSpace=RowReduce[Chop[input, 10^(-400)]];

(*Print[Dimensions[TanSpace]];*)
CurrentJacob = JacobianQuad/.pt;

SingularRelations = Join[SingularRelations, TanSpace.CurrentJacob];

, {pt, PointsZ1}];

Dimensions[SingularRelations];

(* Relations to vanish on z0, z1, and z6 *)
VanishRelations = {};
AllPoints = Join[PointsZ0,PointsZ1, PointsZ6];
Dimensions[AllPoints];

Do [
AppendTo[VanishRelations, Quadratic/.pt];
, {pt, AllPoints}];

Dimensions[VanishRelations];

(* Here we solve for the quadratic *)
Relations = Join[VanishRelations, SingularRelations];
coeffMatrix = Table[Coefficient[Relations[[k]],ListCoeffs[[1;;Length[WeightSixMonomials]]]], {k, 1, Length[Relations]}];

(*Dimensions[coeffMatrix]
MatrixRank[coeffMatrix]
Accuracy[coeffMatrix]*)

reducedCoeff =RowReduce[coeffMatrix,Tolerance-> 10^(-300)];
Accuracy[reducedCoeff];
exactCoeff=RootApproximant[Chop[reducedCoeff,10^(-300)], 2];

(* We can normalize c15 = 1 *)
sol=Solve[Append[exactCoeff.ListCoeffs[[1;;Length[WeightSixMonomials]]],c15-1] == 0];
r1 =(Quadratic/. RootReduce[sol[[1]]])/(s3b*P0)

3.2 - We can compute r2 = z0 z1 z2 z3 / s3’’ P0

(* Now we solve the coefficients from the relations *)
ListCoeffs=Table[ToExpression[StringJoin["c",ToString[i]]],{i,1,100}];
Quadratic = ListCoeffs[[1;;Length[WeightZeroMonomials]]].WeightZeroMonomials;

SingularRelations = {};
(* Relations to vanish on z0, z1, z2 and z3 *)
VanishRelations = {};
AllPoints = Join[PointsZ0,PointsZ1, PointsZ2, PointsZ3];
Dimensions[AllPoints];

Do [
AppendTo[VanishRelations, Quadratic/.pt];
, {pt, AllPoints}];

Dimensions[VanishRelations];

(* Here we solve for the quadratic *)
Relations = Join[VanishRelations, SingularRelations];
coeffMatrix = Table[Coefficient[Relations[[k]],ListCoeffs[[1;;Length[WeightZeroMonomials]]]], {k, 1, Length[Relations]}];

Dimensions[coeffMatrix];
MatrixRank[coeffMatrix];
Accuracy[coeffMatrix];

reducedCoeff =RowReduce[Chop[coeffMatrix, 10^(-300)]];
Accuracy[reducedCoeff];
exactCoeff=RootApproximant[Chop[reducedCoeff,10^(-300)], 2];

(* We set c1 = 1 *)
sol=Solve[Append[exactCoeff.ListCoeffs[[1;;Length[WeightZeroMonomials]]],c1-1] == 0];
r2=(Quadratic/.RootReduce[sol[[1]]])/(s3b*P0)

3.3 - We can compute r3 = s3/P0

r3 = s3/P0;

3.4 - We can compute r4 = z2^2 z4 z0 / s3’ P0

(* Now we solve the coefficients from the relations *)
ListCoeffs=Table[ToExpression[StringJoin["c",ToString[i]]],{i,1,100}];
Quadratic = ListCoeffs[[1;;Length[WeightFiveMonomials]]].WeightFiveMonomials;

(* Singular relations for z2^2 *)
JacobianQuad = Table[D[Quadratic, Ps[[k]]], {k, 1, Length[Ps]}];
SingularRelations = {};
Do [
input =NullSpace[JacobianW/.pt];
TanSpace=RowReduce[Chop[input, 10^(-400)]];

(*Print[Dimensions[TanSpace]];*)
CurrentJacob = JacobianQuad/.pt;

SingularRelations = Join[SingularRelations, TanSpace.CurrentJacob];

, {pt, PointsZ2}];

Dimensions[SingularRelations];

(* Relations to vanish on z0, z2, and z4 *)
VanishRelations = {};
AllPoints = Join[PointsZ0,PointsZ2, PointsZ4];
Dimensions[AllPoints];

Do [
AppendTo[VanishRelations, Quadratic/.pt];
, {pt, AllPoints}];

Dimensions[VanishRelations];

(* Here we solve for the quadratic *)
Relations = Join[VanishRelations, SingularRelations];
coeffMatrix = Table[Coefficient[Relations[[k]],ListCoeffs[[1;;Length[WeightFiveMonomials]]]], {k, 1, Length[Relations]}];

Dimensions[coeffMatrix];
MatrixRank[coeffMatrix];
Accuracy[coeffMatrix];

reducedCoeff =RowReduce[Chop[coeffMatrix, 10^(-300)]];
Accuracy[reducedCoeff];
exactCoeff=RootApproximant[Chop[reducedCoeff,10^(-300)], 2];

(* add c15 = 1 *)
sol=Solve[Append[exactCoeff.ListCoeffs[[1;;Length[WeightFiveMonomials]]],c15-1] == 0];
r4=(Quadratic/. RootReduce[sol[[1]]])/(s3a*P0)

3.5 - We can compute r5 = z1 z5^2 z0 / s3’ P0

(* Now we solve the coefficients from the relations *)
ListCoeffs=Table[ToExpression[StringJoin["c",ToString[i]]],{i,1,100}];
Quadratic = ListCoeffs[[1;;Length[WeightFiveMonomials]]].WeightFiveMonomials;

(* Singular relations for z5^2 *)
JacobianQuad = Table[D[Quadratic, Ps[[k]]], {k, 1, Length[Ps]}];
SingularRelations = {};
Do [
input =NullSpace[JacobianW/.pt];
TanSpace=RowReduce[Chop[input, 10^(-400)]];

(*Print[Dimensions[TanSpace]];*)
CurrentJacob = JacobianQuad/.pt;

SingularRelations = Join[SingularRelations, TanSpace.CurrentJacob];

, {pt, PointsZ5}];

Dimensions[SingularRelations];

(* Relations to vanish on z0, z1, and z5 *)
VanishRelations = {};
AllPoints = Join[PointsZ0,PointsZ1, PointsZ5];
Dimensions[AllPoints];

Do [
AppendTo[VanishRelations, Quadratic/.pt];
, {pt, AllPoints}];

Dimensions[VanishRelations];

(* Here we solve for the quadratic *)
Relations = Join[VanishRelations, SingularRelations];
coeffMatrix = Table[Coefficient[Relations[[k]],ListCoeffs[[1;;Length[WeightFiveMonomials]]]], {k, 1, Length[Relations]}];

Dimensions[coeffMatrix];
MatrixRank[coeffMatrix];
Accuracy[coeffMatrix];

reducedCoeff =RowReduce[Chop[coeffMatrix, 10^(-300)]];
Accuracy[reducedCoeff];
exactCoeff=RootApproximant[Chop[reducedCoeff,10^(-300)], 2];

(* c15 = 1 *)
sol=Solve[Append[exactCoeff.ListCoeffs[[1;;Length[WeightFiveMonomials]]],c15-1] == 0];
r5= (Quadratic/.RootReduce[sol[[1]]])/(s3a*P0)

3.6 - We can compute r6 = z4^2 z3 z0 / s3’’ P0

(* Now we solve the coefficients from the relations *)
ListCoeffs=Table[ToExpression[StringJoin["c",ToString[i]]],{i,1,100}];
Quadratic = ListCoeffs[[1;;Length[WeightSixMonomials]]].WeightSixMonomials;

(* Singular relations for z4^2 *)
JacobianQuad = Table[D[Quadratic, Ps[[k]]], {k, 1, Length[Ps]}];
SingularRelations = {};
Do [
input =NullSpace[JacobianW/.pt];
TanSpace=RowReduce[Chop[input, 10^(-400)]];

(*Print[Dimensions[TanSpace]];*)
CurrentJacob = JacobianQuad/.pt;

SingularRelations = Join[SingularRelations, TanSpace.CurrentJacob];

, {pt, PointsZ4}];

Dimensions[SingularRelations];

(* Relations to vanish on z0, z3, and z4 *)
VanishRelations = {};
AllPoints = Join[PointsZ0,PointsZ3, PointsZ4];
Dimensions[AllPoints];

Do [
AppendTo[VanishRelations, Quadratic/.pt];
, {pt, AllPoints}];

Dimensions[VanishRelations];

(* Here we solve for the quadratic *)
Relations = Join[VanishRelations, SingularRelations];
coeffMatrix = Table[Coefficient[Relations[[k]],ListCoeffs[[1;;Length[WeightSixMonomials]]]], {k, 1, Length[Relations]}];

Dimensions[coeffMatrix];
MatrixRank[coeffMatrix];
Accuracy[coeffMatrix];

reducedCoeff =RowReduce[Chop[coeffMatrix, 10^(-300)]];
Accuracy[reducedCoeff];
exactCoeff=RootApproximant[Chop[reducedCoeff,10^(-300)], 2];

(* c15 = 1 *)
sol=Solve[Append[exactCoeff.ListCoeffs[[1;;Length[WeightSixMonomials]]],c15-1] == 0];
r6= (Quadratic/.RootReduce[sol[[1]]])/(s3b*P0)

3.7 - We can compute r7 = z0 z4 z5 z6 / s3’’ P0

(* Now we solve the coefficients from the relations *)
ListCoeffs=Table[ToExpression[StringJoin["c",ToString[i]]],{i,1,100}];
Quadratic = ListCoeffs[[1;;Length[WeightZeroMonomials]]].WeightZeroMonomials;

SingularRelations = {};
(* Relations to vanish on z0, z4, z5 and z6 *)
VanishRelations = {};
AllPoints = Join[PointsZ0,PointsZ4, PointsZ5, PointsZ6];
Dimensions[AllPoints];

Do [
AppendTo[VanishRelations, Quadratic/.pt];
, {pt, AllPoints}];

Dimensions[VanishRelations];

(* Here we solve for the quadratic *)
Relations = Join[VanishRelations, SingularRelations];
coeffMatrix = Table[Coefficient[Relations[[k]],ListCoeffs[[1;;Length[WeightZeroMonomials]]]], {k, 1, Length[Relations]}];

Dimensions[coeffMatrix];
MatrixRank[coeffMatrix];
Accuracy[coeffMatrix];

reducedCoeff =RowReduce[Chop[coeffMatrix, 10^(-300)]];
Accuracy[reducedCoeff];
exactCoeff=RootApproximant[Chop[reducedCoeff,10^(-300)], 2];

(* c1 = 1 *)
sol=Solve[Append[exactCoeff.ListCoeffs[[1;;Length[WeightZeroMonomials]]],c1-1] == 0];
r7=(Quadratic/.RootReduce[sol[[1]]])/(s3b*P0)

Step 4: Here we compute the multiplication table.

(* Here we build the multiplication table between Z0, ..., Z6. Note there are Z7, ..., Z12, but we can compute that later using a C2 action. *)
Zs = {Z0, Z1, Z2, Z3, Z4, Z5, Z6};
Order3Zs = {Z1 -> Z2, Z2->Z3, Z3->Z1, Z4->Z5, Z5->Z6, Z6->Z4};

(* We also want random points on W, which can get from random points of the FPP*)
WPoints = {};
Do [
currentPoint  = Upoints[[i]];


Do [
AppendTo[currentPoint,Us[[j]]-> ( Us[[j]]/.conversion)/.currentPoint];
, {j, 11, 20}];


newPoint = {};
Do [
AppendTo[newPoint, Ps[[j]]->Us[[j]]/.currentPoint];
, {j, 1, 10}];

Do [
AppendTo[newPoint, Ps[[j]]->(Us[[j]]/(U10)^(1/2))/.currentPoint];
, {j, 11, 20}];

AppendTo[WPoints, newPoint];
, {i, 1, Length[Upoints]}];

Dimensions[WPoints];
Accuracy[WPoints];

(* We need translates of the 4 sections s1, s2, s3, and s4 *)
s1a=s1/.Order3Ps;
s1b=s1a/.Order3Ps;

s2a=s2/.Order3Ps;
s2b=s2a/.Order3Ps;

s3a=s3/.Order3Ps;
s3b=s3a/.Order3Ps;

s4a = s4/.Order3Ps;
s4b=s4a/.Order3Ps;

(* Here we set up the rational functions *)
RFunctions = {R1, R2, R3, R4, R5, R6, R7};
As = {a1, a2, a3, a4, a5, a6, a7};
RatToPs = {R1 -> r1, R2 -> r2, R3-> r3, R4->r4, R5->r5, R6->r6, R7->r7};

(* Here we begin to construct the multiplication table, we enter the 7 manual entries correspondent to r1, ..., r7 and the trivial multiplications first. *)
(* The rest can be populated using a C3 action *)

(* The (i, j)-entry of the matrix will represent z_i/z_0 * z_j/z_0 *)
MultiplicationTable = ConstantArray[0, {7, 7}];

Do [
MultiplicationTable[[1]][[i]] = Zs[[i]]/Z0;
MultiplicationTable[[i]][[1]] = Zs[[i]]/Z0;
, {i, 1, 7}];

(* Populate Calculated Entries *)

index = {{1, 1}, {1, 2}, {1, 4}, {2, 4}, {1, 5}, {4, 4}, {4, 5}};
FilledIndex = {};
Do [
AppendTo[FilledIndex, i];
AppendTo[FilledIndex, Reverse[i]]; 
, {i, index}];
FilledIndex (* see which entries are filled *)

MultiplicationTable [[1+1]][[1+1]] = a1*R1*Z3/Z0;

MultiplicationTable [[1+1]][[2+1]] = a2*R2*Z6/Z0;
MultiplicationTable [[2+1]][[1+1]] = a2*R2*Z6/Z0;

MultiplicationTable [[1+1]][[4+1]] = a3*R3;
MultiplicationTable [[4+1]][[1+1]] = a3*R3;

MultiplicationTable [[2+1]][[4+1]] = a4*R4*Z5/Z0;
MultiplicationTable [[4+1]][[2+1]] = a4*R4*Z5/Z0;

MultiplicationTable [[1+1]][[5+1]] = a5*R5*Z2/Z0;
MultiplicationTable [[5+1]][[1+1]] = a5*R5*Z2/Z0;

MultiplicationTable [[4+1]][[4+1]] = a6*R6*Z6/Z0;

MultiplicationTable [[4+1]][[5+1]] = a7*R7*Z3/Z0;
MultiplicationTable [[5+1]][[4+1]] = a7*R7*Z3/Z0;


 (* Symmetry *)
SymmetricMatrixQ[MultiplicationTable]

MultiplicationTable = MultiplicationTable/.RatToPs;

(* Here we impose the C3 action to fill up the display and actual multiplication table.*)
Order3Index = {2, 3, 1, 5, 6, 4};

Do [
i = pair[[1]];
j = pair[[2]];
entry = MultiplicationTable[[i+1]][[j+1]];

(*Print[{i, j}];*)

i1 = Order3Index[[i]];
j1 = Order3Index[[j]];
MultiplicationTable[[i1+1]][[j1+1]] = entry/.Order3Ps/.Order3Zs;

entry1 = MultiplicationTable[[i1+1]][[j1+1]];

i2 = Order3Index[[i1]];
j2 = Order3Index[[j1]];
MultiplicationTable[[i2+1]][[j2+1]] =  entry1/.Order3Ps/.Order3Zs;


entry2 = MultiplicationTable[[i2+1]][[j2+1]];


(*Print[{i2, j2}];*)

entry3 = entry2/.Order3Ps/.Order3Zs;

(* Checking that the entires are C3 invariant *)
Print[{i1, j1}," ",FullSimplify[entry==entry3]];
, {pair, FilledIndex}];

(* Checking the table is still symmetric *)
SymmetricMatrixQ[MultiplicationTable]

(* Here we build a helper function CustomMultiply that translates what "multiplication means" between z_i/z_0 and z_j/z_0 *)
TranslationTable = {};
Do [
Do [
translation = {Zs[[i]]*Zs[[j]]/Z0^2 -> MultiplicationTable[[i]][[j]]};
AppendTo[TranslationTable, translation];
, {j, 2, i}];
, {i, 2, 7}];

Length[TranslationTable];

CustomMultiply[input_] := Module[{output},
output = Select[input/.TranslationTable,!FullSimplify[input===#] &];
If[Length[output] > 0, output[[1]],  input/.TranslationTable[[1]]]
];

(* There are two free variables in the multiplication table right now. We can change one of them so that the multiplication is associative. The idea is that we will build the multiplications (z_i*z_j / z0^2)* (z_k/z_0) and (z_j*z_k / z0^2)* (z_i/z_0) artifically and solve for the condition that their quotient is 1 *)

(* https://mathematica.stackexchange.com/questions/197590/how-to-generate-all-the-combinations-with-repetition - for generating triplets of indices *)
AssociativeIndex = With[{n=6,k=3},Join@@Table[IntegerPartitions[s,{k},Range[n]],{s,k,n k}]]


Relations = {};
Do [
(*Print[triplet];*)
i = triplet[[1]];
j = triplet[[2]];
k = triplet[[3]];

LHS = CustomMultiply[CustomMultiply[(Zs[[i+1]]*Zs[[j+1]]/Z0^2)]*Zs[[k+1]]/Z0];
RHS = CustomMultiply[CustomMultiply[(Zs[[j+1]]*Zs[[k+1]]/Z0^2)]*Zs[[i+1]]/Z0];
Relations = Join[Relations, (LHS /RHS)/.WPoints[[20;;30]]];

, {triplet, AssociativeIndex}];

Length[Relations];

(* Each relation is just some constant times a term with the coefficients in As. Here we want to extract the coefficients from each relation and make them exact *)
fn[x_,{var__}]:=List@@Pick[x,x,Alternatives[var]^_.];
RemoveAs = {a1->1, a2->1, a3->1, a4->1, a5->1, a6->1, a7->1};
RelationsNoA = Relations/.RemoveAs;
exactCoeff = RootApproximant[Chop[RelationsNoA, 10^(-500)], 2];

Length[exactCoeff] == Length[Relations];
Equations = {};

Do [
rel = Relations[[i]];

(* We only care about non-trivial relations (some of them are just 1) *)
If [! rel === 1, 
term = Times@@ fn[rel, As];
(*Print[exactCoeff[[i]]*term];*)
AppendTo[Equations, exactCoeff[[i]]*term];
];
, {i, 1, Length[Relations]}];

Now we solve for the associative relation. We also want to add the additional conditions a1 = -a6. This is because we want to take into account the C2 action.

The quotient C2(z_i*z_j) / z_i z_j should be +1 when i,j are both non-zero or both zero and -1 otherwise. The condition that a1 = -a6 ensures this (essentially by looking at the quotient of C2(z4^2) / z1^2, the quotient is -a6/a1. The other quotients don’t add any additional relations)

See below:

C2AutPs = {P10->-P10, P11->-P11, P12->-P12,P13->-P13,P14->-P14,P15->-P15,P16->-P16,P17->-P17,P18->-P18, P19->-P19};
C2AutZs = {Z0->-Z0,
Z1->Z4, Z2->Z5, Z3->Z6, 
Z4->Z1, Z5->Z2, Z6->Z3,
Z7->Z10, Z8->Z11, Z9->Z12,
Z10->Z7, Z11->Z8, Z12->Z9};

Simplify[(MultiplicationTable[[5]][[5]]/.Join[C2AutZs,C2AutPs])/MultiplicationTable[[2]][[2]]]

sols = Solve[{Equations==1, a1==-a6}][[1]]
(* We see that there's a free variable a4 left, which we will set to 1 *)
ans = Join[RootReduce[Simplify[sols]/.a4->1], {a4->1}]

AssociativeMultiTable = MultiplicationTable/.ans;

(* We build a similar custom function for the case of this new multiplication table *)
TranslationTable = {};

Do [
Do [
(* Without the division by Z0^2 should be fine, we are just replacing formal symbols *)
translation = {Zs[[i]]*Zs[[j]]/Z0^2 -> AssociativeMultiTable[[i]][[j]]};
AppendTo[TranslationTable, translation];
, {j, 2, i}];
, {i, 2, 7}];

Length[TranslationTable]

AssocMultiply[input_] := Module[{output},
output = Select[input/.TranslationTable,!FullSimplify[input===#] &];
If[Length[output] > 0, output[[1]],  input/.TranslationTable[[1]]]
];

Now we extend the multiplication to Z0,...,Z12 and solve for cubic relations on them.

ExtendedZs = {Z0, Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8, Z9, Z10, Z11, Z12};
ZWeights = {{0, 0}, 
{1, 1}, {4, 2}, {2, 4},
{-1, 1}, {-4, 2}, {-2, 4},
{1, 3}, {4, 6}, {2, 5},
{-1, 3}, {-4, 6}, {-2, 5}
};

LaterZsToEarly = {Z7->Z1*s2/s3, Z8->Z2*s2a/s3a, Z9->Z3*s2b/s3b, 
Z10->Z4*s1/s3, Z11->Z5*s1a/s3a, Z12->Z6*s1b/s3b};

(* Sanity Check this is C3 invariant *)
Print[Sort[LaterZsToEarly/.Order3Ps/.Join[Order3Zs, {Z7->Z8, Z8->Z9, Z9->Z7, Z10->Z11, Z11->Z12, Z12->Z10}]] === Sort[LaterZsToEarly ]];

(* We solve for the image of Z in CP^12 *)
deg=3;

Powers={b1,b2,b3,b4,b5,b6,b7,b8,b9,b10, b11, b12, b13};
AllPowers=Solve[Join[Table[Powers[[k]]≥0,{k,1,Length[Powers]}],{Sum[Powers[[k]],{k,1,Length[Powers]}]==deg}],Integers];

count = 0;
Equations = {};
Equationsrr = {};

(* We only do this on certain weights, the others can be generated as C3 translates of these equations *)
UsefulWeights = {{0,0},{0,1},{0,3},{1,0},{1,1},{1,2},{1,3},{1,4},{1,5},{1,6},{3,0},{3,1},{3,2},{3,3},{3,4},{3,5},{3,6}};

NewOrder3Zs = Join[Order3Zs, {Z7->Z8, Z8->Z9, Z9->Z7, Z10->Z11, Z11->Z12, Z12->Z10}];


RankTable = ConstantArray[0, {7, 7}];

Do [

mod7A = currentWeight[[1]];
mod7B = currentWeight[[2]];

Print["Weight=",{mod7A, mod7B}];

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

(*Print[Values[CorrectPowers]];*)
Monomials = {};
Do [
(*Creates a list of the variables whose product is the desired monomial *)
bsCount ={b1,b2,b3,b4,b5,b6,b7,b8,b9,b10, b11, b12, b13}/.power;

monList ={};
Do [
If [!bsCount[[i]] === 0,
AppendTo[monList, ConstantArray[ExtendedZs[[i]], bsCount[[i]] ]];
];
, {i, 1, Length[bsCount]}];
monList = Flatten[monList];
(*Print[monList];*)

(* Create the product after multiplication *)
monList = monList/.LaterZsToEarly;
rationalFunc = AssocMultiply[AssocMultiply[monList[[1]]*monList[[2]]/Z0^2]*monList[[3]]/Z0];

AppendTo[Monomials, rationalFunc/.{Z0->1, Z1->1, Z2->1, Z3->1, Z4->1, Z5->1, Z6->1}];

, {power, CorrectPowers}];

(* All Monomials are now in Ps *)
TableOfValues=Table[Monomials[[k]]/.WPoints[[l]],{l,1,30},{k,1,Length[Monomials]}];

(* Standard operations to solve for nullspace *)
(* Print["matrix done"];
Print[Dimensions[TableOfValues]];
Print[MatrixRank[TableOfValues]];
Print[Accuracy[TableOfValues]];*)

result = NullSpace[TableOfValues];
(* Print["nullspace done"];
Print[Dimensions[result]];
Print[Accuracy[result]];*)

(* Here we populate the rank table *)
p2 = Mod[{4*mod7A, 2*mod7B}, 7];
p3 = Mod[{16*mod7A, 4*mod7B}, 7];
RankTable[[mod7A+1]][[mod7B+1]] = Dimensions[result][[1]];
RankTable[[p2[[1]]+1]][[p2[[2]]+1]] = Dimensions[result][[1]];
RankTable[[p3[[1]]+1]][[p3[[2]]+1]] = Dimensions[result][[1]];

If [!Dimensions[result][[1]]===0,

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

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

MonomialsInZs=Table[Product[ExtendedZs[[n]]^(Powers[[n]]/.CorrectPowers[[k]]),{n,1,Length[ExtendedZs]}],{k,1,Length[CorrectPowers]}];

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

AppendTo[Equations, newEquations];
AppendTo[Equationsrr, newEquationsrr];

If [!mod7A===0 || !mod7B=== 0,

AppendTo[Equations, newEquations/.NewOrder3Zs];
AppendTo[Equations, newEquations/.NewOrder3Zs/.NewOrder3Zs];

AppendTo[Equationsrr, newEquationsrr/.NewOrder3Zs];
AppendTo[Equationsrr, newEquationsrr/.NewOrder3Zs/.NewOrder3Zs];

];
count = count + Length[Monomials];
];
, {currentWeight, UsefulWeights}];

Print[Length[Flatten[Equationsrr]]];
(* Showing the number of equations from each weight *)
Style[MatrixForm[RankTable], 20]

fname=OpenWrite["Equations/Z_singular_equations.txt"];
Write[fname,Flatten[Equationsrr]];
Close[fname];

Finally we also want to compute random points on B2/Z with high accuracy. This is fine as long as we avoid singular points of the calculation.

ZSingularEquations = Flatten[Equationsrr/.rr->(Sqrt[-7])];
(*

iter =10;
PointsOnZ={};
For[k=1,k≤iter,k++,
Print["Iteration ", k];

current = NSolve[Join[{Z0==RandomChoice[Join[-Range[100],Range[100]]],Z1==RandomChoice[Join[-Range[100],Range[100]]],Z2==RandomChoice[Join[-Range[100],Range[100]]]},Table[0==ZSingularEquations[[i]],{i,1,Length[ZSingularEquations]}]],1100];

Print[Accuracy[current[[1;;10]]]];

If [Length[current] ≠ 0,
PointsOnZ = Join[PointsOnZ,current[[1;;10]]];
];
];
*)

(* We already pre-computed them, so we will just check here that they lie in the nonsingular locus. *)
filename = "Dependency/points_on_Z.txt";
PointsOnZ=ToExpression[Import[filename]];

{Dimensions[PointsOnZ],Accuracy[PointsOnZ]}

JacobianFPP=Factor[Table[D[ZSingularEquations[[i]],ExtendedZs[[k]]],{i,1,Length[ZSingularEquations]},{k,1,Length[ExtendedZs]}]];
Dimensions[JacobianFPP];

Do [
currentJacob = JacobianFPP/.pt;
If [!MatrixRank[currentJacob] === 10, Print[MatrixRank[currentJacob] ]];
, {pt, PointsOnZ}];