This notebook is a supporting file for Section 6 of the paper. It contains the calculation of the twin fake projective planes
(d^2D)_3X_3 and (dD)_3X_3.

We start with the double cover of d3D3 FPP and use the C3 action which is a combination of covering and lift actions.

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

(* this defines the mixed C3 action. The other twin will have w replaced by w^2 *)
PC3Sub={P1->P2,P2->P3,P3->P1,P11->P12,P12->P13,P13->P11,P4->w P5,P5->w P6,P6->w P4,P14->w P15,P15->w P16,P16->w P14,P7->w^2 P8,P8->w^2 P9,P9->w^2 P7,P17->w^2 P18,P18->w^2 P19,P19->w^2 P17};
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. We will write points in the 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];
];
];

Length[SetOfPointsP]
fname="InputSurfacePoints";
OpenWrite[fname];
Write[fname,SetOfPointsP];
Close[fname];

fname="InputSurfacePoints";
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. The normalization line 
Sections4=Simplify[(D[UnivLinearSolved,c14]/(D[UnivLinearSolved,c14] needed to
be modified to get a sensible ratio, and same for the next line.

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

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,39}];
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->0,c10->1,c11->0}))/.CtoBasis];
Sections3=Simplify[(D[UnivLinearSolved,c10]/(D[UnivLinearSolved,c10]/.{c12->0,c13->0,c14->1}))/.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];

Sections1
Sections2
Sections3
Sections4

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

As a trick, we find a set of equations that are algebraically independent once and then use them. Points for r2 are obtained by C2 from points of r1.

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

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

I am now looking for tangency conditions for r1 and r2, including double tangency. I 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. I 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(2K) on double cover of FPP which vanish twice on r2 and also has 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}];
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]

NewSym={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};
MatrixRank[Join[SecWt1,SecWt1/.NewSym]/.SetOfPointsP]

(* as a nice bonus, the SecWt1 elements are an eigenbasis for the NewSym action *)
For[k=1,k≤6,k++,Print[FullSimplify[(SecWt1[[k]]/.NewSym)/SecWt1[[k]]]];
];

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 now compute the equations of the double cover of a twin FPP.

VarsQ={Q0,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9,Q10,Q11,Q12,Q13,Q14,Q15,Q16,Q17,Q18,Q19};
SetOfPointsQ={};
For[numpt=1,numpt<Length[SetOfPointsP],numpt++,
If[Mod[numpt,40]==0,Print[numpt]];
newpt=Join[BasisWt00,BasisWt10]/.SetOfPointsP[[numpt]];
zvalue=( (UnivPolySolved/UnivPolySolved2)/.SetOfPointsP[[numpt]])^(1/3);
For[i=1,i≤6,i++,
AppendTo[newpt,(SecWt1[[i]] zvalue Sections3^(-1))/.SetOfPointsP[[numpt]]];
AppendTo[newpt,((SecWt1[[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"];
];

(* there seemed to be some precisions issues, so we will do a sanity check *)
Sum[Sum[Abs[EqsFPPnewDouble[[k]]/.SetOfPointsQ[[i]]],{i,1,Length[SetOfPointsQ]}],{k,1,100}]

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

(* now getting ready to save it in Magma friendly format *)
RationalPart[x_]:=
(minp=MinimalPolynomial[x,t];
If[Coefficient[minp,t^2]≠0,bb=-Coefficient[minp,t]/Coefficient[minp,t^2]/2];
If[Coefficient[minp,t^2]==0,bb=-(minp/.t->0)/Coefficient[minp,t]];
bb)
QuadRoot5To2[x_]:={RationalPart[x],RationalPart[(x-RationalPart[x])/5^(1/2)]}
ToFour[x_]:=Join[QuadRoot5To2[Re[x]],QuadRoot5To2[Im[x]/3^(1/2)]]

(* testing ToFour *)
cf=Coefficient[EqsFPPnewDouble[[1]],Q15 Q16]
cf- ToFour[cf].{1, (Sqrt[5]),I (Sqrt[3]),I (Sqrt[15])}

(* it didn't work Mod 79, so we will use another Modulus: 661 *) 
For[i=1,i≤661,i++,If[Mod[i^2+3,661]==0,Print[i]]];
For[i=1,i≤661,i++,If[Mod[i^2-5,661]==0,Print[i]]];
Mod[68 115,661]

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

EqsFPPTwinDoubleWithrr={};
For[i=1,i≤100,i++,
newrel=0;
For[j=1,j≤Length[AllMons],j++,
cf=Coefficient[EqsFPPnewDouble[[i]] ,AllMons[[j]]];
newrel+=(ToFour[cf].{1,rr1,rr2,rr3}) AllMons[[j]] ;
];
AppendTo[EqsFPPTwinDoubleWithrr,newrel];
];
EqsFPPTwinDoubleMod661=Table[Sum[Mod[Coefficient[EqsFPPTwinDoubleWithrr[[i]]/.{rr1->115,rr2->68,rr3->549},AllMons[[j]]],661] AllMons[[j]],{j,1,Length[AllMons]}],{i,1,100}];

fname="EqsFPPTwinDoubleMagma";
OpenWrite[fname];
Write[fname,EqsFPPTwinDoubleWithrr];
Close[fname];

fname="EqsFPPTwinDoublemod661";
OpenWrite[fname];
Write[fname,EqsFPPTwinDoubleMod661];
Close[fname];

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

Now we will look for the equations of the FPP.

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

(* now compute the equations *)
(* Beautify command takes a surprisingly long time, so we do something else that was reasonably fast. 
The assumption is that we are in the expected field. There may well be a less ridiculous way to do it. *)
SetOfPointsTwin={};
For[numpt=1,numpt<Length[SetOfPointsP],numpt++,
If[Mod[numpt,40]==0,Print[numpt]];
newpt=BasisWt00/.SetOfPointsP[[numpt]];
zvalue=( (UnivPolySolved/UnivPolySolved2)/.SetOfPointsP[[numpt]])^(1/3);
For[i=1,i≤6,i++,
new1=(SecWt1[[i]] zvalue Sections3^(-1))/.SetOfPointsP[[numpt]];
new2=((SecWt1[[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[SetOfPointsTwin,Table[VarsQ[[i]]->newpt[[i]],{i,1,10}]];
];
Accuracy[SetOfPointsTwin]

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

EqsFPPTwin={};
EqsFPPTwinWithrr={};

UnivPolyQ=Sum[AllMons[[i]] CoeffList[[i]],{i,1,Length[AllMons]}];
Eqs=Table[0==UnivPolyQ/.SetOfPointsTwin[[k]],{k,1,Length[SetOfPointsTwin]}];
Print["solving"];
sols=Solve[Eqs][[1]];
Print["solved"];
UnivEqsNew=UnivPolyQ/.sols;
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++,
ToMakeDeg4=1+(-3)^(1/2) + (-15)^(1/2) + 5^(1/2);
cf = Chop[Coefficient[maybeneweq,AllMons[[k]]],10^(-1800)] + 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;
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[EqsFPPTwin,neweq];
AppendTo[EqsFPPTwinWithrr,neweqwithrr];
];
];
Print[Length[EqsFPPTwin]," equations"];

ByteCount[EqsFPPTwin]
ByteCount[EqsFPPTwinWithrr]

For[i=1,i≤84,i++,
If[Mod[i,20]==5,Print[i]];
EqsFPPTwin[[i]]=Numerator[Factor[EqsFPPTwin[[i]]]];
EqsFPPTwinWithrr[[i]]=Numerator[Factor[EqsFPPTwinWithrr[[i]]]];
];

ByteCount[EqsFPPTwin]
ByteCount[EqsFPPTwinWithrr]

EqsFPPTwinMod661=Table[Sum[Mod[Coefficient[EqsFPPTwinWithrr[[i]]/.{rr1->115,rr2->68,rr3->549},AllMons[[j]]],661] AllMons[[j]],{j,1,Length[AllMons]}],{i,1,84}];

ByteCount[EqsFPPTwinMod661]

fname="EqsFPPTwinMagma";
OpenWrite[fname];
Write[fname,EqsFPPTwinWithrr];
Close[fname];

fname="EqsFPPTwinmod661";
OpenWrite[fname];
Write[fname,EqsFPPTwinMod661];
Close[fname];

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