This notebook constructs a surface FPP/C3 with expected \pi_1=C_{13} and finds its singular points.
It is a supporting file for the Section 2 of the joint paper of L.B., A.B., and E.F.

We use the thesis of Jacob Lurie and define the Cartan cubic in 27 variables.

The variables will be lines on the Fermat cubic x0^3+...+x3^3=0.
These lines are the line through (1:a:0:0), (0:0:1:b) with a^3=b^3=-1
and its permutations.

There are 27 collections of triples are given by translates of  (x0+x1+x2+x3=0)
[(1:-1:0:0),(0:0:-1:1)], [(1:0:0:-1),(0:1:-1:0)],[(1:0:-1:0),(0:1:0:-1)].
Specifically, we can multiply the variables by cube roots of 1 but the first factor can be taken 
as 1. There are additional 18 collections which are translates of (x0+x1=0).
(multiply x1 by a cube root, and also 6 pairs of elements).

The 9 triples that give positive signs are x_i + a x_j = 0 with 1<i<j.

(* setting up the variables *)
VarsP={P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14,P15,P16,P17,P18,P19,P20,P21,P22,P23,P24,P25,P26,P27};

(* setting up the lines *)
w=Exp[2 Pi I/3];
Lines={};
For[i=1,i≤3,i++,
For[j=1,j≤3,j++,
AppendTo[Lines,{{1,-w^i,0,0},{0,0,1,-w^j}}];
AppendTo[Lines,{{1,0,-w^i,0},{0,1,0,-w^j}}];
AppendTo[Lines,{{1,0,0,-w^i},{0,1,-w^j,0}}];
]];
Length[Lines]

(* setting up the extra C3 action 
  we tried permuting last three variables, it didn't work
  now just scaling the first coordinate *)
C3action={};
Pairs={};
For[i=1,i≤27,i++,
For[j=1,j≤27,j++,
line1=Lines[[i]];
line2=Lines[[j]];
line2[[1]][[1]]=w  line2[[1]][[1]];
line2[[2]][[1]]=w  line2[[2]][[1]];
If[MatrixRank[Join[line1,line2]]==2,
AppendTo[C3action,VarsP[[i]]->VarsP[[j]]];
AppendTo[Pairs,{i,j}];
];
];];
{C3action,Length[C3action]}


(* setting up planes *)
PlusPlanes={};
For[i=1,i≤3,i++,
AppendTo[PlusPlanes,{0,1,w^i,0}];
AppendTo[PlusPlanes,{0,1,0,w^i}];
AppendTo[PlusPlanes,{0,0,1,w^i}];
];

Length[PlusPlanes]

MinusPlanes={};
For[i=1,i≤3,i++,
AppendTo[MinusPlanes,{1,w^i,0,0}];
AppendTo[MinusPlanes,{1,0,w^i,0}];
AppendTo[MinusPlanes,{1,0,0,w^i}];
];

For[i=1,i≤3,i++,
For[j=1,j≤3,j++,
For[k=1,k≤3,k++,
AppendTo[MinusPlanes,{1,w^i,w^j,w^k}];
]]];
Length[MinusPlanes]

CartanCubic=0;
(* get the polynomial *)
For[i=1,i≤9,i++,
plane=PlusPlanes[[i]];
mon=1;
For[j=1,j≤27,j++,
If[(Simplify[plane.Lines[[j]][[1]]]==0)&&(Simplify[plane.Lines[[j]][[2]]]==0),
mon=mon  VarsP[[j]];
];
];
CartanCubic+=mon;
];

CartanCubic


For[i=1,i≤36,i++,
plane=MinusPlanes[[i]];
mon=1;
For[j=1,j≤27,j++,
If[(plane.Lines[[j]][[1]]==0)&&(plane.Lines[[j]][[2]]==0),
mon=mon VarsP[[j]];
];
];
CartanCubic-=mon;
];

CartanCubic

Eqs27 = Table[D[CartanCubic,VarsP[[i]]],{i,1,27}];

(* we will find the Cartan subgroup of the E_6 that will be used to find an element of order 13 *)
grading={};
For[i=1,i≤Length[MonomialList[CartanCubic]],i++,
mon=MonomialList[CartanCubic][[i]];
AppendTo[grading,Table[Exponent[mon,VarsP[[k]]],{k,1,27}]];
];
degrel=NullSpace[grading]

count=0;
For[i1=0,i1≤12,i1++,
Print[i1];
For[i2=0,i2≤12,i2++,
For[i3=0,i3≤12,i3++,
For[i4=0,i4≤12,i4++,
For[i5=0,i5≤12,i5++,
For[i6=0,i6≤12,i6++,

wt=Mod[degrel[[1]] i1 + degrel[[2]] i2+ degrel[[3]] i3+ degrel[[4]] i4
+ degrel[[5]] i5+ degrel[[6]] i6,13];
If[Sort[wt]=={0,0,0,1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,10,10,11,11,12,12},
isOK=1;
For[m=1,m≤27,m++,
If[Mod[wt[[Pairs[[m]][[1]]]]-9 wt[[Pairs[[m]][[2]]]],13]≠0,
isOK=0;m=28;
];
];
If[isOK==1,
count++;
Print[{i1,i2,i3,i4,i5,i6}];
i1=13;i2=13;i3=13;i4=13;i5=13;i6=13;
];];
]]]]]];
count

deggood=Mod[{0,1,2,1,12,5}.degrel,13]
Sort[deggood]

Eqs27[[9]]
Eqs27[[18]]
Eqs27[[27]]

Below is an affine chart of the OP^2 in some other coordinates, obtained from another spreadsheet by modifying the description in the paper of Gross and Elkies.

VarsQ={Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,Q10,Q11,Q12,Q13,Q14,Q15,Q16,Q17,Q18,Q19,Q20,Q21,Q22,Q23,Q24,Q25,Q26,Q27}
CubicQ=Q10 Q18 Q2-Q11 Q2 Q20-Q16 Q2 Q23-Q14 Q2 Q25-Q11 Q12 Q26-Q10 Q13 Q26-Q15 Q16 Q26-Q14 Q17 Q26-Q19 Q20 Q27-Q18 Q21 Q27-Q23 Q24 Q27-Q22 Q25 Q27+Q1 Q26 Q27-Q12 Q18 Q3-Q13 Q20 Q3+Q16 Q22 Q3-Q14 Q24 Q3-Q10 Q19 Q4-Q11 Q21 Q4-Q17 Q23 Q4+Q15 Q25 Q4-Q1 Q3 Q4-Q12 Q19 Q5-Q1 Q2 Q5+Q13 Q21 Q5-Q17 Q22 Q5-Q15 Q24 Q5-Q17 Q18 Q6-Q16 Q19 Q6+Q11 Q24 Q6-Q13 Q25 Q6+Q17 Q20 Q7-Q16 Q21 Q7-Q10 Q24 Q7-Q12 Q25 Q7-Q15 Q18 Q8+Q14 Q19 Q8-Q11 Q22 Q8-Q13 Q23 Q8-Q1 Q7 Q8-Q15 Q20 Q9-Q14 Q21 Q9-Q10 Q22 Q9+Q12 Q23 Q9-Q1 Q6 Q9;
VarsU={U1,U2,U3,U4,U5,U6,U7,U8,U9,U10,U11,U12,U13,U14,U15,U16}
QtoU={Q1->(U10^2)+(U11^2)+(U12^2)+(U13^2)+(U14^2)+(U15^2)+(U16^2)+(U9^2),Q26->(U1^2)+(U2^2)+(U3^2)+(U4^2)+(U5^2)+(U6^2)+(U7^2)+(U8^2),Q27->1,Q2->U1-I U4,Q5->U1+I U4,Q10->-I U12+U9,Q13->I U12+U9,Q18->-U10 U2-U11 U3-U12 U4-U13 U5-U14 U6-U15 U7-U16 U8+U1 U9-I (-U1 U12-U11 U2+U10 U3-U16 U5+U15 U6-U14 U7+U13 U8-U4 U9),Q21->-U10 U2-U11 U3-U12 U4-U13 U5-U14 U6-U15 U7-U16 U8+U1 U9+I (-U1 U12-U11 U2+U10 U3-U16 U5+U15 U6-U14 U7+U13 U8-U4 U9),Q3->U2-I U3,Q4->U2+I U3,Q11->U10-I U11,Q12->U10+I U11,Q19->-U1 U10-U12 U3+U11 U4-U14 U5+U13 U6+U16 U7-U15 U8-U2 U9-I (-U1 U11+U12 U2-U10 U4-U15 U5-U16 U6+U13 U7+U14 U8-U3 U9),Q20->-U1 U10-U12 U3+U11 U4-U14 U5+U13 U6+U16 U7-U15 U8-U2 U9+I (-U1 U11+U12 U2-U10 U4-U15 U5-U16 U6+U13 U7+U14 U8-U3 U9),Q6->U5-I U8,Q9->U5+I U8,Q14->U13-I U16,Q17->U13+I U16,Q22->-U1 U13+U14 U2+U15 U3+U16 U4-U10 U6-U11 U7-U12 U8-U5 U9-I (-U1 U16+U15 U2-U14 U3-U13 U4+U12 U5+U11 U6-U10 U7-U8 U9),Q25->-U1 U13+U14 U2+U15 U3+U16 U4-U10 U6-U11 U7-U12 U8-U5 U9+I (-U1 U16+U15 U2-U14 U3-U13 U4+U12 U5+U11 U6-U10 U7-U8 U9),Q7->U6-I U7,Q8->U6+I U7,Q15->U14-I U15,Q16->U14+I U15,Q23->-U1 U14-U13 U2+U16 U3-U15 U4+U10 U5+U12 U7-U11 U8-U6 U9-I (-U1 U15-U16 U2-U13 U3+U14 U4+U11 U5-U12 U6+U10 U8-U7 U9),Q24->-U1 U14-U13 U2+U16 U3-U15 U4+U10 U5+U12 U7-U11 U8-U6 U9+I (-U1 U15-U16 U2-U13 U3+U14 U4+U11 U5-U12 U6+U10 U8-U7 U9)};
Expand[Table[D[CubicQ,VarsQ[[i]]]/.QtoU,{i,1,27}]]

(* trying to match *)
ThirdP[a_,b_]:=
Ret=Sum[k Exponent[1+D[D[CartanCubic,VarsP[[a]]],VarsP[[b]]],VarsP[[k]]],{k,1,27}];
ThirdQ[a_,b_]:=
Ret=Sum[k Exponent[1+D[D[CubicQ,VarsQ[[a]]],VarsQ[[b]]],VarsQ[[k]]],{k,1,27}];

TriplesP={};
For[i=1,i≤27,i++,For[j=i,j≤27,j++,For[k=j,k≤27,k++,
If[Coefficient[CartanCubic,VarsP[[i]] VarsP[[j]] VarsP[[k]]]==1,AppendTo[TriplesP,{i,j,k}]];
]]];
TriplesP
TriplesQ={};
For[i=1,i≤27,i++,For[j=i,j≤27,j++,For[k=j,k≤27,k++,
If[Coefficient[CubicQ,VarsQ[[i]] VarsQ[[j]] VarsQ[[k]]]==1,AppendTo[TriplesQ,{i,j,k}]];
]]];
TriplesQ

(* this was found by starting with the {9,10},{18,18},{27,2} and then trial and error *)
Matches={{1,3},{2,12},{3,7},{4,1},{5,5},{6,24},{7,4},{8,19},{9,10},{10,16},{11,23},{12,17},{13,26},{14,13},{15,6},{16,15},{17,8},{18,18},{19,22},{20,9},{21,20},{22,27},{23,21},{24,11},{25,25},{26,14},{27,2}};
For[numsteps=1,numsteps≤27,numsteps++,
lm=Length[Matches];
For[i=1,i≤lm,i++,
For[j=1,j≤i,j++,
nextp=ThirdP[Matches[[i]][[1]],Matches[[j]][[1]]];
nextq=ThirdQ[Matches[[i]][[2]],Matches[[j]][[2]]];
If[nextp>0,
If[nextq==0,Print["Warning"];numsteps=100];
(* check signs as well *)
If[Coefficient[CartanCubic,VarsP[[nextp]]VarsP[[Matches[[i]][[1]]]] VarsP[[Matches[[j]][[1]]]]]≠
Coefficient[CubicQ,VarsQ[[nextq]]VarsQ[[Matches[[i]][[2]]]] VarsQ[[Matches[[j]][[2]]]]],
Print["Warning: signs"];numsteps=100];

AppendTo[Matches,{nextp,nextq}];
];
]];
Matches=DeleteDuplicates[Matches];
];
Sort[Matches]
Length[Matches]

(* checking substitutions *)
QtoP=Table[VarsQ[[Matches[[i]][[2]]]]->VarsP[[Matches[[i]][[1]]]],{i,1,27}];
PtoQ=Table[VarsP[[Matches[[i]][[1]]]]-> VarsQ[[Matches[[i]][[2]]]],{i,1,27}];
(CubicQ/.QtoP)-CartanCubic
(CartanCubic/.PtoQ)-CubicQ

We will now write down the affine chart on OP^2 in P-variables.

PtoU=PtoQ/.QtoU
Expand[Eqs27/.PtoU] (* yes indeed *)
(* Let us now try to change to variables V with more torus symmetry of the map *)
VarsV={V1,V2,V3,V4,V5,V6,V7,V8,V9,V10,V11,V12,V13,V14,V15,V16};
UtoV=
{U2->(V2+V3)/2,U3->I(V2-V3)/2,
U10->(V10+V11)/2,U11->I(V10-V11)/2,
U6->(V6+V7)/2,U7->I(V6-V7)/2,
U1->(V1+V4)/2,U4->I(V1-V4)/2,
U5->(V5+V8)/2,U8->I(V5-V8)/2,
U9->(V9+V12)/2,U12->I(V9-V12)/2,
U14->(V14+V15)/2,U15->I(V14-V15)/2,
U13->(V13+V16)/2,U16->I(V13-V16)/2
};
PtoV=Simplify[PtoU/.UtoV]
Expand[Eqs27/.PtoV](* yes indeed *)

(* this is one way of seeing the same chart *)
Solve[Table[0==Eqs27[[k]],{k,1,27}],Complement[VarsP,{P22,P1,P2,P3,P5,P7,P9,P10,P12,P14,P15,P16,P17,P20,P24,P26,P27}]]

We will set up a C13 equivariant cut which, for some strange reason, will also admit a C3 action.

Sort[Table[{deggood[[i]],i},{i,1,27}]]
CoeffToFind={d1,d2,d3,d4,d5,d6,d7,d8,d9,d10,d11,d12};
Lineqs={P9,P18,P27,P23+d1 P26,P19 +d2 P24, P5+d3 P8, P13+d4 P16, P10+d5 P15,P1 + d6 P6, P2 + d7 P3, P11+d8 P12, P14+ d9 P17, P4+ d10 P7, P20+ d11 P21, P22+d12 P25};
FullSimplify[(Eqs27[[9]]+Eqs27[[18]]+Eqs27[[27]])/.Solve[Table[0==Lineqs[[k]],{k,1,Length[Lineqs]}],VarsP][[1]]]/.{d12->-1/d1,d11->-1/d2,d10->-1/d3,d9->-1/d4,d8->-1/d5,d7->-1/d6};
Solve[Table[0==Lineqs[[i]],{i,1,15}],{P1,P10,P12,P13,P17,P18,P19,P3,P21,P25,P23,P27,P7,P5,P9}][[1]]/.{d12->-1/d1,d11->-1/d2,d10->-1/d3,d9->-1/d4,d8->-1/d5,d7->-1/d6}
EqsCut=(Eqs27/.Solve[Table[0==Lineqs[[i]],{i,1,15}],{P1,P10,P12,P13,P17,P18,P19,P3,P21,P25,P23,P27,P7,P5,P9}][[1]]/.{d12->-1/d1,d11->-1/d2,d10->-1/d3,d9->-1/d4,d8->-1/d5,d7->-1/d6})
EqsAffine=(Lineqs/.{d12->-1/d1,d11->-1/d2,d10->-1/d3,d9->-1/d4,d8->-1/d5,d7->-1/d6})/.PtoV;
PtoT=Sort[{P6->t6,P14->t9,P24->t2,P26->t1,P4->t10,P15-> t5,P20->t11,P22->t12,P8->t3,P16->t4,P2->t7,P11->t8}];
EqsT=Table[0==EqsCut[[k]]/.PtoT,{k,1,26}]/.{d12->-1/d1,d11->-1/d2,d10->-1/d3,d9->-1/d4,d8->-1/d5,d7->-1/d6};
VarsT={t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12};
(* We will now set all but d1,d2 to 1, which we can do by torus symmetry *)
EqsAffine=EqsAffine/.{d3->1,d4->1,d5->1,d6->1}
EqsT=EqsT/.{d3->1,d4->1,d5->1,d6->1}

Length[EqsAffine];
EqsAffine1=Simplify[EqsAffine/.d1->1];
EqsSome=Complement[EqsAffine1/.{V9->0,V1->0},{0}]
Length[EqsSome];
sls=Solve[{0==EqsSome[[2]],0== EqsSome[[4]],0==EqsSome[[6]]},{V15,V11,V7}][[1]];
EqsSome=Complement[EqsSome/.sls,{0}]
Length[EqsSome]
NewVarsV=Complement[VarsV,{V1,V9,V7,V15,V11}]
(* Now let us try to set up singularity condition, as a specific kind of relation on the equations at a point *)
VarsAll=Join[NewVarsV,{a1,a2,a3,a4,a5,a6,a7,a8,d2}];
AllEqsSing=EqsSome;
For[k=1,k≤11,k++,
AppendTo[AllEqsSing,Sum[({1,1,a1,a2,a3,a4,a5,a6,a7,a8}[[j]])D[EqsSome[[j]],NewVarsV[[k]]],{j,1,10}]];
];
ByteCount[AllEqsSing]
{Length[AllEqsSing],Length[VarsAll]}

(* We will try to randomly look around in the space of solutions using the FindRoot command. 
This may need to be rerun if not lucky the first time around. *)
flag=0;
NumTries=100;
For[m=1,m≤NumTries,m++,
fi=Table[VarsAll[[k]]->RandomReal[1]+I RandomReal[1],{k,1,20}];
If[Mod[m,20]==0,Print[m]];
RemoveIndex=RandomInteger[20]+1;
Eqs20=Complement[AllEqsSing,{AllEqsSing[[RemoveIndex]]}];
fi=FindRoot[Table[0==Eqs20[[k]],{k,1,20}],Table[{VarsAll[[n]],(VarsAll[[n]]/.fi)},{n,1,20}],AccuracyGoal->500,PrecisionGoal->500,MaxIterations->200,WorkingPrecision->1000];
JacMat=Table[D[EqsSome[[i]],NewVarsV[[j]]]/.fi,{i,1,10},{j,1,11}];
If[MatrixRank[JacMat]==8,m=NumTries+1;flag=1;];
];
If[flag==1,d2good=RootApproximant[d2/.fi,48]];
d2good

It appears that there are d2 that are roots of some degree 7 polynomial that fit the bill. We will use the starting data to get to a nearby point.

figood=fi;

Length[EqsAffine];
EqsAffine1=Simplify[EqsAffine];
EqsSome=Complement[EqsAffine1/.{V9->0,V1->0},{0}];
Length[EqsSome];
sls=Solve[{0==EqsSome[[2]],0== EqsSome[[4]],0==EqsSome[[6]]},{V15,V11,V7}][[1]];
EqsSome=Complement[Simplify[EqsSome/.sls],{0}];
Length[EqsSome];
NewVarsV=Complement[VarsV,{V1,V9,V7,V15,V11}];
VarsAll=Join[NewVarsV,{a1,a2,a3,a4,a5,a6,a7,a8,d2}];
AllEqsSing=Simplify[EqsSome];
For[k=1,k≤11,k++,
AppendTo[AllEqsSing,Sum[({1,1,a1,a2,a3,a4,a5,a6,a7,a8}[[j]])D[EqsSome[[j]],NewVarsV[[k]]],{j,1,10}]];
];
ByteCount[AllEqsSing]
{Length[AllEqsSing],Length[VarsAll]}

PtsSingCuts={};
For[numpt=0,numpt≤1,numpt++,
AllEqsSingc1=AllEqsSing/.d1->(1+numpt/10^20);
(* we will use the same RemoveIndex as before *)
Eqs20=Complement[AllEqsSingc1,{AllEqsSing[[RemoveIndex]]}];
fi=FindRoot[Table[0==Eqs20[[k]],{k,1,20}],Table[{VarsAll[[n]],(VarsAll[[n]]/.figood)},{n,1,20}],AccuracyGoal->50000,PrecisionGoal->50000,MaxIterations->300,WorkingPrecision->100000];
JacMat=Table[D[EqsSome[[i]]/.d1->1+numpt/10^20,NewVarsV[[j]]]/.fi,{i,1,10},{j,1,11}];
If[MatrixRank[JacMat]==8,AppendTo[PtsSingCuts,{1+numpt/10^20,d2/.fi}]];
];

(* finding minimal polynomial *)
mp=MinimalPolynomial[RootApproximant[N[PtsSingCuts[[2]][[2]],10000],10],t]

(* we want to see the above coefficients as polynomials in 10^20 with small coefficients *)
Poly=0;
For[i=0,i≤8,i++,
cf=SeriesCoefficient[mp,{t,0,i}];
nl=NullSpace[{Join[Table[10^(20 k),{k,0,24}],{cf}]}];
nl=LatticeReduce[nl][[1]];
nl=Table[-nl[[k]]/nl[[26]],{k,1,26}];
Poly+=Sum[t^i r^k nl[[k+1]],{k,0,24}];
];
(Poly/.r->10^20) - mp

LocusOfSingularCuts=Numerator[Factor[Poly/.{t->d2,r->1/(d1-1)}]]
(* now we find the potential A2 singularities *)
LocusA2=Solve[{LocusOfSingularCuts==0,D[LocusOfSingularCuts,d1]==0,D[LocusOfSingularCuts,d2]==0}];
Length[DeleteDuplicates[RootReduce[{d1,d2}/.LocusA2]]]
LocusA2values=Sort[RootReduce[{d1,d2}/.LocusA2]]

Length[DeleteDuplicates[MinimalPolynomial[LocusA2values]]]
(* There seem to be 4 different cases: (0,0), \pm (i,i) and two size 8 orbits of the Galois group. 
And (0,0) is not viable, so really just three cases. *)
LocusA2valuesReduced={LocusA2values[[3]],LocusA2values[[4]],LocusA2values[[19]]}

RootReduce[(D[D[LocusOfSingularCuts,d1],d1]D[D[LocusOfSingularCuts,d2],d2]-D[D[LocusOfSingularCuts,d1],d2]^2)/.{d1-> LocusA2valuesReduced[[2]][[1]],d2-> LocusA2valuesReduced[[2]][[2]]}]
RootReduce[(D[D[LocusOfSingularCuts,d1],d1]D[D[LocusOfSingularCuts,d2],d2]-D[D[LocusOfSingularCuts,d1],d2]^2)/.{d1-> LocusA2valuesReduced[[3]][[1]],d2-> LocusA2valuesReduced[[3]][[2]]}]
(* These results show that the two cases are different: we guess is that the 2nd one give 78 nodes while the third one gives 39 A_2
singularities *)

Just playing around with the field extensions.

(* Let us inderstand the last solution *)
gal=ToNumberField[LocusA2valuesReduced[[3]][[2]],LocusA2valuesReduced[[3]][[1]]][[2]]

roots=t/.Solve[MinimalPolynomial[LocusA2valuesReduced[[3]][[1]],t]==0]
perm={};
For[i=1,i≤8,i++,
For[j=1,j≤8,j++,
(* For some reason it needs to have powers of 27 *)
If[RootReduce[Sum[gal[[k]] roots[[i]]^(k-1)(27)^(k-1),{k,1,8}]]==roots[[j]],
AppendTo[perm,{i,j}];
];
]];
perm
perm={};
For[i=1,i≤8,i++,
For[j=1,j≤8,j++,
(* For some reason it needs to have powers of 27 *)
If[RootReduce[-1/roots[[i]]]==roots[[j]],
AppendTo[perm,{i,j}];
];
]];
perm
RootReduce[roots[[8]]+roots[[5]]+roots[[2]]+roots[[3]]]
RootReduce[roots[[8]]^2+roots[[5]]^2+roots[[2]]^2+roots[[3]]^2]

(* Can we have (d1,d2)->(-1/d2,d1) ? This would send (8,5)->(3,8)->(2,3)->(5,2)->(8,5). We will use it.
*)

Before continuing further, let us look at the C3 action on these surfaces. The action looks like multiplication of an index by 3 mod 13,
however, there may also be some scaling. We will experiment with it to make it work. We will be scaling variables to make a good match. 
We were successful in it, the new relations are called NewRels.

Unset[{r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12}];
RelsT=(Table[EqsCut[[k]]/.PtoT,{k,1,26}]/.{d12->-1/d1,d11->-1/d2,d10->-1/d3,d9->-1/d4,d8->-1/d5,d7->-1/d6})/.{d3->1,d4->1,d5->1,d6->1}
ActionC3={t1->t3,t3->t9,t9->t1,t2->t6,t6->t5,t5->t2 ,t4->t12,t12->t10,t10->t4,t7->t8,t8->t11,t11->t7}
ScaleT={t1->t1 r1,t2->t2 r2,t3->t3 r3,t4->t4 r4,t5->t5 r5,t6->t6 r6,t7->t7 r7,t8->t8 r8,t9->t9 r9,t10->t10 r10,t11->t11 r11,t12->t12 r12};

AllMonsT={};For[i=1,i≤12,i++,For[j=1,j≤i,j++,AppendTo[AllMonsT,VarsT[[i]]VarsT[[j]]]]];

RelsT=RelsT/.ScaleT;
NonZeroRelsT=Complement[RelsT,{RelsT[[9]],RelsT[[18]]}];
Length[NonZeroRelsT]
EqsSym={};
For[i=1,i≤24,i++,
NewRel={NonZeroRelsT[[i]]/.ActionC3};
grad13=Mod[Sum[Exponent[MonomialList[NewRel[[1]]][[1]],VarsT[[j]]] j ,{j,1,12}],13];
For[n=1,n≤24,n++,
If[Mod[Sum[Exponent[MonomialList[NonZeroRelsT[[n]]][[1]],VarsT[[j]]] j ,{j,1,12}],13]==grad13,
AppendTo[NewRel,NonZeroRelsT[[n]]];
];
];
MatMon=Table[Coefficient[NewRel[[k]],AllMonsT[[l]]],{k,1,3},{l,1,Length[AllMonsT]}];
EqsSym=Union[EqsSym,DeleteDuplicates[Flatten[Minors[MatMon,3]]]];
];

Length[EqsSym]

VarsR={r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12};
Solve[Join[Table[0==EqsSym[[20 k]],{k,1,17}],Table[VarsR[[i]]≠0,{i,1,12}]],VarsR];

DeleteDuplicates[EqsSym/.{r9->(((d2^(1/13)) r3)/(d1^(2/13))),r1->(((d2^(4/13)) r3)/(d1^(8/13))),r5->-(((d2^(7/13)) r2)/(d1^(1/13))),r10->(d1^(5/13)) (d2^(4/13)) r12,r6->(d1^(3/13)) (d2^(5/13)) r2,r4->(d1^(7/13)) (d2^(3/13)) r12,r7->-(((d2^(2/13)) r8)/(d1^(4/13))),r11->-(r8/((d1^(1/13)) (d2^(6/13))))}]

Simplify[Numerator[Factor[RelsT/.{r9->(((d2^(1/13)) r3)/(d1^(2/13))),r1->(((d2^(4/13)) r3)/(d1^(8/13))),r5->-(((d2^(7/13)) r2)/(d1^(1/13))),r10->(d1^(5/13)) (d2^(4/13)) r12,r6->(d1^(3/13)) (d2^(5/13)) r2,r4->(d1^(7/13)) (d2^(3/13)) r12,r7->-(((d2^(2/13)) r8)/(d1^(4/13))),r11->-(r8/((d1^(1/13)) (d2^(6/13))))}/.r2-> r8 d1^(2/13)d2^(-1/13)/.r3-> r8  (d1^(6/13)) (d2^(10/13))/.r12->r8 (d1^(-10/13)) (d2^(-8/13)) /.r8->1 ]]]

NewRels={-(t10^2)-d1 (d2^2) (t2 t5+t1 t6-t11 t9),d1 (d2^2) (t3^2)+t2 t4+t12 t7-t11 t8,d1 (d2^2) t1 t5+t12 t7-d2 (t11 t8+t10 t9),-t12 t4+d1 d2 (-t11 t5+t10 t6+d2 t7 t9),t4 t6+d2 (-t3 t7-t2 t8+d1 d2 t1 t9),d1 d2 (t3 t4-t2 t5+d2 t1 t6)-t12 t8,d1 (d2^2) t1 t2+t12 t4+d2 t10 t6-d2 (t8^2),t11 t12+d1 d2 ((t5^2)-t3 t7+d2 t1 t9),-t11 t2-t10 t3+t6 t7+t4 t9,(t4^2)+d1 (d2^2) (t3 t5+t2 t6-t1 t7),-t12 t6+t11 t7-t10 t8-d1 (d2^2) (t9^2),-d1 (d2^2) t2 t3+d2 t1 t4+d2 t11 t7-t10 t8,t10 t12-d1 d2 (t4 t5-t2 t7+d2 t1 t8),d1 (d2^2) t1 t3+t12 t5-d2 (t11 t6+t8 t9),t10 t11-d1 d2 (d2 t3 t5-t2 t6+t12 t9),-d2 (t11^2)+t10 t12+d2 t4 t5+d1 (d2^2) t3 t6,t10 t7+d1 d2 ((t2^2)+d2 t1 t3-t8 t9),t1 t12-t6 t7+t5 t8-t4 t9,-(t12^2)-d1 (d2^2) (t5 t6-t3 t8+t2 t9),d1 (d2^2) (t1^2)+t11 t4+t10 t5-t7 t8,t11 t4-d2 (t12 t3+t7 t8)+d1 (d2^2) t6 t9,-t10 t4+d1 d2 (t12 t2+d2 t11 t3-t6 t8),d2 t1 t11-t10 t2+d2 t5 t7-d1 (d2^2) t3 t9,-t4 t7+d1 d2 (t1 t10-t5 t6+d2 t2 t9),d2 t12 t2+t10 t4-d2 (t7^2)+d1 (d2^2) t5 t9,t4 t8+d1 d2 (-t1 t11+(t6^2)+d2 t3 t9)};
(* checking that the action is OK for the C13 invariant equations *)
{(NewRels[[9]]/.ActionC3)-NewRels[[18]],(NewRels[[18]]/.ActionC3) + NewRels[[9]]+NewRels[[18]]}

Now we want to understand the symmetry of i -> 13-i for the indices. There may be a substitution that sends d1->-1/d2, d2->-1/d1 or something like that.
Some trial and error gives the following, which works.

NewRelsTwist=Numerator[Factor[NewRels/.{t1->t12,t2->-t11,t3->t10,t4->-t9,t5->-t8,t6->-t7,t7->t6,t8->t5,t9->t4,t10->-t3,t11->t2,t12->-t1}/.{d1->-1/d1,d2->-1/d2}]];
count=0;
For[i=1,i≤26,i++,
For[j=1,j≤26,j++,
If[Simplify[NewRels[[i]]+NewRelsTwist[[j]]]==0,count++];
If[Simplify[NewRels[[i]]-NewRelsTwist[[j]]]==0,count++];
]];
count

Let us now try to understand i -> Mod[5 i ,13].

(* The idea is to match one equation and have the square of the transformation match the previous one. *)
VarsT/.{t1->-t12,t2->t11,t3->t10,t4->-t9,t5->-t8,t6->-t7,t7->t6,t8->t5,t9->t4,t10->-t3,t11->t2,t12->-t1}/.{d1->-1/d1,d2->-1/d2}
ChangeVars={t1->t5  d1^a1 d2^b1,t2->t10 d1^a2 d2^b2,t3-> t2 d1^a3 d2^b3  ,t4->t7  d1^a4 d2^b4 ,t5->t12  d1^a5 d2^b5 ,t6->t4  d1^a6 d2^b6 ,t7->t9  d1^a7 d2^b7,t8->t1  d1^a8 d2^b8,t9->t6 d1^a9 d2^b9 ,t10-> t11  d1^a10 d2^b10 ,t11->t3 d1^a11 d2^b11,t12->t8  d1^a12 d2^b12,d2->d1,d1->-1/d2};
Tchch=Expand[VarsT/.ChangeVars/.ChangeVars]
EqsToSolve=Join[Table[Exponent[Tchch[[k]],d1]==Exponent[Tchch[[1]],d1],{k,1,12}],Table[Exponent[Tchch[[k]],d2]==Exponent[Tchch[[1]],d2],{k,1,12}]]

neweq=Numerator[Factor[NewRels[[17]]/.ChangeVars]];
cmp1=Expand[-neweq/Coefficient[neweq,t10^2]]
cmp2=Expand[NewRels[[1]]]
ListMons={t2 t5,t1 t6, t11 t9};
For[i=1,i≤3,i++,AppendTo[EqsToSolve,Exponent[Coefficient[cmp1,ListMons[[i]]],d1]==Exponent[Coefficient[cmp2,ListMons[[i]]],d1]]];
For[i=1,i≤3,i++,AppendTo[EqsToSolve,Exponent[Coefficient[cmp1,ListMons[[i]]],d2]==Exponent[Coefficient[cmp2,ListMons[[i]]],d2]]];

neweq=Numerator[Factor[NewRels[[16]]/.ChangeVars]];
cmp1=Expand[d1 (d2^2)neweq/Coefficient[neweq,t3^2]]
cmp2=Expand[NewRels[[2]]]
ListMons={t2 t4,t12 t7, t11 t8};
For[i=1,i≤3,i++,AppendTo[EqsToSolve,Exponent[Coefficient[cmp1,ListMons[[i]]],d1]==Exponent[Coefficient[cmp2,ListMons[[i]]],d1]]];
For[i=1,i≤3,i++,AppendTo[EqsToSolve,Exponent[Coefficient[cmp1,ListMons[[i]]],d2]==Exponent[Coefficient[cmp2,ListMons[[i]]],d2]]];

neweq=Numerator[Factor[NewRels[[19]]/.ChangeVars]];
cmp1=Expand[-d2 neweq/Coefficient[neweq,t8^2]]
cmp2=Expand[NewRels[[7]]]
ListMons={t1 t2,t12 t4, t10 t6};
For[i=1,i≤3,i++,AppendTo[EqsToSolve,Exponent[Coefficient[cmp1,ListMons[[i]]],d1]==Exponent[Coefficient[cmp2,ListMons[[i]]],d1]]];
For[i=1,i≤3,i++,AppendTo[EqsToSolve,Exponent[Coefficient[cmp1,ListMons[[i]]],d2]==Exponent[Coefficient[cmp2,ListMons[[i]]],d2]]];

sls=Solve[EqsToSolve][[1]]
Simplify[Tchch/.sls/.{a1->0,b1->0}]

(* Now let's match signs *)
ChangeVars={t1->t5  d1^a1 d2^b1,t2->t10 d1^a2 d2^b2,t3-> t2 d1^a3 d2^b3  ,t4->t7  d1^a4 d2^b4 ,t5->t12  d1^a5 d2^b5 ,t6->t4  d1^a6 d2^b6 ,t7->t9  d1^a7 d2^b7,t8->t1  d1^a8 d2^b8,t9->t6 d1^a9 d2^b9 ,t10-> t11  d1^a10 d2^b10 ,t11->t3 d1^a11 d2^b11,t12->t8  d1^a12 d2^b12,d2->d1,d1->-1/d2}/.sls/.{a1->0,b1->1}
VarsT/.{t1->t12,t2->-t11,t3->t10,t4->-t9,t5->-t8,t6->-t7,t7->t6,t8->t5,t9->t4,t10->-t3,t11->t2,t12->-t1}/.{d1->-1/d1,d2->-1/d2}

ChangeVars={t1-> t5,t2->(t10/d2),t3->-t2,t4->((d1 t7)/d2),t5->-(t12/d2),t6->(t4/d2),t7->d1 t9,t8->  d1 t1,t9->t6,t10->((d1 t11)/d2),t11->d1 t3,t12->((d1 t8)/d2),d2->d1,d1->-(1/d2)};
ChangeVars={t1->d2 t5 c1,t2->t10 c2,t3->d2 t2 c3,t4->d1 t7 c4,t5->t12 c5 ,t6->t4 c6,t7-> d1 d2 t9 c7,t8->d1 d2 t1 c8,t9->d2 t6 c9,t10->d1 t11 c10,t11->d1 d2 t3 c11,t12->d1 t8 c12,d2->d1,d1->-(1/d2)}/.{c1->c2,c10->-c2,c11->-c2,c12->-c2,c3->c2,c4->-c2,c5->c2,c6->c2,c7->-c2,c8->-c2,c9->c2}/.c2->1
Tchch=(VarsT/.ChangeVars/.ChangeVars)

EqsToSolve={}

neweq=Numerator[Factor[NewRels[[17]]/.ChangeVars]];
cmp1=Expand[-neweq/Coefficient[neweq,t10^2]]
cmp2=Expand[NewRels[[1]]]
ListMons={t2 t5,t1 t6, t11 t9};
For[i=1,i≤3,i++,AppendTo[EqsToSolve,Coefficient[cmp1,ListMons[[i]]]==Coefficient[cmp2,ListMons[[i]]]]];

neweq=Numerator[Factor[NewRels[[16]]/.ChangeVars]];
cmp1=Expand[d1 (d2^2)neweq/Coefficient[neweq,t3^2]]
cmp2=Expand[NewRels[[2]]]
ListMons={t2 t4,t12 t7, t11 t8};
For[i=1,i≤3,i++,AppendTo[EqsToSolve,Coefficient[cmp1,ListMons[[i]]]==Coefficient[cmp2,ListMons[[i]]]]];

neweq=Numerator[Factor[NewRels[[19]]/.ChangeVars]];
cmp1=Expand[-d2 neweq/Coefficient[neweq,t8^2]]
cmp2=Expand[NewRels[[7]]]
ListMons={t1 t2,t12 t4, t10 t6};
For[i=1,i≤3,i++,AppendTo[EqsToSolve,Coefficient[cmp1,ListMons[[i]]]==Coefficient[cmp2,ListMons[[i]]]]];

neweq=Numerator[Factor[NewRels[[20]]/.ChangeVars]];
cmp1=Expand[d1 d2 neweq/Coefficient[neweq,t5^2]]
cmp2=Expand[NewRels[[8]]]
ListMons={t11 t12,t3 t7, t1 t9};
For[i=1,i≤3,i++,AppendTo[EqsToSolve,Coefficient[cmp1,ListMons[[i]]]==Coefficient[cmp2,ListMons[[i]]]]];


EqsToSolve=EqsToSolve/.{d1->1,d2->1}
Solve[EqsToSolve]

neweq=Numerator[Factor[NewRels[[17]]/.ChangeVars]];
cmp1=Expand[-neweq/Coefficient[neweq,t10^2]];
cmp2=Expand[NewRels[[1]]];
cmp1-cmp2

neweq=Numerator[Factor[NewRels[[16]]/.ChangeVars]];
cmp1=Expand[d1 (d2^2)neweq/Coefficient[neweq,t3^2]];
cmp2=Expand[NewRels[[2]]];
cmp1-cmp2

neweq=Numerator[Factor[NewRels[[19]]/.ChangeVars]];
cmp1=Expand[-d2 neweq/Coefficient[neweq,t8^2]];
cmp2=Expand[NewRels[[7]]];
cmp1-cmp2

neweq=Numerator[Factor[NewRels[[20]]/.ChangeVars]];
cmp1=Expand[d1 d2 neweq/Coefficient[neweq,t5^2]];
cmp2=Expand[NewRels[[8]]];
cmp1-cmp2

(* Now making sure it works *)
ChangeVars={t1->d2 t5,t2->t10,t3->d2 t2,t4->-d1 t7,t5->t12,t6->t4,t7->-d1 d2 t9,t8->-d1 d2 t1,t9->d2 t6,t10->-d1 t11,t11->-d1 d2 t3,t12->-d1 t8,d2->d1,d1->-(1/d2)};
Tchch=(VarsT/.ChangeVars/.ChangeVars)/d1
count=0;
For[i=1,i≤26,i++,
For[j=1,j≤26,j++,
cmp1=Numerator[Factor[NewRels[[j]]/.ChangeVars]];
cmp2=NewRels[[i]];
If[ByteCount[Simplify[cmp1/cmp2]]<200,Print[{i,j, Simplify[cmp1/cmp2]}];
count++];
]];
count
(* it works *)

ChangeVars={t1->d2 t5,t2->t10,t3->d2 t2,t4->-d1 t7,t5->t12,t6->t4,t7->-d1 d2 t9,t8->-d1 d2 t1,t9->d2 t6,t10->-d1 t11,t11->-d1 d2 t3,t12->-d1 t8,d2->d1,d1->-(1/d2)};
VarsT/.ChangeVars/.ChangeVars
VarsT/.ChangeVars/.ChangeVars/.ChangeVars
Join[VarsT,{d1,d2}]/.ChangeVars/.ChangeVars/.ChangeVars/.ChangeVars

We will now construct quotients by C13. The (i,i) was found to have a curve of singular points (by a separate Magma calculation), so we can ignore it. We will focus on the third point. We will look at degree 2 and degree 3 invariants.
It is important to apply the symmetry to get the smaller field of definition.

VarsT={t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12};
NumCut=3;
EqsCover=NewRels;
GB=GroebnerBasis[EqsCover,VarsT,MonomialOrder->DegreeReverseLexicographic,CoefficientDomain->RationalFunctions];
VarsS={S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12,S13,S14};
degS={2,2,2,2,3,3,3,3,3,3,3,3,3,3};

StoT={S1->t1 t12 + t3 t10 + t9  t4, S2-> t2 t11, S3->t6 t7, S4-> t5 t8,
S5->t1 t2 t10, S6 -> t3 t6 t4, S7->t9 t5 t12, S8->t1 t4 t8,S9->t3 t12 t11, S10->t9 t10 t7,S11->t1 t5 t7,S12->t3 t2 t8, S13->t9 t6 t11,
S14->t1 t3 t9};

ChangeVars={t1->d2 t5,t2->t10,t3->d2 t2,t4->-d1 t7,t5->t12,t6->t4,t7->-d1 d2 t9,t8->-d1 d2 t1,t9->d2 t6,t10->-d1 t11,t11->-d1 d2 t3,t12->-d1 t8,d2->d1,d1->-(1/d2)};
(* symmetrize *)
StoT=Simplify[Table[VarsS[[k]]->(VarsS[[k]]/.StoT)+(VarsS[[k]]/.StoT/.ChangeVars)+(VarsS[[k]]/.StoT/.ChangeVars/.ChangeVars)+(VarsS[[k]]/.StoT/.ChangeVars/.ChangeVars/.ChangeVars),{k,1,Length[VarsS]}]];

EqsFPPC3withrr={};

VarsR={R1,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11,R12,R13,R14};

For[deg=5,deg≤6,deg++,

solvemons=Solve[Append[Table[VarsR[[i]]≥ 0,{i,1,14}],Sum[VarsR[[i]] degS[[i]],{i,1,14}]==deg],VarsR,Integers];
AllMons=Table[Product[VarsS[[i]]^(VarsR[[i]]/.solvemons[[k]]),{i,1,14}],{k,1,Length[solvemons]}];

AllMonsT={};
AllRem={};
For[i=1,i≤Length[AllMons],i++,
If[Mod[i,20]==0,Print[i];];
AppendTo[AllRem,PolynomialReduce[AllMons[[i]]/.StoT,GB,VarsT,MonomialOrder->DegreeReverseLexicographic][[2]] ];
newrem=AllRem[[i]];
monlist=MonomialList[newrem/.{d1->LocusA2valuesReduced[[NumCut]][[1]],d2->LocusA2valuesReduced[[NumCut]][[2]]}];
For[k=1,k≤Length[monlist],k++,
monlist[[k]]=Product[VarsT[[j]]^Exponent[monlist[[k]],VarsT[[j]]],{j,1,12}];
];
AllMonsT=Union[AllMonsT,monlist];
];
RemMatrix=Table[Coefficient[AllRem[[i]],AllMonsT[[j]]],{j,1,Length[AllMonsT]},{i,1,Length[AllMons]}];
RtRm=RootReduce[RemMatrix/.{d1-> LocusA2valuesReduced[[3]][[1]],d2-> LocusA2valuesReduced[[3]][[2]]}];
blah=RowReduce[NullSpace[N[RtRm,10000]]];
Print[Accuracy[blah]];
goodnull=RootApproximant[Chop[blah],3];
For[i=1,i≤Length[goodnull],i++,For[j=1,j≤Length[goodnull[[i]]],j++,
If[Exponent[MinimalPolynomial[goodnull[[i]][[j]],t],t]==3,Print["Warning: not enough accuracy"]];
]];
goodnullwithrr=Table[Rationalize[Re[goodnull[[i]][[j]]]]+Rationalize[Im[goodnull[[i]][[j]]]/(2^(1/2))]  rr,{i,1,Length[goodnull]},{j,1,Length[goodnull[[1]]]}];
Print[DeleteDuplicates[Simplify[(goodnullwithrr/.rr-> I (Sqrt[2]))-goodnull]]];(* should be zeros *)
EqsFPPC3withrr=Join[EqsFPPC3withrr,goodnullwithrr.AllMons];
EqsFPPC3withrr=Table[Numerator[Factor[EqsFPPC3withrr[[i]]]],{i,1,Length[EqsFPPC3withrr]}];
Print[Table[ByteCount[EqsFPPC3withrr[[i]]],{i,1,Length[EqsFPPC3withrr]}]];

];

(* The name of the file needs to be adjusted according to your situation. *)
fname="Desktop/FPP2020/EqsFPPC3withrr";
OpenWrite[fname];
Write[fname,{VarsS,EqsFPPC3withrr}];
Close[fname];

The process of finding the A2 singular points was rather difficult. Specifically, we computed a degree 12 polynomial equation in S1,...,S4, then found singular points on hyperplane cuts, then excluded them to find isolated singular points. If one doesn’t want to repeat these calculation, it is a good idea to scroll all the way down.

(* Let us find points on this surface *)
pts={};
Binomial[3+12,3]
For[i=1,i≤300,i++,
If[Mod[i,20]==5,Print[i];];
sls=NSolve[Join[Table[0==EqsFPPC3withrr[[k]]/.rr->I (Sqrt[2]),{k,1,38}],{1==S1,S2==i , S4==RandomInteger[1000]/1000}],500];
AppendTo[pts,{S1,S2,S3,S4}/.sls[[1]]];
AppendTo[pts,{S1,S2,S3,S4}/.sls[[7]]];
];

AllMons12={};
For[i1=0,i1≤12,i1++,For[i2=0,i2≤12-i1,i2++,For[i3=0,i3≤12-i2-i1,i3++,
i4=12-i1-i2-i3;
AppendTo[AllMons12,S1^i1 S2^i2 S3^i3 S4^i4];
]]];
Length[AllMons12]

ns12=NullSpace[Table[AllMons12[[i]]/.Table[VarsS[[k]]->pts[[j]][[k]],{k,1,4}],{j,1,Length[pts]},{i,1,455}]];

Accuracy[ns12]
Length[ns12]
ns12=RowReduce[ns12];

Eq12=Numerator[Factor[Block[{$MaxExtraPrecision=400},RootApproximant[Chop[ns12],2][[1]].AllMons12]]];

(* The name of the file needs to be adjusted according to your situation. *)
fname="Desktop/FPP2020/Eq12";
OpenWrite[fname];
Write[fname,Numerator[Factor[Eq12]]];
Close[fname];

(* finding singularities on curves *)
ptssing={};
Eq12=Numerator[Factor[Eq12]];
For[mm=1,mm≤60,mm++,
If[Mod[mm,20]==0,Print[mm]];

sls=NSolve[{D[Eq12,S1]==0,D[Eq12,S2]==0,D[Eq12,S3]==0,D[Eq12,S4]==0,S1==1,RandomInteger[100] S2+RandomInteger[100]S3==1},600];
AppendTo[ptssing,sls[[1]]];
AppendTo[ptssing,sls[[2]]];
AppendTo[ptssing,sls[[3]]];
AppendTo[ptssing,sls[[4]]];
AppendTo[ptssing,sls[[5]]];
AppendTo[ptssing,sls[[6]]];
];

deg=9;
Mons=Flatten[Table[S1^i1 S2^i2 S3^i3 S4^(deg-i1-i2-i3),{i1,0,deg},{i2,0,deg-i1},{i3,0,deg-i1-i2}]];
{Length[Mons],Length[ptssing]}
MatrixRank[Table[Mons[[k]]/.ptssing[[l]],{k,1,Length[Mons]},{l,1,Length[ptssing]}]]

(* finding an equation through the these curves *)
ns=RowReduce[NullSpace[Table[Mons[[k]]/.ptssing[[l]],{l,1,Length[ptssing]},{k,1,Length[Mons]}]]];
Accuracy[ns]

EqExtra=RootApproximant[Chop[ns[[1]],10^(-400)],2].Mons;

For[numtry=1,numtry≤1000,numtry++,
If[Mod[numtry,50]==25,Print[numtry]];
fi=FindRoot[{D[Eq12,S1],D[Eq12,S2],D[Eq12,S3],EqExtra-1},{{S1,RandomReal[10]+I RandomReal[10]},{S2,RandomReal[10]+I RandomReal[10]},{S3,RandomReal[10]+I RandomReal[10]},{S4,RandomReal[10]+I RandomReal[10]}},WorkingPrecision->800,AccuracyGoal->300,MaxIterations->200,PrecisionGoal->300,DampingFactor->2];
If[Abs[N[D[Eq12,S4]/.fi,200]]<10^(-100),numtry=5000;
]
];

(* we may have gotten lucky here. *)
(* we can find the symmetric functions nicely *)
RootApproximant[N[((S2+S3+S4)/S1)/.fi,100],2]
RootApproximant[N[((S2^2+S3^2+S4^2)/S1^2)/.fi,100],2]
RootApproximant[N[((S2 S3 S4)/S1^3)/.fi,100],2]

The code below will find these singularities, which are presumably A2.

slsA2=NSolve[{D[Eq12,S1]==0,D[Eq12,S2]==0,D[Eq12,S3]==0,D[Eq12,S4]==0,S2+S3+S4 ==(((-6634-1690 I (Sqrt[2]))/4689))S1,S2^2+S3^2+S4^2==(((19047874+9436960 I (Sqrt[2]))/21986721))S1^2,(S2 S3 S4)/S1^3==((-31580371658-55821877592 I (Sqrt[2]))/927861612921),S1==1},10000];

Length[slsA2]
RootApproximant[{S1,S2,S3,S4}/.slsA2[[1]],6]

(* Let us try to see what Eq12 looks like in the neighborhood *)
(* The Hessian has the correct rank, so it's at least A2 singularity *)
Block[{$MaxExtraPrecision=400},MatrixRank[D[Eq12/.S1->1,{{S2,S3,S4},2}]/.slsA2[[1]]]]

slsA2all=NSolve[Union[Table[0==EqsFPPC3withrr[[k]]/.rr->(-2)^(1/2),{k,1,38}],{S1==1,S2==(NumericalApproximation[Root, (RowDefault["-0.946", FEExpressionsNumericalApproximationElider])-(RowDefault["0.287", FEExpressionsNumericalApproximationElider]) I, Root[7791554052+93589206276 #1+509485577737 ((#1)^2)+1607527101208 ((#1)^3)+2910688975638 ((#1)^4)+2625478328052 ((#1)^5)+927861612921 ((#1)^6)&,1], ]),S3==(NumericalApproximation[Root, (RowDefault["-0.305", FEExpressionsNumericalApproximationElider])+(RowDefault["0.0312", FEExpressionsNumericalApproximationElider]) I, Root[7791554052+93589206276 #1+509485577737 ((#1)^2)+1607527101208 ((#1)^3)+2910688975638 ((#1)^4)+2625478328052 ((#1)^5)+927861612921 ((#1)^6)&,4], ]),S4==(NumericalApproximation[Root, (RowDefault["-0.164", FEExpressionsNumericalApproximationElider])-(RowDefault["0.254", FEExpressionsNumericalApproximationElider]) I, Root[7791554052+93589206276 #1+509485577737 ((#1)^2)+1607527101208 ((#1)^3)+2910688975638 ((#1)^4)+2625478328052 ((#1)^5)+927861612921 ((#1)^6)&,5], ])}],10000];

Length[slsA2all]
(* the two points are really the same in the weighted projective space, i.e. the degree 3 coordinates differ by sign *)
OneSingPt=RootApproximant[VarsS/.slsA2all[[1]],12]

(* Verifying the symmetries *)
EqsFPPC3=EqsFPPC3withrr/.rr->(-2)^(1/2);
SymC3={S2->S3,S3->S4,S4->S2,S5->S6,S6->S7,S7->S5,S8->S9,S9->S10,S10->S8,S11->S12,S12->S13,S13->S11};
AllMons5={};For[i=1,i≤4,i++,For[j=5,j≤14,j++,AppendTo[AllMons5,VarsS[[i]] VarsS[[j]]];]];
MatrixRank[N[Join[Table[Coefficient[EqsFPPC3[[k]],AllMons5[[j]]],{k,1,9},{j,1,Length[AllMons5]}],Table[Coefficient[EqsFPPC3[[k]]/.SymC3,AllMons5[[j]]],{k,1,9},{j,1,Length[AllMons5]}]],400]]
AllMons6={};For[i=5,i≤14,i++,For[j=5,j≤i,j++,AppendTo[AllMons6,VarsS[[i]] VarsS[[j]]];]];
For[i=1,i≤4,i++,For[j=1,j≤i,j++,For[k=1,k≤j,k++,AppendTo[AllMons6,VarsS[[i]] VarsS[[j]] VarsS[[k]]];]]];
MatrixRank[N[Join[Table[Coefficient[EqsFPPC3[[k]],AllMons6[[j]]],{k,10,38},{j,1,Length[AllMons6]}],Table[Coefficient[EqsFPPC3[[k]]/.SymC3,AllMons6[[j]]],{k,10,38},{j,1,Length[AllMons6]}]],400]]

(* this is one of the singular points *)
SingPt={1,(NumericalApproximation[Root, (RowDefault["-0.946", FEExpressionsNumericalApproximationElider])-(RowDefault["0.287", FEExpressionsNumericalApproximationElider]) I, Root[7791554052+93589206276 #1+509485577737 ((#1)^2)+1607527101208 ((#1)^3)+2910688975638 ((#1)^4)+2625478328052 ((#1)^5)+927861612921 ((#1)^6)&,1], ]),(NumericalApproximation[Root, (RowDefault["-0.305", FEExpressionsNumericalApproximationElider])+(RowDefault["0.0312", FEExpressionsNumericalApproximationElider]) I, Root[7791554052+93589206276 #1+509485577737 ((#1)^2)+1607527101208 ((#1)^3)+2910688975638 ((#1)^4)+2625478328052 ((#1)^5)+927861612921 ((#1)^6)&,4], ]),(NumericalApproximation[Root, (RowDefault["-0.164", FEExpressionsNumericalApproximationElider])-(RowDefault["0.254", FEExpressionsNumericalApproximationElider]) I, Root[7791554052+93589206276 #1+509485577737 ((#1)^2)+1607527101208 ((#1)^3)+2910688975638 ((#1)^4)+2625478328052 ((#1)^5)+927861612921 ((#1)^6)&,5], ]),(NumericalApproximation[Root, (RowDefault["-0.0449", FEExpressionsNumericalApproximationElider])+(RowDefault["0.0861", FEExpressionsNumericalApproximationElider]) I, Root[60543101014034339302359321897216+12429199753843940249542310869880000 ((#1)^2)+1407528888440538494286034776825979876 ((#1)^4)+71521352845168767631662499939019863764 ((#1)^6)+1396108467620609455806310632775609087605 ((#1)^8)+15497619198439764302725564534095602352990 ((#1)^10)+424526377259868304341133068224411342619801 ((#1)^12)&,4], ]),(NumericalApproximation[Root, (RowDefault["-0.182", FEExpressionsNumericalApproximationElider])-(RowDefault["0.154", FEExpressionsNumericalApproximationElider]) I, Root[60543101014034339302359321897216+12429199753843940249542310869880000 ((#1)^2)+1407528888440538494286034776825979876 ((#1)^4)+71521352845168767631662499939019863764 ((#1)^6)+1396108467620609455806310632775609087605 ((#1)^8)+15497619198439764302725564534095602352990 ((#1)^10)+424526377259868304341133068224411342619801 ((#1)^12)&,1], ]),(NumericalApproximation[Root, (RowDefault["-6.32", FEExpressionsNumericalApproximationElider, " × ", (Superscript["10", "-3"])])+(RowDefault["0.149", FEExpressionsNumericalApproximationElider]) I, Root[60543101014034339302359321897216+12429199753843940249542310869880000 ((#1)^2)+1407528888440538494286034776825979876 ((#1)^4)+71521352845168767631662499939019863764 ((#1)^6)+1396108467620609455806310632775609087605 ((#1)^8)+15497619198439764302725564534095602352990 ((#1)^10)+424526377259868304341133068224411342619801 ((#1)^12)&,6], ]),(NumericalApproximation[Root, (RowDefault["0.137", FEExpressionsNumericalApproximationElider])+(RowDefault["0.0533", FEExpressionsNumericalApproximationElider]) I, Root[8774003351756980932483681+672833255435424282364104654 ((#1)^2)-47969767545943553537000401121 ((#1)^4)-741157855976810274095630917276 ((#1)^6)+226394812295088348972266045484255 ((#1)^8)-7297605804433029837148382978293938 ((#1)^10)+88757919455440397866081312302967329 ((#1)^12)&,10], ]),(NumericalApproximation[Root, (RowDefault["0.198", FEExpressionsNumericalApproximationElider])+(RowDefault["0.0608", FEExpressionsNumericalApproximationElider]) I, Root[8774003351756980932483681+672833255435424282364104654 ((#1)^2)-47969767545943553537000401121 ((#1)^4)-741157855976810274095630917276 ((#1)^6)+226394812295088348972266045484255 ((#1)^8)-7297605804433029837148382978293938 ((#1)^10)+88757919455440397866081312302967329 ((#1)^12)&,12], ]),(NumericalApproximation[Root, (RowDefault["0.0108", FEExpressionsNumericalApproximationElider])-(RowDefault["0.103", FEExpressionsNumericalApproximationElider]) I, Root[8774003351756980932483681+672833255435424282364104654 ((#1)^2)-47969767545943553537000401121 ((#1)^4)-741157855976810274095630917276 ((#1)^6)+226394812295088348972266045484255 ((#1)^8)-7297605804433029837148382978293938 ((#1)^10)+88757919455440397866081312302967329 ((#1)^12)&,7], ]),(NumericalApproximation[Root, (RowDefault["-0.0910", FEExpressionsNumericalApproximationElider])-(RowDefault["0.191", FEExpressionsNumericalApproximationElider]) I, Root[164443232352822997829932876275122176+1195802406701844764334137174595338240 ((#1)^2)-15215628937860468439384337043256672704 ((#1)^4)-1534227413301414894026705303432679582240 ((#1)^6)+6538739076896138007723797020037997671700 ((#1)^8)+132026059514896430456644419098449783766916 ((#1)^10)+424526377259868304341133068224411342619801 ((#1)^12)&,5], ]),(NumericalApproximation[Root, (RowDefault["0.0967", FEExpressionsNumericalApproximationElider])-(RowDefault["0.449", FEExpressionsNumericalApproximationElider]) I, Root[164443232352822997829932876275122176+1195802406701844764334137174595338240 ((#1)^2)-15215628937860468439384337043256672704 ((#1)^4)-1534227413301414894026705303432679582240 ((#1)^6)+6538739076896138007723797020037997671700 ((#1)^8)+132026059514896430456644419098449783766916 ((#1)^10)+424526377259868304341133068224411342619801 ((#1)^12)&,9], ]),(NumericalApproximation[Root, (RowDefault["0.256", FEExpressionsNumericalApproximationElider])-(RowDefault["0.0269", FEExpressionsNumericalApproximationElider]) I, Root[164443232352822997829932876275122176+1195802406701844764334137174595338240 ((#1)^2)-15215628937860468439384337043256672704 ((#1)^4)-1534227413301414894026705303432679582240 ((#1)^6)+6538739076896138007723797020037997671700 ((#1)^8)+132026059514896430456644419098449783766916 ((#1)^10)+424526377259868304341133068224411342619801 ((#1)^12)&,11], ]),(NumericalApproximation[Root, (RowDefault["0.780", FEExpressionsNumericalApproximationElider])-(RowDefault["0.739", FEExpressionsNumericalApproximationElider]) I, Root[53258647842318380176-4964124530200959864 ((#1)^2)+39941399978020282041 ((#1)^4)&,3], ])};
Block[{$MaxExtraPrecision=10000},Chop[N[EqsFPPC3/.Table[VarsS[[l]]->SingPt[[l]],{l,1,14}],1000]]]

Let us find a set of points on this surface, which we will use to change to new coordinates.

PtsOnXinS={};
For[numpt=1,numpt≤50,numpt++,
If[Mod[numpt,10]==2,Print[numpt]];
sls=NSolve[Join[Table[0==EqsFPPC3[[k]],{k,1,38}],{S1==1,S2==RandomInteger[100]/100+I RandomInteger[100]/100,
S3==RandomInteger[100]/100+I RandomInteger[100]/100}],2000];
AppendTo[PtsOnXinS,sls[[1]]];
AppendTo[PtsOnXinS,sls[[5]]];
];

MatrixRank[{{SingPt[[2]],SingPt[[3]],SingPt[[4]]},{SingPt[[3]],SingPt[[4]],SingPt[[2]]},{SingPt[[4]],SingPt[[2]],SingPt[[3]]}}]

(* making a variable change *)
VarsZ={Z1,Z2,Z3,Z4,Z5,Z6,Z7,Z8,Z9,Z10,Z11,Z12,Z13,Z14};
StoZ4={S1->Z1+Z2+Z3+Z4, S2->SingPt[[2]]Z2 + SingPt[[4]]Z3 + SingPt[[3]]Z4,S3->SingPt[[2]]Z3 + SingPt[[4]]Z4 + SingPt[[3]]Z2, S4->SingPt[[2]]Z4 + SingPt[[4]]Z2 + SingPt[[3]]Z3};

k=3;l=3;
StoZ=Join[StoZ4, 
{VarsS[[2+k]]-> SingPt[[2+l]] VarsZ[[2+k]]+SingPt[[4+l]] VarsZ[[3+k]]+SingPt[[3+l]] VarsZ[[4+k]],
VarsS[[3+k]]-> SingPt[[2+l]] VarsZ[[3+k]]+SingPt[[4+l]] VarsZ[[4+k]]+SingPt[[3+l]] VarsZ[[2+k]],
VarsS[[4+k]]-> SingPt[[2+l]] VarsZ[[4+k]]+SingPt[[4+l]] VarsZ[[2+k]]+SingPt[[3+l]] VarsZ[[3+k]]}];

k=6;l=6;
StoZ=Join[StoZ, 
{VarsS[[2+k]]-> SingPt[[2+l]] VarsZ[[2+k]]+SingPt[[4+l]] VarsZ[[3+k]]+SingPt[[3+l]] VarsZ[[4+k]],
VarsS[[3+k]]-> SingPt[[2+l]] VarsZ[[3+k]]+SingPt[[4+l]] VarsZ[[4+k]]+SingPt[[3+l]] VarsZ[[2+k]],
VarsS[[4+k]]-> SingPt[[2+l]] VarsZ[[4+k]]+SingPt[[4+l]] VarsZ[[2+k]]+SingPt[[3+l]] VarsZ[[3+k]]}];

k=9;l=9;
StoZ=Join[StoZ, 
{VarsS[[2+k]]-> SingPt[[2+l]] VarsZ[[2+k]]+SingPt[[4+l]] VarsZ[[3+k]]+SingPt[[3+l]] VarsZ[[4+k]],
VarsS[[3+k]]-> SingPt[[2+l]] VarsZ[[3+k]]+SingPt[[4+l]] VarsZ[[4+k]]+SingPt[[3+l]] VarsZ[[2+k]],
VarsS[[4+k]]-> SingPt[[2+l]] VarsZ[[4+k]]+SingPt[[4+l]] VarsZ[[2+k]]+SingPt[[3+l]] VarsZ[[3+k]]}];

AppendTo[StoZ, S14->Z14 SingPt[[14]]];
ZtoS = Solve[Table[VarsS[[k]]==(VarsS[[k]]/.StoZ),{k,1,14}],VarsZ][[1]];
PtsOnXinZ={};
For[j=1,j≤Length[PtsOnXinS],j++,
AppendTo[PtsOnXinZ,N[Table[VarsZ[[k]]->(VarsZ[[k]]/.ZtoS)/.PtsOnXinS[[j]],{k,1,14}],2000]];
];

(* finding new equations *)
AllMons5={};For[i=1,i≤4,i++,For[j=5,j≤14,j++,AppendTo[AllMons5,VarsZ[[i]] VarsZ[[j]]];]];
AllMons6={};For[i=5,i≤14,i++,For[j=5,j≤i,j++,AppendTo[AllMons6,VarsZ[[i]] VarsZ[[j]]];]];
For[i=1,i≤4,i++,For[j=1,j≤i,j++,For[k=1,k≤j,k++,AppendTo[AllMons6,VarsZ[[i]] VarsZ[[j]] VarsZ[[k]]];]]];

ns=RowReduce[NullSpace[Table[AllMons5[[k]]/.PtsOnXinZ[[l]],{l,1,Length[PtsOnXinZ]},{k,1,Length[AllMons5]}]]];
ns=Chop[ns,10^(-Round[Accuracy[ns]]+5)];
Accuracy[ns]
ns=Block[{$MaxExtraPrecision=400},RootApproximant[ns,2]];
EqsFPPC3Zdeg5=ns.AllMons5;
ns=RowReduce[NullSpace[Table[AllMons6[[k]]/.PtsOnXinZ[[l]],{l,1,Length[PtsOnXinZ]},{k,1,Length[AllMons6]}]]];
ns=Chop[ns,10^(-Round[Accuracy[ns]]+5)];
Accuracy[ns]
ns=Block[{$MaxExtraPrecision=400},RootApproximant[ns,2]];
EqsFPPC3Zdeg6=ns.AllMons6;

EqsFPPC3Z=Join[EqsFPPC3Zdeg5,EqsFPPC3Zdeg6];
Length[EqsFPPC3Z]
EqsFPPC3Z=Table[Numerator[Factor[EqsFPPC3Z[[k]]]],{k,1,38}];
(* The name of the file needs to be adjusted according to your situation. *)
fname="Desktop/FPP2020/EqsFPPC3Z";
OpenWrite[fname];
Write[fname,EqsFPPC3Z];
Close[fname];

sls=NSolve[Join[Table[0==EqsFPPC3Z[[k]],{k,1,38}],{Z2==1,Z3==0,Z4==0,Z5==1,Z6==0,Z7==0,Z8==1,Z9==0,
Z10==0,Z11==1,Z12==0}],600];
Length[sls]
SingPtZ=RootApproximant[VarsZ/.sls[[1]]]

We will try to make the variables nicer. Our first goal is to have good first four variables.

CyclicZ={Z2->Z3,Z3->Z4,Z4->Z2,Z5->Z6,Z6->Z7,Z7->Z5,Z8->Z9,Z9->Z10,Z10->Z8,Z11->Z12,Z12->Z13,Z13->Z12};
fixedpts=NSolve[Join[Join[Table[0==EqsFPPC3Z[[k]],{k,1,38}],Table[VarsZ[[k]]==(VarsZ[[k]]/.CyclicZ),{k,1,14}]],{Z1==1}],600];
Length[fixedpts]
FirstFourFixed=Table[RootApproximant[{Z1,Z2,Z3,Z4}/.fixedpts[[k+1]],12],{k,1,3}]
w=Exp[2Pi I/3]
RootReduce[(NumericalApproximation[Root, (RowDefault["0.346", FEExpressionsNumericalApproximationElider])-(RowDefault["0.313", FEExpressionsNumericalApproximationElider]) I, Root[95976324187162929-327346819240559508 #1+575126053567930852 ((#1)^2)-302805817555985640 ((#1)^3)+292145372160311475 ((#1)^4)&,3], ]) w +(NumericalApproximation[Root, (RowDefault["0.173", FEExpressionsNumericalApproximationElider])+(RowDefault["1.22", FEExpressionsNumericalApproximationElider]) I, Root[95976324187162929-327346819240559508 #1+575126053567930852 ((#1)^2)-302805817555985640 ((#1)^3)+292145372160311475 ((#1)^4)&,2], ]) w^2]
RootReduce[(NumericalApproximation[Root, (RowDefault["0.346", FEExpressionsNumericalApproximationElider])-(RowDefault["0.313", FEExpressionsNumericalApproximationElider]) I, Root[95976324187162929-327346819240559508 #1+575126053567930852 ((#1)^2)-302805817555985640 ((#1)^3)+292145372160311475 ((#1)^4)&,3], ]) w^2 +(NumericalApproximation[Root, (RowDefault["0.173", FEExpressionsNumericalApproximationElider])+(RowDefault["1.22", FEExpressionsNumericalApproximationElider]) I, Root[95976324187162929-327346819240559508 #1+575126053567930852 ((#1)^2)-302805817555985640 ((#1)^3)+292145372160311475 ((#1)^4)&,2], ]) w]
FullSimplify[Solve[{a1 FirstFourFixed[[1]][[1]] + a2 FirstFourFixed[[1]][[2]]==1}][[1]]]

(* We can not make these value too nice, but we can try to make FirstFourFixed[[1]][[2]] and the w-traces of the others as small as possible. 
Specifically, we want to use a1, a2 in Q(I (Sqrt[2])) to get the following three values smaller. *)
wantsmall=FullSimplify[{a1 + ((-7232-6409 I (Sqrt[2]))/22815)a2,-a1 +((478046747-96008522 I (Sqrt[2]))/448337565)a2,-a1+((-710395495-190928201 I (Sqrt[2]))/448337565)a2 }]
AllOptions=FullSimplify[{wantsmall/.{a1->1,a2->0},wantsmall/.{a1->I (Sqrt[2]),a2->0},wantsmall/.{a2->1,a1->0},wantsmall/.{a2->I (Sqrt[2]),a1->0}}]

PossibleLinComb={{a1->1,a2->0},{a1->I (Sqrt[2]),a2->0},{a2->1,a1->0},{a2->I (Sqrt[2]),a1->0}};
ToBeSmall={};
For[i=1,i≤4,i++,
row=wantsmall/.PossibleLinComb[[i]];
newrow={};
For[j=1,j≤3,j++,
AppendTo[newrow,(-1/2)SeriesCoefficient[MinimalPolynomial[row[[j]]+100+100I (Sqrt[2]),t],{t,0,1}]/SeriesCoefficient[MinimalPolynomial[row[[j]]+100+100I (Sqrt[2]),t],{t,0,2}]-100];
AppendTo[newrow,(-1/2)SeriesCoefficient[MinimalPolynomial[row[[j]]/I (Sqrt[2])+100+100I (Sqrt[2]),t],{t,0,1}]/SeriesCoefficient[MinimalPolynomial[row[[j]]+100+100I (Sqrt[2]),t],{t,0,2}]-100];
];
AppendTo[ToBeSmall,newrow];
];
ToBeSmall

(* we want to find a linear combination of these four vectors in dimension 6 space that would have small integer values *)
MatrixRank[ToBeSmall]
{u,r,v}=SmithDecomposition[ToBeSmall];
vt=Transpose[v];
SmallBestRels={};
For[i=1,i≤Length[ToBeSmall],i++,
If[Mod[i,10]==0,Print[i];];
basisvec=Table[0,{k,1,Length[ToBeSmall[[1]]]}];
basisvec[[i]]=1;
AppendTo[SmallBestRels,LinearSolve[vt,basisvec]];
];
ByteCount[SmallBestRels] 
SmallerBestRels=LatticeReduce[SmallBestRels];
ByteCount[SmallerBestRels]

SmallerBestRels
({b1,b2,b3,b4}/.
Solve[{b1,b2,b3,b4}.ToBeSmall=={13,-35,27,54,-14,-90}][[1]])
{a1,a2}/.{{a1->1,a2->0},{a1->I (Sqrt[2]),a2->0},{a2->1,a1->0},{a2->I (Sqrt[2]),a1->0}};
FullSimplify[wantsmall/.{a1-> (114916/23859)-(185699/47718) I(Sqrt[2]),a2-> (243097205/12430539)+(636441325/24861078)I (Sqrt[2])  }];

Table[RootApproximant[{Z1,a1 Z1 + a2 Z2,a1 Z1 + a2 Z3,a1 Z1 + a2 Z4}/.{a1->(( 114916)/23859)-(185699/47718)I (Sqrt[2]),a2-> (243097205/12430539)+(636441325/24861078)I (Sqrt[2])  }/.fixedpts[[k+1]],12],{k,1,3}]
(* So this looks pretty good. How do we make a change of variables so that new Z-s are related to old Z-s this way *)
ZtoNewZ=FullSimplify[Solve[{Z1n,Z2n,Z3n,Z4n}==({Z1,a1 Z1 + a2 Z2,a1 Z1 + a2 Z3,a1 Z1 + a2 Z4}/.{a1->(( 114916)/23859)-(185699/47718)I (Sqrt[2]),a2-> (243097205/12430539)+(636441325/24861078)I (Sqrt[2])  }),{Z1,Z2,Z3,Z4}]][[1]];
NewZtoZ={Z1n->Z1,Z2n-> a1 Z1 + a2 Z2,Z3n-> a1 Z1 + a2 Z3,Z4n-> a1 Z1 + a2 Z4}/.{a1->(( 114916)/23859)-(185699/47718)I (Sqrt[2]),a2-> (243097205/12430539)+(636441325/24861078)I (Sqrt[2])  };

(* Let us see if the equation of degree 12 looks good in these Z-s *)
(* The name of the file needs to be adjusted according to your situation. *)
fname="Desktop/FPP2020/Eq12";
OpenRead[fname];
Eq12S=Read[fname];
Close[fname];

ptsZnew={};
For[numpt=1,numpt≤500,numpt++,
If[Mod[numpt,100]==20,Print[numpt]];
pt={S1,S2,S3,S4}/.NSolve[{0==Eq12S,S1==1,S2==RandomInteger[1000]/1000+I RandomInteger[1000]/1000, S3==RandomInteger[1000]/1000+I RandomInteger[1000]/1000},600][[1]];
AppendTo[ptsZnew,NewZtoZ/.NSolve[({S1,S2,S3,S4}/.StoZ)==pt][[1]]];
];

MonsZ12n={};
For[i1=0,i1≤12,i1++,
For[i2=0,i2≤12-i1,i2++,
For[i3=0,i3≤12-i1-i2,i3++,
AppendTo[MonsZ12n,Z1n^i1 Z2n^i2 Z3n^i3 Z4n^(12-i1-i2-i3)];
]]];
BigMatrix=MonsZ12n/.ptsZnew;

ns=NullSpace[BigMatrix];
Length[ns]
Accuracy[ns]
nsr=Block[{$MaxExtraPrecision=1000},RootApproximant[RowReduce[Chop[ns,10^(-500)]],2]];
Eq12Znew=Numerator[Factor[nsr[[1]].MonsZ12n]];
ByteCount[Eq12S]
ByteCount[Eq12Znew]

Let us just write equations in new variables W and see if these are at all shorter.

VarsW={W1,W2,W3,W4,W5,W6,W7,W8,W9,W10,W11,W12,W13,W14};
WtoZ=Join[{W1->Z1n,W2->Z2n,W3->Z3n,W4->Z4n}/.NewZtoZ,Table[VarsW[[k]]->VarsZ[[k]]/ (NumericalApproximation[Root, (RowDefault["1.61", FEExpressionsNumericalApproximationElider, " × ", (Superscript["10", "-4"])])+(RowDefault["3.79", FEExpressionsNumericalApproximationElider, " × ", (Superscript["10", "-3"])]) I, Root[15365966018055877446552+2133382522369818357486433000 ((#1)^2)+74586114717210050853770127140625 ((#1)^4)&,4], ]),{k,5,14}]];
PtsOnXinW={};

For[j=1,j≤Length[PtsOnXinZ],j++,
AppendTo[PtsOnXinW,N[Table[VarsW[[k]]->(VarsW[[k]]/.WtoZ)/.PtsOnXinZ[[j]],{k,1,14}],2000]];
];

AllMons5={};For[i=1,i≤4,i++,For[j=5,j≤14,j++,AppendTo[AllMons5,VarsW[[i]] VarsW[[j]]];]];
AllMons6={};For[i=5,i≤14,i++,For[j=5,j≤i,j++,AppendTo[AllMons6,VarsW[[i]] VarsW[[j]]];]];
For[i=1,i≤4,i++,For[j=1,j≤i,j++,For[k=1,k≤j,k++,AppendTo[AllMons6,VarsW[[i]] VarsW[[j]] VarsW[[k]]];]]];

ns=RowReduce[NullSpace[Table[AllMons5[[k]]/.PtsOnXinW[[l]],{l,1,Length[PtsOnXinW]},{k,1,Length[AllMons5]}]]];
ns=Chop[ns,10^(-Round[Accuracy[ns]]+5)];
Accuracy[ns]
ns=Block[{$MaxExtraPrecision=400},RootApproximant[ns,2]];
EqsFPPC3Wdeg5=ns.AllMons5;
ns=RowReduce[NullSpace[Table[AllMons6[[k]]/.PtsOnXinW[[l]],{l,1,Length[PtsOnXinW]},{k,1,Length[AllMons6]}]]];
ns=Chop[ns,10^(-Round[Accuracy[ns]]+5)];
Accuracy[ns]
ns=Block[{$MaxExtraPrecision=400},RootApproximant[ns,2]];
EqsFPPC3Wdeg6=ns.AllMons6;
EqsFPPC3W=Join[Table[Numerator[Factor[EqsFPPC3Wdeg5[[k]]]],{k,1,9}],Table[Numerator[Factor[EqsFPPC3Wdeg6[[k]]]],{k,1,29}]];

sls=NSolve[Join[Table[0==EqsFPPC3W[[k]],{k,1,38}],{W1==0,W2==1,W3==0,W4==0}],500];
RootApproximant[VarsW/.sls]

Let us not try to make this nice but just go with what we have. Let us find a bunch of points on this surface with high accuracy.

SetOfPointsW = {};
prec=6000;
(* to be skipped the second time around *)
numpts=200;
For[i=1,i≤numpts,i++,
If[Mod[i,numpts/10]==0,Print[i]];
rpt1=(RandomInteger[{1,100000}]+I RandomInteger[{1,100000}])/100000;
rpt2=(RandomInteger[{1,100000}]+I RandomInteger[{1,100000}])/100000;
Eqs=Join[{W1==1,W2+W3+W4==rpt1,W5+W6+W7==rpt2},Table[0==EqsFPPC3W[[j]],{j,1,38}]];
sols=NSolve[Eqs,200];
fr=FindRoot[Table[Eqs[[k]],{k,1,14}],Table[{VarsW[[j]],VarsW[[j]]/.sols[[1]]},{j,1,14}],MaxIterations->100,WorkingPrecision->prec,AccuracyGoal->prec];
AppendTo[SetOfPointsW,fr];
];
(* quality control *)
For[i=1,i≤numpts,i++, If[Chop[Sum[Abs[EqsFPPC3W[[j]]]/.SetOfPointsW[[i]],{j,1,38}]]≠0,Print["warning, point ",i," looks fishy"];];];

(* The name of the file needs to be adjusted according to your situation. *)fname="Desktop/FPP2020/SetOfPointsW";
OpenWrite[fname];
Write[fname,SetOfPointsW];
Close[fname];

fname="Desktop/FPP2020/EqsFPPC3W";
OpenWrite[fname];
Write[fname,EqsFPPC3W];
Close[fname];