This notebook is a supporting file for Section 5 of the paper.

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

fname="Sectionsfandd";
OpenRead[fname];
fandd=Read[fname];
Close[fname];

Sectionf=fandd[[1]];Sectiond=fandd[[2]];

fname="EquationsOfXintermsofW";
OpenRead[fname];
AllEqs23Wgood=Read[fname];
Close[fname];

(* we start by finding some random points on X with high accuracy *)
VarsW={W0,W1,W2,W3,W4,W5,W6,W7};
prec=6000;
SetOfPointsW = {};
numpts=100;
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[{W0==1,W1+W2+W3==rpt1,W4+W5+W6==rpt2},Table[0==AllEqs23Wgood[[j]],{j,1,9}]];
sols=NSolve[Eqs,200];
fr=FindRoot[Table[Eqs[[k]],{k,1,8}],Table[{VarsW[[j]],VarsW[[j]]/.sols[[1]]},{j,1,8}],MaxIterations->100,WorkingPrecision->prec,AccuracyGoal->prec];
AppendTo[SetOfPointsW,fr];
];

(* we will find a nicer formula for f *)

Numf=Sectionf (W4+W5+W6);

maxlength=500;
CoeffList=Table[ToExpression[StringJoin["A",ToString[i]]],{i,1,maxlength}];

AllMons=MonomialList[(W1+W2+W3+W0)^2(W4+W5+W6+W7) + (W4+W5+W6+W7)^3];
count=Length[AllMons]
AllMons=Table[AllMons[[i]]/(AllMons[[i]]/.Table[VarsW[[k]]->1,{k,1,8}]),{i,1,count}]
UnivPoly=Sum[AllMons[[k]] CoeffList[[k]],{k,1,count}];
Eqs=Table[0==(Numf-UnivPoly)/.SetOfPointsW[[k]],{k,1,count+20}];
sols=Chop[ Solve[Eqs][[1]]];
up=N[(UnivPoly/.sols)/.Table[CoeffList[[k]]->0,{k,1,maxlength}],2000];
upexact=Sum[RootApproximant[Coefficient[up,AllMons[[k]]],2]AllMons[[k]],{k,1,count}]

Table[Chop[N[(Sectionf -upexact/(W4+W5+W6) ) /.SetOfPointsW[[k]]]],{k,1,100}]  (* quality control *)

I will find multiple points on C, sigma(C) and sigma^2(C).

(* finding points on C1 *)
SetOfPointsC={};
SetOfPointsSigmaC={};
SetOfPointsSigma2C={};
prec=6000;

For[numpt=1,numpt≤100,numpt++,
If[Mod[numpt,20]==0,Print[numpt]];
rpt1=(RandomInteger[100000]+1)/10000 + I(RandomInteger[100000]+1)/10000;
RelsToSolve=Join[Table[AllEqs23Wgood[[i]]==0,{i,1,9}],{0==upexact,0==Sectiond,W2==rpt1,W4+W5+W6==1}]; (* so that we don't pick up solutions of upexact with W4+W5+W6=0 *)
sols=NSolve[RelsToSolve,600][[1]];
If[Abs[upexact/.sols]<10^(-40), (* to avoid some spurrious results *)
fr= FindRoot[Join[Table[AllEqs23Wgood[[i]]==0,{i,1,5}],{0==Sectiond,W2==rpt1,W4+W5+W6==1}],Table[{VarsW[[k]],(VarsW[[k]]/.sols)},{k,1,8}],AccuracyGoal->prec,WorkingPrecision->(3/2)prec,PrecisionGoal->prec];(* it is better to use Sectiond which doesn't have the multiple root *)
If[Sum[Abs[AllEqs23Wgood[[k]]/.fr],{k,1,9}]≤10^(-prec+1),
AppendTo[SetOfPointsC,N[fr,prec]];
AppendTo[SetOfPointsSigma2C,N[fr,prec]/.{W1->W2,W2->W3,W3->W1,W4->W5,W5->W6,W6->W4}];
AppendTo[SetOfPointsSigmaC,N[fr,prec]/.{W1->W3,W2->W1,W3->W2,W4->W6,W5->W4,W6->W5}];
];
];
];
Print[Length[SetOfPointsC]," points"];

(* sanity check *)
Sectionf= upexact/(W4+W5+W6);
SectionSigmaf=Sectionf/.{W1->W2,W2->W3,W3->W1,W4->W5,W5->W6,W6->W4};
SectionSigma2f=SectionSigmaf/.{W1->W2,W2->W3,W3->W1,W4->W5,W5->W6,W6->W4};
Table[N[SectionSigmaf/.SetOfPointsSigma2C[[k]]],{k,1,100}]
Table[N[SectionSigma2f/.SetOfPointsSigmaC[[k]]],{k,1,100}]
Table[N[(Sectionf SectionSigmaf SectionSigma2f - Sectiond^3)/.SetOfPointsW[[k]]],{k,1,100}]

Defining sigma_lift action and the cube of z.

SubSigma={W1->W2,W2->W3,W3->W1,W4->W5,W5->W6,W6->W4,z->z Sectiond /SectionSigmaf };
zcube = SectionSigmaf /Sectionf;

N[(((z^3/.SubSigma )/.z^3->SectionSigmaf /Sectionf) - (zcube/.SubSigma))/.SetOfPointsW[[45]]]

Finding conditions of tangency at C, sigma(C), sigma^2(C) in the form of a tangent vector to the surface which is normal
to the curve.

TangentsAtC={};
For[i=1,i≤100,i++,
If[Mod[i,20]==0,Print[i]];
pt=(VarsW/.SetOfPointsC[[i]])+ t{0,0,1,dW3,dW4,dW5,dW6,dW7};
pts=Table[VarsW[[k]]->pt[[k]],{k,1,8}];
EqsPlus=Table[AllEqs23Wgood[[i]],{i,1,5}];
sols=Solve[Table[0== SeriesCoefficient[EqsPlus[[j]]/.pts,{t,0,1}],{j,1,5}]][[1]];
AppendTo[TangentsAtC, pts/.sols];
];
TangentsAtSigma2C=TangentsAtC/.{W1->W2,W2->W3,W3->W1,W4->W5,W5->W6,W6->W4};
TangentsAtSigmaC=TangentsAtC/.{W1->W3,W2->W1,W3->W2,W4->W6,W5->W4,W6->W5};

(* sanity check *)
N[Series[Sectionf/.TangentsAtC[[10]],{t,0,2}]]
N[Series[SectionSigma2f/.TangentsAtSigmaC[[10]],{t,0,2}]]
N[Series[SectionSigmaf/.TangentsAtSigma2C[[10]],{t,0,2}]]

AllMons=MonomialList[(W1+W2+W3+W0)^1(W4+W5+W6)^2 + (W0+W1+W2+W3)^3]; (* experimentally found to give one solution up to scaling *)
count=Length[AllMons];
AllMons=Table[AllMons[[i]]/(AllMons[[i]]/.Table[VarsW[[k]]->1,{k,1,8}]),{i,1,count}];

UnivPoly=Sum[CoeffList[[i]] AllMons[[i]],{i,1,Length[AllMons]}];
Eqs = Join[Table[0==UnivPoly/.SetOfPointsC[[i]],{i,1,81}],Table[0==UnivPoly/.SetOfPointsSigmaC[[i]],{i,1,81}]];
Eqs=Join[Eqs,Table[0==SeriesCoefficient[UnivPoly/.TangentsAtC[[i]],{t,0,1}],{i,1,81}]];
sols=Chop[Solve[Eqs]][[1]];
UnivPolySolved =UnivPoly/.sols;

PossibleEqs={};
For[j=1,j≤Length[AllMons],j++,
UPs=Coefficient[UnivPolySolved,CoeffList[[j]]];
If[Abs[UPs/.SetOfPointsW[[1]]]>.1,
Print[j];
AppendTo[PossibleEqs,Sum[RootApproximant[Coefficient[UPs,AllMons[[i]]],2]AllMons[[i]],{i,1,Length[AllMons]}]];
];
];
U4= Factor[z PossibleEqs[[1]]/Sectiond];
U5= Factor[U4/.SubSigma];
U6=Factor[U5/.SubSigma];

AllMons=MonomialList[(W1+W2+W3+W0)^1(W4+W5+W6)^2 + (W0+W1+W2+W3)^3]; (* experimentally found to give one solution up to scaling *)
count=Length[AllMons];
AllMons=Table[AllMons[[i]]/(AllMons[[i]]/.Table[VarsW[[k]]->1,{k,1,8}]),{i,1,count}];

UnivPoly=Sum[CoeffList[[i]] AllMons[[i]],{i,1,Length[AllMons]}];
Eqs = Join[Table[0==UnivPoly/.SetOfPointsSigma2C[[i]],{i,1,81}],Table[0==UnivPoly/.SetOfPointsSigmaC[[i]],{i,1,81}]];
Eqs=Join[Eqs,Table[0==SeriesCoefficient[UnivPoly/.TangentsAtSigma2C[[i]],{t,0,1}],{i,1,81}]];
sols=Chop[Solve[Eqs]][[1]];
UnivPolySolved =UnivPoly/.sols;

PossibleEqs={};
For[j=1,j≤Length[AllMons],j++,
UPs=Coefficient[UnivPolySolved,CoeffList[[j]]];
If[Abs[UPs/.SetOfPointsW[[1]]]>.1,
Print[j];
AppendTo[PossibleEqs,Sum[RootApproximant[Coefficient[UPs,AllMons[[i]]],2]AllMons[[i]],{i,1,Length[AllMons]}]];
];
];
U7= Factor[z^(-1) PossibleEqs[[1]]/ Sectiond] ;
U8=Factor[ U7/.SubSigma];
U9=Factor[U8/.SubSigma];

We will make a set of points for the triple cover Z (or rather on its double cover).

SetOfPointsZ=Table[Append[SetOfPointsW[[i]],z->N[(zcube/.SetOfPointsW[[i]])^(1/3),6000]],{i,1,Length[SetOfPointsW]}];
SetOfPointsP = Table[{P0->W0/.SetOfPointsZ[[i]],P1->W1/.SetOfPointsZ[[i]],P2->W2/.SetOfPointsZ[[i]],P3->W3/.SetOfPointsZ[[i]],
P4->U4/.SetOfPointsZ[[i]],P5->U5/.SetOfPointsZ[[i]],P6->U6/.SetOfPointsZ[[i]],P7->U7/.SetOfPointsZ[[i]],P8->U8/.SetOfPointsZ[[i]],P9->U9/.SetOfPointsZ[[i]]},{i,1,Length[SetOfPointsZ]}];

Now we find the cubic equations. We will call the variables P0,...,P9.

maxlength=500;
CoeffList=Table[ToExpression[StringJoin["A",ToString[i]]],{i,1,maxlength}];

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

deg=3;
VarsP={P0,P1,P2,P3,P4,P5,P6,P7,P8,P9};
AllMons=MonomialList[(P0+P1+P2+P3+P4+P5+P6+P7+P8+P9)^deg];
count=Length[AllMons]
AllMons=Table[AllMons[[i]]/(AllMons[[i]]/.Table[VarsP[[m]]->1,{m,1,10}]),{i,1,count}];
WtinZ={0,0,0,0,1,1,1,-1,-1,-1};

EqsFPP={};
For[k=0,k≤2,k++,
UnivPolyP=0;
Print[k];
count1=0;
For[j=1,j≤count, j++,
expall=Sum[Exponent[AllMons[[j]],VarsP[[m]]] WtinZ[[m]],{m,1,10}];
If[Mod[expall,3]==k,
UnivPolyP+=CoeffList[[j]] AllMons[[j]];
count1++;
];
];

Eqs=Table[0==UnivPolyP/.SetOfPointsP[[l]],{l,1,count1+10}];
sols=NSolve[Eqs][[1]];

UpSolved=UnivPolyP/.sols;
For[j=1,j≤count, j++,
neweq=Coefficient[UpSolved,CoeffList[[j]]];
If[ByteCount[neweq]>100,
AppendTo[EqsFPP,Sum[Beautify[Coefficient[neweq,AllMons[[l]]]]AllMons[[l]],{l,1,count}]];
];
];

];

Print[Length[EqsFPP]];

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

Find a nicer basis of the space of sections. We will look at nicer values at the fixed points of the group action.
I want to find three numbers a1, a2, a3 so that the values of a1 P4 + a2 P5 + a6 P6 on the three kinds of fixed points 
is small. The numbers a1, a2, a3 must be in Q[sqrt(-15)], and we will want the overall size to be small.
So we will have 6 coefficients.

solsfixedlift=NSolve[Join[Table[0==EqsFPP[[i]],{i,1,84}],{P0==1,P1==P2,P2==P3,P4==P5,P5==P6,P7==P8,P8==P9}],400];
fixed1=Table[Table[VarsP[[i]]->Beautify[VarsP[[i]]/.solsfixedlift[[j]]],{i,1,10}],{j,1,3}];

w=Exp[2 Pi I /3]

solsfixedliftw=NSolve[Join[Table[0==EqsFPP[[i]],{i,1,84}],{P0==1,P1== P2,P2== P3,P4==w P5,P5==w P6,P7==w^2 P8,P8==w^2 P9}],400];
fixed2=Table[Table[VarsP[[i]]->Beautify[VarsP[[i]]/.solsfixedliftw[[j]]],{i,1,10}],{j,1,3}];

solsfixedliftw2=NSolve[Join[Table[0==EqsFPP[[i]],{i,1,84}],{P0==1,P1== P2,P2== P3,P4==w ^2P5,P5==w^2 P6,P7==w P8,P8==w P9}],400];
fixed3 = Table[Table[VarsP[[i]]->Beautify[VarsP[[i]]/.solsfixedliftw2[[j]]],{i,1,10}],{j,1,3}];

ToBasis[x_]:=
Ret={AlgebraicNumberTrace[x,Extension->w +I (Sqrt[15])],
        AlgebraicNumberTrace[x w ,Extension->w +I (Sqrt[15])],
     AlgebraicNumberTrace[x I (Sqrt[15]),Extension->w +I (Sqrt[15])],
     AlgebraicNumberTrace[x w I (Sqrt[15]),Extension->w +I (Sqrt[15])]
};

ValuesAtPts={};
as={P4,P5,P6, I (Sqrt[15]) P4, I (Sqrt[15]) P5, I (Sqrt[15]) P6};  
For[k=1,k≤6,k++,
nextrow={};
For[l=1,l≤3,l++,
nextrow=Join[nextrow,ToBasis[as[[k]]/.fixed1[[l]]]];
];
For[l=1,l≤3,l++,
nextrow=Join[nextrow,ToBasis[as[[k]]/.fixed2[[l]]]];
];
For[l=1,l≤3,l++,
nextrow=Join[nextrow,ToBasis[as[[k]]/.fixed3[[l]]]];
];
AppendTo[ValuesAtPts,nextrow];
];

{u,r, v} = SmithDecomposition[ValuesAtPts]; 
vt=Transpose[v];
SmallValues={};
For[i=1,i≤6,i++,
basisvec=Table[0,{k,1,36}];
basisvec[[i]]=1;
AppendTo[SmallValues,LinearSolve[vt,basisvec]];
];
ValuesSmall=LatticeReduce[SmallValues]

goodb456=Solve[{b1,b2,b3,b4,b5,b6}.ValuesAtPts==6 ValuesSmall[[1]]+2 ValuesSmall[[6]]][[1]];
goodP =( P4(b1 + I (Sqrt[15]) b4) +P5(b2 + I (Sqrt[15]) b5) +P6(b3 + I (Sqrt[15]) b6)) /.goodb456;

Beautify[goodP/.fixed1]
Beautify[N[goodP/.fixed2,4000]]
Beautify[N[goodP/.fixed3,4000]]

(* now for the other three coords *)

ValuesAtPts={};
as={P7,P8,P9, I (Sqrt[15]) P7, I (Sqrt[15]) P8, I (Sqrt[15]) P9};  
For[k=1,k≤6,k++,
nextrow={};
For[l=1,l≤3,l++,
nextrow=Join[nextrow,ToBasis[as[[k]]/.fixed1[[l]]]];
];
For[l=1,l≤3,l++,
nextrow=Join[nextrow,ToBasis[as[[k]]/.fixed2[[l]]]];
];
For[l=1,l≤3,l++,
nextrow=Join[nextrow,ToBasis[as[[k]]/.fixed3[[l]]]];
];
AppendTo[ValuesAtPts,nextrow];
];
ValuesAtPts;

{u,r, v} = SmithDecomposition[ValuesAtPts]; 
vt=Transpose[v];
SmallValues={};
For[i=1,i≤6,i++,
basisvec=Table[0,{k,1,36}];
basisvec[[i]]=1;
AppendTo[SmallValues,LinearSolve[vt,basisvec]];
];
ValuesSmall=LatticeReduce[SmallValues]

goodb789=Solve[{b1,b2,b3,b4,b5,b6}.ValuesAtPts==6ValuesSmall[[2]]+2ValuesSmall[[6]]][[1]]
goodP =( P7(b1 + I (Sqrt[15]) b4) +P8(b2 + I (Sqrt[15]) b5) +P9(b3 + I (Sqrt[15]) b6)) /.goodb789;
Det[{{(b1 + I (Sqrt[15]) b4),(b2 + I (Sqrt[15]) b5) ,(b3 + I (Sqrt[15]) b6)},{(b2 + I (Sqrt[15]) b5) ,(b3 + I (Sqrt[15]) b6),(b1 + I (Sqrt[15]) b4)},{(b3 + I (Sqrt[15]) b6),(b1 + I (Sqrt[15]) b4),(b2 + I (Sqrt[15]) b5) }}/.goodb789]
Beautify[N[goodP/.fixed1,4000]]
Beautify[N[goodP/.fixed2,4000]]
Beautify[N[goodP/.fixed3,4000]]

ValuesAtPts={};
as={P0,P1,P0 I (Sqrt[15]),P1 I (Sqrt[15])};  
For[k=1,k≤4,k++,
nextrow={};
nextrow=Join[nextrow,ToBasis[as[[k]]/.fixed1[[1]]]];
nextrow=Join[nextrow,ToBasis[as[[k]]/.fixed2[[1]]]];
nextrow=Join[nextrow,ToBasis[as[[k]]/.fixed3[[1]]]];
AppendTo[ValuesAtPts,nextrow];
];
ValuesAtPts

{u,r, v} = SmithDecomposition[ValuesAtPts]; 
vt=Transpose[v];
SmallValues={};
For[i=1,i≤4,i++,
basisvec=Table[0,{k,1,12}];
basisvec[[i]]=1;
AppendTo[SmallValues,LinearSolve[vt,basisvec]];
];
ValuesSmall=LatticeReduce[SmallValues]

goodb123=Solve[{b1,b2,b3,b4}.ValuesAtPts==ValuesSmall[[3]]][[1]]
goodP =( P0(b1 + I (Sqrt[15]) b3) +P1(b2 + I (Sqrt[15]) b4) ) /.goodb123

Beautify[N[goodP/.fixed1,4000]]
Beautify[N[goodP/.fixed2,4000]]
Beautify[N[goodP/.fixed3,4000]]

sols0={a1->b1+b3 I (Sqrt[15]),a2->b2+b4 I (Sqrt[15])}/.goodb123;
sols1={a1->b1+b4 I (Sqrt[15]),a2->b2+b5 I (Sqrt[15]),a3->b3+b6 I (Sqrt[15])}/.goodb456;
sols2={a1->b1+b4 I (Sqrt[15]),a2->b2+b5 I (Sqrt[15]),a3->b3+b6 I (Sqrt[15])}/.goodb789;

NewPtoW={P0->W0,
P1->(a2 W1+ a1 W0)/.sols0,
P2->(a2 W2+ a1 W0)/.sols0,
P3->(a2 W3+ a1 W0)/.sols0,
P4->((a1 U4 + a2 U5 + a3 U6)/.sols1),
P5->((a3 U4 + a1 U5 + a2 U6)/.sols1),
P6->((a2 U4 + a3 U5 + a1 U6)/.sols1),
P7->((a1 U7 + a2 U8 + a3 U9)/.sols2),
P8->((a3 U7 + a1 U8 + a2 U9)/.sols2),
P9->((a2 U7 + a3 U8 + a1 U9)/.sols2)};
SetOfPointsP =Table[Table[VarsP[[k]]->((VarsP[[k]]/.NewPtoW)/.SetOfPointsZ[[numpt]]),{k,1,10}],{numpt,1,Length[SetOfPointsZ]}];
deg=3;
VarsP={P0,P1,P2,P3,P4,P5,P6,P7,P8,P9};
AllMons=MonomialList[(P0+P1+P2+P3+P4+P5+P6+P7+P8+P9)^deg];
count=Length[AllMons]
AllMons=Table[AllMons[[i]]/(AllMons[[i]]/.Table[VarsP[[m]]->1,{m,1,10}]),{i,1,count}];
WtinZ={0,0,0,0,1,1,1,-1,-1,-1};

EqsFPPnew={};
For[k=0,k≤2,k++,
UnivPolyP=0;
Print[k];
count1=0;
For[j=1,j≤count, j++,
expall=Sum[Exponent[AllMons[[j]],VarsP[[m]]] WtinZ[[m]],{m,1,10}];
If[Mod[expall,3]==k,
UnivPolyP+=CoeffList[[j]] AllMons[[j]];
count1++;
];
];

Eqs=Table[0==UnivPolyP/.SetOfPointsP[[l]],{l,1,count1+10}];
sols=NSolve[Eqs][[1]];
Print[count1-Length[sols]];

UpSolved=UnivPolyP/.sols;
For[j=1,j≤count, j++,
neweq=Coefficient[UpSolved,CoeffList[[j]]];
If[ByteCount[neweq]>100,
AppendTo[EqsFPPnew,Sum[Beautify[Coefficient[neweq,AllMons[[l]]]]AllMons[[l]],{l,1,count}]];
];
];

];

Print[Length[EqsFPPnew]];

ptP={P0->0,P1->1,P2->0,P3->0,P4->0,P5->0,P6->0,P7->0,P8->0,P9->0};
Simplify[EqsFPPnew/.ptP]
JacMat=Table[D[EqsFPPnew[[k]],VarsP[[i]]]/.ptP,{k,1,84},{i,1,10}];
RootReduce[NullSpace[JacMat]]
Table[ByteCount[EqsFPPnew[[i]]],{i,1,84}]

Table[ByteCount[Numerator[Factor[EqsFPPnew[[i]]]]],{i,1,84}]
Table[ByteCount[Numerator[Factor[EqsFPP[[i]]]]],{i,1,84}]

solsfixedlift=NSolve[Join[Table[0==EqsFPPnew[[i]],{i,1,84}],{P0==1,P1==P2,P2==P3,P4==P5,P5==P6,P7==P8,P8==P9}],400];
fixed1=Table[Table[VarsP[[i]]->Beautify[VarsP[[i]]/.solsfixedlift[[j]]],{i,1,10}],{j,1,3}]
w=Exp[2 Pi I /3]

solsfixedliftw=NSolve[Join[Table[0==EqsFPPnew[[i]],{i,1,84}],{P0==1,P1== P2,P2== P3,P4==w P5,P5==w P6,P7==w^2 P8,P8==w^2 P9}],400];
fixed2=Table[Table[VarsP[[i]]->Beautify[VarsP[[i]]/.solsfixedliftw[[j]]],{i,1,10}],{j,1,3}]

solsfixedliftw2=NSolve[Join[Table[0==EqsFPPnew[[i]],{i,1,84}],{P0==1,P1== P2,P2== P3,P4==w ^2P5,P5==w^2 P6,P7==w P8,P8==w P9}],400];
fixed3 = Table[Table[VarsP[[i]]->Beautify[VarsP[[i]]/.solsfixedliftw2[[j]]],{i,1,10}],{j,1,3}]

Save, normally and in Magma format. Will also save it Mod 79.

For[i=1,i≤84,i++,
EqsFPPnew[[i]]=Numerator[Factor[EqsFPPnew[[i]]]];
];
EqsFPPnewWithrr={};
For[i=1,i≤84,i++,
newrel=0;
For[j=1,j≤Length[AllMons],j++,
cf=Coefficient[EqsFPPnew[[i]] ,AllMons[[j]]];
newrel+=(Re[cf]+Im[cf] (rr/Sqrt[15])) AllMons[[j]] ;
];
AppendTo[EqsFPPnewWithrr,newrel];
];


EqsFPPnewMod79=Table[Sum[Mod[Coefficient[EqsFPPnewWithrr[[i]]/.rr->8,AllMons[[j]]],79]
AllMons[[j]],{j,1,Length[AllMons]}],{i,1,84}];

fname="EqsFPPd3D3";
OpenWrite[fname];
Write[fname,EqsFPPnew];
Close[fname];

fname="EqsFPPd3D3Magma";
OpenWrite[fname];
Write[fname,EqsFPPnewWithrr];
Close[fname];

fname="EqsFPPd3D3mod79";
OpenWrite[fname];
Write[fname,EqsFPPnewMod79];
Close[fname];

We will now look for equations for the double cover of FPP.

(* will find odd sections of 2K on the double cover of FPP *)
AllMons=MonomialList[(W1+W2+W3+W0)^2(W4+W5+W6)^1 + (W1^2+W2^2+W3^2) W7+ (W4+W5+W6+W7)^3]; (* experimentally found to give one solution up to scaling *)
count=Length[AllMons];
AllMons=Table[AllMons[[i]]/(AllMons[[i]]/.Table[VarsW[[k]]->1,{k,1,8}]),{i,1,count}];

UnivPoly=Sum[CoeffList[[i]] AllMons[[i]],{i,1,Length[AllMons]}];
Eqs = Join[Table[0==UnivPoly/.SetOfPointsC[[i]],{i,1,81}],Table[0==UnivPoly/.SetOfPointsSigmaC[[i]],{i,1,81}]];
Eqs=Join[Eqs,Table[0==SeriesCoefficient[UnivPoly/.TangentsAtC[[i]],{t,0,1}],{i,1,81}]];
sols=Chop[Solve[Eqs]][[1]];
UnivPolySolved =UnivPoly/.sols;

PossibleEqs={};
For[j=1,j≤Length[AllMons],j++,
UPs=Coefficient[UnivPolySolved,CoeffList[[j]]];
If[Abs[UPs/.SetOfPointsW[[1]]]>.1,
Print[j];
AppendTo[PossibleEqs,Sum[RootApproximant[Coefficient[UPs,AllMons[[i]]],2]AllMons[[i]],{i,1,Length[AllMons]}]];
];
];
U14= Factor[z PossibleEqs[[1]]/Sectiond];
U15= Factor[U14/.SubSigma];
U16=Factor[U15/.SubSigma];

AllMons=MonomialList[(W1+W2+W3+W0)^2(W4+W5+W6)^1 + (W1^2+W2^2+W3^2) W7+ (W4+W5+W6+W7)^3]; (* experimentally found to give one solution up to scaling *)
count=Length[AllMons];
AllMons=Table[AllMons[[i]]/(AllMons[[i]]/.Table[VarsW[[k]]->1,{k,1,8}]),{i,1,count}];

UnivPoly=Sum[CoeffList[[i]] AllMons[[i]],{i,1,Length[AllMons]}];
Eqs = Join[Table[0==UnivPoly/.SetOfPointsSigma2C[[i]],{i,1,81}],Table[0==UnivPoly/.SetOfPointsSigmaC[[i]],{i,1,81}]];
Eqs=Join[Eqs,Table[0==SeriesCoefficient[UnivPoly/.TangentsAtSigma2C[[i]],{t,0,1}],{i,1,81}]];
sols=Chop[Solve[Eqs]][[1]];
UnivPolySolved =UnivPoly/.sols;

PossibleEqs={};
For[j=1,j≤Length[AllMons],j++,
UPs=Coefficient[UnivPolySolved,CoeffList[[j]]];
If[Abs[UPs/.SetOfPointsW[[1]]]>.1,
Print[j];
AppendTo[PossibleEqs,Sum[RootApproximant[Coefficient[UPs,AllMons[[i]]],2]AllMons[[i]],{i,1,Length[AllMons]}]];
];
];
U17= Factor[z^(-1) PossibleEqs[[1]]/ Sectiond] ;
U18=Factor[ U17/.SubSigma];
U19=Factor[U18/.SubSigma];

VarsP20={P0,P1,P2,P3,P4,P5,P6,P7,P8,P9,P10,P11,P12,P13,P14,P15,P16,P17,P18,P19};
P20toW={P0->W0,
P1->(a2 W1+ a1 W0)/.sols0,
P2->(a2 W2+ a1 W0)/.sols0,
P3->(a2 W3+ a1 W0)/.sols0,
P4->((a1 U4 + a2 U5 + a3 U6)/.sols1),
P5->((a3 U4 + a1 U5 + a2 U6)/.sols1),
P6->((a2 U4 + a3 U5 + a1 U6)/.sols1),
P7->((a1 U7 + a2 U8 + a3 U9)/.sols2),
P8->((a3 U7 + a1 U8 + a2 U9)/.sols2),
P9->((a2 U7 + a3 U8 + a1 U9)/.sols2),
P10->W7,
P11->W4,
P12->W5,
P13->W6,
P14->U14,
P15->U15,
P16->U16,
P17->U17,
P18->U18,
P19->U19
};
SetOfPointsP20 =Table[Table[VarsP20[[k]]->((VarsP20[[k]]/.P20toW)/.SetOfPointsZ[[numpt]]),{k,1,20}],{numpt,1,Length[SetOfPointsZ]}];
deg=2;

AllMons=MonomialList[Sum[VarsP20[[k]],{k,1,20}]^deg];
count=Length[AllMons]
AllMons=Table[AllMons[[i]]/(AllMons[[i]]/.Table[VarsP20[[m]]->1,{m,1,20}]),{i,1,count}];
WtinZ={{0,0},{0,0},{0,0},{0,0},{0,1},{0,1},{0,1},{0,-1},{0,-1},{0,-1},
{1,0},{1,0},{1,0},{1,0},{1,1},{1,1},{1,1},{1,-1},{1,-1},{1,-1}
};

EqsFPPdouble={};
For[parity=0,parity≤1,parity++,
For[k=0,k≤2,k++,
UnivPolyP=0;
Print[{parity,k}];
count1=0;
For[j=1,j≤count, j++,
expall=Sum[Exponent[AllMons[[j]],VarsP20[[m]]] WtinZ[[m]],{m,1,20}];
If[Mod[expall[[2]],3]==k,
If[Mod[expall[[1]],2]==parity,
UnivPolyP+=CoeffList[[j]] AllMons[[j]];
count1++;
];];
];


Eqs=Table[0==UnivPolyP/.SetOfPointsP20[[l]],{l,1,count1+10}];
sols=NSolve[Eqs][[1]];
Print[count1-Length[sols]];

UpSolved=Chop[UnivPolyP/.sols];
For[j=1,j≤count, j++,
neweq=Coefficient[UpSolved,CoeffList[[j]]];
If[ByteCount[neweq]>100,
AppendTo[EqsFPPdouble,Sum[Beautify[Coefficient[neweq,AllMons[[l]]]]AllMons[[l]],{l,1,count}]];
];
];

];];

Print[Length[EqsFPPdouble]];

Table[ByteCount[EqsFPPdouble[[i]]],{i,1,100}]

(* now we will make the other variables nicer, as before *)
solsfixedlift=NSolve[Join[Table[0==EqsFPPdouble[[i]],{i,1,100}],{P0==1,P1==P2,P2==P3,P4==P5,P5==P6,P7==P8,P8==P9,P11==P12,P12==P13,P14==P15,P15==P16,P17==P18,P18==P19}],400];
fixed1=Table[Table[VarsP20[[i]]->Beautify[VarsP20[[i]]/.solsfixedlift[[j]]],{i,1,20}],{j,1,3}];
Length[fixed1]

solsfixedlift=NSolve[Join[Table[0==EqsFPPdouble[[i]],{i,1,100}],{P0==1,P1==P2,P2==P3,P4== w P5,P5== w P6,P7==w^2  P8,P8==w^2 P9,P11== P12,P12==  P13,P14==w P15,P15==w P16,P17==w^2 P18,P18==w^2 P19}],400];
fixed2=Table[Table[VarsP20[[i]]->Beautify[VarsP20[[i]]/.solsfixedlift[[j]]],{i,1,20}],{j,1,3}];
Length[fixed2]

solsfixedlift=NSolve[Join[Table[0==EqsFPPdouble[[i]],{i,1,100}],{P0==1,P1==P2,P2==P3,P4== w^2 P5,P5== w^2 P6,P7==w  P8,P8==w P9,P11== P12,P12==  P13,P14==w^2 P15,P15==w^2 P16,P17==w P18,P18==w P19}],400];
fixed3=Table[Table[VarsP20[[i]]->Beautify[VarsP20[[i]]/.solsfixedlift[[j]]],{i,1,20}],{j,1,3}];
Length[fixed3]

(* making P14, P15, P16 nicer *)
ValuesAtPts={};
as={P14,P15,P16, I (Sqrt[15]) P14, I (Sqrt[15]) P15, I (Sqrt[15]) P16};  
For[k=1,k≤6,k++,
nextrow={};
For[l=1,l≤3,l++,
nextrow=Join[nextrow,ToBasis[as[[k]]/.fixed1[[l]]]];
];
For[l=1,l≤3,l++,
nextrow=Join[nextrow,ToBasis[as[[k]]/.fixed2[[l]]]];
];
For[l=1,l≤3,l++,
nextrow=Join[nextrow,ToBasis[as[[k]]/.fixed3[[l]]]];
];
AppendTo[ValuesAtPts,nextrow];
];

{u,r, v} = SmithDecomposition[ValuesAtPts]; 
vt=Transpose[v];
SmallValues={};
For[i=1,i≤6,i++,
basisvec=Table[0,{k,1,36}];
basisvec[[i]]=1;
AppendTo[SmallValues,LinearSolve[vt,basisvec]];
];
ValuesSmall=LatticeReduce[SmallValues]

goodb141516=Solve[{b1,b2,b3,b4,b5,b6}.ValuesAtPts==6 ValuesSmall[[1]]+2 ValuesSmall[[6]]][[1]];
goodP =( P14(b1 + I (Sqrt[15]) b4) +P15(b2 + I (Sqrt[15]) b5) +P16(b3 + I (Sqrt[15]) b6)) /.goodb141516;

Beautify[goodP/.fixed1]
Beautify[N[goodP/.fixed2,4000]]
Beautify[N[goodP/.fixed3,4000]]

(* making P17, P18, P19 nicer *)
ValuesAtPts={};
as={P17,P18,P19, I (Sqrt[15]) P17, I (Sqrt[15]) P18, I (Sqrt[15]) P19};  
For[k=1,k≤6,k++,
nextrow={};
For[l=1,l≤3,l++,
nextrow=Join[nextrow,ToBasis[as[[k]]/.fixed1[[l]]]];
];
For[l=1,l≤3,l++,
nextrow=Join[nextrow,ToBasis[as[[k]]/.fixed2[[l]]]];
];
For[l=1,l≤3,l++,
nextrow=Join[nextrow,ToBasis[as[[k]]/.fixed3[[l]]]];
];
AppendTo[ValuesAtPts,nextrow];
];

{u,r, v} = SmithDecomposition[ValuesAtPts]; 
vt=Transpose[v];
SmallValues={};
For[i=1,i≤6,i++,
basisvec=Table[0,{k,1,36}];
basisvec[[i]]=1;
AppendTo[SmallValues,LinearSolve[vt,basisvec]];
];
ValuesSmall=LatticeReduce[SmallValues]

goodb171819=Solve[{b1,b2,b3,b4,b5,b6}.ValuesAtPts==6 ValuesSmall[[1]]+2 ValuesSmall[[6]]][[1]];
goodP =( P17(b1 + I (Sqrt[15]) b4) +P18(b2 + I (Sqrt[15]) b5) +P19(b3 + I (Sqrt[15]) b6)) /.goodb171819;

(* it is not perfect, but it will be OK *)
Beautify[goodP/.fixed1]
Beautify[N[goodP/.fixed2,4000]]
Beautify[N[goodP/.fixed3,4000]]

ValuesAtPts={};
as={P10,P11,P10 I (Sqrt[15]),P11 I (Sqrt[15])};  
For[k=1,k≤4,k++,
nextrow={};
nextrow=Join[nextrow,ToBasis[as[[k]]/.fixed1[[1]]]];
nextrow=Join[nextrow,ToBasis[as[[k]]/.fixed2[[1]]]];
nextrow=Join[nextrow,ToBasis[as[[k]]/.fixed3[[1]]]];
AppendTo[ValuesAtPts,nextrow];
];
ValuesAtPts

{u,r, v} = SmithDecomposition[ValuesAtPts]; 
vt=Transpose[v];
SmallValues={};
For[i=1,i≤4,i++,
basisvec=Table[0,{k,1,12}];
basisvec[[i]]=1;
AppendTo[SmallValues,LinearSolve[vt,basisvec]];
];
ValuesSmall=LatticeReduce[SmallValues]

goodb111213=Solve[{b1,b2,b3,b4}.ValuesAtPts==ValuesSmall[[3]]][[1]]
goodP =( P10(b1 + I (Sqrt[15]) b3) +P11(b2 + I (Sqrt[15]) b4) ) /.goodb111213

Beautify[N[goodP/.fixed1,4000]]
Beautify[N[goodP/.fixed2,4000]]
Beautify[N[goodP/.fixed3,4000]]

sols10={a1->b1+b3 I (Sqrt[15]),a2->b2+b4 I (Sqrt[15])}/.goodb111213;
sols11={a1->b1+b4 I (Sqrt[15]),a2->b2+b5 I (Sqrt[15]),a3->b3+b6 I (Sqrt[15])}/.goodb141516;
sols12={a1->b1+b4 I (Sqrt[15]),a2->b2+b5 I (Sqrt[15]),a3->b3+b6 I (Sqrt[15])}/.goodb171819;

P20toW={P0->W0,
P1->(a2 W1+ a1 W0)/.sols0,
P2->(a2 W2+ a1 W0)/.sols0,
P3->(a2 W3+ a1 W0)/.sols0,
P4->((a1 U4 + a2 U5 + a3 U6)/.sols1),
P5->((a3 U4 + a1 U5 + a2 U6)/.sols1),
P6->((a2 U4 + a3 U5 + a1 U6)/.sols1),
P7->((a1 U7 + a2 U8 + a3 U9)/.sols2),
P8->((a3 U7 + a1 U8 + a2 U9)/.sols2),
P9->((a2 U7 + a3 U8 + a1 U9)/.sols2),
P10->W7,
P11->(a2 W4+ a1 W7)/.sols10, 
P12->(a2 W5+ a1 W7)/.sols10,
P13->(a2 W6+ a1 W7)/.sols10,
P14->((a1 U14 + a2 U15 + a3 U16)/.sols11),
P15->((a3 U14 + a1 U15 + a2 U16)/.sols11),
P16->((a2 U14 + a3 U15 + a1 U16)/.sols11),
P17->((a1 U17 + a2 U18 + a3 U19)/.sols12),
P18->((a3 U17 + a1 U18 + a2 U19)/.sols12),
P19->((a2 U17 + a3 U18 + a1 U19)/.sols12)
};

SetOfPointsP20 =Table[Table[VarsP20[[k]]->((VarsP20[[k]]/.P20toW)/.SetOfPointsZ[[numpt]]),{k,1,20}],{numpt,1,Length[SetOfPointsZ]}];
deg=2;

AllMons=MonomialList[Sum[VarsP20[[k]],{k,1,20}]^deg];
count=Length[AllMons]
AllMons=Table[AllMons[[i]]/(AllMons[[i]]/.Table[VarsP20[[m]]->1,{m,1,20}]),{i,1,count}];
WtinZ={{0,0},{0,0},{0,0},{0,0},{0,1},{0,1},{0,1},{0,-1},{0,-1},{0,-1},
{1,0},{1,0},{1,0},{1,0},{1,1},{1,1},{1,1},{1,-1},{1,-1},{1,-1}
};

EqsFPPnewdouble={};
For[parity=0,parity≤1,parity++,
For[k=0,k≤2,k++,
UnivPolyP=0;
Print[{parity,k}];
count1=0;
For[j=1,j≤count, j++,
expall=Sum[Exponent[AllMons[[j]],VarsP20[[m]]] WtinZ[[m]],{m,1,20}];
If[Mod[expall[[2]],3]==k,
If[Mod[expall[[1]],2]==parity,
UnivPolyP+=CoeffList[[j]] AllMons[[j]];
count1++;
];];
];


Eqs=Table[0==UnivPolyP/.SetOfPointsP20[[l]],{l,1,count1+10}];
sols=NSolve[Eqs][[1]];
Print[count1-Length[sols]];

UpSolved=Chop[UnivPolyP/.sols];
For[j=1,j≤count, j++,
neweq=Coefficient[UpSolved,CoeffList[[j]]];
If[ByteCount[neweq]>100,
AppendTo[EqsFPPnewdouble,Sum[Beautify[Coefficient[neweq,AllMons[[l]]]]AllMons[[l]],{l,1,count}]];
];
];

];];

Print[Length[EqsFPPnewdouble]];

We will write these equations into the files, as usual.

Table[ByteCount[Numerator[Factor[EqsFPPnewdouble[[i]]]]],{i,1,100}]
For[i=1,i≤100,i++,
EqsFPPnewdouble[[i]]=Numerator[Factor[EqsFPPnewdouble[[i]]]];
];
EqsFPPnewdoubleWithrr={};
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[EqsFPPnewdoubleWithrr,newrel];
];
EqsFPPnewdoubleMod79=Table[Sum[Mod[Coefficient[EqsFPPnewdoubleWithrr[[i]]/.rr->8,AllMons[[j]]],79]
AllMons[[j]],{j,1,Length[AllMons]}],{i,1,100}];

fname="EqsFPPd3D3double";
OpenWrite[fname];
Write[fname,EqsFPPnewdouble];
Close[fname];

fname="EqsFPPd3D3doubleMagma";
OpenWrite[fname];
Write[fname,EqsFPPnewdoubleWithrr];
Close[fname];

fname="EqsFPPd3D3doublemod79";
OpenWrite[fname];
Write[fname,EqsFPPnewdoubleMod79];
Close[fname];

We will compute the quadratic and cubic equations on invariant variables in this basis.

P20toW={P0->W0,
P1->(a2 W1+ a1 W0)/.sols0,
P2->(a2 W2+ a1 W0)/.sols0,
P3->(a2 W3+ a1 W0)/.sols0,
P10->W7,
P11->(a2 W4+ a1 W7)/.sols0,
P12->(a2 W5+ a1 W7)/.sols0,
P13->(a2 W6+ a1 W7)/.sols0,
};

EqsFPPmodC3double={};
WtinZ={{0,0},{0,0},{0,0},{0,0},{0,1},{0,1},{0,1},{0,-1},{0,-1},{0,-1},
{1,0},{1,0},{1,0},{1,0},{1,1},{1,1},{1,1},{1,-1},{1,-1},{1,-1}
};
For[deg=2,deg≤3,deg++,
For[parity=0,parity≤1,parity++,
AllMons=MonomialList[(P0+P1+P2+P3+P10+P11+P12+P13)^deg];
count=Length[AllMons];
AllMons=Table[AllMons[[i]]/(AllMons[[i]]/.Table[VarsP20[[m]]->1,{m,1,20}]),{i,1,count}];

UnivPolyP=0;
Print["parity=",parity];
count1=0;
For[j=1,j≤count, j++,
expall=Sum[Exponent[AllMons[[j]],VarsP20[[m]]] WtinZ[[m]],{m,1,20}];
If[Mod[expall[[1]],2]==parity,
UnivPolyP+=CoeffList[[j]] AllMons[[j]];
count1++;
];
];


Eqs=Table[0==UnivPolyP/.SetOfPointsP20[[l]],{l,1,count1+10}];
sols=NSolve[Eqs][[1]];
Print[count1-Length[sols]];

UpSolved=Chop[UnivPolyP/.sols];
For[j=1,j≤count, j++,
neweq=Coefficient[UpSolved,CoeffList[[j]]];
If[ByteCount[neweq]>100,
AppendTo[EqsFPPmodC3double,Sum[Beautify[Coefficient[neweq,AllMons[[l]]]]AllMons[[l]],{l,1,count}]];
];
];
];];

EqsFPPmodC3double=Table[Numerator[Factor[EqsFPPmodC3double[[i]]]],{i,1,Length[EqsFPPmodC3double]}];
Table[ByteCount[EqsFPPmodC3double[[i]]],{i,1,Length[EqsFPPmodC3double]}]