Map of Keum’s FPP to P5 Q

In this notebook, we will present the code to map Keum’s FPP to P5 via |5H|

(*To simplify file I/O, we will use relative file paths. You should set the below working directory to the file folder this notebook is in, along with the relevant supporting files.*)
SetDirectory[""]
FileNames[]

(*Loading equations for Keum's FPP from Borisov21*)
Qs={Q0,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9};

fname="Dependencies\\BeautifiedFPPwithrrn.txt";
OpenRead[fname];
EqsFPPKeum=Read[fname]/.rr->(Sqrt[-7]);
Close[fname];

(*Based on the holomorphic Lefschetz fixed point formula we have three sections in H^0(FPP,4H) of weight 3, 5, and 6, which we will call R3, R5, R6 respectively.
Now let S3 = R3^3 \in H^0(FPP, 12H). Then it has degree 2 (mod 7) and hence is formed by quadratic expressions in the P_i. Our goal is to compute this
*)
(*We know that S3 passes through the two C7 fixed points defined below. Let's also compute the Jacobian of the FPP*)
p1=Table[Qs[[i]]->KroneckerDelta[i-8],{i,1,10}];
p2=Table[Qs[[i]]->KroneckerDelta[i-9],{i,1,10}];
p3=Table[Qs[[i]]->KroneckerDelta[i-10],{i,1,10}];

JacobianFPP=Factor[Table[D[EqsFPPKeum[[i]],Qs[[k]]],{i,1,Length[EqsFPPKeum]},{k,1,Length[Qs]}]];

(*Now let's get the basis for the tangent space at p1. There will be 3 since we are working over the cone. Let's express the tangent space in u-v coordinates*)
TanSpacep1=RowReduce[NullSpace[JacobianFPP/.p1]]
(*Tells us that P7=1, P5=u, P6=v*)
(*We will compute the order 2 space for p1*)

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

(*Now, let's compute the order 1 tangent space*)
ListCoeffs=Table[ToExpression[StringJoin["a",ToString[i]]],{i,1,1000}];
count=1;
ParamTanSpaceOrd1={};

Do[
AppendTo[ParamTanSpaceOrd1, Which[i==8, Q7->1,i==6, Q5->u,i==7,Q6->v,True,	Qs[[i]]->ListCoeffs[[count++]]u+ListCoeffs[[count++]]v]];
,{i,Length[Qs]}];
tbz={};
Do[
Do[
AppendTo[tbz,SeriesCoefficient[EqListp1/.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]}]]]
Ord1Eqsp1=Simplify[ParamTanSpaceOrd1/.FirstOrderCoeffs[[1]]]

(*Now let's compute the order 2 tangent space*)
ParamTanSpaceOrd2={};
Do[
AppendTo[ParamTanSpaceOrd2, Which[i==8, Q7->1,i==6, Q5->u,i==7,Q6->v,True,	Qs[[i]]->(Qs[[i]]/.Ord1Eqsp1[[i]])+ListCoeffs[[count++]]u^2+ListCoeffs[[count++]] u v + ListCoeffs[[count++]]v^2]];
,{i,Length[Qs]}];

tbz={};
Do[
Do[
AppendTo[tbz,SeriesCoefficient[EqListp1/.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]}]]]
Ord2Eqsp1=Simplify[ParamTanSpaceOrd2/.SecondOrderCoeffs[[1]]]

Now we compute the order 2 neighborhood for p2

(*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*)
rk=0;
EqListp2={};
rows={};
Do[
newRow=JacobianFPP[[i]]/.p2;
If[MatrixRank[Append[rows,newRow]]>rk,AppendTo[EqListp2,EqsFPPKeum[[i]]]; rk++;AppendTo[rows,newRow]];
,{i,Length[EqsFPPKeum]}];

(*Now, let's compute the order 1 tangent space*)
ListCoeffs=Table[ToExpression[StringJoin["a",ToString[i]]],{i,1,1000}];
count=1;
ParamTanSpaceOrd1={};
Do[
AppendTo[ParamTanSpaceOrd1, Which[i==9, Q8->1,i==7, Q6->v,i==5,Q4->u,True,	Qs[[i]]->ListCoeffs[[count++]]u+ListCoeffs[[count++]]v]];
,{i,Length[Qs]}];

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 tangent space*)
ParamTanSpaceOrd2={};
Do[
AppendTo[ParamTanSpaceOrd2, Which[i==9, Q8->1,i==7, Q6->v,i==5,Q4->u,True,	Qs[[i]]->(Qs[[i]]/.Ord1Eqsp2[[i]])+ListCoeffs[[count++]]u^2+ListCoeffs[[count++]] u v + ListCoeffs[[count++]]v^2]];
,{i,Length[Qs]}];
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=Simplify[ParamTanSpaceOrd2/.SecondOrderCoeffs[[1]]]

We can now generate S3

Mod7Weights={0,-1,-4,-2,-3 ,-5,-6,-3,-5,-6}; (*C_7 weights on the variables, as in the paper*)
deg=2;
Powers={d1,d2,d3,d4,d5,d6,d7,d8,d9,d10};
AllPowers=
Solve[Join[Table[Powers[[k]]≥0,{k,1,10}],{Sum[Powers[[k]],{k,1,10}]==deg}],Integers];

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]}];
(*Generates list of suitable monomials*)
Monomials=Table[Product[Qs[[n]]^(Powers[[n]]/.CorrectPowers[[k]]),{n,1,10}],{k,1,Length[CorrectPowers]}];

ListCoeffs=Table[ToExpression[StringJoin["b",ToString[i]]],{i,1,Length[Monomials]}];

S3=Sum[ListCoeffs[[i]] Monomials[[i]],{i,1,Length[Monomials]}];

tbz={};
Do[
Do[
AppendTo[tbz,SeriesCoefficient[S3/.Ord2Eqsp1,{u,0,i},{v,0,j}]];
AppendTo[tbz,SeriesCoefficient[S3/.Ord2Eqsp2,{u,0,i},{v,0,j}]];
,{j,0,2-i}];
,{i,0,2}];
CoeffsS3=FullSimplify[Solve[Table[tbz[[i]]==0,{i,1,Length[tbz]}]]]
S3=Together[S3/.CoeffsS3]

bsS3=Intersection[Variables[S3], ListCoeffs];

We will now compute random points on the FPP

(*Qpoints=Table[NSolve[Join[{Q0==RandomInteger[100]+1,Q1==RandomInteger[100]+1,Q2==RandomInteger[100]+1},Table[0==EqsFPPKeum[[i]],{i,Length[EqsFPPKeum]}]],1000][[1]],{j,100}];*)
(*To save time, we will use the pre-computed points in "FPPPoints.txt". If you wish, you may recompute points using the code above.*)
fname="Dependencies\\FPPPoints.txt"; OpenRead[fname]; Qpoints=Read[fname];Close[fname];

Using the order 3 automorphism on the surface, we will use S3 to obtain S5 and S6 (quadratics in degree 6). Then, we will get D3=S3 S5 S6 (a sextic in degree 0). We will use the points we just generated to generate a list of values, in the coefficients, of S3 S5 S6 - D3 = 0.

Order3Aut={Q1->Q2, Q2->Q3,Q3->Q1,Q4->Q5,Q5->Q6,Q6->Q4,Q7->Q8,Q8->Q9,Q9->Q7}; (*Order 3 automorphism*)
S5=Simplify[S3 /. Order3Aut];
S6=Simplify[S5/.Order3Aut];

(*Since D3 is of degree 0 and is C_7 invariant, we have by degree and invariant conditions the following form for D3*)
es={e1,e2,e3};
D3=(e1 Q0^2+e2(Q1 Q6+Q2 Q4+Q3 Q5)+e3(Q1 Q9+Q2 Q7+Q3 Q8))^3;
eAndbs=Join[es,bsS3];

(*tries to compute the matrix with the ceofficients e_i and b_i as the independent variables*)

Monomialseb=MonomialList[(S3 S5 S6 - D3),eAndbs][[1]]/(MonomialList[(S3 S5 S6 - D3),eAndbs][[1]]/.Table[eAndbs[[i]]->1,{i,1,Length[eAndbs]}]);

ListCoeffsAndPs=N[MonomialList[(S3 S5 S6 - D3),eAndbs][[1]]/Monomialseb,1000];


ListCoeffsAndPsValues=Table[ListCoeffsAndPs/.Qpoints[[i]],{i,1,Length[Qpoints]}];
MatrixRank[ListCoeffsAndPsValues]

ListCoeffsAndPsReducedValues=Block[{$MaxExtraPrecision=1000},RootApproximant[Chop[RowReduce[ListCoeffsAndPsValues],10^(-900)]]];

(*We will export this to magma to reduce mod 37 and take linear cuts to find a solution*)
(*Since we need to take it mod p, want to make the coefficients integers*)

DenominatorOfAll=Apply[LCM,Join[Table[Apply[LCM,Denominator[ListCoeffsAndPsReducedValues[[i]]]],{i,1,Length[ListCoeffsAndPsReducedValues]}]]];
ListCoeffsAndPsReducedValues=Cancel[ListCoeffsAndPsReducedValues*DenominatorOfAll];

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

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

fname=OpenWrite["Coeffs Vanishing.txt"];
Write[fname, tbzrr];
Close[fname];

Our magma calculations find one of the solutions mod 37. Here we will attempt to lift it back up.

varseb={b1,b2,b3,b4,b5,b8,e1,e2,e3};
(*We will make the further assumption that e1==1. So we will eventually find D3 up to scaling*)
Point1mod37=Solve[Table[varseb[[i]]=={19,1,34,7,9,0,1,16,17}[[i]],{i,Length[varseb]}]][[1]]
Mod[tbzrr/.rr->17/.Point1mod37,37]

(*Let's also make a record of the two other points mod 37*)
Point2mod37=Solve[Table[varseb[[i]]=={13,26,33,34,12,0,1,16,17}[[i]],{i,Length[varseb]}]][[1]]
Point3mod37=Solve[Table[varseb[[i]]=={5,10,7,33,16,0,1,16,17}[[i]],{i,Length[varseb]}]][[1]]
Mod[tbzrr/.rr->17/.Point2mod37,37]
Mod[tbzrr/.rr->17/.Point3mod37,37]

We will now try to lift the solutions to the p-adics, and then to the complex numbers

NewPt={v1,v2,v3,v4,v5,v6,v7,v8,v9};
Pt=varseb/.Point3mod37
rvalue=17;

maxdep=200;

For[dep=1,dep≤maxdep,dep++,
If[Mod[dep,20]==0,Print[dep," "];];

(* computing the square root of (-7) *)
For[i=0,i≤36,i++,
If[Mod[(rvalue+37^dep i)^2+7,37^(dep+1)]==0,
rvalue+=37^dep i;
Break;]  (*Found the new rvalue, and stop the loop*)
];

(* trying on new point *)
Ptplus=Table[varseb[[i]]->(Pt[[i]]+37^dep  NewPt[[i]]),{i,1,Length[NewPt]}];
Eqs={0==v7}; (* we want to normalize e1=1 in the solution *)

For[i=1,i≤Length[tbzrr],i++,
AppendTo[Eqs,0==Simplify[((tbzrr[[i]]/.rr->rvalue)/.Ptplus)/37^dep]] (*Problem is here*)
];

sols=Solve[Eqs,Modulus->37];
Pt=(varseb/.Ptplus)/.sols[[1]]
];
Mod[(tbzrr/.rr->rvalue)/.Table[varseb[[i]]->Pt[[i]],{i,1,Length[Pt]}],37^(dep)]
IntegerExponent[(tbzrr/.rr->rvalue)/.Table[varseb[[i]]->Pt[[i]],{i,1,Length[Pt]}], 37]

(* we will now try to translate from approximate p-adic to complex *)
smallsols={};
rootscx={};
For[i=1,i≤Length[Pt],i++,
FourNumbers={Pt[[i]],1, Mod[rvalue ,37^(dep-1)], 37^(dep-1)};
vec1={f1,f2,f3,f4}/.Solve[{f1==1,f3==0,f4==0,{f1,f2,f3,f4}.FourNumbers==0}][[1]];
vec2={f1,f2,f3,f4}/.Solve[{f1==0,f3==1,f4==0,{f1,f2,f3,f4}.FourNumbers==0}][[1]];
vec3={f1,f2,f3,f4}/.Solve[{f1==0,f3==0,f4==1,{f1,f2,f3,f4}.FourNumbers==0}][[1]];
AppendTo[smallsols,LatticeReduce[{vec1,vec2,vec3}][[1]]]; (* finding a small linear combination *)
AppendTo[rootscx,-(smallsols[[i]][[2]]+Sqrt[-7] smallsols[[i]][[3]])/(smallsols[[i]][[1]])];
];

rootscx
relations=tbzrr/.rr->Sqrt[-7];
FullSimplify[relations/.Table[varseb[[i]]->rootscx[[i]],{i,1,Length[varseb]}]]

(*The second point Point3mod37 worked! *)

Let’s plug in the solutions to S3, S5, S6, D3.

NewS3=FullSimplify[S3/.Table[varseb[[i]]->rootscx[[i]],{i,1,Length[varseb]}]][[1]] 
NewS5=FullSimplify[S5/.Table[varseb[[i]]->rootscx[[i]],{i,1,Length[varseb]}]][[1]];
NewS6=FullSimplify[S6/.Table[varseb[[i]]->rootscx[[i]],{i,1,Length[varseb]}]][[1]];
NewD3=FullSimplify[D3/.Table[varseb[[i]]->rootscx[[i]],{i,1,Length[varseb]}]]
NewD =FullSimplify[( e1 Q0^2+e2(Q1 Q6+Q2 Q4+Q3 Q5)+e3(Q1 Q9+Q2 Q7+Q3 Q8))/.Table[varseb[[i]]->rootscx[[i]],{i,1,Length[varseb]}]]

Now, we will construct the double fiber (|9H|)

(*We need points on the double fiber, but they are given in Borisov21 in old P coordinates. We will make the coordinate change to Q coordinates.*)

(*Files below are from Borisov21*)
fname="Dependencies\\TenFunctionsKeum";  
OpenRead[fname];
TenFunctions=Read[fname];
Close[fname];

QuadricsW=Flatten[ToExpression[Import["Dependencies\\QuadricsWforKeum.txt", "Data"]]];
(*This is the seventh root of SeventhPower that we need to attach to the function field of Y_0 to get the FPP*)
g7=(1/(121 (u0^3-u1^2)^7 w1^3))32 (-43 I+17 Sqrt[7]) u0^7 (4 I u0^4 w1-4 I u0 u1^2 w1+18 I u0 u1 w2+6 Sqrt[7] u0 u1 w2-3 I u0^2 w3+Sqrt[7] u0^2 w3+5 I u0^2 w4-3 Sqrt[7] u0^2 w4-I u1 w5+3 Sqrt[7] u1 w5-19 I u1 w6-7 Sqrt[7] u1 w6) (66 I u0^2 u1 w2-8 I u0^3 w3+4 Sqrt[7] u0^3 w3+2 I u1^2 w3+10 Sqrt[7] u1^2 w3+31 I u0^3 w4+Sqrt[7] u0^3 w4+34 I u1^2 w4-6 Sqrt[7] u1^2 w4+17 I u0 u1 w5-3 Sqrt[7] u0 u1 w5-142 I u0 u1 w6-6 Sqrt[7] u0 u1 w6)^2;

(*We will get some random points on the double fiber*)
PointsOnDoubleFiber={};
While[Length[PointsOnDoubleFiber]<80,
PointsOnDoubleFiber=Join[PointsOnDoubleFiber, Take[NSolve[
Join[Table[QuadricsW[[i]]==0,{i,1,Length[QuadricsW]}], {u1==0,w1==RandomInteger[{1,100}],w2==RandomInteger[{1,100}],w3==RandomInteger[{1,100}]}]
,1000], 3]];

If[Mod[Length[PointsOnDoubleFiber],10]==0,Print[Length[PointsOnDoubleFiber]]];
];
(*We will now input the random points we obtained from the double fiber. In particular, we will pick a seventh root of g7, which we will call g. *)
PointsOnDoubleFiberWithg={};
For[i=1,i<=Length[PointsOnDoubleFiber],i++,
If[Mod[i,20]==0,Print[i]];
gvalue=Solve[g^7==g7/.PointsOnDoubleFiber[[i]]][[1]];
AppendTo[PointsOnDoubleFiberWithg,Join[PointsOnDoubleFiber[[i]],gvalue]];
]

Ps={P0,P1,P2,P3,P4,P5,P6,P7,P8,P9};
PPointsOnDoubleFiber=Table[Ps[[i]]->(TenFunctions[[i]]/.PointsOnDoubleFiberWithg[[j]]),{j,1,Length[PointsOnDoubleFiberWithg]},{i,1,Length[Ps]}];
QtoP={Q0->P0,Q1->((4 I P1)/(-I+(Sqrt[7]))),Q2->((4 I P2)/(-I+(Sqrt[7]))),Q3->((4 I P3)/(-I+(Sqrt[7]))),Q4->((1856 I P4)/(15 I+(Sqrt[7]))),Q5->((1856 I P5)/(15 I+(Sqrt[7]))),Q6->((1856 I P6)/(15 I+(Sqrt[7]))),Q7->-((464 ((67 I+47 (Sqrt[7])) P4+4 (15 I+(Sqrt[7])) P7))/(7 (-27 I+4 (Sqrt[7])))),Q8->-((464 ((67 I+47 (Sqrt[7])) P5+4 (15 I+(Sqrt[7])) P8))/(7 (-27 I+4 (Sqrt[7])))),Q9->-((464 ((67 I+47 (Sqrt[7])) P6+4 (15 I+(Sqrt[7])) P9))/(7 (-27 I+4 (Sqrt[7]))))};
QPointsOnDoubleFiber=QtoP/.PPointsOnDoubleFiber;

(*We have now successfully constructed points on the double fiber in Q coordinates. Now we want to look for homogenous polynomials that vanish at these points.*)
(*We will try degree 3 polynomials. We have tried degree 2 polynomials, but the null space is empty and they do not exist*)
Mod7Weights={0,-1,-4,-2,-3 ,-5,-6,-3,-5,-6}; (*C_7 weights on the variables, as in the paper*)

deg=3;
Deg3EqsDoubleFiber={};

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

(*Find the monomials that respect the C7 weighting*)
For[mod7=0,mod7≤6,mod7++,
Print["mod7=",mod7];

CorrectPowers={};
For[k=1,k≤Length[AllPowers],k++,
If[Mod[Sum[Mod7Weights[[k]] Powers[[k]],{k,1,10}]/.AllPowers[[k]],7]==mod7,
AppendTo[CorrectPowers,AllPowers[[k]]];
];
];
(*Generates list of suitable monomials*)
Monomials=Table[Product[Qs[[n]]^(Powers[[n]]/.CorrectPowers[[k]]),{n,1,10}],{k,1,Length[CorrectPowers]}];

TableOfValues=Table[Monomials[[k]]/.QPointsOnDoubleFiber[[l]],{l,1,Length[QPointsOnDoubleFiber]},{k,1,Length[Monomials]}];
Print[Dimensions[TableOfValues]];
nullspace=RootApproximant[Chop[RowReduce[NullSpace[TableOfValues]],10^(-900)]];
Deg3EqsDoubleFiber=Join[Deg3EqsDoubleFiber,Table[Monomials.nullspace[[i]],{i,1,Length[nullspace]}]];
];
Length[Deg3EqsDoubleFiber]
Table[ByteCount[Deg3EqsDoubleFiber[[k]]],{k,1,Length[Deg3EqsDoubleFiber]}]

fname="Deg3EqsDoubleFiberQ.txt";
OpenWrite[fname];
Write[fname, Deg3EqsDoubleFiber];
Close[fname];

We now have all the ingredients necessary to compute |5H|

DoubleFiber=Deg3EqsDoubleFiber;

PointsS3={};
For[k=1,k<=10,k++,
AppendTo[PointsS3,NSolve[Join[{0==NewS3, 0==NewD, 1==Q0, RandomInteger[30]==Q1},Table[0==EqsFPPKeum[[i]],{i,1,Length[EqsFPPKeum]}]],1000][[1]]];
];

(*Let's compute the Z's*)
Relations16=Partition[DoubleFiber,16]; (*Splitting it into the 112/7=16 eigenspaces of C_7 action*)

DoubleFiberValues=Table[Table[Relations16[[k]]/.PointsS3[[i]],{i,1,Length[PointsS3]}], {k,1,Length[Relations16]}];
Table[MatrixRank[DoubleFiberValues[[i]]],{i,1,Length[DoubleFiberValues]}]

(*Now for each space we take the nullspace. We will have to ignore the fourth index since the weights are shifted down*)
SevenNewBasis=Table[RootApproximant[Chop[RowReduce[NullSpace[DoubleFiberValues[[i]]]],10^(-900)],2], {i,1,Length[DoubleFiberValues]}];
(*Weights are shifted down three. We will have Z1,...,Z6, but Z1 really corresponds to index 4 in DoubleFiberValues and SevenNewBasis*)
ZFunctions={Z1,Z2,Z3,Z4,Z5,Z6};
For[i=1,i<=6,i++,
Print["i: ", i];

ShiftUp3Index=Mod[i+3,7]+1; (*adding +1 to account for 1-indexing*)
FoundBasis=False;
k=1;

While[FoundBasis==False,
If[k>Length[SevenNewBasis[[ShiftUp3Index]]], Print["WTF"];];

ZFunctions[[i]]=SevenNewBasis[[ShiftUp3Index]][[k]].Relations16[[ShiftUp3Index]];

If[!(AllTrue[ZFunctions[[i]]/.Qpoints, PossibleZeroQ]),
FoundBasis=True
];
k++;
];
];

fname="ZsIn5H.txt";
OpenWrite[fname];
Write[fname, ZFunctions];
Close[fname];

The construction of the six sections (in the code, “ZFunctions”) is not C3-invariant. We will try to express the sections in a nicer basis so as to reconstruct the C3 group action. To do this, we compute some C3-fixed points.

FixedPts=Solve[Join[{Q0==1,Q1==Q2,Q2==Q3,Q4==Q5,Q5==Q6,Q7==Q8,Q8==Q9},Table[EqsFPPKeum[[i]]==0,{i,Length[EqsFPPKeum]}]]];
Length[FixedPts]
FixedPt=FixedPts[[1]]

We consider the value of the ratio between each section and the first section at the fixed point. We try to simplify this value by multiplying each section by a scalar.

ZFunctionNew ={};
AppendTo[ZFunctionNew,ZFunctions[[1]]];

Fixing second section

FixedVal2=RootReduce[(ZFunctions[[2]]/ZFunctions[[1]])/.FixedPt]
MinimalPolynomial[RootReduce[FixedVal2]]

FactorInteger[1927561216 ]

denom=2*7;
MinimalPolynomial[FixedVal2*denom]

FactorInteger[37256083456]

sols2=Solve[a^2+7b^2==2^3 263, Integers];
minpolyarr2=Table[MinimalPolynomial[RootReduce[FixedVal2*denom/((a+b(Sqrt[-7]))/.sols2[[i]])]],{i,1,Length[sols2]}];
Sort[Table[{ByteCount[minpolyarr2[[i]]],i},{i,1,Length[minpolyarr2]}]]

minpolyarr2[[1]]
sols2[[1]]

AppendTo[ZFunctionNew,RootReduce[ZFunctions[[2]]*denom/(a+b(Sqrt[-7]))/.sols2[[1]]]];

Fixing third section

FixedVal3=RootReduce[(ZFunctions[[3]]/ZFunctions[[1]])/.FixedPt]
MinimalPolynomial[RootReduce[FixedVal3]]

FactorInteger[729 ]

denom=3;
MinimalPolynomial[FixedVal3*denom]

FactorInteger[780448]

sols3=Solve[a^2+7b^2==2^3 29, Integers];
minpolyarr3=Table[MinimalPolynomial[RootReduce[FixedVal3*denom/((a+b(Sqrt[-7]))/.sols3[[i]])]],{i,1,Length[sols3]}];
Sort[Table[{ByteCount[minpolyarr3[[i]]],i},{i,1,Length[minpolyarr3]}]]

minpolyarr3[[1]]
sols3[[1]]

AppendTo[ZFunctionNew,RootReduce[ZFunctions[[3]]*denom/(a+b(Sqrt[-7]))/.sols3[[1]]]];

Fixing fourth section

FixedVal4=RootReduce[(ZFunctions[[4]]/ZFunctions[[1]])/.FixedPt]
MinimalPolynomial[RootReduce[FixedVal4]]

FactorInteger[194672 ]

denom=2 23;
MinimalPolynomial[FixedVal4*denom]

FactorInteger[34187517952]

sols4=Solve[a^2+7b^2==2^4 7 23, Integers];
minpolyarr4=Table[MinimalPolynomial[RootReduce[FixedVal4*denom/((a+b(Sqrt[-7]))/.sols4[[i]])]],{i,1,Length[sols4]}];
Sort[Table[{ByteCount[minpolyarr4[[i]]],i},{i,1,Length[minpolyarr4]}]]

minpolyarr4[[1]]
sols4[[1]]

AppendTo[ZFunctionNew,RootReduce[ZFunctions[[4]]*denom/(a+b(Sqrt[-7]))/.sols4[[1]]]];

Fixing fifth section

FixedVal5=RootReduce[(ZFunctions[[5]]/ZFunctions[[1]])/.FixedPt]
MinimalPolynomial[RootReduce[FixedVal5]]

FactorInteger[5488 ]

denom=2 7;
MinimalPolynomial[FixedVal5*denom]

FactorInteger[524386566144]

sols5=Solve[a^2+7b^2==2^7 3^2 7, Integers];
minpolyarr5=Table[MinimalPolynomial[RootReduce[FixedVal5*denom/((a+b(Sqrt[-7]))/.sols5[[i]])]],{i,1,Length[sols5]}];
Sort[Table[{ByteCount[minpolyarr5[[i]]],i},{i,1,Length[minpolyarr5]}]]

minpolyarr5[[1]]
sols5[[1]]

AppendTo[ZFunctionNew,RootReduce[ZFunctions[[5]]*denom/(a+b(Sqrt[-7]))/.sols5[[1]]]];

Fixing sixth section

FixedVal6=RootReduce[(ZFunctions[[6]]/ZFunctions[[1]])/.FixedPt]
MinimalPolynomial[RootReduce[FixedVal6]]

FactorInteger[2916 ]

denom=2 3;
MinimalPolynomial[FixedVal6*denom]

FactorInteger[232202043392]

sols6=Solve[a^2+7b^2==2^5 11^2, Integers];
minpolyarr6=Table[MinimalPolynomial[RootReduce[FixedVal6*denom/((a+b(Sqrt[-7]))/.sols6[[i]])]],{i,1,Length[sols6]}];
Sort[Table[{ByteCount[minpolyarr6[[i]]],i},{i,1,Length[minpolyarr6]}]]

minpolyarr6[[1]]
sols6[[1]]

AppendTo[ZFunctionNew,RootReduce[ZFunctions[[6]]*denom/(a+b(Sqrt[-7]))/.sols6[[1]]]];

We’re done fixing the group action.  Let’s see what the new values look like.

Table[RootReduce[(ZFunctionNew[[i]]/.FixedPt)/(ZFunctionNew[[1]]/.FixedPt)],{i,6}]

Let’s compute the 56 degree 6 equations

(*Now we generate points in the image of the Zs, with the goal of finding degree 6 polynomials in the Zs that vanish on it*) 
Zs={Z1,Z2,Z3,Z4,Z5,Z6};

Zpoints=Factor[Table[Zs[[j]]->ZFunctionNew[[j]]/.Qpoints[[i]],{i,1,Length[Qpoints]},{j,1,Length[Zs]}]];

(*Now we will generate the polynomials*)
EqsFPPInP5={};


Zweights={1,2,3,4,5,6};
(*We will split the calculation up by the weights of the polynomials based on C_7 action. 56/7=8*)
deg=6;
powers={a1,a2,a3,a4,a5,a6};
allsols=
Solve[Join[Table[powers[[k]]≥0,{k,1,Length[powers]}],{Sum[powers[[k]],{k,1,Length[powers]}]==deg}],Integers];
(*Weights mod 7*)

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

(*Generate general form of the degree 6 polynomials*)
sols={};
For[k=1,k≤Length[allsols],k++,
If[Mod[Sum[Zweights[[k]] powers[[k]],{k,1,Length[powers]}]/.allsols[[k]],7]==mod7,
AppendTo[sols,allsols[[k]]];
];
];
Print["sols done"];

(*Generate table of values of the polynomials, with Qpoints filled in*)
Monomials=Table[Product[Zs[[n]]^(powers[[n]]/.sols[[k]]),{n,1,Length[powers]}],{k,1,Length[sols]}];
Print["monoms done"];
TableOfValues=Table[Monomials[[k]]/.Zpoints[[l]],{l,1,Length[Zpoints]},{k,1,Length[Monomials]}];
(*Approximate the zero set of polynomials with RootApproximant*)
Print["values done"];
ns=Block[{$MaxExtraPrecision=100000},RootApproximant[Chop[RowReduce[NullSpace[TableOfValues],Tolerance->10^(-500)],10^(-900)],2]];
(*The final equations*)
Print["ns done"];
EqsFPPInP5=Join[EqsFPPInP5,Table[Sum[Monomials[[k]]ns[[i]][[k]],{k,1,Length[Monomials]}],{i,1,Length[ns]}]];
Print[Length[EqsFPPInP5]];
];
Length[EqsFPPInP5]
ByteCount[EqsFPPInP5]

Let’s check the group action

cs={c1,c2,c3,c4,c5,c6};
NoScaleC3Aut={Z1->  Z2, Z2->   Z4,Z3-> Z6,Z4-> Z1,Z5-> Z3,Z6-> Z5 };
EqsFPPInP5NoScaleC3 = Together[EqsFPPInP5/.NoScaleC3Aut];

GetMonoms[poly_,vars_]:=MonomialList[poly,vars]/(MonomialList[poly,vars]/.Table[vars[[i]]->1,{i,1,Length[vars]}]);

(*Take the first equation and look at its monomials*)
Eq1Monoms=GetMonoms[EqsFPPInP5NoScaleC3[[1]],Zs];
(*Let's try to solve for coefficients that give us zero*)
as=Table[ToExpression[StringJoin["a",ToString[i]]],{i,8}];
CurrentMonoms=GetMonoms[Sum[EqsFPPInP5[[i]],{i,1,8}],Zs];
MonomsTBZ=Complement[CurrentMonoms,Eq1Monoms];
CoeffsTBZ=RowReduce[Transpose[Table[Coefficient[EqsFPPInP5[[i]],MonomsTBZ[[j]]],{i,8},{j,1,Length[MonomsTBZ]}]]];
Dimensions[CoeffsTBZ]
MatrixRank[CoeffsTBZ]
sols=RootReduce[Solve[Table[CoeffsTBZ[[i]].as==0,{i,1,Length[CoeffsTBZ]}]][[1]]]

LinComb=RootReduce[Table[EqsFPPInP5[[i]],{i,8}].(as/.sols)];
LinCombMonoms=GetMonoms[LinComb,Zs];
LinCombMonoms==Eq1Monoms

C3Aut={Z1->  c1 Z2, Z2->c2   Z4,Z3->c3 Z6,Z4->c4 Z1,Z5-> c5 Z3,Z6->c6 Z5 };
EqsFPPInP5C3 = Together[EqsFPPInP5/.C3Aut];
CoeffConditions=Join[{c1 !=0,c2!=0,c3!=0,c4!=0,c5!=0,c6!=0},
Table[Coefficient[EqsFPPInP5C3[[1]],LinCombMonoms[[i]]]==Coefficient[LinComb,LinCombMonoms[[i]]],{i,1,Length[LinCombMonoms]}]];
csols=Simplify[Solve[CoeffConditions], Assumptions->c2>0];
Length[csols]

(*Out of the 7 solutions, there is one such that c1 * c2 * c4 is in the correct field. We can just check each one to see that the 5th solution is defined over the field *)
c124=RootReduce[((c1 c2 c4)/.csols[[5]])/.c2->1]

RootReduce[Solve[(1+3(Sqrt[-7]))/8==c2^3 c124  ]]
c2val =c2/. RootReduce[Solve[(1+3(Sqrt[-7]))/8==c2^3 c124  ]][[1]]

(*Similarly, we pick the value of c2 that gives the correct field. We chose (1+3(Sqrt[-7]))/8 through experimentation*)
csvals =RootApproximant[ N[cs/.csols[[5]]/.c2 ->c2val,1000]]

(*Now we can define new sections and recompute the equations*)
SixSections=ZFunctionNew;
SixSections[[2]]=csvals[[1]]ZFunctionNew[[2]];
SixSections[[3]]=ZFunctionNew[[3]];
SixSections[[4]]=1/csvals[[4]] ZFunctionNew[[4]];
SixSections[[5]]=1/csvals[[5]]ZFunctionNew[[5]];
SixSections[[6]]=csvals[[3]]ZFunctionNew[[6]];

fname=OpenWrite["SixSections.txt"]
Write[fname,SixSections];
Close[fname];

NewZpoints=Factor[Table[Zs[[j]]->SixSections[[j]]/.Qpoints[[i]],{i,1,Length[Qpoints]},{j,1,Length[Zs]}]];

(*Now we will generate the polynomials*)
EqsFPPInP5={};


Zweights={1,2,3,4,5,6};
(*We will split the calculation up by the weights of the polynomials based on C_7 action. 56/7=8*)
deg=6;
powers={a1,a2,a3,a4,a5,a6};
allsols=
Solve[Join[Table[powers[[k]]≥0,{k,1,Length[powers]}],{Sum[powers[[k]],{k,1,Length[powers]}]==deg}],Integers];
(*Weights mod 7*)

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

(*Generate general form of the degree 6 polynomials*)
sols={};
For[k=1,k≤Length[allsols],k++,
If[Mod[Sum[Zweights[[k]] powers[[k]],{k,1,Length[powers]}]/.allsols[[k]],7]==mod7,
AppendTo[sols,allsols[[k]]];
];
];
Print["sols done"];

(*Generate table of values of the polynomials, with Qpoints filled in*)
Monomials=Table[Product[Zs[[n]]^(powers[[n]]/.sols[[k]]),{n,1,Length[powers]}],{k,1,Length[sols]}];
Print["monoms done"];
TableOfValues=Table[Monomials[[k]]/.NewZpoints[[l]],{l,1,Length[NewZpoints]},{k,1,Length[Monomials]}];
(*Approximate the zero set of polynomials with RootApproximant*)
Print["values done"];
ns=Block[{$MaxExtraPrecision=100000},RootApproximant[Chop[RowReduce[NullSpace[TableOfValues],Tolerance->10^(-500)],10^(-900)],2]];
(*The final equations*)
Print["ns done"];
EqsFPPInP5=Join[EqsFPPInP5,Table[Sum[Monomials[[k]]ns[[i]][[k]],{k,1,Length[Monomials]}],{i,1,Length[ns]}]];
Print[Length[EqsFPPInP5]];
];
Length[EqsFPPInP5]
ByteCount[EqsFPPInP5]

fname="EqsFPPInP5WithZCoords.txt";
OpenWrite[fname];
Write[fname, EqsFPPInP5];
Close[fname];

(*Let's also store it in magma format with rr to check the Hilbert polynomial*)
EqsFPPInP5withrr={};
AllMonomials=Table[Product[Zs[[n]]^(powers[[n]]/.allsols[[k]]),{n,1,Length[Zs]}],{k,1,Length[allsols]}];
For[numeq=1,numeq≤Length[EqsFPPInP5],numeq++,
eq=Expand[EqsFPPInP5[[numeq]]];
eqwithrr=0;
For[k=1,k≤Length[AllMonomials],k++,
coef=Coefficient[eq,AllMonomials[[k]]];
eqwithrr+=(Re[coef]  + rr Im[coef/(Sqrt[7])]) AllMonomials[[k]];
];
AppendTo[EqsFPPInP5withrr,eqwithrr];
];
(* just checking *)
Simplify[(EqsFPPInP5withrr/.rr->((Sqrt[-7])))-EqsFPPInP5]

fname="EqsFPPInP5WithZCoordsrr.txt";
OpenWrite[fname];
Write[fname, EqsFPPInP5withrr];
Close[fname];

EqsFPPInP5withrrn=Numerator[Factor[EqsFPPInP5withrr]];
fname="EqsFPPInP5WithZCoordsrrn.txt";
OpenWrite[fname];
Write[fname, EqsFPPInP5withrrn];
Close[fname];

(*Recomputing the action*)
C3Action={Z1->Z2,Z2->(1+3(Sqrt[-7]))/8 Z4,Z4->Z1,Z3->Z6,Z6->(1+3(Sqrt[-7]))/8 Z5, Z5->Z3}
(*Verify that it leaves the rank unchanged*)
AllEqsUnderAction=Join[EqsFPPInP5/.C3Action,EqsFPPInP5];
AllMonoms=GetMonoms[Sum[AllEqsUnderAction[[i]],{i,1,Length[AllEqsUnderAction]}],Zs];
coeffs=Table[Coefficient[AllEqsUnderAction,AllMonoms[[i]]],{i,1,Length[AllMonoms]}];
Dimensions[coeffs]
MatrixRank[coeffs]

Finally, we further simplify the equations by taking random linear combinations of the equations in each weight and keeping those that have smallest bytecount and are linearly independent from those previously chosen.

GetMonoms[poly_,vars_]:=MonomialList[poly,vars]/(MonomialList[poly,vars]/.Table[vars[[i]]->1,{i,1,Length[vars]}]);
Zs={Z1,Z2,Z3,Z4,Z5,Z6}

BeautifiedEqsFPPInP5={};
EqsFPPInP5ByWeight=Partition[EqsFPPInP5,8];
DesiredMaxMonoms=55;
For[mod7=1, mod7<=7,mod7++,
Print["Graded part: ",mod7];
GradeEqs=EqsFPPInP5ByWeight[[mod7]];
cs=Table[ToExpression[StringJoin["c",ToString[i]]],{i,1,Length[GradeEqs]}];
GradeMonoms={};
Do[
GradeMonoms=Union[GradeMonoms,GetMonoms[GradeEqs[[i]],Zs]];
,{i,Length[GradeEqs]}];

NumGradeMonoms=Length[GradeMonoms];
(*Each column represents a polynomial. Rows have same monomial*)
CoeffMatrix=Table[Coefficient[GradeEqs,GradeMonoms[[i]]],{i,1,Length[GradeMonoms]}];
GenLinComb=GradeEqs.cs;
	AllSolutions={};

For[iteration=1,iteration<=5000,iteration++,
If[Mod[iteration,1000]==0, Print["Iteration: ",iteration, " #Solutions: ", Length[AllSolutions]]];
RandomToBeOne=RandomInteger[{1,Length[GradeMonoms]}];
IndicesInUse={RandomToBeOne};
MatrixCoeffsInUse={CoeffMatrix[[RandomToBeOne]]};

EquationsToSolve={CoeffMatrix[[RandomToBeOne]].cs == 1};

If[EquationsToSolve[[1]]==False,Print["WTF"];
Print[RandomToBeOne];Print[CoeffMatrix[[RandomToBeOne]].cs];
];

(*Now we set monomials to zero*)
For[i=1, MatrixRank[N[MatrixCoeffsInUse]]<Length[GradeEqs], i++,
RandomToBeZero=RandomChoice[Complement[Range[Length[GradeMonoms]],IndicesInUse]];
AppendTo[IndicesInUse,RandomToBeZero];

TestMatrixCoeffsInUse=Append[MatrixCoeffsInUse,CoeffMatrix[[RandomToBeZero]]];

If[MatrixRank[N[TestMatrixCoeffsInUse]]>MatrixRank[N[MatrixCoeffsInUse]],
AppendTo[EquationsToSolve, 0==CoeffMatrix[[RandomToBeZero]].cs];
MatrixCoeffsInUse=TestMatrixCoeffsInUse;
];
];

SolutionsC=Chop[NSolve[EquationsToSolve,cs,50],10^(-900)];
NewEquation=Chop[Expand[GenLinComb/.SolutionsC[[1]]]];


NumMonoms=Length[MonomialList[NewEquation]];
If[NumMonoms<=DesiredMaxMonoms,
(*Generate new polynomial with RootApproximant.*)
NewEquation=Sum[RootApproximant[Coefficient[NewEquation,GradeMonoms[[i]]],2] GradeMonoms[[i]],{i,1,Length[GradeMonoms]}];
AppendTo[AllSolutions,NewEquation];
];
];

(*Sort by ByteCount*)

SortedSolutionsIndex=Sort[Table[{ByteCount[AllSolutions[[i]]],i},{i,1,Length[AllSolutions]}]];
BestEqs={};
currentrank=0;
For[i=1, i<=Length[SortedSolutionsIndex],i++,
If[Mod[i,500]==0,Print["Iteration: ",i]];
index=SortedSolutionsIndex[[i]][[2]];
BestEqsTry=Append[BestEqs,AllSolutions[[index]]];
MatrixCoeffsBest=Table[Coefficient[BestEqsTry,GradeMonoms[[j]]],{j,1,Length[GradeMonoms]}];
If[MatrixRank[N[MatrixCoeffsBest]]>currentrank,
BestEqs=BestEqsTry;
currentrank++;
];
];


Sort[Table[Length[MonomialList[BestEqs[[i]]]],{i,1,Length[BestEqs]}]];
ByteCount[BestEqs];
AppendTo[BeautifiedEqsFPPInP5,BestEqs];
];
Length[BeautifiedEqsFPPInP5]
ByteCount[BeautifiedEqsFPPInP5]

BeautifiedEqsFPPInP5=Flatten[BeautifiedEqsFPPInP5];

fname="BeautifiedEqsFPPInP5WithZCoords.txt";
OpenWrite[fname];
Write[fname, BeautifiedEqsFPPInP5];
Close[fname];

(*Let's also store it in magma format with rr to check the Hilbert polynomial*)
BeautifiedEqsFPPInP5withrr={};
AllMonomials=Table[Product[Zs[[n]]^(powers[[n]]/.allsols[[k]]),{n,1,Length[Zs]}],{k,1,Length[allsols]}];
For[numeq=1,numeq≤Length[EqsFPPInP5],numeq++,
eq=Expand[BeautifiedEqsFPPInP5[[numeq]]];
eqwithrr=0;
For[k=1,k≤Length[AllMonomials],k++,
coef=Coefficient[eq,AllMonomials[[k]]];
eqwithrr+=(Re[coef]  + rr Im[coef/(Sqrt[7])]) AllMonomials[[k]];
];
AppendTo[BeautifiedEqsFPPInP5withrr,eqwithrr];
];
(* just checking *)
Simplify[(BeautifiedEqsFPPInP5withrr/.rr->((Sqrt[-7])))-BeautifiedEqsFPPInP5]

fname="BeautifiedEqsFPPInP5WithZCoordsrr.txt";
OpenWrite[fname];
Write[fname, BeautifiedEqsFPPInP5withrr];
Close[fname];

BeautifiedEqsFPPInP5withrrn=Numerator[Factor[BeautifiedEqsFPPInP5withrr]];
fname="BeautifiedEqsFPPInP5WithZCoordsrrn.txt";
OpenWrite[fname];
Write[fname, BeautifiedEqsFPPInP5withrrn];
Close[fname];