In this file we will do a finite field search with the goal of reducing the computations to roughly p^2 points for which we have singularities at pt1, pt2, pt3

fname ="Desktop/MumfordFPP2021/KeumFPP/Quadricscontainsline2";
OpenRead[fname];
Quadricscontainsline2=Read[fname];
Close[fname];

fname ="Desktop/MumfordFPP2021/KeumFPP/transf";
OpenRead[fname];
transf=Read[fname];
Close[fname];

fname ="Desktop/MumfordFPP2021/KeumFPP/ThreepointsPoly";
OpenRead[fname];
{pt1,pt2,pt4}=Read[fname];
Close[fname];

We will run this for different primes and save the results in GoodRulesPrime. The error messages come from dividing by 0 mod p, it is easier to live with them than to write various checks.

vs={v1,v2,v3,v4,v5,v6};
(* we can use a shorter version of transf *)
qcl2num=Numerator[Factor[Quadricscontainsline2]];
transfshort={transf[[3]],transf[[5]],transf[[9]],transf[[10]],transf[[11]],transf[[12]],transf[[13]],transf[[14]]};


varqcl2={c110,c4,r1,r2,r4,r5,r6};
uvs={u0,u1,v1,v2,v3,v4,v5,v6};

(* this is from the previous part of the project - a necessary condition for having worse than a node *)
H={{h11,h12,h13,h14},{h12,h22,h23,h24},{h13,h23,h33,h34},{h14,h24,h34,h44}};
W={W1,W2,W3,W4};
tbz=Factor[Coefficient[Expand[Det[H+t Table[W[[i]]W[[j]],{i,1,4},{j,1,4}]]],t]]

p=53;
fname =StringJoin[StringJoin["Desktop/MumfordFPP2021/KeumFPP/out",ToString[p]],".txt"]
OpenRead[fname];
outp=Read[fname];
Close[fname];
Length[outp]

fivevars={s1,s2,s3,s4,b2};
ds={ds1,ds2,ds3,ds4,db2};


goodindices={};
For[m=1,m≤Length[outp],m++,

If[Mod[m,500]==20,Print["m=",m]];
news=outp[[m]];

rules0=Table[fivevars[[k]]->news[[k]],{k,1,5}];
rules=rules0;
For[i=1,i≤8,i++,
rules=Join[rules,
Table[transfshort[[9-i]][[k]][[1]]-> PolynomialMod[transfshort[[9-i]][[k]][[2]]/.rules,p]
,{k,1,Length[transfshort[[9-i]]]}]];
];
rulesqcl2=Table[varqcl2[[k]]->PolynomialMod[varqcl2[[k]]/.rules,p],{k,1,7}];
quadmod=PolynomialMod[qcl2num/.rulesqcl2,p];
pt1mod=Table[PolynomialMod[uvs[[k]]/.pt1/.rules0,p],{k,1,8}];
pt2mod=Table[PolynomialMod[uvs[[k]]/.pt2/.rules0,p],{k,1,8}];
pt4mod=Table[PolynomialMod[uvs[[k]]/.pt4/.rules0,p],{k,1,8}];

(* Let's put conditions on being singular at each point *)

sols=Solve[{quadmod[[4]]==0,quadmod[[6]]==0},{v5,v6}][[1]];
Eqs={PolynomialMod[v3 v4 quadmod[[1]]/.sols,p],PolynomialMod[v3 v4 quadmod[[2]]/.sols,p]};

Eqs1=PolynomialMod[Eqs/.{u0->1,v1->pt4mod[[3]]},p]; (* taking a slice *)
Wptspoint=Mod[{D[Eqs1[[1]],u1],D[Eqs1[[1]],v2],D[Eqs1[[1]],v3],D[Eqs1[[1]],v4]}/.Table[uvs[[k]]-> pt4mod[[k]],{k,1,8}],p] ;
(* finding W *)
EqsSing= PolynomialMod[Wptspoint[[1]]Eqs1[[2]]-(D[Eqs1[[2]],u1]/.Table[uvs[[k]]-> pt4mod[[k]],{k,1,8}])Eqs1[[1]],p];
Hptspoint=Mod[Table[D[D[EqsSing,{u1,v2,v3,v4}[[k]]],{u1,v2,v3,v4}[[l]]],{k,1,4},{l,1,4}]/.Table[uvs[[k]]-> pt4mod[[k]],{k,1,8}],p];
SubHW=Join[Flatten[Table[H[[i]][[j]]->Hptspoint[[i]][[j]],{j,1,4},{i,1,j}]],Table[W[[k]]-> Wptspoint[[k]],{k,1,4}]];
tbz1spoint4=Mod[tbz/.SubHW,p];

Eqs1=PolynomialMod[Eqs/.{u0->1,v1-> pt2mod[[3]]},p]; (* taking a slice *)
Wptspoint=Mod[{D[Eqs1[[1]],u1],D[Eqs1[[1]],v2],D[Eqs1[[1]],v3],D[Eqs1[[1]],v4]}/.Table[uvs[[k]]->  pt2mod[[k]],{k,1,8}],p] ;
(* finding W *)
EqsSing= PolynomialMod[Wptspoint[[1]]Eqs1[[2]]-(D[Eqs1[[2]],u1]/.Table[uvs[[k]]->  pt2mod[[k]],{k,1,8}])Eqs1[[1]],p];
Hptspoint=Mod[Table[D[D[EqsSing,{u1,v2,v3,v4}[[k]]],{u1,v2,v3,v4}[[l]]],{k,1,4},{l,1,4}]/.Table[uvs[[k]]->  pt2mod[[k]],{k,1,8}],p];
SubHW=Join[Flatten[Table[H[[i]][[j]]->Hptspoint[[i]][[j]],{j,1,4},{i,1,j}]],Table[W[[k]]-> Wptspoint[[k]],{k,1,4}]];
tbz1spoint2=Mod[tbz/.SubHW,p];

If[tbz1spoint4== 0, If[tbz1spoint2== 0,
(* Print[m," ",news];*)
AppendTo[goodindices,m];
]];

];

(* trying to do the lifting *)
fivevars={s1,s2,s3,s4,b2};
ds={ds1,ds2,ds3,ds4,db2};
maxdeg=200;
Print["Length[goodindices]=",Length[goodindices]];
For[goodindex=1,goodindex≤Length[goodindices],goodindex++,
sold = outp[[goodindices[[goodindex]]]];
deg=1;

For[deg=1,deg≤maxdeg,deg++,
If[Mod[deg,40]==3,Print["goodindex=",goodindex," deg=",deg]];

quadmod={};
pt1mod={};
pt2mod={};
pt4mod={};
detpt1={};
detpt2={};
detpt4={};
A2atpt4={};
A2atpt2={};

For[m=0,m≤5,m++,

news=sold;
If[m>0,news[[m]]+=p^deg]; (* we will compute all partial derivatives manually, the m=0 term is the original value *)
rules0=Table[fivevars[[k]]->news[[k]],{k,1,5}];
rules=rules0;
For[i=1,i≤8,i++,
rules=Join[rules,
Table[transfshort[[9-i]][[k]][[1]]-> PolynomialMod[transfshort[[9-i]][[k]][[2]]/.rules,p^(deg+1)],{k,1,Length[transfshort[[9-i]]]}]];
];
rulesqcl2=Table[varqcl2[[k]]->PolynomialMod[varqcl2[[k]]/.rules,p^(deg+1)],{k,1,7}];
AppendTo[quadmod,PolynomialMod[qcl2num/.rulesqcl2,p^(deg+1)]];

AppendTo[pt1mod,Table[PolynomialMod[uvs[[k]]/.pt1/.rules0,p^(deg+1)],{k,1,8}]];
AppendTo[pt2mod,Table[PolynomialMod[uvs[[k]]/.pt2/.rules0,p^(deg+1)],{k,1,8}]];
AppendTo[pt4mod,Table[PolynomialMod[uvs[[k]]/.pt4/.rules0,p^(deg+1)],{k,1,8}]];
(* Let's put conditions on being singular at each point *)
matrixpt=Table[D[quadmod[[m+1]][[k]],uvs[[l]]]/.Table[uvs[[j]]-> pt1mod[[m+1]][[j]],{j,1,8}],{k,1,4},{l,1,4}];
AppendTo[detpt1,Mod[Det[matrixpt] /p^deg,p]];
matrixpt=Table[D[quadmod[[m+1]][[k]],uvs[[l]]]/.Table[uvs[[j]]-> pt2mod[[m+1]][[j]],{j,1,8}],{k,1,4},{l,1,4}];
AppendTo[detpt2,Mod[Det[matrixpt] /p^deg,p]];
matrixpt=Table[D[quadmod[[m+1]][[k]],uvs[[l]]]/.Table[uvs[[j]]-> pt4mod[[m+1]][[j]],{j,1,8}],{k,1,4},{l,1,4}];
AppendTo[detpt4,Mod[Det[matrixpt] /p^deg,p]];

sols=Solve[{quadmod[[m+1]][[4]]==0,quadmod[[m+1]][[6]]==0},{v5,v6}][[1]];
Eqs={PolynomialMod[v3 v4 quadmod[[m+1]][[1]]/.sols,p^(deg+1)],PolynomialMod[v3 v4 quadmod[[m+1]][[2]]/.sols,p^(deg+1)]};

Eqs1=PolynomialMod[Eqs/.{u0->1,v1->pt4mod[[m+1]][[3]]},p^(deg+1)]; (* taking a slice *)
Wptspoint=Mod[{D[Eqs1[[1]],u1],D[Eqs1[[1]],v2],D[Eqs1[[1]],v3],D[Eqs1[[1]],v4]}/.Table[uvs[[k]]-> pt4mod[[m+1]][[k]],{k,1,8}],p^(deg+1)] ;
(* finding W *)
EqsSing= PolynomialMod[Wptspoint[[1]]Eqs1[[2]]-(D[Eqs1[[2]],u1]/.Table[uvs[[k]]-> pt4mod[[m+1]][[k]],{k,1,8}])Eqs1[[1]],p^(deg+1)];
Hptspoint=Mod[Table[D[D[EqsSing,{u1,v2,v3,v4}[[k]]],{u1,v2,v3,v4}[[l]]],{k,1,4},{l,1,4}]/.Table[uvs[[k]]-> pt4mod[[m+1]][[k]],{k,1,8}],p^(deg+1)];
SubHW=Join[Flatten[Table[H[[i]][[j]]->Hptspoint[[i]][[j]],{j,1,4},{i,1,j}]],Table[W[[k]]-> Wptspoint[[k]],{k,1,4}]];
tbz1spoint=Mod[tbz/.SubHW,p^(deg+1)];
If[Mod[tbz1spoint,p^deg]≠0,Print["pt1warning"]];
AppendTo[A2atpt4,tbz1spoint/p^deg];

Eqs1=PolynomialMod[Eqs/.{u0->1,v1-> pt2mod[[m+1]][[3]]},p^(deg+1)]; (* taking a slice *)
Wptspoint=Mod[{D[Eqs1[[1]],u1],D[Eqs1[[1]],v2],D[Eqs1[[1]],v3],D[Eqs1[[1]],v4]}/.Table[uvs[[k]]->  pt2mod[[m+1]][[k]],{k,1,8}],p^(deg+1)] ;
(* finding W *)
EqsSing= PolynomialMod[Wptspoint[[1]]Eqs1[[2]]-(D[Eqs1[[2]],u1]/.Table[uvs[[k]]->  pt2mod[[m+1]][[k]],{k,1,8}])Eqs1[[1]],p^(deg+1)];
Hptspoint=Mod[Table[D[D[EqsSing,{u1,v2,v3,v4}[[k]]],{u1,v2,v3,v4}[[l]]],{k,1,4},{l,1,4}]/.Table[uvs[[k]]->  pt2mod[[m+1]][[k]],{k,1,8}],p^(deg+1)];
SubHW=Join[Flatten[Table[H[[i]][[j]]->Hptspoint[[i]][[j]],{j,1,4},{i,1,j}]],Table[W[[k]]-> Wptspoint[[k]],{k,1,4}]];
tbz1spoint=Mod[tbz/.SubHW,p^(deg+1)];
If[Mod[tbz1spoint,p^deg]≠0,Print["pt2warning"]];
AppendTo[A2atpt2,tbz1spoint/p^deg];

];


(* now let us solve to get the shift *)
eqs={0== {c,ds1,ds2,ds3,ds4,db2}.detpt1,0== {c,ds1,ds2,ds3,ds4,db2}.detpt2,0== {c,ds1,ds2,ds3,ds4,db2}.detpt4,
0== {c,ds1,ds2,ds3,ds4,db2}.A2atpt4,0== {c,ds1,ds2,ds3,ds4,db2}.A2atpt2,1==c+ds1+ds2+ds3+ds4+db2};

sols=Solve[eqs,Modulus->p][[1]];
If[Length[Variables[Table[sols[[k]][[2]],{k,1,Length[sols]}]]]>0,(* Print[goodindex," nope"];*)deg=maxdeg+100];
sold =sold +   (ds/.sols)p^deg;

];
If[deg==maxdeg+1,
acc=maxdeg+1; 
smallpolys={};
degpoly=12;
For[i=1,i≤5,i++,
whatis=sold[[i]];
thebasis=Append[Table[Mod[whatis^k,p^acc],{k,0,degpoly}],p^acc]; 
smalllincomb=LatticeReduce[NullSpace[{thebasis}]][[1]];
smallpoly=Sum[smalllincomb[[k]] t^(k-1),{k,1,degpoly+1}];
AppendTo[smallpolys,smallpoly];
];
(* Print[sold]; *)
Print[smallpolys];
goodindex=Length[goodindices]+1; (* stop it *)
];
];

It seems to have worked! We will now calculate the coefficients in the number field rather than p-adically.

acc=maxdeg+1; 
smallpolys={};
degpoly=12;
For[i=1,i≤5,i++,
If[i==5,degpoly=6];(* the last one turns out to be degree 6 *)
whatis=sold[[i]];
thebasis=Append[Table[Mod[whatis^k,p^acc],{k,0,degpoly}],p^acc]; 
smalllincomb=LatticeReduce[NullSpace[{thebasis}]][[1]];
smallpoly=Sum[smalllincomb[[k]] t^(k-1),{k,1,degpoly+1}];
AppendTo[smallpolys,smallpoly];
];
Print[smallpolys];

We are trying to figure out precisely which roots of the defining polynomials correspond to the p-adic solution. After running this for different choices of the roots, we found specific combination of roots compatible with the p-adic data. Here we just present the outcome. The ToNumberField command takes surprisingly long time.

strue={t/.Solve[smallpolys[[1]]==0][[1]],t/.Solve[smallpolys[[2]]==0][[4]],RootReduce[t/.Solve[smallpolys[[3]]==0][[9]]],t/.Solve[smallpolys[[4]]==0][[1]],t/.Solve[smallpolys[[5]]==0][[2]]};

(* verification *)
For[m=1,m≤5,m++,
Print[m];
stry=strue[[m]];
tostudy=ToNumberField[stry,strue[[1]]];
deg=12;
(* the generator for some reason is 2072 strue[[1]] *)
Print[
{N[tostudy[[2]].Table[(2072 strue[[1]])^(k),{k,0,deg-1}],200]-N[stry,200],PolynomialMod[tostudy[[2]].Table[(2072 sold[[1]])^(k),{k,0,deg-1}] - sold[[m]],p^10]}];
];

Thus we know what the values are. Let’s compute the quadrics (with some accuracy).

strue={t/.Solve[smallpolys[[1]]==0][[1]],t/.Solve[smallpolys[[2]]==0][[4]],RootReduce[t/.Solve[smallpolys[[3]]==0][[9]]],t/.Solve[smallpolys[[4]]==0][[1]],t/.Solve[smallpolys[[5]]==0][[2]]}

prec=2000;
rules0=Table[fivevars[[k]]->N[strue[[k]],prec],{k,1,5}];
rules=rules0;
For[i=1,i≤8,i++,
rules=Join[rules,
Table[transfshort[[9-i]][[k]][[1]]-> (transfshort[[9-i]][[k]][[2]]/.rules)
,{k,1,Length[transfshort[[9-i]]]}]];
];
rulesqcl2=Table[varqcl2[[k]]->(varqcl2[[k]]/.rules),{k,1,7}];
quadtrue=qcl2num/.rulesqcl2;
pt1true=Table[uvs[[k]]/.pt1/.rules0,{k,1,8}];
pt2true=Table[uvs[[k]]/.pt2/.rules0,{k,1,8}];
pt4true=Table[uvs[[k]]/.pt4/.rules0,{k,1,8}];

(* reality check *)
JacMat=Table[D[quadtrue[[k]],uvs[[l]]],{k,1,9},{l,1,8}]/.{u0->1,u1->1};
JacMatFiber=Table[D[quadtrue[[k]],uvs[[l]]],{k,1,9},{l,3,8}]/.{u0->1,u1->1};
MatrixRank[JacMat/.Table[uvs[[k]]-> pt1true[[k]],{k,1,8}]]
MatrixRank[JacMat/.Table[uvs[[k]]-> pt2true[[k]],{k,1,8}]]
MatrixRank[JacMat/.Table[uvs[[k]]-> pt4true[[k]],{k,1,8}]]
MatrixRank[JacMatFiber/.Table[uvs[[k]]-> pt1true[[k]],{k,1,8}]]
MatrixRank[JacMatFiber/.Table[uvs[[k]]-> pt2true[[k]],{k,1,8}]]
MatrixRank[JacMatFiber/.Table[uvs[[k]]-> pt4true[[k]],{k,1,8}]]

fname ="Desktop/MumfordFPP2021/KeumFPP/LastPoint";
OpenRead[fname];
pt3true=uvs/.(Read[fname]/.rules);
Close[fname];
MatrixRank[JacMatFiber/.Table[uvs[[k]]-> pt3true[[k]],{k,1,8}]]
MatrixRank[JacMat/.Table[uvs[[k]]-> pt3true[[k]],{k,1,8}]]
(* perfect! *)

Now we will transport the four points to nice positions. This will hopefully reduce the field of definition to a quadratic one.

ws={w1,w2,w3,w4,w5,w6};
uws={u0,u1,w1,w2,w3,w4,w5,w6};

pt1sm=Table[pt1true[[k]]/pt1true[[3]],{k,3,8}];
pt2sm=Table[pt2true[[k]]/pt2true[[3]],{k,3,8}];
pt3sm=Table[pt3true[[k]]/pt3true[[3]],{k,3,8}];
pt4sm=Table[pt4true[[k]]/pt4true[[3]],{k,3,8}];

VarChange={v1->w1,v2-> l1 w2+l2 u1 w1, v3-> l3 w3 + l4 w4 + l5 u0^2 w1, v4-> l6 w3 + l7 w4 + l8 u0^2 w1, v5-> l9 w5 + l10 w6 + l11 u0 w2 + l12 u0 u1 w1,v6-> l13 w5 + l14 w6 + l15 u0 w2 + l16 u0 u1 w1}/.{u0->1,u1->1};
prec=1800;
MatrixVarChange=N[Table[Coefficient[vs[[k]]/.VarChange,ws[[l]]],{k,1,6},{l,1,6}],prec];
(*
standardpoint2={1,0,0,0,0,0};
standardpoint4={1,1,1,1,1,1};
standardpoint1={1,*,1,0,1,0};
standardpoint3={1,*,*,*,0,0};
*)
sls1=NSolve[Table[0==( MatrixVarChange.N[pt2sm,prec])[[k]],{k,2,6}]];
MatrixVarChange1=(MatrixVarChange/.sls1[[1]]);
newpt=MatrixVarChange1.N[pt4sm,prec];
sls2=NSolve[{newpt[[1]]==newpt[[3]],newpt[[1]]==newpt[[4]],newpt[[1]]==newpt[[5]],newpt[[1]]==newpt[[6]],newpt[[1]]==newpt[[2]]}];
MatrixVarChange2=(MatrixVarChange1/.sls2[[1]]);
newpt=MatrixVarChange2.N[pt1sm,prec];
sls3=NSolve[{newpt[[1]]==newpt[[3]],0==newpt[[4]],newpt[[1]]==newpt[[5]],0==newpt[[6]]}];
MatrixVarChange3=(MatrixVarChange2/.sls3[[1]]);
newpt=MatrixVarChange3.N[pt3sm,prec];
sls4=NSolve[{newpt[[5]]==0,0==newpt[[6]]}];
MatrixVarChange4=(MatrixVarChange3/.sls4[[1]]);

RootApproximant[Chop[MatrixVarChange4.pt1sm,10^(-1000)]]
RootApproximant[Chop[MatrixVarChange4.pt2sm,10^(-1000)]]
RootApproximant[Chop[MatrixVarChange4.pt3sm,10^(-1000)]]
RootApproximant[Chop[MatrixVarChange4.pt4sm,10^(-1000)]]



InvVC=Inverse[MatrixVarChange4];
VtoW={v1->w1,v2-> l1 w2+l2 u1 w1, v3-> l3 w3 + l4 w4 + l5 u0^2 w1, v4-> l6 w3 + l7 w4 + l8 u0^2 w1, v5-> l9 w5 + l10 w6 + l11 u0 w2 + l12 u0 u1 w1,v6-> l13 w5 + l14 w6 + l15 u0 w2 + l16 u0 u1 w1}/.{l1->InvVC[[2]][[2]],l2->InvVC[[2]][[1]],l3->InvVC[[3]][[3]],l4->InvVC[[3]][[4]],l5->InvVC[[3]][[1]],l6->InvVC[[4]][[3]],l7->InvVC[[4]][[4]],l8->InvVC[[4]][[1]],l9->InvVC[[5]][[5]],l10->InvVC[[5]][[6]],l11->InvVC[[5]][[2]],l12->InvVC[[5]][[1]],l13->InvVC[[6]][[5]],l14->InvVC[[6]][[6]],l15->InvVC[[6]][[2]],l16->InvVC[[6]][[1]]};

QuadricsW={}; (* let's compute new quadrics *)
Quadrics=quadtrue;

degs={{2,0},{3,0},{0,1},{3,1},{4,1},{4,1},{5,1},{5,1}};
as={a1,a2,a3,a4,a5,a6,a7,a8};

dega=as.degs;
count=0;
deg28=Solve[Join[{dega[[1]]==8,dega[[2]]==2},{a1≥0,a2≥0,a3≥0,a4≥0,a5≥0,a6≥0,a7≥0,a8≥0}],Integers];
deg29=Solve[Join[{dega[[1]]==9,dega[[2]]==2},{a1≥0,a2≥0,a3≥0,a4≥0,a5≥0,a6≥0,a7≥0,a8≥0}],Integers];
deg210=Solve[Join[{dega[[1]]==10,dega[[2]]==2},{a1≥0,a2≥0,a3≥0,a4≥0,a5≥0,a6≥0,a7≥0,a8≥0}],Integers];


qtry=(a1 Quadrics[[1]] + a2 Quadrics[[2]] + a3 Quadrics[[3]])/.VtoW;
sls=Solve[{1==Coefficient[qtry,w3^2],0==Coefficient[qtry,w3 w4],0==Coefficient[qtry,w4^2]}][[1]];
AppendTo[QuadricsW,
Sum[RootApproximant[Chop[Coefficient[qtry/.sls,Product[uws[[k]]^(as[[k]]/.deg28[[l]]),{k,1,8}]]]] Product[uws[[k]]^(as[[k]]/.deg28[[l]]),{k,1,8}],{l,1,Length[deg28]}]];
qtry=(a1 Quadrics[[1]] + a2 Quadrics[[2]] + a3 Quadrics[[3]])/.VtoW;
sls=Solve[{0==Coefficient[qtry,w3^2],1==Coefficient[qtry,w3 w4],0==Coefficient[qtry,w4^2]}][[1]];
AppendTo[QuadricsW,
Sum[RootApproximant[Chop[Coefficient[qtry/.sls,Product[uws[[k]]^(as[[k]]/.deg28[[l]]),{k,1,8}]]]] Product[uws[[k]]^(as[[k]]/.deg28[[l]]),{k,1,8}],{l,1,Length[deg28]}]];
qtry=(a1 Quadrics[[1]] + a2 Quadrics[[2]] + a3 Quadrics[[3]])/.VtoW;
sls=Solve[{0==Coefficient[qtry,w3^2],0==Coefficient[qtry,w3 w4],1==Coefficient[qtry,w4^2]}][[1]];
AppendTo[QuadricsW,
Sum[RootApproximant[Chop[Coefficient[qtry/.sls,Product[uws[[k]]^(as[[k]]/.deg28[[l]]),{k,1,8}]]]] Product[uws[[k]]^(as[[k]]/.deg28[[l]]),{k,1,8}],{l,1,Length[deg28]}]];

qtry=(a1 Quadrics[[4]] + a2 Quadrics[[5]] + a3 Quadrics[[6]])/.VtoW;
sls=Solve[{1==Coefficient[qtry,w3 w6],0==Coefficient[qtry,w3 w5],0==Coefficient[qtry,w4 w5]}][[1]];
AppendTo[QuadricsW,
Sum[RootApproximant[Chop[Coefficient[qtry/.sls,Product[uws[[k]]^(as[[k]]/.deg29[[l]]),{k,1,8}]]]] Product[uws[[k]]^(as[[k]]/.deg29[[l]]),{k,1,8}],{l,1,Length[deg29]}]];
qtry=(a1 Quadrics[[4]] + a2 Quadrics[[5]] + a3 Quadrics[[6]])/.VtoW;
sls=Solve[{0==Coefficient[qtry,w3 w6],1==Coefficient[qtry,w3 w5],0==Coefficient[qtry,w4 w5]}][[1]];
AppendTo[QuadricsW,
Sum[RootApproximant[Chop[Coefficient[qtry/.sls,Product[uws[[k]]^(as[[k]]/.deg29[[l]]),{k,1,8}]]]] Product[uws[[k]]^(as[[k]]/.deg29[[l]]),{k,1,8}],{l,1,Length[deg29]}]];
qtry=(a1 Quadrics[[4]] + a2 Quadrics[[5]] + a3 Quadrics[[6]])/.VtoW;
sls=Solve[{0==Coefficient[qtry,w3 w6],0==Coefficient[qtry,w3 w5],1==Coefficient[qtry,w4 w5]}][[1]];
AppendTo[QuadricsW,
Sum[RootApproximant[Chop[Coefficient[qtry/.sls,Product[uws[[k]]^(as[[k]]/.deg29[[l]]),{k,1,8}]]]] Product[uws[[k]]^(as[[k]]/.deg29[[l]]),{k,1,8}],{l,1,Length[deg29]}]];

qtry=(a1 Quadrics[[7]] + a2 Quadrics[[8]] + a3 Quadrics[[9]])/.VtoW;
sls=Solve[{1==Coefficient[qtry,w5^2],0==Coefficient[qtry,w5 w6],0==Coefficient[qtry,w6^2]}][[1]];
AppendTo[QuadricsW,
Sum[RootApproximant[Chop[Coefficient[qtry/.sls,Product[uws[[k]]^(as[[k]]/.deg210[[l]]),{k,1,8}]]]] Product[uws[[k]]^(as[[k]]/.deg210[[l]]),{k,1,8}],{l,1,Length[deg210]}]];
qtry=(a1 Quadrics[[7]] + a2 Quadrics[[8]] + a3 Quadrics[[9]])/.VtoW;
sls=Solve[{0==Coefficient[qtry,w5^2],1==Coefficient[qtry,w5 w6],0==Coefficient[qtry,w6^2]}][[1]];
AppendTo[QuadricsW,
Sum[RootApproximant[Chop[Coefficient[qtry/.sls,Product[uws[[k]]^(as[[k]]/.deg210[[l]]),{k,1,8}]]]] Product[uws[[k]]^(as[[k]]/.deg210[[l]]),{k,1,8}],{l,1,Length[deg210]}]];
qtry=(a1 Quadrics[[7]] + a2 Quadrics[[8]] + a3 Quadrics[[9]])/.VtoW;
sls=Solve[{0==Coefficient[qtry,w5^2],0==Coefficient[qtry,w5 w6],1==Coefficient[qtry,w6^2]}][[1]];
AppendTo[QuadricsW,
Sum[RootApproximant[Chop[Coefficient[qtry/.sls,Product[uws[[k]]^(as[[k]]/.deg210[[l]]),{k,1,8}]]]] Product[uws[[k]]^(as[[k]]/.deg210[[l]]),{k,1,8}],{l,1,Length[deg210]}]];

QuadricsW

fname ="Desktop/MumfordFPP2021/KeumFPP/QuadricsWforKeum";
OpenWrite[fname];
Write[fname, QuadricsW];
Close[fname];

Now we will find S1 and S2 and then the 7-fold cover. The method is the same as before.

(* start by finding lines in the fiber *)
pt1w=RootApproximant[Chop[MatrixVarChange4.pt1sm,10^(-1000)]]
pt2w=RootApproximant[Chop[MatrixVarChange4.pt2sm,10^(-1000)]]
pt3w=RootApproximant[Chop[MatrixVarChange4.pt3sm,10^(-1000)]]
pt4w=RootApproximant[Chop[MatrixVarChange4.pt4sm,10^(-1000)]]

line12=Join[{u0->1,u1->1},Table[ws[[k]]-> pt1w[[k]]+ t pt2w[[k]],{k,1,6}]];
line34=Join[{u0->1,u1->1},Table[ws[[k]]-> pt3w[[k]]+ t pt4w[[k]],{k,1,6}]];
line13=Join[{u0->1,u1->1},Table[ws[[k]]-> pt1w[[k]]+ t pt3w[[k]],{k,1,6}]];
{Expand[QuadricsW/.line12],Expand[QuadricsW/.line34],Expand[QuadricsW/.line13]} (* sanity check *)

Let us now parameterize S and compute the intersection points of S with the singular fiber. We readily see two points on lines and two pairs of points on conics. Then we will compute a conic through the A1 and A2 singularity.

ParamS={u0-> a1 + a2 t + a3 t^2, u1-> a4 + a5 t + a6 t^2 + a7 t^3, w1-> 0, w2-> 1, w3-> 1, w4-> t, w5-> a8 + a9 t + a10 t^2, w6-> a11 + a12 t + a13 t^2};
eqs={};
For[i=1,i≤9,i++,
tbz=QuadricsW[[i]]/.ParamS;
For[k=0,k≤Exponent[tbz,t],k++,
AppendTo[eqs,0== SeriesCoefficient[tbz,{t,0,k}]];
];
];
ParamSsolved=ParamS/.RootReduce[Solve[eqs]][[1]]

solSfiber=RootReduce[Solve[(u0^3==u1^2)/.ParamSsolved]]

ptsS={};
For[i=1,i≤6,i++,
pt=RootReduce[uws/.ParamSsolved/.solSfiber[[i]]];
lambda=(pt[[2]]/pt[[1]]);
lambda2=pt[[3]];If[lambda2==0,lambda2=pt[[4]]/lambda^3];
pt=RootReduce[{pt[[1]]/lambda^2,pt[[2]]/lambda^3,pt[[3]]/lambda2,pt[[4]]/lambda^3/lambda2,pt[[5]]/lambda^4/lambda2,pt[[6]]/lambda^4/lambda2,pt[[7]]/lambda^5/lambda2,pt[[8]]/lambda^5/lambda2}];
AppendTo[ptsS,Table[pt[[k]],{k,3,8}]];
];

{MatrixRank[{ptsS[[1]],pt1w,pt2w}],MatrixRank[{ptsS[[2]],pt1w,pt3w}],MatrixRank[{ptsS[[3]],pt3w,pt4w}]} (* this tells us where the points lie *)

We will now look for the other curves. We know from the intersection data that S2 should be cut out by a (1,8) divisor which vanishes three times on line12, and once on line13 and  line 34. Let us figure out at least the double vanishing condition. This means that the equation vanishes on the tangent space. We can ensure it by looking at several points on the line.

Equation18=a1 u0^4 w1 + a2 u0 u1^2 w1 + a3 u0 u1 w2 + a4 u0^2 w3 + a5 u0^2 w4 + a6 u1 w5 + a7 u1 w6;
Equation18=RootReduce[Equation18/.Solve[{1==a1, 0== Equation18/.line12/.t->0, 0==  Equation18/.line12/.t->1
, 0== Equation18/.line13/.t->0, 0==  Equation18/.line13/.t->1
,0== Equation18/.line34/.t->0, 0==  Equation18/.line34/.t->1}][[1]]];

eqs={};
For[m=2,m≤4,m++,
Print[m];
pt=RootReduce[line12/.t->m];
TanSp=RootReduce[NullSpace[Table[D[QuadricsW[[k]],uws[[l]]]/.pt,{k,1,9},{l,1,8}]]];
For[i=1,i≤4,i++,
AppendTo[eqs,0== Sum[(D[Equation18,uws[[k]]]/.pt) TanSp[[i]][[k]],{k,1,8}]];
];
];
Equation18=RootReduce[Equation18/.Solve[eqs][[1]]] (* we get the desired equation *)

We will compute the parametric equation of S2 following previous paper. For a given value of t we can compute the corresponding point on S and six points of S2 in that fiber. If we knew which of the six points is the right one, we could use several such calculations to find the coefficients of the parametrization.

One way to determine which points are right is by looking at field extensions.  We expect that the correct point on S2 will be defined over the quadratic imaginary field.

correctpoints={};
For[m=2,m≤10,m++,
Print[m];
samplet=t->m;
sixpoints=RootApproximant[Chop[uws/.NSolve[Join[Table[0== QuadricsW[[k]],{k,1,9}],{0== Equation18,0== w1-1,0== u0-(u0/.ParamSsolved/.samplet), 0== u1- (u1/.ParamSsolved/.samplet)}],2000]],12];
For[k=1,k≤6,k++,
If[Exponent[MinimalPolynomial[sixpoints[[k]][[5]],s],s]≤2,
AppendTo[correctpoints,sixpoints[[k]]];
];
];
]

ParamS2solved={}; (* now we can easily solve for S2 and verify *)
For[k=1,k≤8,k++,
solS2=Solve[Table[correctpoints[[m-1]][[k]] == a1 + a2 m + a3 m^2 + a4 m^3 + a5 m^4 + a6 m^5,{m,2,10}]];
AppendTo[ParamS2solved,uws[[k]]->  (a1 + a2 t + a3 t^2 + a4 t^3 + a5 t^4 + a6 t^5)/.solS2[[1]]];
];
ParamS2solved
Factor[QuadricsW/.ParamS2solved]

We will now construct the curve S1 as well as an order three automorphism of the surface. The curve S1 could be constructed as follows. For a given value of t we have points on S and S2. We can employ the group law on the elliptic curve to construct the corresponding point on S1. The idea is to pick a hyperplane that’s tangent to the point on S, then compute four extra intersection points, then take hyperplane through these points and the point from S2. It’s best to do it numerically, with subsequent RootApproximant command.

correctpointsS1={};
For[m=2,m≤10,m++,
Print[m];
samplet=t->m;
ptS=uws/.ParamSsolved/.samplet;
ns=RootReduce[NullSpace[Table[D[QuadricsW[[l]]/.{u0-> ptS[[1]],u1->ptS[[2]]},ws[[k]]]/.ParamSsolved/.samplet,{l,1,9},{k,1,6}]]];
tanpl=a1 w2 +a2 w3 + w4;
tanplsolved=tanpl/.Solve[{Sum[(D[tanpl,ws[[k]]]/.ParamSsolved/.samplet)ns[[1]][[k]],{k,1,6}]==0,0== tanpl/.ParamSsolved/.samplet}][[1]];
sols=NSolve[Join[Table[0== QuadricsW[[k]],{k,1,9}],{0==tanplsolved,0== w2-1,0== u0-(u0/.ParamSsolved/.samplet), 0== u1- (u1/.ParamSsolved/.samplet)}],2000];
fourextra={};
For[k=1,k≤6,k++,
If[Chop[w1/.sols[[k]]]≠0,
AppendTo[fourextra,sols[[k]]];
];
];
fivepts=Append[fourextra,ParamS2solved/.samplet];
planew=a1 w1 + a2 w2 +a3 w3 + a4 w4 + a5 w5 + w6;
planewsolved=planew/.Solve[Table[0==planew/.fivepts[[k]],{k,1,5}]][[1]];
pluss1 = 
NSolve[Join[Table[0== QuadricsW[[k]],{k,1,9}],{0==planewsolved,0== w1-1,0== u0-ptS[[1]], 0== u1- ptS[[2]]}],1900];
For[k=1,k≤6,k++,
w2k=RootApproximant[Chop[w2/.pluss1[[k]]],12];
If[Exponent[MinimalPolynomial[w2k,s],s]≤2,
If[w2k≠ (w2/.ParamS2solved/.samplet),
AppendTo[correctpointsS1,RootApproximant[Chop[uws/.pluss1[[k]]]]];
]];
];
];

ParamS1solved={}; (* now we can easily solve for S1 and verify *)
For[k=1,k≤8,k++,
solS1=Solve[Table[correctpointsS1[[m-1]][[k]] == a1 + a2 m + a3 m^2 + a4 m^3 + a5 m^4 + a6 m^5,{m,2,10}]];
AppendTo[ParamS1solved,uws[[k]]->  (a1 + a2 t + a3 t^2 + a4 t^3 + a5 t^4 + a6 t^5)/.solS1[[1]]];
];
ParamS1solved
Factor[QuadricsW/.ParamS1solved]

Let us now compute the degree (1,10) equation that vanishes on S1.

Equation110=a1 u0^5 w1 + a2 u0^2 u1^2 w1 +  u0^2 u1 w2 + a4 u0^3 w3 + a5 u1^2 w3 + a6 u0^3 w4 + a7 u1^2 w4 + a8 u0 u1 w5 + a9 u0 u1 w6; (* scale to have a3 =1 *)
tbz=Factor[Equation110/.ParamS1solved];
Equation110=Equation110/.RootReduce[Solve[Table[0==SeriesCoefficient[tbz,{t,0,k}],{k,0,10}]]][[1]]
Factor[Equation110/.ParamS1solved]

We would like to now compute the order three automorphism. Perhaps the best way of thinking about it is that we have a curve over C[u0,u1] and we would like to find its automorphism. Moreover, we know that, up to scaling, w1 -> Equation18 -> Equation110->w1.

Order3Matrices={};
For[m=2,m≤30,m++,
If[Mod[m,5]==0,Print[m]];
samplet=t->m;
sixt=Solve[((u1^2/u0^3)/.ParamSsolved )==( (u1^2/u0^3)/.ParamSsolved /.samplet)];
For[k=1,k≤6,k++,
If[m== t/.sixt[[k]], goodindex=k];
];
sixt=N[sixt,2000];

ptsS=Table[uws/.ParamSsolved/.sixt[[k]],{k,1,6}];
ptsS1=Table[uws/.ParamS1solved/.sixt[[k]],{k,1,6}];
ptsS2=Table[uws/.ParamS2solved/.sixt[[k]],{k,1,6}];

(* normalize to the good point *)
lambda0=( ptsS[[goodindex]][[1]]/ptsS[[goodindex]][[2]]);

For[k=1,k≤6,k++,
pt=ptsS[[k]];
lambda=(pt[[1]]/pt[[2]] )/lambda0;
ptsS[[k]]=Table[uws[[l]]-> {pt[[1]] lambda^2, pt[[2]] lambda^3,pt[[3]],pt[[4]]lambda^3, pt[[5]] lambda^4, pt[[6]] lambda^4, pt[[7]] lambda^5, pt[[8]] lambda^5}[[l]],{l,1,8}];
pt=ptsS1[[k]];
lambda=pt[[1]]/pt[[2]]/lambda0;
ptsS1[[k]]=Table[uws[[l]]-> {pt[[1]] lambda^2, pt[[2]] lambda^3,pt[[3]],pt[[4]]lambda^3, pt[[5]] lambda^4, pt[[6]] lambda^4, pt[[7]] lambda^5, pt[[8]] lambda^5}[[l]],{l,1,8}];
pt=ptsS2[[k]];
lambda=pt[[1]]/pt[[2]]/lambda0;
ptsS2[[k]]=Table[uws[[l]]-> {pt[[1]] lambda^2, pt[[2]] lambda^3,pt[[3]],pt[[4]]lambda^3, pt[[5]] lambda^4, pt[[6]] lambda^4, pt[[7]] lambda^5, pt[[8]] lambda^5}[[l]],{l,1,8}];
];

AutoOrder3 = Table[ws[[k]]-> Sum[ToExpression[StringJoin["a",ToString[6 k + l-6]]]ws[[l]],{l,1,6}],{k,1,6}];

Eqs={};
For[k=1,k≤6,k++,
(* w2 will be nonzero *)
pt1=ptsS[[k]];
pt2=ptsS1[[k]];
pt3=ptsS2[[k]];
newpt1=Table[ws[[l]]->(ws[[l]]/. AutoOrder3/.pt1),{l,1,6}];
newpt2=Table[ws[[l]]->(ws[[l]]/. AutoOrder3/.pt2),{l,1,6}];
newpt3=Table[ws[[l]]->(ws[[l]]/. AutoOrder3/.pt3),{l,1,6}];
Eqs=Join[Eqs,Flatten[Minors[ws/.{pt1,newpt3},2]]];
Eqs=Join[Eqs,Flatten[Minors[ws/.{pt2,newpt1},2]]];
Eqs=Join[Eqs,Flatten[Minors[ws/.{pt3,newpt2},2]]];
];
sols=Solve[Table[0==Eqs[[k]],{k,1,Length[Eqs]}]];
AutoOrder3=N[AutoOrder3/.sols[[1]]/.(a1-> Coefficient[Equation18,w1]/.{u0->(u0/.ParamSsolved/.samplet),u1->(u1/.ParamSsolved/.samplet)}),1900];
AppendTo[Order3Matrices,AutoOrder3];
];

wts={0,3,4,4,5,5};
Order3Computed={u0->u0,u1->u1};
For[i=1,i≤6,i++,
wito=0;
For[j=1,j≤6,j++,
diff=wts[[i]]-wts[[j]]+8;
powers=Solve[{2a1+3a2==diff,a1≥0,a2≥0},Integers];
term=Sum[ToExpression[StringJoin["c",ToString[k]]]u0^(a1/.powers[[k]])u1^(a2/.powers[[k]]), {k,1,Length[powers]}];
sls=Solve[Table[(term/.ParamSsolved/.t->m+1)==Coefficient[(ws[[i]]/.Order3Matrices[[m]] ),ws[[j]]],{m,2,29}]];
wito+=ws[[j]] Sum[RootApproximant[Chop[ToExpression[StringJoin["c",ToString[k]]]/.sls[[1]]]]u0^(a1/.powers[[k]])u1^(a2/.powers[[k]]), {k,1,Length[powers]}];
];
AppendTo[Order3Computed,ws[[i]]->wito];
];
(* we scale it so that the cube is simple *)
Order3Computed =RootReduce[ Join[{Order3Computed[[1]],Order3Computed[[2]]},Table[ws[[k]]-> (ws[[k]]/.Order3Computed) /(((3-I (Sqrt[7]))/8)(u0^3-u1^2)),{k,1,6}]]];
FullSimplify[((ws/.Order3Computed)/.Order3Computed)/.Order3Computed]

fname ="Desktop/MumfordFPP2021/KeumFPP/Order3ComputedKeum";
OpenWrite[fname];
Write[fname,Order3Computed];
Close[fname];

checking={};
For[k=1,k≤9,k++,
blah=QuadricsW[[k]]/.Order3Computed;
AppendTo[checking,Expand[blah-Coefficient[blah,w6^2]QuadricsW[[9]]-Coefficient[blah,w6 w5]QuadricsW[[8]] -Coefficient[blah,w5^2]QuadricsW[[7]] - Coefficient[blah, w4 w5]QuadricsW[[6]] - Coefficient[blah,w3 w5] QuadricsW[[5]] - Coefficient[blah,w3 w6]QuadricsW[[4]] - Coefficient[blah, w4^2]QuadricsW[[3]] -Coefficient[blah, w4 w3]QuadricsW[[2]]-Coefficient[blah, w3^2]QuadricsW[[1]]]]
];
checking (* no surprises *)

fname ="Desktop/MumfordFPP2021/KeumFPP/Order3ComputedKeum";
OpenRead[fname];
Order3Computed=Read[fname];
Close[fname];

The following rational function has orders that are divisible by 7 at the three lines in I9 fiber that do not get contracted to (1,3)/7 singular points as well as the u0=0 divisor. The divisor of it is
(-2B -C  -B1 -4C1 -7A2 - 4 B2 - 2 C2 - 3S +2 S1 +S2  + 7 F2).
We also made some adjustments and choices to have the Z_3 action extend to the seventh root of the function.

SeventhPower=  Factor[(-198144-78336 I (Sqrt[7]))(Equation18/(w1 u0^4)) (Equation110/(w1 u0^5))^2  (u0^21)/(u0^3-u1^2)^7 ];
SeventhPowerSigma=SeventhPower/.Order3Computed;
Factor[SeventhPowerSigma/ SeventhPower^4,Extension->I (Sqrt[7])] ;
RootReduce[Factor[(SeventhPowerSigma/SeventhPower^4) /((1/32) (3+7 I (Sqrt[7]))u0^(-3)((u0^3)-(u1^2))^4 w1^2/((4 I (u0^4) w1-4 I u0 (u1^2) w1+(18 I+6 (Sqrt[7])) u0 u1 w2+(-3 I+(Sqrt[7])) (u0^2) w3+(5 I-3 (Sqrt[7])) (u0^2) w4+(-I+3 (Sqrt[7])) u1 w5+(-19 I-7 (Sqrt[7])) u1 w6)(66 I (u0^2) u1 w2+(-8 I+4 (Sqrt[7])) (u0^3) w3+(2 I+10 (Sqrt[7])) (u1^2) w3+(31 I+(Sqrt[7])) (u0^3) w4+(34 I-6 (Sqrt[7])) (u1^2) w4+(17 I-3 (Sqrt[7])) u0 u1 w5+(-142 I-6 (Sqrt[7])) u0 u1 w6)))^7,Extension->I (Sqrt[7])] ]
(* this tells us how to change g *)
Order3ComputedG=Append[Order3Computed,g->g^4((1/32) (3+7 I (Sqrt[7]))u0^(-3)((u0^3)-(u1^2))^4w1^2/((4 I (u0^4) w1-4 I u0 (u1^2) w1+(18 I+6 (Sqrt[7])) u0 u1 w2+(-3 I+(Sqrt[7])) (u0^2) w3+(5 I-3 (Sqrt[7])) (u0^2) w4+(-I+3 (Sqrt[7])) u1 w5+(-19 I-7 (Sqrt[7])) u1 w6)(66 I (u0^2) u1 w2+(-8 I+4 (Sqrt[7])) (u0^3) w3+(2 I+10 (Sqrt[7])) (u1^2) w3+(31 I+(Sqrt[7])) (u0^3) w4+(34 I-6 (Sqrt[7])) (u1^2) w4+(17 I-3 (Sqrt[7])) u0 u1 w5+(-142 I-6 (Sqrt[7])) u0 u1 w6)))];
Factor[((g^7/SeventhPower)/.Order3ComputedG)/ (g^28/SeventhPower^4),Extension->I (Sqrt[7])]

Here is the intersection matrix in the basis A,B,C,...,C2, S, S1, S2, F3, F2.
We have relations that correspond to Equation18/(w1 u0^4), Equation110/(w1 u0^5) and SeventhPower.

M1={{-2,1,0,0,0,0,0,0,1,1,3,1,0,0},{1,-2,1,0,0,0,0,0,0,1,0,0,0,0},{0,1,-2,1,0,0,0,0,0,0,0,0,0,0},{0,0,1,-2,1,0,0,0,0,1,1,3,0,0},{0,0,0,1,-2,1,0,0,0,0,1,0,0,0},{0,0,0,0,1,-2,1,0,0,0,0,0,0,0},{0,0,0,0,0,1,-2,1,0,3,1,1,0,0},{0,0,0,0,0,0,1,-2,1,0,0,1,0,0},{1,0,0,0,0,0,0,1,-2,0,0,0,0,0},{1,1,0,1,0,0,3,0,0,-3,0,0,3,2},{3,0,0,1,1,0,1,0,0,0,-3,0,3,2},{1,0,0,3,0,0,1,1,0,0,0,-3,3,2},{0,0,0,0,0,0,0,0,0,3,3,3,0,0},{0,0,0,0,0,0,0,0,0,2,2,2,0,0}};
nm1=NullSpace[M1];
nm1//MatrixForm

M1[[10]]
relS2=-(nm1[[3]]/3 -4  nm1[[1]]/3 )
relS1=-(nm1[[4]]/3 - 5 nm1[[1]]/3 )

relg7=-(relS2 + 2 relS1 -7 nm1[[1]])

(* this is the canonical divisor of the singular model *)
KSing={0,2/7,1/7,0,2/7,1/7,0,2/7,1/7,3/7,3/7,3/7,0,1/2};

The canonical divisor on the FPP cover is given as the pullback of 
F  + 2/7 B + 1/7 C + 2/7 B1 _ 1/7 C1 + 2/7 B2 + 1/7 C2  + 3/7 S + 3/7 S1 + 3/7 S2. 
We can look at the global sections of 2K on the blowup of the FPP and push them forward.

(* here we look to compute the conditions on the pushforwards of the sections of the 
2K from the fake projective plane *)
KSing={0,2/7,1/7,0,2/7,1/7,0,2/7,1/7,3/7,3/7,3/7,0,1/2};
Divg=relg7/7;
For[m=-6,m≤6,m++,
d=-Ceiling[-2KSing- m Divg];
(* print RR data *)
Print[{d,m,(d.M1.d-d.M1.{0,0,0,0,0,0,0,0,0,0,0,0,0,1/2})/2+1}];
]

Let us focus on {{0,0,0,0,0,0,1,1,0,1,0,0,0,0},-1,1}. This means that g^(-1) can be used for this  variable. 
Eigenspace of wt 0:   1.

FunctionWt0=1;
FunctionWtminus1=g^(-1);

The {{0,1,0,0,1,2,3,2,1,2,0,0,0,-2},-3,2} is also interesting. It looks like we should be looking for a function f g^(-3) u0^9/ (u0^3 - u1^2)^3 with 
f = s/(w1^2 u0^7) with section s of weight (2,14) that vanishes at 3A, 2B and 3A1. This will hopefully be easy to find, although the dimension of (2,14) space is somewhat large.

{0,1,0,0,1,2,3,2,1,2,0,0,0,-2} + 3 nm1[[1]]
monsdeg2=Complement[DeleteDuplicates[Flatten[Table[ws[[i]]ws[[j]],{i,1,6},{j,1,6}]]],{w3^2,w3 w4, w4^2, w4 w5, w3 w5, w3 w6, w5^2,w5 w6, w6^2}]; (* removing the ones we can get rid of via equations *)
monsuniv={};
count=0;
UnivPoly=0;
For[k=1,k≤Length[monsdeg2],k++,
deg=14-Sum[wts[[l]] Exponent[monsdeg2[[k]],ws[[l]]],{l,1,6}];
sols=Solve[{2 a1 + 3 a2 ==deg,a1≥0,a2≥ 0},Integers]; 
For[l=1,l≤Length[sols],l++,
UnivPoly+=monsdeg2[[k]] u0^(a1/.sols[[l]]) u1^(a2/.sols[[l]]) ToExpression[StringJoin["c",ToString[count]]];
AppendTo[monsuniv,monsdeg2[[k]] u0^(a1/.sols[[l]]) u1^(a2/.sols[[l]]) ];
count++;
];
];
UnivPoly

(* the A1 is line12, B is line13 and A is line24 *)
(* Looking for double tangency at all three lines for starters *)

eqs={};
For[m=2,m≤6,m++,
pt=RootReduce[line13/.t->m];
AppendTo[eqs,0== UnivPoly/.pt];
TanSp=RootReduce[NullSpace[Table[D[QuadricsW[[k]],uws[[l]]]/.pt,{k,1,9},{l,1,8}]]];
For[i=1,i≤4,i++,
AppendTo[eqs,0== Sum[(D[UnivPoly,uws[[k]]]/.pt) TanSp[[i]][[k]],{k,1,8}]];
];
];
UnivPoly=RootReduce[UnivPoly/.Solve[eqs][[1]]];
Print[Length[Complement[Variables[UnivPoly],uws]]" coeffs"];


eqs={};
For[m=2,m≤6,m++,
pt=RootReduce[line12/.t->m];
AppendTo[eqs,0== UnivPoly/.pt];
TanSp=RootReduce[NullSpace[Table[D[QuadricsW[[k]],uws[[l]]]/.pt,{k,1,9},{l,1,8}]]];
For[i=1,i≤4,i++,
AppendTo[eqs,0== Sum[(D[UnivPoly,uws[[k]]]/.pt) TanSp[[i]][[k]],{k,1,8}]];
];
];
UnivPoly=RootReduce[UnivPoly/.Solve[eqs][[1]]];
Print[Length[Complement[Variables[UnivPoly],uws]]" coeffs"];

eqs={};
For[m=2,m≤6,m++,
pt=RootReduce[line34/.t->m];
AppendTo[eqs,0== UnivPoly/.pt];
TanSp=RootReduce[NullSpace[Table[D[QuadricsW[[k]],uws[[l]]]/.pt,{k,1,9},{l,1,8}]]];
For[i=1,i≤4,i++,
AppendTo[eqs,0== Sum[(D[UnivPoly,uws[[k]]]/.pt) TanSp[[i]][[k]],{k,1,8}]];
];
];
UnivPoly=RootReduce[UnivPoly/.Solve[eqs][[1]]];
Print[Length[Complement[Variables[UnivPoly],uws]]" coeffs"];

Variables[UnivPoly]

(* now we will deal with the order three vanishing along Line12 and line34*)
Eqs={};
coeffs=Table[ToExpression[StringJoin["r",ToString[k]]],{k,1,50}]
For[m=3,m≤9,m++,
Print[m];
pt=line12/.t->m;
curve=Join[{u0->1,w1->(w1/.pt), u1->1+eps,w2-> (w2/.pt)},Table[ws[[k]]->(ws[[k]]/.pt)+ Sum[coeffs[[6 k+l-18]] eps^l,{l,1,6}],{k,3,6}]];
curve=curve/.NSolve[Flatten[Table[0== SeriesCoefficient[QuadricsW[[k]]/.curve,{eps,0,l}],{k,1,9},{l,0,6}]],2000][[1]];
AppendTo[ Eqs,SeriesCoefficient[UnivPoly/.curve,{eps,0,0}]];
AppendTo[ Eqs,SeriesCoefficient[UnivPoly/.curve,{eps,0,1}]];
AppendTo[ Eqs,SeriesCoefficient[UnivPoly/.curve,{eps,0,2}]];
];
sls=Solve[Table[0==Eqs[[k]],{k,1,Length[Eqs]}]][[1]];
UnivPoly=UnivPoly/.sls;

Eqs={};
For[m=3,m≤9,m++,
Print[m];
pt=line34/.t->m;
curve=Join[{u0->1,w1->(w1/.pt), u1->1+eps,w2-> (w2/.pt)},Table[ws[[k]]->(ws[[k]]/.pt)+ Sum[coeffs[[6 k+l-18]] eps^l,{l,1,6}],{k,3,6}]];
curve=curve/.NSolve[Flatten[Table[0== SeriesCoefficient[QuadricsW[[k]]/.curve,{eps,0,l}],{k,1,9},{l,0,6}]],2000][[1]];
AppendTo[ Eqs,SeriesCoefficient[UnivPoly/.curve,{eps,0,0}]];
AppendTo[ Eqs,SeriesCoefficient[UnivPoly/.curve,{eps,0,1}]];
AppendTo[ Eqs,SeriesCoefficient[UnivPoly/.curve,{eps,0,2}]];
];
sls=Solve[Table[0==Eqs[[k]],{k,1,Length[Eqs]}]][[1]];
UnivPolySolved=UnivPoly/.sls;
Complement[Variables[UnivPolySolved],uws]

section1=Factor[Sum[RootApproximant[Chop[Coefficient[UnivPolySolved,c0 monsuniv[[k]]]]]  monsuniv[[k]],{k,1,Length[monsuniv]}]];
section2=Factor[Sum[RootApproximant[Chop[Coefficient[UnivPolySolved,c10  monsuniv[[k]]]]]  monsuniv[[k]],{k,1,Length[monsuniv]}]];

FunctionWtminus31=section1 g^(-3) u0^2/(u0^3-u1^2)^3/ w1^2;
FunctionWtminus32=section2 g^(-3) u0^2/(u0^3-u1^2)^3/ w1^2;

Order3ComputedG=Factor[Order3ComputedG];
Print[ByteCount[%]];
FunctionWtminus1n=Factor[FunctionWtminus1/.Order3ComputedG];
Print[ByteCount[%]];
FunctionWtminus1nn=Factor[FunctionWtminus1n/.Order3ComputedG];
Print[ByteCount[%]];
FunctionWtminus31n=Factor[FunctionWtminus31/.Order3ComputedG];
Print[ByteCount[%]];
FunctionWtminus31nn=Factor[FunctionWtminus31n/.Order3ComputedG];
Print[ByteCount[%]];
FunctionWtminus32n=Factor[FunctionWtminus32/.Order3ComputedG];
Print[ByteCount[%]];
FunctionWtminus32nn=Factor[FunctionWtminus32n/.Order3ComputedG];
TenFunctions={1,FunctionWtminus1,FunctionWtminus1n,FunctionWtminus1nn,
FunctionWtminus31,FunctionWtminus31n,FunctionWtminus31nn,
FunctionWtminus32,FunctionWtminus32n,FunctionWtminus32nn};
Mod7Weights={0,-1,-4,-2,-3,-5,-6,-3,-5,-6};

We are now ready to compute the equations. The process is simple: we will compute a lot of values of uw-s and g to high accuracy, which will then give us values of ten functions to high accuracy. We will then compute cubic relations. If the run gives an error, redo it.

UWGpts={};

Ps={P0,P1,P2,P3,P4,P5,P6,P7,P8,P9};

prec=2000;
For[n=1,n≤50,n++,
If[Mod[n,10]==0, Print[n]];
sols=NSolve[Join[{u0==1,w1==1,u1==RandomInteger[30]+ I RandomInteger[30],w2==RandomInteger[30]+ I RandomInteger[30]},Table[0==QuadricsW[[k]],{k,1,9}]],prec];
For[i=1,i≤6,i++,
gvalue=g->(g/.Solve[g^7== SeventhPower/.sols[[i]]][[1]]);
AppendTo[UWGpts,Append[sols[[i]],gvalue]];
];
];
Ppoints=Table[Table[Ps[[n]]-> (TenFunctions[[n]]/.UWGpts[[k]]),{n,1,10}],{k,1,Length[UWGpts]}];

deg=3;
EqsFPP={};

powers={a1,a2,a3,a4,a5,a6,a7,a8,a9,a10};
allsols=
Solve[Join[Table[powers[[k]]≥0,{k,1,10}],{Sum[powers[[k]],{k,1,10}]==deg}],Integers];

For[mod7=0,mod7≤6,mod7++,
Print["mod7=",mod7];
sols={};
For[k=1,k≤Length[allsols],k++,
If[Mod[Sum[Mod7Weights[[k]] powers[[k]],{k,1,10}]/.allsols[[k]],7]==mod7,
AppendTo[sols,allsols[[k]]];
];
];
Monomials=Table[Product[Ps[[n]]^(powers[[n]]/.sols[[k]]),{n,1,10}],{k,1,Length[sols]}];
TableOfValues=Table[Monomials[[k]]/.Ppoints[[l]],{l,1,Length[Ppoints]},{k,1,Length[Monomials]}];
ns=RootApproximant[Chop[RowReduce[NullSpace[TableOfValues]],10^(-900)]];
EqsFPP=Join[EqsFPP,Table[Sum[Monomials[[k]]ns[[i]][[k]],{k,1,Length[Monomials]}],{i,1,Length[ns]}]];
];
Length[EqsFPP]
Table[ByteCount[EqsFPP[[k]]],{k,1,Length[EqsFPP]}]

EqsFPPn=Numerator[Factor[EqsFPP]];
Table[ByteCount[EqsFPPn[[k]]],{k,1,Length[EqsFPP]}]

fname ="Desktop/MumfordFPP2021/KeumFPP/EqsFPPnumKeum";
OpenWrite[fname];
Write[fname,EqsFPPn];
Close[fname];

ToRR[x_]:=
Ret=Rationalize[Chop[N[Re[x],1900],10^(-1800)]]+rr Rationalize[Chop[N[Im[x]/(Sqrt[7]),1900],10^(-1800)]];

AllMonomials=Table[Product[Ps[[n]]^(powers[[n]]/.allsols[[k]]),{n,1,10}],{k,1,Length[allsols]}];
EqsFPPwithrr=Table[Sum[ToRR[Coefficient[EqsFPPn[[k]],AllMonomials[[l]]]]AllMonomials[[l]],{l,1,Length[AllMonomials]}],{k,1,Length[EqsFPPn]}];
fname ="Desktop/MumfordFPP2021/KeumFPP/EqsFPPwithrrKeum";
OpenWrite[fname];
Write[fname,EqsFPPwithrr];
Close[fname];