This notebook is a supporting file for Section 6 of the paper. It contains the calculation of the X9 fake projective plane.

We first read the double cover of the surface D3X3. We need to rewrite it in good coordinates. We will also need to use the equations of the FPP itself as an intermediate set of equations to generate points. 


The action is:
C3EqsFPPnewSub={Q1->Q2, Q2->Q3, Q3->Q1, Q5->Q6,Q6->Q7,Q7->Q5, Q8->Q10,Q10->Q12,Q12->Q8, Q9->Q11,Q11->Q13,Q13->Q9, Q14->Q16,Q16->Q18,Q18->Q14,Q15->Q17,Q17->Q19,Q19->Q15};
C2actionQ={Q0->Q0,Q1->Q1,Q2->Q2,Q3->Q3,Q4->-Q4,Q5->-Q5,Q6->-Q6,Q7->-Q7,Q8->Q9, Q9->Q8,Q10->Q11,Q11->Q10,Q12->Q13,Q13->Q12,Q14->Q15,Q15->Q14,Q16->Q17,Q17->Q16,Q18->Q19,Q19->Q18};

SetDirectory[$InitialDirectory];
SetDirectory["DeskTop/EnricoFPP/paperfiles"];

InputSurfaceFileName="EqsFPPD3X3double";

OpenRead[InputSurfaceFileName];
InputEqs=Read[InputSurfaceFileName];
Close[InputSurfaceFileName];

fname="EqsFPPD3X3";
OpenRead[fname];
EqsFPPD3X3=Read [fname];
Close[fname];

VarsP = {P0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14,P15,P16,P17,P18,P19};
VarsQ = {Q0,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,Q10,Q11,Q12,Q13,Q14,Q15,Q16,Q17,Q18,Q19};
InputEqs=InputEqs/.{Q0->P0,Q1->P1, Q2->P2,Q3->P3,Q4->P10,Q5->P11,Q6->P12,Q7->P13,
Q8->P4+P14,Q9->P4-P14,Q10->P5+P15,Q11->P5-P15,Q12->P6+P16,Q13->P6-P16,
Q14->P7+P17,Q15->P7-P17,Q16->P8+P18,Q17->P8-P18,Q18->P9+P19,Q19->P9-P19
};

EqsFPPD3X3=EqsFPPD3X3/.{Q0->P0,Q1->P1, Q2->P2,Q3->P3,Q4->P4,Q5->P5,Q6->P6,Q7->P7,
Q8->P8,Q9->P9};

PC3Sub={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};
C2actionP=Table[VarsP[[i]]->-VarsP[[i]],{i,11,20}];
Length[InputEqs]
Length[EqsFPPD3X3]

(* Let's verify that these are symmetries *)
AllMons={};
For[i=1,i≤20,i++,For[j=1,j≤i,j++, AppendTo[AllMons,VarsP[[i]] VarsP[[j]]]]];

MatrixRank[Table[N[Coefficient[InputEqs[[i]],AllMons[[j]]],500],{i,1,100},{j,1,210}]]
MatrixRank[Join[Table[N[Coefficient[InputEqs[[i]],AllMons[[j]]],500],{i,1,100},{j,1,210}],
Table[N[Coefficient[InputEqs[[i]]/.PC3Sub,AllMons[[j]]],500],{i,1,100},{j,1,210}]]]

MatrixRank[Join[Table[N[Coefficient[InputEqs[[i]],AllMons[[j]]],500],{i,1,100},{j,1,210}],
Table[N[Coefficient[InputEqs[[i]]/.C2actionP,AllMons[[j]]],500],{i,1,100},{j,1,210}]]]

Then we compute a number of points on this surface. This is complicated. The trick is to compute the points for the FPP first and then go to the double cover.

(* trial run, and finding independent equations *)
rpt1=(RandomInteger[1000]+I RandomInteger[1000])/100;
rpt2=(RandomInteger[1000]+I RandomInteger[1000])/100;
sub3={P0->1,P1->rpt1,P2->rpt2};
Eqs=Table[0==EqsFPPD3X3[[k]]/.sub3,{k,1,84}];
sols=FindInstance[N[Eqs,1200],Table[VarsP[[i]],{i,4,10}]];
Eqs1=Table[0==N[InputEqs[[k]],1200]/.Join[sub3,sols[[1]]],{k,1,100}];
sols1=FindInstance[N[Eqs1,1200],Table[VarsP[[i]],{i,11,20}]];
sols17=Join[sols1[[1]],sols[[1]]];
GoodIndices={};
rk=0;
EqsP0=InputEqs/.sub3;
JacMatAtp={};
For[k=1,k≤100,k++,
JacMatAtpnew=Append[JacMatAtp,Table[D[N[EqsP0[[k]],1200],VarsP[[l]]]/.sols17,{l,4,20}]];
If[MatrixRank[JacMatAtpnew]==rk+1,
JacMatAtp=JacMatAtpnew;rk++;AppendTo[GoodIndices,k];];
If[rk==17,k=200];
];
GoodIndices

Warning: finding points is time-consuming.  One may want to write points into a file afterwards, just in case the code needs to be restarted.

NumPoints=150;
SetOfPointsP={};
For[i=1,i≤NumPoints,i++,
rpt1=(RandomInteger[1000]+I RandomInteger[1000])/100;
rpt2=(RandomInteger[1000]+I RandomInteger[1000])/100;
sub3=N[{P0->1,P1->rpt1,P2->rpt2},5500];
Eqs=Table[0==EqsFPPD3X3[[k]]/.sub3,{k,1,84}];
sols=FindInstance[N[Eqs,1200],Table[VarsP[[i]],{i,4,10}]];
Eqs1=Table[0==N[InputEqs[[k]],1200]/.Join[sub3,sols[[1]]],{k,1,100}];
sols1=FindInstance[N[Eqs1,1200],Table[VarsP[[i]],{i,11,20}]];
sols17=Join[sols1[[1]],sols[[1]]];
rk=0;
GoodEqs=Table[InputEqs[[GoodIndices[[k]]]]/.sub3,{k,1,17}];
fr= FindRoot[GoodEqs,Table[{VarsP[[k]],VarsP[[k]]/.sols17},{k,4,20}],MaxIterations->200,WorkingPrecision->5000,AccuracyGoal->5000];
fr=Join[fr,sub3];
If[Sum[Abs[InputEqs[[k]]/.fr],{k,1,100}]<10^(-1400), (* quality control *)
AppendTo[SetOfPointsP,fr];
];
If[Mod[i,30]==0,Print[{i, Length[SetOfPointsP]}]];
];

(* saving to a file *)
Length[SetOfPointsP]
fname="FPPD3X3DoublePoints";
OpenWrite[fname];
Write[fname,SetOfPointsP];
Close[fname];

(* reading from a file *)
fname="FPPD3X3DoublePoints";
OpenRead[fname];
SetOfPointsP=Read[fname];
Close[fname];
Length[SetOfPointsP]

The next step is to solve the equation that determines the divisors for the new FPP.
We look at Sym^2(Wt_{0,0}+Sym^2(Wt_{1,0}) + Wt_{0,1} \otimes Wt_{0,2}.
We first find the four spaces in question.

Beautify[x_]:=
Ret=Expand[RootApproximant[Chop[N[Re[x],4900],10^(-4800)],2]+RootApproximant[Chop[N[Im[x]/(Sqrt[15]),4900],10^(-4800)],2] I(Sqrt[15])];

w=Exp[2 Pi I /3];
BasisWt00={};
rk=0;
For[i=1,i≤10,i++,
NewBasis=Append[BasisWt00,Simplify[(1/3)(VarsP[[i]]+(VarsP[[i]]/.PC3Sub )+ ((VarsP[[i]]/.PC3Sub )/.PC3Sub))]];
If[MatrixRank[Table[NewBasis/.SetOfPointsP[[k]],{k,1,10}]]==rk+1, rk++;BasisWt00=NewBasis;];
];
If[Length[BasisWt00]≠4, Print["Warning: unexpected dimension of the space"];];

BasisWt10={};
rk=0;
For[i=11,i≤20,i++,
NewBasis=Append[BasisWt10,Simplify[(1/3)(VarsP[[i]]+(VarsP[[i]]/.PC3Sub )+ ((VarsP[[i]]/.PC3Sub )/.PC3Sub))]];
If[MatrixRank[Table[NewBasis/.SetOfPointsP[[k]],{k,1,10}]]==rk+1, rk++;BasisWt10=NewBasis;];
];
If[Length[BasisWt10]≠4, Print["Warning: unexpected dimension of the space"];];

BasisWt01={};
rk=0;
For[i=1,i≤10,i++,
NewBasis=Append[BasisWt01,Simplify[(1/3)(VarsP[[i]]+ w^2 (VarsP[[i]]/.PC3Sub )+ w ((VarsP[[i]]/.PC3Sub )/.PC3Sub))]];
If[MatrixRank[Table[NewBasis/.SetOfPointsP[[k]],{k,1,10}]]==rk+1, rk++;BasisWt01=NewBasis;];
];
If[Length[BasisWt01]≠3, Print["Warning: unexpected dimension of the space"];];

BasisWt02={};
rk=0;
For[i=1,i≤10,i++,
NewBasis=Append[BasisWt02,Simplify[(1/3)(VarsP[[i]]+ w (VarsP[[i]]/.PC3Sub )+ w^2 ((VarsP[[i]]/.PC3Sub )/.PC3Sub))]];
If[MatrixRank[Table[NewBasis/.SetOfPointsP[[k]],{k,1,10}]]==rk+1, rk++;BasisWt02=NewBasis;];
];
If[Length[BasisWt02]≠3, Print["Warning: unexpected dimension of the space"];];

VarsC={c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14};
CtoBasis=Join[Join[Table[VarsC[[i]]->BasisWt00[[i]],{i,1,4}],Table[VarsC[[i+4]]->BasisWt10[[i]],{i,1,4}]],Join[Table[VarsC[[i+8]]->BasisWt01[[i]],{i,1,3}],Table[VarsC[[i+11]]->BasisWt02[[i]],{i,1,3}]]];
(* looking for quadratic equations *)
AllMons={c1^2,c1 c2, c1 c3, c1 c4, c2^2, c2 c3, c2 c4, c3^2, c3 c4, c4^2,c5^2,c5 c6, c5 c7, c5 c8, c6^2, c6 c7, c6 c8, c7^2, c7 c8, c8^2, c9 c12, c9 c13, c9 c14,c10 c12, c10 c13, c10 c14,c11 c12, c11 c13, c11 c14};
CoeffList=Table[ToExpression[StringJoin["A",ToString[i]]],{i,1,500}];
UnivLinear=Sum[CoeffList[[i]] AllMons[[i]] ,{i,1,29}];

Eqs=Table[0==UnivLinear/.CtoBasis/.SetOfPointsP[[i]],{i,1,Length[SetOfPointsP]}];
FirstSolve=Chop[Solve[Eqs][[1]]];
FirstSolve=Table[CoeffList[[i]]->
Sum[RootApproximant[Coefficient[CoeffList[[i]]/.FirstSolve,CoeffList[[k]]]] CoeffList[[k]],{k,1,29}],{i,1,29}];

Eqs={A1==1};
MatrixDer=Table[D[D[UnivLinear,VarsC[[i+8]]],VarsC[[j+11]]],{i,1,3},{j,1,3}]/.FirstSolve;
MinDer=Minors[MatrixDer,2];
For[i=1,i≤3,i++,For[j=1,j≤3,j++,AppendTo[Eqs,0==MinDer[[i]][[j]]];]];

MatrixDer=Table[D[D[UnivLinear,VarsC[[i]]],VarsC[[j]]],{i,1,4},{j,1,4}]/.FirstSolve;
MinDer=Minors[MatrixDer,2];
For[i=1,i≤6,i++,For[j=1,j≤6,j++,AppendTo[Eqs,0==MinDer[[i]][[j]]];]];

MatrixDer=Table[D[D[UnivLinear,VarsC[[i+4]]],VarsC[[j+4]]],{i,1,4},{j,1,4}]/.FirstSolve;
MinDer=Minors[MatrixDer,2];
For[i=1,i≤6,i++,For[j=1,j≤6,j++,AppendTo[Eqs,0==MinDer[[i]][[j]]];]];

Length[Eqs]
sols=NSolve[Eqs,1500];
Accuracy[sols]
Length[sols]

UnivLinearSolved=UnivLinear/.FirstSolve/.Chop[sols[[1]]];
UnivLinearSolved = Sum[Beautify[Coefficient[UnivLinearSolved,AllMons[[i]]]]AllMons[[i]],{i,1,29}];

Sections4=Simplify[(D[UnivLinearSolved,c14]/(D[UnivLinearSolved,c14]/.{c9->1,c10->0,c11->0}))/.CtoBasis]
Sections3=Simplify[(D[UnivLinearSolved,c10]/(D[UnivLinearSolved,c10]/.{c12->1,c13->0,c14->0}))/.CtoBasis]
prod=Factor[UnivLinearSolved/.{c9->0,c10->0,c11->0,c12->0,c13->0,c14->0},Extension->Automatic];
Sections1=Simplify[I prod[[2]]/.CtoBasis]
Sections2=Simplify[I prod[[3]]/.CtoBasis]

We have thus computed the sections s1, s2, s3 and s4 up to scaling. The next step is to compute the curves which these sections are composed of. We have s3  = r1 r2 and s4 = r1’ r2’. We can look at r1 as the zero set of <s3, s1> and r2 as of <s3, s2>.

Then we want to look at sections on input doubleFPP/C3 of the twist of 2K by the line bundles that are [r1]-[r2] and [r2]-[r1]. Indeed, we expect these to give the triple cover. It is not particularly important which one is r1 and which one is r2, as this only changes w to w^2. To make it work, we will need to have a rational function which realizes 3[r1]-3[r2]. Indeed, if we are looking for quadratic relations, then we must be able to go from O(2[r1]-2[r2]) to O([r2]-[r1]), so we need it. 

We will need to find H^0(2.FPP/C3, O(2K+[r2]-[r1])) and H^0(2.FPP/C3, O(2K+[r1]-[r2])). Each of these should be of dimension 6, and the action of C2 should interchange them. We can have O(2K+[r2]-[r1]) embedded into O(4K) by a section of 2K-[r1]+[r2], namely s3. In fact, this section vanished twice at [r2]. So we can try to see the sections of 2K+[r2]-[r1] as sections of O(2K) on double cover of FPP which vanish twice on r2 and also has the weight w (same as s3) with respect to C3 action. 

To construct 3[r1]-3[r2] we will look for all quadratic C3 invariants in P that are zero on 3[r1]. This should be 3[r1]+[theta] for the unique antiinvariant [theta], and likewise for 3[r2].

SetOfPointsr1={};
For[numpt=1,numpt≤30,numpt++,
If[Mod[numpt,5]==0,Print[numpt," ",Length[SetOfPointsr1]]];

rpt1=(RandomInteger[1000]+I RandomInteger[1000])/1000;
RelsToSolve=Join[Table[EqsFPPD3X3[[i]]==0,{i,1,84}],{0==Sections3,P1==rpt1,P0==1}];
sols=FindInstance[N[RelsToSolve,900],Table[VarsP[[i]],{i,1,10}]][[1]];
solslast=FindInstance[Table[0==InputEqs[[k]]/.sols,{k,1,100}],Table[VarsP[[i]],{i,11,20}]][[1]];
GoodEqs=Join[{Sections3,P1-rpt1, P0-1},Table[InputEqs[[GoodIndices[[k]]]],{k,1,17}]];
fr= FindRoot[GoodEqs,Table[{VarsP[[k]],VarsP[[k]]/.Join[sols,solslast]},{k,1,20}],MaxIterations->200,WorkingPrecision->2000,AccuracyGoal->2000];

If[Abs[N[Sections1/.fr]]<10^(-1000),AppendTo[SetOfPointsr1,fr],
frnew=Solve[Join[Table[VarsP[[k]]==(VarsP[[k]]/.fr),{k,1,10}],Table[VarsP[[k]]==-(VarsP[[k]]/.fr),{k,11,20}]]];
If[Abs[N[Sections1/.frnew[[1]]]]<10^(-1000),AppendTo[SetOfPointsr1,frnew[[1]]]]];
];

SetOfPointsr2={};
For[i=1,i≤Length[SetOfPointsr1],i++,
newpoint={};
For[l=1,l≤10,l++,
AppendTo[newpoint,VarsP[[l]]->(VarsP[[l]]/.SetOfPointsr1[[i]])];
];
For[l=11,l≤20,l++,
AppendTo[newpoint,VarsP[[l]]->-(VarsP[[l]]/.SetOfPointsr1[[i]])];
];
AppendTo[SetOfPointsr2,newpoint];
];

We are now looking for tangency conditions for r1 and r2, including double tangency. We want to have a parametric equation to order 2 of normal curves at their computed points. The accuracy is 1500.

ParamAtr1={};
Block[{$MaxExtraPrecision=200},
For[i=1,i≤30,i++,
If[Mod[i,30]==0,Print[i]];
pt=(VarsP/.SetOfPointsr1[[i]])+ t{0,0,1,dP3,dP4,dP5,dP6,dP7,dP8,dP9,dP10,dP11,dP12, dP13,dP14, dP15, dP16, dP17, dP18, dP19}
+t^2{0,0,0,ddP3,ddP4,ddP5,ddP6,ddP7,ddP8,ddP9,ddP10,ddP11,ddP12, ddP13,ddP14, ddP15, ddP16, ddP17, ddP18, ddP19};
VarsdP={dP3,dP4,dP5,dP6,dP7,dP8,dP9,dP10,dP11,dP12, dP13,dP14, dP15, dP16, dP17, dP18, dP19,ddP3,ddP4,ddP5,ddP6,ddP7,ddP8,ddP9,ddP10,ddP11,ddP12, ddP13,ddP14, ddP15, ddP16, ddP17, ddP18, ddP19};
pts=Table[VarsP[[k]]->pt[[k]],{k,2,20}];
Eqs=Table[InputEqs[[k]]/.P0->1,{k,1,100}];
ToSolve=Join[Table[ SeriesCoefficient[Eqs[[GoodIndices[[k-1]]]]/.pts,{t,0,1}],{k,2,18}],
Table[SeriesCoefficient[Eqs[[GoodIndices[[k-1]]]]/.pts,{t,0,2}],{k,2,18}]];
sols=NSolve[ToSolve][[1]];
fr=FindRoot[ToSolve,Table[{VarsdP[[k]],VarsdP[[k]]/.sols},{k,1,34}],MaxIterations->200,WorkingPrecision->1500,AccuracyGoal->2000];
AppendTo[ParamAtr1, Join[pts/.fr,{P0->1}]];
];
];

ParamAtr2={};
Block[{$MaxExtraPrecision=200},
For[i=1,i≤30,i++,
If[Mod[i,30]==0,Print[i]];
pt=(VarsP/.SetOfPointsr2[[i]])+ t{0,0,1,dP3,dP4,dP5,dP6,dP7,dP8,dP9,dP10,dP11,dP12, dP13,dP14, dP15, dP16, dP17, dP18, dP19}
+t^2{0,0,0,ddP3,ddP4,ddP5,ddP6,ddP7,ddP8,ddP9,ddP10,ddP11,ddP12, ddP13,ddP14, ddP15, ddP16, ddP17, ddP18, ddP19};
VarsdP={dP3,dP4,dP5,dP6,dP7,dP8,dP9,dP10,dP11,dP12, dP13,dP14, dP15, dP16, dP17, dP18, dP19,ddP3,ddP4,ddP5,ddP6,ddP7,ddP8,ddP9,ddP10,ddP11,ddP12, ddP13,ddP14, ddP15, ddP16, ddP17, ddP18, ddP19};
pts=Table[VarsP[[k]]->pt[[k]],{k,2,20}];
Eqs=Table[InputEqs[[k]]/.P0->1,{k,1,100}];
ToSolve=Join[Table[ SeriesCoefficient[Eqs[[GoodIndices[[k-1]]]]/.pts,{t,0,1}],{k,2,18}],
Table[SeriesCoefficient[Eqs[[GoodIndices[[k-1]]]]/.pts,{t,0,2}],{k,2,18}]];
sols=NSolve[ToSolve][[1]];
fr=FindRoot[ToSolve,Table[{VarsdP[[k]],VarsdP[[k]]/.sols},{k,1,34}],MaxIterations->200,WorkingPrecision->1500,AccuracyGoal->2000];
AppendTo[ParamAtr2, Join[pts/.fr,{P0->1}]];
];
];

Now calculating the section of O(4K) which vanishes three times at r1=0. These should be invariant sections with respect to the C3 action. We will start by writing a generic polynomial which is invariant with respect to C3 action by evaluating at points of SetOfPointsP.

AllMons={};
For[i=1,i≤20,i++,For[j=1,j≤i,j++,AppendTo[AllMons,VarsP[[i]]VarsP[[j]]];]];
CoeffList=Table[ToExpression[StringJoin["A",ToString[i]]],{i,1,500}];
ListPolySym={};
rk=0;
PolySymValues={};
For[i=1,i≤Length[AllMons],i++,
If[Mod[i,30]==0,Print[i," ",Length[ListPolySym]," ",rk]];
newpoly=(1/3)(AllMons[[i]] + (AllMons[[i]]/.PC3Sub) + ((AllMons[[i]]/.PC3Sub)/.PC3Sub));
NewPolySymValues=Append[PolySymValues, newpoly/.SetOfPointsP];
If[MatrixRank[NewPolySymValues]==rk+1,
rk++;
PolySymValues=NewPolySymValues;
AppendTo[ListPolySym,newpoly];
If[rk==38, i=500];
];
];

Block[{$MaxExtraPrecision=200},
UnivPoly=Sum[CoeffList[[k]] ListPolySym[[k]],{k,1,Length[ListPolySym]}];
Eqs={};
For[k=1,k≤20,k++,
AppendTo[Eqs, 0==UnivPoly/.ParamAtr1[[k]]/.t->0];
AppendTo[Eqs, 0==SeriesCoefficient[UnivPoly/.ParamAtr1[[k]],{t,0,1}]];
AppendTo[Eqs, 0==SeriesCoefficient[UnivPoly/.ParamAtr1[[k]],{t,0,2}]];
];
sols=Chop[Solve[Eqs][[1]],10^(-1800)];
];

Block[{$MaxExtraPrecision=200},
UnivPolySolved=UnivPoly/.sols;
For[i=1,i≤38,i++,
If[ByteCount[Coefficient[UnivPolySolved,CoeffList[[i]]]]>100,
Cgood=CoeffList[[i]]; i=100;
];
];
UnivPolySolved=Sum[Beautify[Coefficient[CoeffList[[i]]/.sols,Cgood]] ListPolySym[[i]],{i,1,38}];
];

Block[{$MaxExtraPrecision=200},
UnivPoly=Sum[CoeffList[[k]] ListPolySym[[k]],{k,1,Length[ListPolySym]}];
Eqs={};
For[k=1,k≤20,k++,
AppendTo[Eqs, 0==UnivPoly/.ParamAtr2[[k]]/.t->0];
AppendTo[Eqs, 0==SeriesCoefficient[UnivPoly/.ParamAtr2[[k]],{t,0,1}]];
AppendTo[Eqs, 0==SeriesCoefficient[UnivPoly/.ParamAtr2[[k]],{t,0,2}]];
];
sols=Chop[Solve[Eqs][[1]]];

UnivPolySolved2=UnivPoly/.sols;
For[i=1,i≤38,i++,
If[ByteCount[Coefficient[UnivPolySolved2,CoeffList[[i]]]]>100,
Cgood=CoeffList[[i]]; i=100;
];
];
UnivPolySolved2=Sum[Beautify[Coefficient[CoeffList[[i]]/.sols,Cgood]] ListPolySym[[i]],{i,1,38}];
];

Simplify[UnivPolySolved2-(UnivPolySolved/.Table[VarsP[[i]]->-VarsP[[i]],{i,11,20}])]

Computing 6 sections of 2K+[r2]-[r1] as sections of O(4K) on double cover of FPP which vanish twice on r2 and also have the weight w^2 (same as s3) with respect to C3 action.

AllMons={};
For[i=1,i≤20,i++,For[j=1,j≤i,j++,AppendTo[AllMons,VarsP[[i]]VarsP[[j]]];]];
CoeffList=Table[ToExpression[StringJoin["A",ToString[i]]],{i,1,500}];
(* finding a basis of the C3 eigenspace *)
ListPolySym={};
rk=0;
PolySymValues={};
For[i=1,i≤Length[AllMons],i++,
If[Mod[i,30]==0,Print[i," ",Length[ListPolySym]," ",rk]];
newpoly=(1/3)(AllMons[[i]] + w (AllMons[[i]]/.PC3Sub) + w^2((AllMons[[i]]/.PC3Sub)/.PC3Sub));
NewPolySymValues=Append[PolySymValues, newpoly/.SetOfPointsP];
If[MatrixRank[NewPolySymValues]==rk+1,
rk++;
PolySymValues=NewPolySymValues;
AppendTo[ListPolySym,newpoly];
If[rk==36,i=400];
];
];

Block[{$MaxExtraPrecision=500},
UnivPoly=Sum[CoeffList[[k]] ListPolySym[[k]],{k,1,Length[ListPolySym]}];
Eqs={};
For[k=1,k≤30,k++,
AppendTo[Eqs, 0==UnivPoly/.ParamAtr2[[k]]/.t->0];
AppendTo[Eqs, 0==SeriesCoefficient[UnivPoly/.ParamAtr2[[k]],{t,0,1}]];
];
sols=Chop[Solve[Eqs][[1]]];

UnivMoreSec=UnivPoly/.sols;
SecWt1={};
For[i=1,i≤36,i++,
If[ByteCount[Coefficient[UnivMoreSec,CoeffList[[i]]]]>100,
AppendTo[SecWt1,Sum[Beautify[Coefficient[CoeffList[[k]]/.sols,CoeffList[[i]]]] ListPolySym[[k]],{k,1,36}]];
Print[Length[SecWt1]];
];
];
];

We could find the double cover as well, but since we are interested in the FPP, let’s do just that first.

(* now finding more points *)


NumPoints=100;
For[i=1,i≤NumPoints,i++,
rpt1=(RandomInteger[1000]+I RandomInteger[1000])/100;
rpt2=(RandomInteger[1000]+I RandomInteger[1000])/100;
sub3=N[{P0->1,P1->rpt1,P2->rpt2},5500];
Eqs=Table[0==EqsFPPD3X3[[k]]/.sub3,{k,1,84}];
sols=FindInstance[N[Eqs,1200],Table[VarsP[[i]],{i,4,10}]];
Eqs1=Table[0==N[InputEqs[[k]],1200]/.Join[sub3,sols[[1]]],{k,1,100}];
sols1=FindInstance[N[Eqs1,1200],Table[VarsP[[i]],{i,11,20}]];
sols17=Join[sols1[[1]],sols[[1]]];
rk=0;
GoodEqs=Table[InputEqs[[GoodIndices[[k]]]]/.sub3,{k,1,17}];
fr= FindRoot[GoodEqs,Table[{VarsP[[k]],VarsP[[k]]/.sols17},{k,4,20}],MaxIterations->200,WorkingPrecision->5000,AccuracyGoal->5000];
fr=Join[fr,sub3];
If[Sum[Abs[InputEqs[[k]]/.fr],{k,1,100}]<10^(-1400), (* quality control *)
AppendTo[SetOfPointsP,fr];
];
If[Mod[i,30]==0,Print[{i, Length[SetOfPointsP]}]];
];

(* We realized that the accuracy was not sufficient, so we needed to go beyond 5000 digits *)
SetOfPointsPnew={};
prec=20000;
For[numpt=1,numpt≤Length[SetOfPointsP],numpt++,
If[Mod[numpt,30]==0,Print[numpt];];
sub3={P0->N[1,prec+500],P1->N[RootApproximant[P1/.SetOfPointsP[[numpt]]],prec+500],P2->N[RootApproximant[P2/.SetOfPointsP[[numpt]]],prec+500]};
GoodEqs=Table[InputEqs[[GoodIndices[[k]]]]/.sub3,{k,1,17}];
fr= FindRoot[GoodEqs,Table[{VarsP[[k]],VarsP[[k]]/.SetOfPointsP[[numpt]]},{k,4,20}],MaxIterations->200,WorkingPrecision->prec,AccuracyGoal->prec];
fr=Join[fr,sub3];
If[Sum[Abs[InputEqs[[k]]/.fr],{k,1,100}]<10^(-1400), (* quality control *)
AppendTo[SetOfPointsPnew,fr];
];
];

SetOfPointsP=SetOfPointsPnew;
Beautify[x_]:=
Ret=Expand[RootApproximant[Chop[N[Re[x],19000],10^(-19000)],2]+RootApproximant[Chop[N[Im[x]/(Sqrt[15]),19000],10^(-19000)],2] I(Sqrt[15])];

SecWt1new=SecWt1;
VarsQ={Q0,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9};
SetOfPointsX9={};
For[numpt=1,numpt<Length[SetOfPointsP],numpt++,
If[Mod[numpt,40]==0,Print[numpt]];
newpt={P0,P1+P2+P3,P4+P5+P6,P7+P8+P9}/.SetOfPointsP[[numpt]];
zvalue=( (UnivPolySolved/UnivPolySolved2)/.SetOfPointsP[[numpt]])^(1/3);
If[Abs[ zvalue  Sections3/.SetOfPointsP[[numpt]]]>10^(-100),
For[i=1,i≤6,i++,
new1=(SecWt1new[[i]] zvalue Sections3^(-1))/.SetOfPointsP[[numpt]];
new2=((SecWt1new[[i]] zvalue^(-1) Sections3^(-1))/.{P10->-P10,P11->-P11,P12->-P12,P13->-P13,P14->-P14,P15->-P15,P16->-P16,P17->-P17,P18->-P18,P19->-P19})/.SetOfPointsP[[numpt]];
AppendTo[newpt,(new1+new2)/2]
];
AppendTo[SetOfPointsX9,Table[VarsQ[[i]]->newpt[[i]],{i,1,10}]];
];
];
Accuracy[SetOfPointsX9]

AllMons={};
For[i=1,i≤10,i++,For[j=1,j≤i,j++,For[k=1,k≤j,k++,AppendTo[AllMons,VarsQ[[i]]VarsQ[[j]] VarsQ[[k]]];]]];

EqsFPPX9={};

UnivPolyQ=Sum[AllMons[[i]] CoeffList[[i]],{i,1,Length[AllMons]}];
Eqs=Table[0==UnivPolyQ/.SetOfPointsX9[[k]],{k,1,Length[SetOfPointsX9]}];
Print["solving"];
sols=Solve[Eqs][[1]];
Print["solved"];
UnivEqsNew=UnivPolyQ/.sols;

ByteCount[UnivEqsNew]

(* this took a couple of days *)
EqsFPPX9Withrr={};
EqsFPPX9={};
For[i=1,i≤Length[AllMons],i++,
If[Mod[i,20]==0,Print[i]];
maybeneweq=Coefficient[UnivEqsNew,CoeffList[[i]]];
If[ByteCount[maybeneweq]>100,
neweq=0;
neweqwithrr=0;
For[k=1,k≤220,k++,
Block[
{$MaxExtraPrecision=2000},
ToMakeDeg4=1+(-3)^(1/2) + (-15)^(1/2) + 5^(1/2);
cf = Chop[Coefficient[maybeneweq,AllMons[[k]]],10^(-10000)] + ToMakeDeg4;
ft=MinimalPolynomial[RootApproximant[cf,4],t];
aa1=-Coefficient[ft,t^3]/Coefficient[ft,t^4]/4-1;
ft=MinimalPolynomial[RootApproximant[cf/5^(1/2),4],t];
aa2=-Coefficient[ft,t^3]/Coefficient[ft,t^4]/4-1;
ft=MinimalPolynomial[RootApproximant[cf/(-3)^(1/2),4],t];
aa3=-Coefficient[ft,t^3]/Coefficient[ft,t^4]/4-1;
ft=MinimalPolynomial[RootApproximant[cf/(-15)^(1/2),4],t];
aa4=-Coefficient[ft,t^3]/Coefficient[ft,t^4]/4-1;
If[Abs[Coefficient[maybeneweq,AllMons[[k]]]-{aa1,aa2,aa3,aa4}.{1,5^(1/2),(-3)^(1/2),(-15)^(1/2)}]>10^(-10000),Print[{"warning ",i, k}];];
neweq+={aa1,aa2,aa3,aa4}.{1,5^(1/2),(-3)^(1/2),(-15)^(1/2)} AllMons[[k]];
neweqwithrr+={aa1,aa2,aa3,aa4}.{1,rr1,rr2,rr3} AllMons[[k]];
];
];
AppendTo[EqsFPPX9,neweq];
AppendTo[EqsFPPX9Withrr,neweqwithrr];
Print[Length[EqsFPPX9]," equation(s)"];
];
];

EqsFPPX9good=EqsFPPX9;
EqsFPPX9Withrrgood=EqsFPPX9Withrr;
For[i=1,i≤84,i++,
If[Mod[i,20]==5,Print[i]];
EqsFPPX9good[[i]]=Numerator[Factor[EqsFPPX9[[i]]]];
EqsFPPX9Withrrgood[[i]]=Numerator[Factor[EqsFPPX9Withrr[[i]]]];
];

FactorInteger[1231]
For[i=1,i≤1231,i++,
If[Mod[i^2+3,1231]==0,Print[i]];
];
For[i=1,i≤1231,i++,
If[Mod[i^2-5,1231]==0,Print[i]];
];
Mod[253 539,1231]

EqsFPPX9Mod1231=Table[Sum[Mod[Coefficient[EqsFPPX9Withrrgood[[i]]/.{rr1->539,rr2->253,rr3->957},AllMons[[j]]],1231] AllMons[[j]],{j,1,Length[AllMons]}],{i,1,84}];
ByteCount[EqsFPPX9good]
ByteCount[EqsFPPX9Withrrgood]
ByteCount[EqsFPPX9Mod1231]

fname="EqsFPPX9Magma";
OpenWrite[fname];
Write[fname,EqsFPPX9Withrrgood];
Close[fname];

fname="EqsFPPX9mod1231";
OpenWrite[fname];
Write[fname,EqsFPPX9Mod1231];
Close[fname];

fname="EqsFPPX9";
OpenWrite[fname];
Write[fname,EqsFPPX9good];
Close[fname];

Might as well calculate the double cover. We will not do anything with it, so will only save the Mathematica format equations.

(* now finding the double FPP *)
VarsQ={Q0,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,Q10,Q11,Q12,Q13,Q14,Q15,Q16,Q17,Q18,Q19};
C2actionQ={Q0->Q0,Q1->Q1,Q2->Q2,Q3->Q3,Q4->-Q4,Q5->-Q5,Q6->-Q6,Q7->-Q7,Q8->Q9, Q9->Q8,Q10->Q11,Q11->Q10,Q12->Q13,Q13->Q12,Q14->Q15,Q15->Q14,Q16->Q17,Q17->Q16,Q18->Q19,Q19->Q18};

(* we will build a list of points on the new surface *)
SetOfPointsQ={};
For[numpt=1,numpt≤150,numpt++,
If[Mod[numpt,40]==0,Print[numpt]];
newpt={P0,P1+P2+P3,P4+P5+P6,P7+P8+P9,P10,P11+P12+P13,P14+P15+P16,P17+P18+P19}/.SetOfPointsP[[numpt]];
zvalue=( (UnivPolySolved/UnivPolySolved2)/.SetOfPointsP[[numpt]])^(1/3);
For[i=1,i≤6,i++,
AppendTo[newpt,(SecWt1new[[i]] zvalue Sections3^(-1))/.SetOfPointsP[[numpt]]];
AppendTo[newpt, ((SecWt1new[[i]] zvalue^(-1) Sections3^(-1))/.{P10->-P10,P11->-P11,P12->-P12,P13->-P13,P14->-P14,P15->-P15,P16->-P16,P17->-P17,P18->-P18,P19->-P19})/.SetOfPointsP[[numpt]]];
];
AppendTo[SetOfPointsQ,Table[VarsQ[[i]]->newpt[[i]],{i,1,20}]];
];
Accuracy[SetOfPointsQ]

AllMons={};
For[i=1,i≤20,i++,For[j=1,j≤i,j++,AppendTo[AllMons,VarsQ[[i]]VarsQ[[j]]];]];

EqsFPPnewDouble={};
For[mod3=0,mod3≤2,mod3++,

degreesinZ={0,0,0,0,0,0,0,0,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1};

AllMons1={};
For[i=1,i≤Length[AllMons],i++,
If[Mod[Sum[Exponent[AllMons[[i]],VarsQ[[j]]] degreesinZ [[j]],{j,1,20}],3]==mod3,
AppendTo[AllMons1,AllMons[[i]]];
];
];
UnivPolyQ=Sum[AllMons1[[i]] CoeffList[[i]],{i,1,Length[AllMons1]}];
Eqs=Table[0==UnivPolyQ/.SetOfPointsQ[[k]],{k,1,Length[SetOfPointsQ]}];
Print[mod3," ", Length[AllMons1]];
sols=Solve[Eqs][[1]];
Print["solved"];
UnivEqsNew=UnivPolyQ/.sols;
For[i=1,i≤Length[AllMons1],i++,
If[ByteCount[Coefficient[UnivEqsNew,CoeffList[[i]]]]>100,
Block[{$MaxExtraPrecision=2000},AppendTo[EqsFPPnewDouble,Sum[RootApproximant[Chop[Coefficient[CoeffList[[k]]/.sols,CoeffList[[i]]],10^(-10000)],4] AllMons1[[k]],{k,1,Length[AllMons1]}]];
];
If[Mod[Length[EqsFPPnewDouble],10]==0,Print[Length[EqsFPPnewDouble]," equations"];];
];
];
];

Length[EqsFPPnewDouble]
Table[ByteCount[EqsFPPnewDouble[[k]]],{k,1,100}]
Table[Length[MonomialList[EqsFPPnewDouble[[k]]]],{k,1,100}]

(* write to file *)


fname="EqsFPPX9double";
OpenWrite[fname];
Write[fname,EqsFPPnewDouble];
Close[fname];