Us = {U0, U1, U2, U3, U4, U5, U6, U7, U8, U9};
Mod7Weights={0,6,5,3,1 ,2,4,1,2,4}; 

(*Order 3 automorphism*)
Order3Aut={U1->U2,U2->U3,U3->U1,U4->U5,U5->U6,U6->U4,U7->U8,U8->U9,U9->U7};

(*Import file of equations for FPP*)
SetDirectory[NotebookDirectory[]];
FileNames[]

filename = "Dependency/EqsFPP.txt";
EqsFPP=ToExpression[Import[filename]]/.rr->(Sqrt[-7]);

(*
fname="Dependency/PointsOnFPP1.txt";
Upoints=ToExpression[Import[filename]];
*)

Calculating s3 and d.
This part of the code is adapted from the code Map_to_P5.nb in the paper “Realizing a Fake Projective Plane as a Degree 25 Surface in P^5” by Borisov and Lihn.

(*p1, p2, p3 are the 3 fixed points of the C7 action *)
p1=Table[Us[[i]]->KroneckerDelta[i-8],{i,1,10}];
p2=Table[Us[[i]]->KroneckerDelta[i-9],{i,1,10}];
p3=Table[Us[[i]]->KroneckerDelta[i-10],{i,1,10}];

(* We will only be using p2 and p3, they are points in the curve {s3 = 0} but note that it will turn out that p1 is also a point in {s3 = 0} *)
(*First, we will compute some local charts on the FPP to speed up the process. We only care about local behavior so we can take a subset of the polynomials that define the point*)

(* Let us compute order three neighborhood of p2 *)
JacobianFPP=Factor[Table[D[EqsFPP[[i]],Us[[k]]],{i,1,Length[EqsFPP]},{k,1,Length[Us]}]];

EqListp2={};
rows={};
rk=0;
Do[
newRow=JacobianFPP[[i]]/.p2;
If[MatrixRank[Append[rows,newRow]]>rk,AppendTo[EqListp2,EqsFPP[[i]]]; rk++;AppendTo[rows,newRow]];
,{i,Length[EqsFPP]}];

(*Now, let's compute the order 1 neighborhood *)
ListCoeffs=Table[ToExpression[StringJoin["a",ToString[i]]],{i,1,1000}];
count=1;
ParamTanSpaceOrd1={};
Do[
AppendTo[ParamTanSpaceOrd1, Which[i==9,U8->1,i==7,U6->v,i==5,U4->u,True,	Us[[i]]->ListCoeffs[[count++]]u+ListCoeffs[[count++]]v]];
,{i,Length[Us]}];

tbz={};
Do[
Do[
AppendTo[tbz,SeriesCoefficient[EqListp2/.ParamTanSpaceOrd1,{u,0,i},{v,0,j}]];
,{j,0,1-i}];
,{i,0,1}];
FirstOrderCoeffs=FullSimplify[Solve[Table[tbz[[i]]==0,{i,1,Length[tbz]}]]];
Ord1Eqsp2=Simplify[ParamTanSpaceOrd1/.FirstOrderCoeffs[[1]]]


(*Now let's compute the order 2 neighborhood *)
ParamTanSpaceOrd2={};
Do[
AppendTo[ParamTanSpaceOrd2, Which[i==9,U8->1,i==7,U6->v,i==5,U4->u,True,	Us[[i]]->(Us[[i]]/.Ord1Eqsp2[[i]])+ListCoeffs[[count++]]u^2+ListCoeffs[[count++]] u v + ListCoeffs[[count++]]v^2]];
,{i,Length[Us]}];
tbz={};
Do[
Do[
AppendTo[tbz,SeriesCoefficient[EqListp2/.ParamTanSpaceOrd2,{u,0,i},{v,0,j}]];
,{j,0,2-i}];
,{i,0,2}];
SecondOrderCoeffs=FullSimplify[Solve[Table[tbz[[i]]==0,{i,1,Length[tbz]}]]];
Ord2Eqsp2=ParamTanSpaceOrd2/.SecondOrderCoeffs[[1]];
Ord2Eqsp2 = Simplify[Ord2Eqsp2]

(*Now let's compute the order 3 neighborhood *)
ParamTanSpaceOrd3={};
Do[
AppendTo[ParamTanSpaceOrd3, Which[i==9,U8->1,i==7,U6->v,i==5,U4->u,True,	Us[[i]]->(Us[[i]]/.Ord2Eqsp2[[i]])+ListCoeffs[[count++]]u^3+ListCoeffs[[count++]] u^2 v +ListCoeffs[[count++]] v^2 u + ListCoeffs[[count++]]v^3]];
,{i,Length[Us]}];
tbz={};
Do[
Do[
AppendTo[tbz,SeriesCoefficient[EqListp2/.ParamTanSpaceOrd3,{u,0,i},{v,0,j}]];
,{j,0,3-i}];
,{i,0,3}];

ThirdOrderCoeffs=FullSimplify[Solve[Table[tbz[[i]]==0,{i,1,Length[tbz]}]]];
Ord3Eqsp2=ParamTanSpaceOrd3/.ThirdOrderCoeffs[[1]];
Ord3Eqsp2 = Simplify[Ord3Eqsp2]

Print["Order 3 neighborhood of p2 is computed"];

(* Let us compute order three neighborhood of p3 *)

(*First, we will compute some local FPPs to speed up the process. We only care about local behavior so we can take a subset of the polynomials that define the point*)
EqListp3={};
rows={};
rk=0;
Do[
newRow=JacobianFPP[[i]]/.p3;
If[MatrixRank[Append[rows,newRow]]>rk,AppendTo[EqListp3,EqsFPP[[i]]]; rk++;AppendTo[rows,newRow]];
,{i,Length[EqsFPP]}];

(*Now let's compute the order 1 neighborhood *)
ListCoeffs=Table[ToExpression[StringJoin["a",ToString[i]]],{i,1,1000}];
count=1;
ParamTanSpaceOrd1={};
Do[
AppendTo[ParamTanSpaceOrd1, Which[i==10, U9->1,i==6, U5->v,i==5,U4->u,True,Us[[i]]->ListCoeffs[[count++]]u+ListCoeffs[[count++]]v]];
,{i,Length[Us]}];

tbz={};
Do[
Do[
AppendTo[tbz,SeriesCoefficient[EqListp3/.ParamTanSpaceOrd1,{u,0,i},{v,0,j}]];
,{j,0,1-i}];
,{i,0,1}];
FirstOrderCoeffs=FullSimplify[Solve[Table[tbz[[i]]==0,{i,1,Length[tbz]}]]];
Ord1Eqsp3=Simplify[ParamTanSpaceOrd1/.FirstOrderCoeffs[[1]]]


(*Now let's compute the order 2 neighborhood *)
ParamTanSpaceOrd2={};
Do[
AppendTo[ParamTanSpaceOrd2, Which[i==10, U9->1,i==6, U5->v,i==5,U4->u,True,	Us[[i]]->(Us[[i]]/.Ord1Eqsp3[[i]])+ListCoeffs[[count++]]u^2+ListCoeffs[[count++]] u v + ListCoeffs[[count++]]v^2]];
,{i,Length[Us]}];
tbz={};
Do[
Do[
AppendTo[tbz,SeriesCoefficient[EqListp3/.ParamTanSpaceOrd2,{u,0,i},{v,0,j}]];
,{j,0,2-i}];
,{i,0,2}];
SecondOrderCoeffs=FullSimplify[Solve[Table[tbz[[i]]==0,{i,1,Length[tbz]}]]];
Ord2Eqsp3=ParamTanSpaceOrd2/.SecondOrderCoeffs[[1]];
Ord2Eqsp3 = Simplify[Ord2Eqsp3]

(*Now let's compute the order 3 neighborhood *)
ParamTanSpaceOrd3={};
Do[
AppendTo[ParamTanSpaceOrd3, Which[i==10, U9->1,i==6, U5->v,i==5,U4->u,True,Us[[i]]->(Us[[i]]/.Ord2Eqsp3[[i]])+ListCoeffs[[count++]]u^3+ListCoeffs[[count++]] u^2 v +ListCoeffs[[count++]] v^2 u + ListCoeffs[[count++]]v^3]];
,{i,Length[Us]}];

(*Print[ParamTanSpaceOrd3];*)

tbz={};
Do[
Do[
AppendTo[tbz,SeriesCoefficient[EqListp3/.ParamTanSpaceOrd3,{u,0,i},{v,0,j}]];
,{j,0,3-i}];
,{i,0,3}];

ThirdOrderCoeffs=FullSimplify[Solve[Table[tbz[[i]]==0,{i,1,Length[tbz]}]]];
Ord3Eqsp3=ParamTanSpaceOrd3/.ThirdOrderCoeffs[[1]];
Ord3Eqsp3 = Simplify[Ord3Eqsp3]

Print["Order 3 neighborhood of p3 is computed"];

(* sanity check: the equations from EqsFPP should vanish on the neighborhoods to third order *)

ShouldBeZero={};
Do[
Do[
AppendTo[ShouldBeZero,SeriesCoefficient[EqsFPP/.Ord3Eqsp2, {u,0,i},{v,0,j}]];
,{j,0,3-i}];
,{i,0,3}];

Do[
Do[
AppendTo[ShouldBeZero,SeriesCoefficient[EqsFPP/.Ord3Eqsp3, {u,0,i},{v,0,j}]];
,{j,0,3-i}];
,{i,0,3}];
If[DeleteDuplicates[Flatten[ShouldBeZero]]=={0},Print["Sanity check successful."]];

Now we construct s3, s5 and s6. Vanishing of s3 to third order at p2 and p3 will put some linear constraints on the coefficients which we solve for.

deg=2;
Powers={d1,d2,d3,d4,d5,d6,d7,d8,d9,d10};
(* All powers of the variables that sum up to 2 *)
AllPowers=
Solve[Join[Table[Powers[[k]]≥0,{k,1,10}],{Sum[Powers[[k]],{k,1,10}]==deg}],Integers];

(*s3 = r_3^3 has weight 2 = 3 x 3 mod 7*)
mod7=2;
CorrectPowers={};
Do[
If[Mod[Sum[Mod7Weights[[k]] Powers[[k]],{k,1,10}]/.AllPowers[[k]],7]==mod7,
AppendTo[CorrectPowers,AllPowers[[k]]];
];
,{k,1,Length[AllPowers]}];

(* CorrectPowers is the combination of variables in all powers that gives 2 mod 7 weight when summed up *)
(* Print[CorrectPowers];*)

(* Monomials given by the correct powers *)
Monomials=Table[Product[Us[[n]]^(Powers[[n]]/.CorrectPowers[[k]]),{n,1,10}],{k,1,Length[CorrectPowers]}];
(* Print[Monomials]; *)

(* To each monomial we give assign a coefficient b_i *)
ListCoeffs=Table[ToExpression[StringJoin["b",ToString[i]]],{i,1,Length[Monomials]}];

(* Initalize S3 as linear combination of monomials and their b_i *)
S3=Sum[ListCoeffs[[i]] Monomials[[i]],{i,1,Length[Monomials]}]


(* We apply the constraints given by the Order 2 neighborhoods of p2 and p3 in *)
tbz={};
Do[
Do[
AppendTo[tbz,SeriesCoefficient[S3/.Ord2Eqsp2,{u,0,i},{v,0,j}]];
AppendTo[tbz,SeriesCoefficient[S3/.Ord2Eqsp3,{u,0,i},{v,0,j}]];
,{j,0,2-i}];
,{i,0,2}];

CoeffsS3=FullSimplify[Solve[Table[tbz[[i]]==0,{i,1,Length[tbz]}]][[1]]];
S3=Together[S3/.CoeffsS3]

(* Compute s5 and s6 from the C3 action. *)
S5=Simplify[S3/.Order3Aut]
S6=Simplify[S5/.Order3Aut]

(* See which variables survive the linear conditions. *)
bsS3=Intersection[Variables[S3], ListCoeffs]

We now write the general form for d and translate the equation d^3 = s3 s5 s6 into equations on the unknown coefficients.

(* Since d is of weight 0 and is C_7 invariant, we have by degree and invariant conditions the following form for D3 = d^3 *)
es={e1,e2,e3};
D3=(e1 U0^2+e2(U1 U4+U2 U5+U3 U6)+e3(U1 U7+U2 U8+U3 U9))^3;
eAndbs=Join[es,bsS3];

(* we will compute 30 points on the FPP *)
Upoints={};
numpoints=30;
prec=1000;

For[k=1,k≤Round[numpoints/10],k++,
Print[k];
EqsFPPplus= Join[EqsFPP, {U0-1, Sum[Us[[l]] RandomInteger[50],{l,1,10}], Sum[Us[[l]] RandomInteger[50],{l,1,10}]}];
sols=NSolve[EqsFPPplus==0,prec];
Upoints=Join[Upoints,Table[sols[[l]],{l,1,10}]];
];

(* Tries to compute the matrix with the ceofficients e_i and b_i as the independent variables *)
Monomialseb=MonomialList[(S3 S5 S6-D3),eAndbs]/(MonomialList[(S3 S5 S6-D3),eAndbs]/.Table[eAndbs[[i]]->1,{i,1,Length[eAndbs]}]);

ListCoeffsAndPs=N[MonomialList[(S3 S5 S6-D3),eAndbs]/Monomialseb,1000];
ListCoeffsAndPsValues=Table[ListCoeffsAndPs/.Upoints[[i]],{i,1,Length[Upoints]}];

Dimensions[ListCoeffsAndPsValues]
MatrixRank[ListCoeffsAndPsValues]
Accuracy[ListCoeffsAndPsValues]

(*We will export the matrix to magma to find a solution, see the file 4H-Quadratic.txt *)

ListCoeffsAndPsReducedValues=Block[{$MaxExtraPrecision=1000},RootApproximant[Chop[RowReduce[ListCoeffsAndPsValues],10^(-800)]]];
DenominatorOfAll=Apply[LCM,Join[Table[Apply[LCM,Denominator[ListCoeffsAndPsReducedValues[[i]]]],{i,1,Length[ListCoeffsAndPsReducedValues]}]]];
ListCoeffsAndPsReducedValues=Cancel[ListCoeffsAndPsReducedValues*DenominatorOfAll];

(* Here we replace sqrt(-7) as rr *)
ToRR[x_]:=Ret=Rationalize[Chop[N[Re[x],1900],10^(-700)]]+rr Rationalize[Chop[N[Im[x]/(Sqrt[7]),1900],10^(-700)]];
ListCoeffsAndPsReducedValuesrr=Table[ToRR[ListCoeffsAndPsReducedValues[[i]]],{i,1,Length[ListCoeffsAndPsReducedValues]}];

tbzrr=Cases[ListCoeffsAndPsReducedValuesrr.Monomialseb,Except@0];

fname=OpenWrite["Equations/4H_Coeffs_Vanishing.txt"];
Write[fname,tbzrr];
Close[fname];

Magma finds 3 irreducible components defined over the number field Q(\sqrt{-7}):
---------------------------------------------
[
    Scheme over FF defined by
    b5*b8^2 + 1/2*(189*rr + 399)*e2*e3^2 + (203*rr - 511)*e3^3,
    b8^3 + (81*rr - 21)*e2*e3^2 + (18*rr - 490)*e3^3,
    b5^2 + (rr - 2)*b5*b8 + 1/4*(-4*rr - 3)*b8^2,
    b5*e2 + 1/24*(-11*rr + 17)*b5*e3 + 1/48*(63*rr + 259)*b8*e3,
    b8*e2 + 1/12*(-7*rr + 21)*b5*e3 + 1/24*(59*rr + 31)*b8*e3,
    e2^2 + (2*rr + 2)*e2*e3 + (2*rr - 6)*e3^2,
    b1,
    b2 + 1/8*(-rr - 5)*b5 + 1/16*(-3*rr + 1)*b8,
    b3 + 1/4*(-rr - 1)*b5 + 1/8*(-3*rr - 3)*b8,
    b4 + 1/4*(rr + 1)*b5 + 1/8*(-3*rr - 3)*b8,
    e1 + 1/4*(-rr - 1)*e2 + 2*e3,
    Scheme over FF defined by
    b8^2 + (-3*rr + 7)*b8*e3 + (-42*rr - 14)*e3^2,
    b1,
    b2 + 1/8*(rr - 15)*b8,
    b3 + 1/4*(rr - 13)*b8,
    b4 + 1/8*(-rr - 1)*b8,
    b5 + 1/32*(7*rr + 3)*b8,
    e1 + 1/2*(rr + 1)*e3,
    e2 + (rr + 1)*e3,
    Scheme over FF defined by
    b1,
    b2 + 1/2*(13*rr - 21)*e3,
    b3 + 1/2*(23*rr - 35)*e3,
    b4 + 1/2*(-rr - 7)*e3,
    b5 + 1/4*(5*rr + 21)*e3,
    b8 + (3*rr - 7)*e3,
    e1 + 1/2*(rr + 1)*e3,
    e2 + (rr + 1)*e3
]
---------------------------------------------
The last component has solution:
[ <0, 1, 1/58*(-rr + 101), 1/29*(2*rr + 1), 1/232*(-27*rr + 1), 1/29*(rr + 15), 1/812*(-17*rr + 35), 1/406*(-17*rr + 35), 1/812*(13*rr + 21)> ] 
which we will use.

(* Here we check the solution and write the result into the file *)
possibleSol = {b1 -> 0, b2-> 1,b3->1/58*(-rr+101), b4->1/29*(2*rr+1), b5->1/232*(-27*rr+1), b8->1/29*(rr+15),e1->1/812*(-17*rr+35),e2->1/406*(-17*rr+35),e3->1/812*(13*rr+21)};
Simplify[tbzrr/.possibleSol/.rr->(Sqrt[-7])]

NewS3 = S3/.possibleSol
NewD = e1 U0^2+e2(U1 U4+U2 U5+U3 U6)+e3(U1 U7+U2 U8+U3 U9);
NewD = Simplify[NewD/.possibleSol]
Clear[S3]
fname=OpenWrite["Equations/s3_and_D.txt"];
Write[fname, {S3->NewS3, D->NewD}];
Close[fname];

We now compute quadratic equations which are zero on r_3=0. This space splits according to characters of C7, so we will only need 20 points.

(* We already calculated previously the explicit equation of s3 = r3^3 below, where r3 is one of the global sections of H^0(FPP, 4H). We can calculate for equations on s3 and use the C3 action. *)
NewS3 =(((1-27*rr)*U1*U3)/29+(4*(101-rr)*U4^2)/29+(8*(15+rr)*U0*U5)/29+(8*(1+2*rr)*U2*U6)/29+8*U4*U7-4*U0*U8+((101-rr)*U0*U8)/29+((101-rr)*U2*U9)/58+(I/58)*Sqrt[7]*(101-rr)*U2*U9)/8;
NewS3 = NewS3/.rr->(Sqrt[-7]);

NewD = ((35-17*rr)*U0^2+2*(35-17*rr)*(U1*U4+U2*U5+U3*U6)+(21+13*rr)*(U1*U7+U2*U8+U3*U9))/812;
NewD = NewD/.rr->(Sqrt[-7]);

(* We compute 20 points on s3 = d = 0. *)

iter = 2;
PointsOnR3= {};
Do [
RandomCut =RandomSample[Join[-Range[100],Range[100]], 10].Us;
InputEquations = Join[EqsFPP, {NewS3, NewD, U2 - 1,RandomCut}];
current = NSolve[InputEquations==0, 500];
PointsOnR3 = Join[PointsOnR3, current[[1;;10]]];
Print[Length[PointsOnR3]];
, {i, 1, iter}];

(* Because accuracy can drop quite a bit here, we are making sure that it is still decent. *)
Accuracy[PointsOnR3]


(* Now computing the polynomials. *)
deg=2;

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

Equations = {};
Equationsrr = {};

For[mod7=0,mod7≤6,mod7++,
Print["mod7=",mod7];

CorrectPowers={};
Do[
If[Mod[Sum[Mod7Weights[[k]] Powers[[k]],{k,1,10}]/.AllPowers[[k]],7]==mod7,
AppendTo[CorrectPowers,AllPowers[[k]]];
];
,{k,1,Length[AllPowers]}];

Monomials=Table[Product[Us[[n]]^(Powers[[n]]/.CorrectPowers[[k]]),{n,1,10}],{k,1,Length[CorrectPowers]}];

TableOfValues=Table[Monomials[[k]]/.PointsOnR3[[l]],{l,1,Length[PointsOnR3]},{k,1,Length[Monomials]}];

(* Print[Dimensions[TableOfValues]];
Print[MatrixRank[TableOfValues]];
Print[Accuracy[TableOfValues]]; *)

result =NullSpace[Chop[RowReduce[TableOfValues,Tolerance->10^(-270)],10^(-270)]];
(* Print["nullspace done"];
Print[Dimensions[result]]; 
Print[Accuracy[result]]; *)

rowResult = RowReduce[Chop[result, 10^(-270)], Tolerance->10^(-270)];
(* Print["Rowreduce done"]; 
Print[Accuracy[rowResult]]; *)

exactResult = RootApproximant[Chop[rowResult, 10^(-270)], 2];

rowResultrr =Table[ToRR[exactResult[[i]]],{i,1,Length[exactResult]}];
newEquationsrr = rowResultrr.Monomials;
newEquations = exactResult.Monomials;
(*Print[newEquations];*)

AppendTo[Equations, newEquations];
AppendTo[Equationsrr, newEquationsrr];

];

Equations

fname=OpenWrite["Equations/4H_one_section.txt"];
Write[fname,Flatten[Equationsrr]];
Close[fname];

Length[Flatten[Equationsrr]]