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

The main idea is this. Starting from a double cover of FPP with a C6 automorphism group (C2 for the cover and C3 for the FPP), we attempt to construct a double cover of FPP with S3 automorphism group (C2 quotient gives FPP).

The variables are  P0,...,P19, with the first ten even and the last ten odd. The action refers to an order three action. 

We first construct the surface D3X3 and its double cover.

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

InputSurfaceFileName="EqsFPPd3D3double";

The first step is to read the input file and verify that the action is OK. If it is,
then the two ranks would be both 100. Then we compute a number of points on this surface.

VarsP = {P0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14,P15,P16,P17,P18,P19};
OpenRead[InputSurfaceFileName];
InputEqs=Read[InputSurfaceFileName];
Close[InputSurfaceFileName];

(* this is the action of C3 on this surface *)
w=Exp[2 Pi I /3];
PC3Sub={P4->w P4,P5->w P5,P6->w P6,P14->w P14,P15->w P15,P16->w P16,P7->P7/w ,P8->P8/w ,P9->P9/w ,P17->P17/w ,P18->P18/w ,P19->P19/w };

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]]],300],{i,1,100},{j,1,210}]]
MatrixRank[Join[Table[N[Coefficient[InputEqs[[i]],AllMons[[j]]],300],{i,1,100},{j,1,210}],
Table[N[Coefficient[InputEqs[[i]]/.PC3Sub,AllMons[[j]]],300],{i,1,100},{j,1,210}]]]

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++,
If[Mod[i,30]==0,Print[i]];
rpt1=(RandomInteger[1000]+I RandomInteger[1000])/100;
rpt2=(RandomInteger[1000]+I RandomInteger[1000])/100;
sub3={P0->1,P1->rpt1,P2->rpt2};
Eqs=Table[0==InputEqs[[k]]/.sub3,{k,1,100}];
sols=NSolve[Eqs,200];
(* we will look for equations with a nondegenerate Jacobian matrix at sols[1] that will allow us to find a point with high precision *)
GoodEqs={};
rk=0;
EqsP0=InputEqs/.sub3;
JacMatAtp={};
For[k=1,k≤100,k++,
JacMatAtpnew=Append[JacMatAtp,Table[D[EqsP0[[k]],VarsP[[l]]]/.sols[[1]],{l,4,20}]];
If[MatrixRank[JacMatAtpnew]==rk+1,
JacMatAtp=JacMatAtpnew;rk++;AppendTo[GoodEqs,EqsP0[[k]]];
If[rk==17,k=200];
];
];
fr= FindRoot[Table[GoodEqs[[k]],{k,1,17}],Table[{VarsP[[k]],VarsP[[k]]/.sols[[1]]},{k,4,20}],MaxIterations->200,WorkingPrecision->2000,AccuracyGoal->2000];
fr=Join[fr,sub3];
If[Sum[Abs[InputEqs[[k]]/.fr],{k,1,100}]<10^(-1900), (* quality control *)
AppendTo[SetOfPointsP,fr];
];
];

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

(* reading from a file *)
fname="FPPd3D3DoublePoints";
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.

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[RootApproximant[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[prod[[2]]/.CtoBasis]
Sections2=Simplify[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={};
rpt1=(RandomInteger[1000]+I RandomInteger[1000])/1000;
RelsToSolve=Join[Table[InputEqs[[i]]==0,{i,1,100}],{0==Sections1,0==Sections3,P1==rpt1,P0==1}];
sols=NSolve[RelsToSolve,500][[1]];
Print["here"];
rk=0;
EqsP0=Join[{Sections1},InputEqs]/.{P0->1,P1->rpt1};
JacMatAtp={};
GoodEqsIndices={};
For[k=1,k≤100,k++,
JacMatAtpnew=Append[JacMatAtp,Table[D[EqsP0[[k]],VarsP[[l]]]/.sols,{l,3,20}]];
If[MatrixRank[N[JacMatAtpnew,500]]==rk+1,
JacMatAtp=JacMatAtpnew;rk++;AppendTo[GoodEqsIndices,k];If[rk==18,k=200];
];
];

Print[GoodEqsIndices];  (* as a trick, we find set of algebraically independent equations once and then reuse it *)

For[numpt=1,numpt≤100,numpt++,
If[Mod[numpt,10]==0,Print[numpt," ",Length[SetOfPointsr1]]];
rpt1=(RandomInteger[1000]+I RandomInteger[1000])/1000;
RelsToSolve=Join[Table[InputEqs[[i]]==0,{i,1,100}],{0==Sections1,0==Sections3,P1==rpt1,P0==1}];
sols=NSolve[RelsToSolve,100][[1]];
rk=0;
EqsP0=Join[{Sections1},InputEqs]/.{P0->1,P1->rpt1};
JacMatAtp={};
GoodEqs=Table[EqsP0[[GoodEqsIndices[[k]]]],{k,1,18}];
fr= FindRoot[Table[GoodEqs[[k]]/.{P0->1,P1->rpt1},{k,1,18}],Table[{VarsP[[k]],VarsP[[k]]/.sols},{k,3,20}],MaxIterations->200,WorkingPrecision->2000,AccuracyGoal->2000];
AppendTo[fr,P0->1];
AppendTo[fr,P1->rpt1];
AppendTo[SetOfPointsr1,fr];
];

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={};
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[[GoodEqsIndices[[k]]-1]]/.pts,{t,0,1}],{k,2,18}],
Table[SeriesCoefficient[Eqs[[GoodEqsIndices[[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={};
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[[GoodEqsIndices[[k]]-1]]/.pts,{t,0,1}],{k,2,18}],
Table[SeriesCoefficient[Eqs[[GoodEqsIndices[[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];
];
];

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

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

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[RootApproximant[Coefficient[CoeffList[[i]]/.sols,Cgood]] ListPolySym[[i]],{i,1,38}];

Simplify[UnivPolySolved- (UnivPolySolved2/.{P10->-P10,P11->-P11,P12->-P12,P13->-P13,P14->-P14,P15->-P15,P16->-P16,P17->-P17,P18->-P18,P19->-P19})]

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

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[RootApproximant[Coefficient[CoeffList[[k]]/.sols,CoeffList[[i]]]] ListPolySym[[k]],{k,1,36}]];
];
];
Length[SecWt1]

We will look to find nicer elements in the linear span of SecWt1. Since these are quadratic in VarsP, we will look for linear combinations of them and the relations of InputEqs. Minimality is defined in terms of coefficients of the monomials. To simplify our work, we only want to look at the linear combinations of the equations that have no coefficients outside of the ones in SecWt1.

msecall={};
For[n=1,n≤6,n++,
msec=MonomialList[SecWt1[[n]]];
msec=Table[msec[[i]]/(msec[[i]]/.Table[VarsP[[k]]->1,{k,1,20}]),{i,1,Length[msec]}];
msecall=Union[msecall,msec];
];
msecall
Length[msecall]

(msecall/.PC3Sub)(E^((2 I π)/3))-msecall

(* figuring out relevant monomials and equations *)
RelsThatMatter={};
InputEqsWt1= Simplify[(InputEqs +(InputEqs/.PC3Sub) (E^((2 I π)/3)) +((InputEqs/.PC3Sub)/.PC3Sub)(E^(-((2 I π)/3))))/3];
For[i=1,i≤100,i++,
If[ByteCount[InputEqsWt1[[i]]]>100,
meq=MonomialList[InputEqsWt1[[i]]];
meq=Table[meq[[j]]/(meq[[j]]/.Table[VarsP[[k]]->1,{k,1,20}]),{j,1,Length[meq]}];
msecall=Union[msecall,meq];
AppendTo[RelsThatMatter,InputEqsWt1[[i]]];
];
];

Length[RelsThatMatter]
Length[msecall]

(* We will now go through the usual process of finding small relations *)
AllSecWt1=Join[SecWt1,RelsThatMatter];
BigMatrix={};
For[i=1,i≤Length[AllSecWt1],i++,
AppendTo[BigMatrix,Table[RootReduce[Coefficient[AllSecWt1[[i]]  ,msecall[[j]] ]],{j,1,Length[msecall]}]];
];


For[k=1,k≤Length[BigMatrix],k++,
firsteq=BigMatrix[[k]];
lcm=Denominator[firsteq[[1]]];
For[i=2,i≤Length[firsteq],i++,
lcm=LCM[lcm,Denominator[firsteq[[i]]]];
];
BigMatrix[[k]]=Expand[firsteq lcm];
];


BigMatrixReIm={};
For[i=1,i≤Length[BigMatrix],i++,
NewRow1={};
NewRow2={};
For[j=1,j≤Length[BigMatrix[[1]]],j++,
AppendTo[NewRow1,Re[BigMatrix[[i]][[j]]]];
AppendTo[NewRow1,Im[BigMatrix[[i]][[j]]]/Sqrt[15]];
AppendTo[NewRow2,Im[BigMatrix[[i]][[j]]](-Sqrt[15])];
AppendTo[NewRow2,Re[BigMatrix[[i]][[j]]]];
];
AppendTo[BigMatrixReIm,NewRow1];
AppendTo[BigMatrixReIm,NewRow2];
];


BestRels=LatticeReduce[BigMatrixReIm];

{u,r, v} = SmithDecomposition[BestRels];
vt=Transpose[v];
SmallBestRels={};
For[i=1,i≤Length[BigMatrixReIm],i++,
If[Mod[i,50]==1,Print[i]];
basisvec=Table[0,{k,1, Length[BigMatrixReIm[[1]]]}];
basisvec[[i]]=1;
AppendTo[SmallBestRels,LinearSolve[vt,basisvec]];
];

SmallerBestRels=LatticeReduce[SmallBestRels];

BackToCxMatrix={};
For[i=1,i≤Length[SmallerBestRels],i++,
newrow=Table[SmallerBestRels[[i]][[2k-1]]+ I Sqrt[15] SmallerBestRels[[i]][[2k]],{k,1,Length[SmallerBestRels[[1]]]/2}];
AppendTo[BackToCxMatrix,newrow];
];

BackToCxReduced={};
rk=0;
For[i=1,i≤Length[BackToCxMatrix],i++,
newred=Append[BackToCxReduced,BackToCxMatrix[[i]]];
If[MatrixRank[N[newred,1000]]>rk, rk++;Print[i," ",rk];BackToCxReduced=newred;];
];

(* we will now try to find elements of SecWt1 which generate it and have good C3 action *)

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

maybenonzero={};
For[k=1,k≤Length[BackToCxReduced],k++,
If[Chop[Sum[BackToCxReduced[[k]][[i]] msecall[[i]],{i,1,Length[msecall]}]/.SetOfPointsP[[2]]]≠0,AppendTo[maybenonzero,k];];
];
maybenonzero

GoodSec1 = Sum[BackToCxReduced[[16]][[i]] msecall[[i]],{i,1,Length[msecall]}];
GoodSec2=GoodSec1/.CyclicD3X3;
GoodSec3=GoodSec2/.CyclicD3X3;
MatrixRank[{GoodSec1,GoodSec2,GoodSec3}/.SetOfPointsP]

GoodSec4 = Sum[BackToCxReduced[[17]][[i]] msecall[[i]],{i,1,Length[msecall]}];
GoodSec5=GoodSec4/.CyclicD3X3;
GoodSec6=GoodSec5/.CyclicD3X3;
MatrixRank[{GoodSec1,GoodSec2,GoodSec3,GoodSec4,GoodSec5,GoodSec6}/.SetOfPointsP]

(* we got 6, which is what we needed *)

We will now construct the double of the fake projective plane D3X3.

SecWt1new={GoodSec1,GoodSec2,GoodSec3,GoodSec4,GoodSec5,GoodSec6};
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};

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


(* we will build a list of points on the new surface *)
SetOfPointsQ={};
For[numpt=1,numpt<Length[SetOfPointsP],numpt++,
If[Mod[numpt,40]==0,Print[numpt]];
newpt={P0,P1,P2,P3,P10,P11,P12,P13}/.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,
AppendTo[EqsFPPnewDouble,Sum[Beautify[Coefficient[CoeffList[[k]]/.sols,CoeffList[[i]]]] AllMons1[[k]],{k,1,Length[AllMons1]}]];
];
];
Print[Length[EqsFPPnewDouble]," equations"];
];

Length[EqsFPPnewDouble]
Table[ByteCount[EqsFPPnewDouble[[k]]],{k,1,100}]
For[i=1,i≤100,i++,
EqsFPPnewDouble[[i]]=Numerator[Factor[EqsFPPnewDouble[[i]]]];
];

(* write to file *)
AllMons={};
For[i=1,i≤20,i++,For[j=1,j≤i,j++,AppendTo[AllMons,VarsQ[[i]]VarsQ[[j]]];]];

EqsFPPD3X3DoubleWithrr={};
For[i=1,i≤100,i++,
newrel=0;
For[j=1,j≤Length[AllMons],j++,
cf=Coefficient[EqsFPPnewDouble[[i]] ,AllMons[[j]]];
newrel+=(Re[cf]+Im[cf] (rr/Sqrt[15])) AllMons[[j]] ;
];
AppendTo[EqsFPPD3X3DoubleWithrr,newrel];
];
EqsFPPD3X3DoubleMod79=Table[Sum[Mod[Coefficient[EqsFPPD3X3DoubleWithrr[[i]]/.rr->8,AllMons[[j]]],79]
AllMons[[j]],{j,1,Length[AllMons]}],{i,1,84}];


fname="EqsFPPD3X3DoubleMagma";
OpenWrite[fname];
Write[fname,EqsFPPD3X3DoubleWithrr];
Close[fname];

fname="EqsFPPD3X3Doublemod79";
OpenWrite[fname];
Write[fname,EqsFPPD3X3DoubleMod79];
Close[fname];

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

(* making sure the action is correct *)
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};
MatrixRank[Join[Table[N[Coefficient[EqsFPPnewDouble[[l]],AllMons[[k]]],2000],{l,1,100},{k,1,Length[AllMons]}],
Table[N[Coefficient[EqsFPPnewDouble[[l]]/.C2actionQ,AllMons[[k]]],2000],{l,1,100},{k,1,Length[AllMons]}]]]
MatrixRank[Join[Table[N[Coefficient[EqsFPPnewDouble[[l]],AllMons[[k]]],2000],{l,1,100},{k,1,Length[AllMons]}],
Table[N[Coefficient[EqsFPPnewDouble[[l]]/.C3EqsFPPnewSub,AllMons[[k]]],2000],{l,1,100},{k,1,Length[AllMons]}]]]

We will now write the equations of the FPP D3X3

(* we need more points *)
NumPoints=100;
For[i=1,i≤NumPoints,i++,
If[Mod[i,30]==0,Print[i]];
rpt1=(RandomInteger[1000]+I RandomInteger[1000])/100;
rpt2=(RandomInteger[1000]+I RandomInteger[1000])/100;
sub3={P0->1,P1->rpt1,P2->rpt2};
Eqs=Table[0==InputEqs[[k]]/.sub3,{k,1,100}];
sols=NSolve[Eqs,200];
(* we will look for equations with a nondegenerate Jacobian matrix at sols[1] that will allow us to find a point with high precision *)
GoodEqs={};
rk=0;
EqsP0=InputEqs/.sub3;
JacMatAtp={};
For[k=1,k≤100,k++,
JacMatAtpnew=Append[JacMatAtp,Table[D[EqsP0[[k]],VarsP[[l]]]/.sols[[1]],{l,4,20}]];
If[MatrixRank[JacMatAtpnew]==rk+1,
JacMatAtp=JacMatAtpnew;rk++;AppendTo[GoodEqs,EqsP0[[k]]];
If[rk==17,k=200];
];
];
fr= FindRoot[Table[GoodEqs[[k]],{k,1,17}],Table[{VarsP[[k]],VarsP[[k]]/.sols[[1]]},{k,4,20}],MaxIterations->200,WorkingPrecision->2000,AccuracyGoal->2000];
fr=Join[fr,sub3];
If[Sum[Abs[InputEqs[[k]]/.fr],{k,1,100}]<10^(-1900), (* quality control *)
AppendTo[SetOfPointsP,fr];
];
];

Beautify[x_]:=
Ret=Expand[Rationalize[Chop[N[Re[x],1900],10^(-1800)]]+Rationalize[Chop[N[Im[x]/(Sqrt[15]),1900],10^(-1800)]] I(Sqrt[15])];

SetOfPointsD3X3={};
For[numpt=1,numpt<Length[SetOfPointsP],numpt++,
If[Mod[numpt,40]==0,Print[numpt]];
newpt={P0,P1,P2,P3}/.SetOfPointsP[[numpt]];
zvalue=( (UnivPolySolved/UnivPolySolved2)/.SetOfPointsP[[numpt]])^(1/3);
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[SetOfPointsD3X3,Table[VarsQ[[i]]->newpt[[i]],{i,1,10}]];
];
Accuracy[SetOfPointsD3X3]

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

EqsFPPD3X3={};

UnivPolyQ=Sum[AllMons[[i]] CoeffList[[i]],{i,1,Length[AllMons]}];
Eqs=Table[0==UnivPolyQ/.SetOfPointsD3X3[[k]],{k,1,Length[SetOfPointsD3X3]}];
Print["solving"];
sols=Solve[Eqs][[1]];
Print["solved"];
UnivEqsNew=UnivPolyQ/.sols;
For[i=1,i≤Length[AllMons],i++,
If[ByteCount[Coefficient[UnivEqsNew,CoeffList[[i]]]]>100,
AppendTo[EqsFPPD3X3,Sum[Beautify[Coefficient[CoeffList[[k]]/.sols,CoeffList[[i]]]] AllMons[[k]],{k,1,Length[AllMons]}]];
];
];
Print[Length[EqsFPPD3X3]," equations"];

Table[ByteCount[EqsFPPD3X3[[i]]],{i,1,84}]

For[i=1,i≤84,i++,
EqsFPPD3X3[[i]]=Numerator[Factor[EqsFPPD3X3[[i]]]];
];
Table[ByteCount[EqsFPPD3X3[[i]]],{i,1,84}]

(* writing to the file *)
EqsFPPD3X3Withrr={};
For[i=1,i≤84,i++,
newrel=0;
For[j=1,j≤Length[AllMons],j++,
cf=Coefficient[EqsFPPD3X3[[i]] ,AllMons[[j]]];
newrel+=(Re[cf]+Im[cf] (rr/Sqrt[15])) AllMons[[j]] ;
];
AppendTo[EqsFPPD3X3Withrr,newrel];
];
EqsFPPD3X3Mod79=Table[Sum[Mod[Coefficient[EqsFPPD3X3Withrr[[i]]/.rr->8,AllMons[[j]]],79]
AllMons[[j]],{j,1,Length[AllMons]}],{i,1,84}];

fname="EqsFPPD3X3";
OpenWrite[fname];
Write[fname,EqsFPPD3X3];
Close[fname];

fname="EqsFPPD3X3Magma";
OpenWrite[fname];
Write[fname,EqsFPPD3X3Withrr];
Close[fname];

fname="EqsFPPD3X3mod79";
OpenWrite[fname];
Write[fname,EqsFPPD3X3Mod79];
Close[fname];