In this notebook I will study the equations of the special Dolgachev surface and then construct the minimal model of the 7-fold cover that would be a fake projective plane.

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

We will compute the basics, such as the singular points.

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

pt2=Table[uws[[k]]->1,{k,1,8}] ;(* this is a smooth intersection of a conic and a line in the singular fiber *)
Expand[QuadricsW/.pt2] 
MatrixRank[Table[D[QuadricsW[[k]],uws[[l]]]/.pt2,{k,1,9},{l,1,8}]]
pt1={u0->1,u1->1,w1->1,w2->0,w3->0,w4->0,w5->0,w6->0};
Expand[QuadricsW/.pt1] (* this is a singular point *) 
MatrixRank[Table[D[QuadricsW[[k]],uws[[l]]]/.pt1,{k,1,9},{l,1,8}]]

We find the lines in the I9 fiber.

LineInFiber=Join[{u0->1,u1->1},Table[ws[[k]]-> ({1,0,a1,a2,a3,a4}+t{0,1,a5,a6,a7,a8})[[k]],{k,1,6}]];
eqs={};
For[i=1,i≤9,i++,
AppendTo[eqs,0== Coefficient[QuadricsW[[i]]/.LineInFiber,t^2]];
AppendTo[eqs,0== Coefficient[QuadricsW[[i]]/.LineInFiber,t]];
AppendTo[eqs,0== QuadricsW[[i]]/.LineInFiber/.t->0];
];
sols=RootReduce[Solve[eqs]]
Line1=LineInFiber/.sols[[1]];
Line2=LineInFiber/.sols[[2]];

We will now write down the parametric equation of S.

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

Let us now 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.

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

ptsS[[1]](*intersection with the line1 *)

SpanP2=Table[ws[[m]]->RootReduce[a1(ws/.pt1 ) + a2 (ptsS[[5]]+ptsS[[6]])
+ a3 ((ptsS[[5]]-ptsS[[6]])) /(NumericalApproximation[Root, (RowDefault["-0.523", FEExpressionsNumericalApproximationElider])+(RowDefault["0.375", FEExpressionsNumericalApproximationElider]) I, Root[11-17 ((#1)^2)+64 ((#1)^4)&,2], ])][[m]],{m,1,6}] (* this makes it defined over the quadratic field *)
EqP2=Factor[(QuadricsW/.{u0->1,u1->1})/.SpanP2,Extension->I (Sqrt[7])];
(* this verifies that all of these are multiples of one equation *)
Expand[Table[EqP2[[k]] Coefficient[EqP2[[4]],a3^2] - EqP2[[4]]Coefficient[EqP2[[k]],a3^2],{k,1,9}]]
MatrixRank[Table[Coefficient[EqP2[[k]],{a1^2,a1 a2, a1 a3, a2^2, a2 a3, a3^2}[[l]]],{k,1,9},{l,1,6}]]

We will now look for other S.

We know from the intersection data that S2 should be cut out by a (1,8) divisor which vanishes, among other things, at 3 x the other line. 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. We would also like to understand the conics in the singular fiber.

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[{0== Equation18/.Line2/.t->0, 0==  Equation18/.Line2/.t->1}][[1]]]

eqs={a1==1,0== Equation18/.pt2, 0== Equation18/.ParamSsolved/.solSfiber[[2]],0== Equation18/.ParamSsolved/.solSfiber[[3]]};
For[m=2,m≤6,m++,
Print[m];
pt=RootReduce[Line2/.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 can try compute the parametric equation of S2 in multiple ways, but some of the more straightforward approaches failed. We will do the following. 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 take 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=u0^5 w1 + a2 u0^2 u1^2 w1 + a3 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;
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) /((1/596) (383+29 I (Sqrt[7]))(u0^3-u1^2)),{k,1,6}]]];
FullSimplify[((ws/.Order3Computed)/.Order3Computed)/.Order3Computed]

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

Just to be sure, we will verify that this is an automorphism (birationally).

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 *)

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 + 3B1 + 5C1 + 7A2 + 6B2 + 3C2 - 3S + S1 + 2S2 - 7F2.
We also made some adjustments and choices to have the Z_3 action extend to the seventh root of the function.

SeventhPower=  Factor[  (-51216-78896 I (Sqrt[7]))(Equation18/(w1 u0^4))^2 (Equation110/(w1 u0^5)) u0^6/(u0^3-u1^2)^2]
SeventhPowerSigma=SeventhPower/.Order3Computed;
Factor[SeventhPowerSigma/SeventhPower^2] ;
RootReduce[Factor[(SeventhPowerSigma/SeventhPower^2) /(u0 ((u0^3)-(u1^2))    w1 /((1/112) (-217 (u0^4) w1-245 I (Sqrt[7]) (u0^4) w1+77 u0 (u1^2) w1+217 I (Sqrt[7]) u0 (u1^2) w1+119 u0 u1 w2+3 I (Sqrt[7]) u0 u1 w2-455 (u0^2) w3+21 I (Sqrt[7]) (u0^2) w3+560 (u0^2) w4+140 u1 w5-124 I (Sqrt[7]) u1 w5-224 u1 w6+128 I (Sqrt[7]) u1 w6)))^7,Extension->I (Sqrt[7]) ] ]
(* this tells us how to change g *)
Order3ComputedG=Append[Order3Computed,g->g^2 (u0 ((u0^3)-(u1^2)) w1 /((1/112) (-217 (u0^4) w1-245 I (Sqrt[7]) (u0^4) w1+77 u0 (u1^2) w1+217 I (Sqrt[7]) u0 (u1^2) w1+119 u0 u1 w2+3 I (Sqrt[7]) u0 u1 w2-455 (u0^2) w3+21 I (Sqrt[7]) (u0^2) w3+560 (u0^2) w4+140 u1 w5-124 I (Sqrt[7]) u1 w5-224 u1 w6+128 I (Sqrt[7]) u1 w6)))];
Factor[((g^7/SeventhPower)/.Order3ComputedG)/ (g^14/SeventhPower^2),Extension->I (Sqrt[7])]

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. 

Eigenspace of wt 0:   1.

Eigenspace of wt -3: We are looking for g^(3) f where 
div(g^7) =  -2B - C + 3B1 + 5C1 + 7A2 + 6B2 + 3C2 - 3S + S1 + 2S2 - 7F2. 
Then we have
div(g^(-3) f) =  3/7(-2B - C + 3B1 + 5C1 + 7A2 + 6B2 + 3C2 - 3S + S1 + 2S2 - 7F2)+ div(f) >=  
             -2F - 4/7 B -2/7 C  - 4/7 B1 - 2/7 C1 - 4/7 B2 - 2/7 C2 - 6/7 S - 6/7 S1 - 6/7 S2. 
div(f) >= -2F - 4/7 B -2/7 C  - 4/7 B1 - 2/7 C1 - 4/7 B2 - 2/7 C2 - 6/7 S - 6/7 S1 - 6/7 S2 +  3/7(-2B - C + 3B1 + 5C1 + 7A2 + 6B2 + 3C2 - 3S + S1 + 2S2 - 7F2) =   -10/7B  -5/7 C + 5/7B1 + 13/7C1  + 3A2 +2B2 + C2  - 15/7S  - 3/7 S1  -4F2 
div(f) >= -B +B1+2C1+3A2+2B2+C2-2S-4F2.
The -B is a bit of a problem, we cure it by looking at f1 = f (u1^2-u0^3)/u0^3. Then 
div f1 = div f + A+...+C2 - 3F2 >= A + C + A1 + 2B1 + 3C1 + 4A2 + 3B2 + 2C2 -2S-7F2, so we are looking for  s/w1^2u0^7 where s is a section of (2,14) with vanishing at A,  A1, 4A2 (the other vanishing follow from it).

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

(* Mow let us build vanishing conditions.
First it must vanish at the two conics. *)
Eqs={};

(* rebuild the special points of the fiber and points of S in it *)
solSfiber=RootReduce[Solve[(u0^3==u1^2)/.ParamSsolved]];
pt2=Table[uws[[k]]->1,{k,1,8}] ;
pt1={u0->1,u1->1,w1->1,w2->0,w3->0,w4->0,w5->0,w6->0};
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}]];
];


SpanP2=Join[{u0->1,u1->1},Table[ws[[m]]->RootReduce[a1(ws/.pt1 ) + a2 (ptsS[[5]]+ptsS[[6]])
+ a3 ((ptsS[[5]]-ptsS[[6]])) /(NumericalApproximation[Root, (RowDefault["-0.523", FEExpressionsNumericalApproximationElider])+(RowDefault["0.375", FEExpressionsNumericalApproximationElider]) I, Root[11-17 ((#1)^2)+64 ((#1)^4)&,2], ])][[m]],{m,1,6}] ];
EqP2=Factor[(QuadricsW/.{u0->1,u1->1})/.SpanP2,Extension->I (Sqrt[7])];
tbz=Expand[(UnivPoly/.SpanP2) ( Coefficient[EqP2[[4]]/.SpanP2,a3^2] )
- (EqP2[[4]]/.SpanP2) Coefficient[UnivPoly/.SpanP2,a3^2]];
AppendTo[Eqs,Coefficient[tbz,a1^2]];
AppendTo[Eqs,Coefficient[tbz,a2^2]];
AppendTo[Eqs,Coefficient[tbz,a1 a2]];
AppendTo[Eqs,Coefficient[tbz,a1 a3]];
AppendTo[Eqs,Coefficient[tbz,a2 a3]];
(* MatrixRank[Table[Coefficient[EqP2[[k]],{a1^2,a1 a2, a1 a3, a2^2, a2 a3, a3^2}[[l]]],{k,1,9},{l,1,6}]] *)

SpanP2=Join[{u0->1,u1->1},Table[ws[[m]]->RootReduce[a1(ws/.pt2) + a2 (ptsS[[2]]+ptsS[[3]])
+ a3 ((ptsS[[2]]-ptsS[[3]])) /(NumericalApproximation[Root, (RowDefault["-0.287", FEExpressionsNumericalApproximationElider])+(RowDefault["1.15", FEExpressionsNumericalApproximationElider]) I, Root[4+5 ((#1)^2)+2 ((#1)^4)&,2], ])][[m]],{m,1,6}] ];
EqP2=Factor[(QuadricsW/.{u0->1,u1->1})/.SpanP2,Extension->I (Sqrt[7])];
tbz=Expand[(UnivPoly/.SpanP2) ( Coefficient[EqP2[[4]]/.SpanP2,a3^2] )
- (EqP2[[4]]/.SpanP2) Coefficient[UnivPoly/.SpanP2,a3^2]];
AppendTo[Eqs,Coefficient[tbz,a1^2]];
AppendTo[Eqs,Coefficient[tbz,a2^2]];
AppendTo[Eqs,Coefficient[tbz,a1 a2]];
AppendTo[Eqs,Coefficient[tbz,a1 a3]];
AppendTo[Eqs,Coefficient[tbz,a2 a3]];
(* MatrixRank[Table[Coefficient[EqP2[[k]],{a1^2,a1 a2, a1 a3, a2^2, a2 a3, a3^2}[[l]]],{k,1,9},{l,1,6}]] *)

sols1=Chop[NSolve[Table[0== Eqs[[k]],{k,1,10}],2000][[1]]];
UnivPoly=UnivPoly/.sols1;

(* now we will deal with the order four vanishing along Line2 *)
Eqs={};
coeffs=Table[ToExpression[StringJoin["r",ToString[k]]],{k,1,50}
];
For[m=3,m≤9,m++,
Print[m];
pt=Line2/.t->m;
curve=Join[{u0->1,w1->1, 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}]];
 AppendTo[ Eqs,SeriesCoefficient[UnivPoly/.curve,{eps,0,3}]];
];
sls=Solve[Table[0==Eqs[[k]],{k,1,Length[Eqs]}]][[1]];
UnivPolySolved=UnivPoly/.sls;
Variables[UnivPolySolved]
section1=Factor[Sum[RootApproximant[Chop[Coefficient[UnivPolySolved,c1 monsuniv[[k]]]]]  monsuniv[[k]],{k,1,Length[monsuniv]}]];
section2=Factor[Sum[RootApproximant[Chop[Coefficient[UnivPolySolved,c11  monsuniv[[k]]]]]  monsuniv[[k]],{k,1,Length[monsuniv]}]];
(* We modify these sections in hopes of getting nicer relations *)
nicer1=RootReduce[Solve[{1== Coefficient[a1 section1 + a2 section2,u0^7 w1^2],
0==( a1 section1 + a2 section2)/.Join[{u0->1,u1->1},Table[ws[[k]]-> ptsS[[1]][[k]],{k,1,6}]]}]][[1]];
nicer2=RootReduce[Solve[{0== Coefficient[a1 section1 + a2 section2,u0^7 w1^2],
1==( a1 section1 + a2 section2)/.Join[{u0->1,u1->1},Table[ws[[k]]-> ptsS[[1]][[k]],{k,1,6}]]}]][[1]];
FunctionWtminus31=g^(-3)/(w1^2u0^4 (u0^3-u1^2))( a1 section1 + a2 section2)/.nicer1;
FunctionWtminus32=g^(-3)/(w1^2u0^4 (u0^3-u1^2))( a1 section1 + a2 section2)/.nicer2;

Eigenspace of wt = -1 
div(g^(-1)f) >=  -2F - 4/7 B -2/7 C  - 4/7 B1 - 2/7 C1 - 4/7 B2 - 2/7 C2 - 6/7 S - 6/7 S1 - 6/7 S2
div(f)>= -2F- 4/7 B -2/7 C  - 4/7 B1 - 2/7 C1 - 4/7 B2 - 2/7 C2 - 6/7 S - 6/7 S1 - 6/7 S2 -2/7B - 1/7C + 3/7B1 + 5/7C1 + A2 + 6/7B2 + 3/7C2 - 3/7S + 1//7S1 + 2/7S2 - F2. 
div(f) >=  -2F2 - 6/7 B -3/7 C  - 1/7 B1 + 3/7 C1+A2 +2/7 B2 +1/7 C2 - 9/7 S - 5/7 S1 - 4/7 S2 
so div(f)  >= -2F2+C1+A2+B2+C2-S; so we are looking for f1/(w1 u0^2) for an element f1 of weight (1,4) which vanishes on A2 (line2).

(* one of these is related to the section of S+4F which is zero at the line A2 through two cusps *)
EquationWtminus1= u0^2 w1 + a2 w3 + a4 w4;
FunctionWtminus1=g^(-1)(EquationWtminus1/.RootReduce[Solve[{0== EquationWtminus1/.Line2/.t->0,0== EquationWtminus1/.Line2/.t->1}]][[1]])/(w1 u0^2)

The following computation is useful to be able to compute the values easier. In the last two cases it is best not to bother with the field extension to save time.

Order3ComputedG=Factor[Order3ComputedG,Extension->I (Sqrt[7])];
Print[ByteCount[%]];
FunctionWtminus1n=Factor[FunctionWtminus1/.Order3ComputedG,Extension->I (Sqrt[7])];
Print[ByteCount[%]];
FunctionWtminus1nn=Factor[FunctionWtminus1n/.Order3ComputedG,Extension->I (Sqrt[7])];
Print[ByteCount[%]];
FunctionWtminus31n=Factor[FunctionWtminus31/.Order3ComputedG,Extension->I (Sqrt[7])];
Print[ByteCount[%]];
FunctionWtminus31nn=Factor[FunctionWtminus31n/.Order3ComputedG,Extension->I (Sqrt[7])];
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,-2,-4,-3,-6,-5,-3,-6,-5};

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=1000;
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/EqsFPPnum";
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/EqsFPPwithrr";
OpenWrite[fname];
Write[fname,EqsFPPwithrr];
Close[fname];

These equations are not optimal, but we will not try to optimize them at this point.

We will now compute some torsion line bundles on the fake projective plane. We use Magma to find some nonreduced C_3 invariant hyperplane cuts and points on them. We will gradually impose conditions of being nodal at these points to get to higher powers of the prime.

At each step we will look at two tangent vectors at each point, defined over prime^depth, and an equation which vanishes on these points and tangent vectors.

prime=29;
rootvalue=rr->14;
cut1=P0+5(P1+P2+P3)+15(P4+P5+P6)+(P7+P8+P9);
cut2=P0+22(P1+P2+P3)+22(P4+P5+P6)+4(P7+P8+P9);
cut3=P0+24(P1+P2+P3)+0(P4+P5+P6)+28(P7+P8+P9);
cuts={cut1,cut2, cut3};

PointsMod29={
{{P0+P9,P1,P2+10 P9,P3+7 P9,P4+20 P9,P5+6 P9,P6+4 P9,P7+20 P9,P8+25 P9},{P0+28 P9,P1+28 P9,P2+18 P9,P3+19 P9,P4+14 P9,P5+24 P9,P6+28 P9,P7+16 P9,P8+5 P9},{P0+21 P9,P1+21 P9,P2+5 P9,P3+19 P9,P4+22 P9,P5+5 P9,P6+24 P9,P7+11 P9,P8+23 P9},{P0+11 P9,P1+4 P9,P2+18 P9,P3+3 P9,P4+10 P9,P5+17 P9,P6+16 P9,P7+6 P9,P8+26 P9},{P0+6 P9,P1+18 P9,P2+22 P9,P3+6 P9,P4+18 P9,P5+11 P9,P6+27 P9,P7+19 P9,P8+8 P9},{P0+16 P9,P1+22 P9,P2+23 P9,P3+P9,P4+9 P9,P5+6 P9,P6+P9,P7+22 P9,P8+15 P9},{P0+26 P8,P1+11 P8,P2+13 P8,P3+5 P8,P4+4 P8,P5+3 P8,P6+7 P8,P7+26 P8,P9},{P0+22 P9,P1+9 P9,P2,P3+17 P9,P4+P9,P5+5 P9,P6+16 P9,P7+7 P9,P8+5 P9},{P0+26 P9,P1+5 P9,P2+11 P9,P3+13 P9,P4+7 P9,P5+4 P9,P6+3 P9,P7,P8+26 P9}},
{{P0+14 P9,P1+12 P9,P2+P9,P3+5 P9,P4+6 P9,P5+24 P9,P6+14 P9,P7+4 P9,P8+15 P9},{P0+14 P9,P1+11 P9,P2+17 P9,P3+4 P9,P4+17 P9,P5+11 P9,P6+22 P9,P7+10 P9,P8+15 P9},{P0+3 P9,P1+21 P9,P2+4 P9,P3+19 P9,P4+11 P9,P5+16 P9,P6+13 P9,P7+16 P9,P8+8 P9},{P0+25 P9,P1+26 P9,P2+19 P9,P3,P4+13 P9,P5+6 P9,P6+18 P9,P7+17 P9,P8+27 P9},{P0+6 P9,P1+19 P9,P2+6 P9,P3,P4+24 P9,P5+P9,P6+28 P9,P7+13 P9,P8+7 P9}},
{{P0+16 P9,P1,P2+27 P9,P3+15 P9,P4+9 P9,P5+18 P9,P6+17 P9,P7+22 P9,P8+17 P9},{P0+19 P9,P1+19 P9,P2+9 P9,P3+17 P9,P4+17 P9,P5+23 P9,P6+13 P9,P7+P9,P8+26 P9},{P0+10 P9,P1+20 P9,P2+19 P9,P3+6 P9,P4+2 P9,P5+28 P9,P6+26 P9,P7+2 P9,P8+16 P9},{P0+10 P9,P1+20 P9,P2+12 P9,P3+10 P9,P4+6 P9,P5+16 P9,P6+12 P9,P7+3 P9,P8+P9},{P0+22 P9,P1+8 P9,P2+26 P9,P3+20 P9,P4+16 P9,P5+25 P9,P6+5 P9,P7+27 P9,P8+16 P9},{P0+21 P9,P1+5 P9,P2+26 P9,P3+25 P9,P4+9 P9,P5+10 P9,P6+3 P9,P7+6 P9,P8+26 P9},{P0+11 P9,P1+4 P9,P2+18 P9,P3+3 P9,P4+10 P9,P5+17 P9,P6+16 P9,P7+6 P9,P8+26 P9},{P0+P9,P1+3 P9,P2+28 P9,P3+11 P9,P4+27 P9,P5+21 P9,P6+8 P9,P7+26 P9,P8+27 P9}}};

We will now compute non-reduced cuts up to powers of prime.

Pns={Pn0,Pn1,Pn2,Pn3,Pn4,Pn5,Pn6,Pn7,Pn8,Pn9};


Coeffs=Table[ToExpression[StringJoin["c",ToString[m]]],{m,1,2000}];
maxdepth=40;
slsroot=Solve[x^2+7==0,Modulus->prime^maxdepth];
Mod[x/.slsroot,29]

GoodCuts={};
For[numcut=1,numcut≤3,numcut++,

Print["numcut=",numcut];

(* building pairs of tangent vectors at each point *)
pointscut=PointsMod29[[numcut]];
ptssub={};
tanvecs={};
For[k=1,k≤Length[pointscut],k++,
pt=Solve[Join[Table[pointscut[[k]][[l]]==0,{l,1,9}],{P0==1}],Modulus->prime][[1]];
AppendTo[ptssub,pt];

tanvector=Table[Ps[[k]]->( Ps[[k]]/.pt)+ t Pns[[k]] ,{k,1,10}];
Eqs=Join[Coefficient[((EqsFPPwithrr/.rootvalue)/.P0->1)/.tanvector,t],{Pn0}];
(* assuming that these are possible parameters in the tangent space *)
sls1=Solve[Join[Table[0==Eqs[[m]],{m,1,Length[Eqs]}],{Pn1==0,Pn2==1}],Modulus->prime][[1]];
sls2=Solve[Join[Table[0==Eqs[[m]],{m,1,Length[Eqs]}],{Pn1==1,Pn2==0}],Modulus->prime][[1]];
AppendTo[tanvecs,{tanvector/.sls1,tanvector/.sls2}];
];

currentcut=cuts[[numcut]];
currentptssub=ptssub;
currenttanvecs=tanvecs;

For[depth=1,depth≤30,depth++,
If[Mod[depth,10]==5,Print[depth]];
rootvalue=rr-> Mod[x/.slsroot[[1]],prime^(depth+1)];
EqsFPPnewdepth=PolynomialMod[(EqsFPPwithrr/.P0->1)/.rootvalue,prime^(depth+1)];

count=1;
(* make extensions to p^(depth+1) *)
newcut=currentcut + prime^depth (a1 (P1+P2+P3)+a2 (P4+P5+P6)+a3(P7+P8+P9));
newptssub={};
newtanvecs={};
For[k=1,k≤Length[currentptssub],k++,
newpt=(Ps/.currentptssub[[k]]) + prime^depth  
{0,0,Coeffs[[count]] , Coeffs[[count+1]] ,  Coeffs[[count+2]] ,  Coeffs[[count+3]] ,
 Coeffs[[count+4]] ,  Coeffs[[count+5]] ,  Coeffs[[count+6]] ,  Coeffs[[count+7]]  };
count+=8;
AppendTo[newptssub,Table[Ps[[j]]->newpt[[j]],{j,1,10}]];
newvec1=newpt + Coefficient[Ps/.currenttanvecs[[k]][[1]],t]t +  prime^depth  t 
{0,0,0,Coeffs[[count]], Coeffs[[count+1]] ,  Coeffs[[count+2]] ,  Coeffs[[count+3]] ,
 Coeffs[[count+4]] ,  Coeffs[[count+5]] ,  Coeffs[[count+6]] };
count+=7;
newvec2=newpt + Coefficient[Ps/.currenttanvecs[[k]][[2]],t]t +  prime^depth  t 
{0,0,0,Coeffs[[count]], Coeffs[[count+1]] ,  Coeffs[[count+2]] ,  Coeffs[[count+3]] ,
 Coeffs[[count+4]] ,  Coeffs[[count+5]] ,  Coeffs[[count+6]] };
count+=7;
AppendTo[newtanvecs,{Table[Ps[[j]]->newvec1[[j]],{j,1,10}],Table[Ps[[j]]->newvec2[[j]],{j,1,10}]}];
];

(* make equations *)
Eqs={};

For[l=1,l≤Length[tanvecs],l++,
For[k=1,k≤84,k++,
tbz=(EqsFPPnewdepth[[k]]/.newtanvecs[[l]][[1]])/prime^depth;
AppendTo[Eqs,PolynomialMod[Factor[tbz/.t->0],prime]];
AppendTo[Eqs,PolynomialMod[Factor[Coefficient[tbz,t]],prime]];
tbz=(EqsFPPnewdepth[[k]]/.newtanvecs[[l]][[2]])/prime^depth;
AppendTo[Eqs,PolynomialMod[Factor[Coefficient[tbz,t]],prime]];
];
tbz=(newcut/.newtanvecs[[l]][[1]])/prime^depth;
AppendTo[Eqs,PolynomialMod[Factor[tbz/.t->0],prime]];
AppendTo[Eqs,PolynomialMod[Factor[Coefficient[tbz,t]],prime]];
];

sls=Solve[Table[0==Eqs[[k]],{k,1,Length[Eqs]}],Modulus->prime][[1]];

currentcut=newcut/.sls;
currentptssub=newptssub/.sls;
currenttanvecs=newtanvecs/.sls;

];
AppendTo[GoodCuts,currentcut];
];

We now compute minimal polynomials for the coefficients of the cuts. It looks like the second cut is defined over the quadratic field.

For[n=1,n≤3,n++,
Print[n];
r7=Mod[rr/.rootvalue,29^31];

whatis=Coefficient[GoodCuts[[n]],P1];
thebasis=Append[Join[Table[Mod[whatis^k,29^31],{k,0,2}],Table[Mod[whatis^k r7,29^31],{k,0,2}]],29^31];
smalllincomb=LatticeReduce[NullSpace[{thebasis}]][[1]];
smallpoly1=Sum[smalllincomb[[k]] t^(k-1),{k,1,3}]+I (Sqrt[7])Sum[smalllincomb[[k]] t^(k-4),{k,4,6}];


whatis=Coefficient[GoodCuts[[n]],P4];
thebasis=Append[Join[Table[Mod[whatis^k,29^31],{k,0,2}],Table[Mod[whatis^k r7,29^31],{k,0,2}]],29^31];
smalllincomb=LatticeReduce[NullSpace[{thebasis}]][[1]];
smallpoly2=Sum[smalllincomb[[k]] t^(k-1),{k,1,3}]+I (Sqrt[7])Sum[smalllincomb[[k]] t^(k-4),{k,4,6}];


whatis=Coefficient[GoodCuts[[n]],P7];
thebasis=Append[Join[Table[Mod[whatis^k,29^31],{k,0,2}],Table[Mod[whatis^k r7,29^31],{k,0,2}]],29^31];
smalllincomb=LatticeReduce[NullSpace[{thebasis}]][[1]];
smallpoly3=Sum[smalllincomb[[k]] t^(k-1),{k,1,3}]+I (Sqrt[7])Sum[smalllincomb[[k]] t^(k-4),{k,4,6}];
Print[{smallpoly1,smallpoly2,smallpoly3}];
];

We will study the second cut in more detail.

r7=Mod[rr/.rootvalue,29^31];
whatis=Coefficient[GoodCuts[[2]],P1];
thebasis=Append[Join[Table[Mod[whatis^k,29^31],{k,0,1}],Table[Mod[whatis^k r7,29^31],{k,0,1}]],29^31];
smalllincomb=LatticeReduce[NullSpace[{thebasis}]][[1]];
smallpoly1=Sum[smalllincomb[[k]] t^(k-1),{k,1,2}]+I (Sqrt[7])Sum[smalllincomb[[k]] t^(k-3),{k,3,4}];
a1=t/.Solve[smallpoly1==0][[1]];

whatis=Coefficient[GoodCuts[[2]],P4];
thebasis=Append[Join[Table[Mod[whatis^k,29^31],{k,0,1}],Table[Mod[whatis^k r7,29^31],{k,0,1}]],29^31];
smalllincomb=LatticeReduce[NullSpace[{thebasis}]][[1]];
smallpoly2=Sum[smalllincomb[[k]] t^(k-1),{k,1,2}]+I (Sqrt[7])Sum[smalllincomb[[k]] t^(k-3),{k,3,4}];
a2=t/.Solve[smallpoly2==0][[1]];

whatis=Coefficient[GoodCuts[[2]],P7];
thebasis=Append[Join[Table[Mod[whatis^k,29^31],{k,0,1}],Table[Mod[whatis^k r7,29^31],{k,0,1}]],29^31];
smalllincomb=LatticeReduce[NullSpace[{thebasis}]][[1]];
smallpoly3=Sum[smalllincomb[[k]] t^(k-1),{k,1,2}]+I (Sqrt[7])Sum[smalllincomb[[k]] t^(k-3),{k,3,4}];
a3=t/.Solve[smallpoly3==0][[1]];

SecondCut= RootReduce[P0+a1(P1+P2+P3)+a2(P4+P5+P6)+ a3(P7+P8+P9)]

We also want to be able to get algebraic solutions for the first and third cuts. This involves computing a2 and a3 in terms of a1.

r7=Mod[rr/.rootvalue,29^31];
whatis=Coefficient[GoodCuts[[1]],P1];
thebasis=Append[Join[Table[Mod[whatis^k,29^31],{k,0,2}],Table[Mod[whatis^k r7,29^31],{k,0,2}]],29^31];
smalllincomb=LatticeReduce[NullSpace[{thebasis}]][[1]];
smallpoly1=Sum[smalllincomb[[k]] t^(k-1),{k,1,3}]+I (Sqrt[7])Sum[smalllincomb[[k]] t^(k-4),{k,4,6}];
a1=RootReduce[t/.Solve[smallpoly1==0][[1]]];

whatis=Coefficient[GoodCuts[[1]],P4];
thebasis=Append[Mod[{1,whatis, Coefficient[GoodCuts[[1]],P1],whatis  Coefficient[GoodCuts[[1]],P1], r7,  whatis r7 },29^31 ],29^31];
smalllincomb=LatticeReduce[NullSpace[{thebasis}]][[1]];
smallpoly2=smalllincomb.{1,t, a1,t a1, I (Sqrt[7]), t I (Sqrt[7]),0 };
a2=RootReduce[t/.Solve[smallpoly2==0][[1]]];

whatis=Coefficient[GoodCuts[[1]],P7];
thebasis=Append[Mod[{1,whatis, Coefficient[GoodCuts[[1]],P1],whatis  Coefficient[GoodCuts[[1]],P1], r7,  whatis  r7 },29^31 ],29^31];
smalllincomb=LatticeReduce[NullSpace[{thebasis}]][[1]];
smallpoly3=smalllincomb.{1,t, a1,t a1,  I (Sqrt[7]), t  I (Sqrt[7]),0 };
a3=RootReduce[t/.Solve[smallpoly3==0][[1]]];

FirstCut=P0+a1(P1+P2+P3)+a2(P4+P5+P6)+ a3(P7+P8+P9);

whatis=Coefficient[GoodCuts[[3]],P1];
thebasis=Append[Join[Table[Mod[whatis^k,29^31],{k,0,2}],Table[Mod[whatis^k  r7,29^31],{k,0,2}]],29^31];
smalllincomb=LatticeReduce[NullSpace[{thebasis}]][[1]];
smallpoly1=Sum[smalllincomb[[k]] t^(k-1),{k,1,3}]+I (Sqrt[7])Sum[smalllincomb[[k]] t^(k-4),{k,4,6}];
a1=RootReduce[t/.Solve[smallpoly1==0][[1]]];

whatis=Coefficient[GoodCuts[[3]],P4];
thebasis=Append[Mod[{1,whatis, Coefficient[GoodCuts[[3]],P1],whatis  Coefficient[GoodCuts[[3]],P1], r7,  whatis r7 },29^31 ],29^31];
smalllincomb=LatticeReduce[NullSpace[{thebasis}]][[1]];
smallpoly2=smalllincomb.{1,t, a1,t a1, I (Sqrt[7]), t I (Sqrt[7]),0 };
a2=RootReduce[t/.Solve[smallpoly2==0][[1]]];

whatis=Coefficient[GoodCuts[[3]],P7];
thebasis=Append[Mod[{1,whatis, Coefficient[GoodCuts[[3]],P1],whatis  Coefficient[GoodCuts[[3]],P1], r7,  whatis r7 },29^31 ],29^31];
smalllincomb=LatticeReduce[NullSpace[{thebasis}]][[1]];
smallpoly3=smalllincomb.{1,t, a1,t a1, I (Sqrt[7]), t I (Sqrt[7]),0 };
a3=RootReduce[t/.Solve[smallpoly3==0][[1]]];

ThirdCut=P0+a1(P1+P2+P3)+a2(P4+P5+P6)+ a3(P7+P8+P9);

We now verify that these are non-reduced. In lieu of formal verification, we will do it numerically.

sols=NSolve[Join[Table[0==EqsFPP[[k]],{k,1,84}],{0==FirstCut,P0==1,P1==RandomInteger[10]+1}],50];
Length[DeleteDuplicates[sols]]

sols=NSolve[Join[Table[0==EqsFPP[[k]],{k,1,84}],{0==SecondCut,P0==1,P1==RandomInteger[10]+1}],50];
Length[DeleteDuplicates[sols]]

sols=NSolve[Join[Table[0==EqsFPP[[k]],{k,1,84}],{0==ThirdCut,P0==1,P1==RandomInteger[10]+1}],50];
Length[DeleteDuplicates[sols]]

Finally, we will use these non-reduced cuts to try to get a nicer set of coordinates Q0,...,Q9 out of P0,...,P9. The idea is to have the second cut to be nice. We also try to have somewhat smaller coefficients in some of the equations.

Ps={P0,P1,P2,P3,P4,P5,P6,P7,P8,P9};
Qs={Q0,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9};
PtoQ={P0->Q0,P1->c1 Q1,P2->c1 Q2, P3-> c1 Q3, P4-> c2 Q4, P5 ->  c2 Q5, P6 ->  c2 Q6, P7-> c3 Q4 + c4 Q7, P8->c3 Q5 +c4 Q8, P9-> c3 Q6 + c4 Q9};
(* the choice of c2 is arbitrary at first and is fixed later *)
solcs=RootReduce[Solve[{c2==((-415-387 I (Sqrt[7]))/(2 152576)),Coefficient[SecondCut/.PtoQ,Q1]==1,Coefficient[SecondCut/.PtoQ,Q4]==0,Coefficient[SecondCut/.PtoQ,Q9]==1}]][[1]]
PtoQ = PtoQ/.solcs;
QtoP =Solve[Table[Ps[[k]]==(Ps[[k]]/.PtoQ),{k,1,10}],Qs][[1]];

(* change P-points *)
Qpoints={};
For[k=1,k≤Length[Ppoints],k++,
AppendTo[Qpoints,Table[Qs[[m]]->(( Qs[[m]]/.QtoP)/.Ppoints[[k]]),{m,1,10}]];
];
Accuracy[Ppoints]
Accuracy[Qpoints]

(* compute new equations *)
deg=3;
EqsFPPq={};

Unset[{a1,a2,a3,a4,a5,a6,a7,a8,a9,a10}];
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[Qs[[n]]^(powers[[n]]/.sols[[k]]),{n,1,10}],{k,1,Length[sols]}];
TableOfValues=Table[Monomials[[k]]/.Qpoints[[l]],{l,1,Length[Ppoints]},{k,1,Length[Monomials]}];
ns=RootApproximant[Chop[RowReduce[NullSpace[TableOfValues]],10^(-900)]];
EqsFPPq=Join[EqsFPPq,Table[Sum[Monomials[[k]]ns[[i]][[k]],{k,1,Length[Monomials]}],{i,1,Length[ns]}]];
];
Length[EqsFPPq]
Table[ByteCount[EqsFPPq[[k]]],{k,1,Length[EqsFPPq]}]
EqsFPPqn=Numerator[Factor[EqsFPPq]];
Table[ByteCount[EqsFPPqn[[k]]],{k,1,Length[EqsFPPqn]}]

We will now write the equations to a file, after replacing root(-7) with rr.

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