Simplifying Keum’s FPP

In this notebook, we will present the code to simplify the Keum fake projective plane.

(*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*)
Ps={P0,P1,P2,P3,P4,P5,P6,P7,P8,P9};

fname="EqsFPPwithrrKeum.txt";
OpenRead[fname];
EqsFPPwithrrKeum=Read[fname];
Close[fname];

prime=43;
rootvalue=rr ->6;
cut=P0+24(P1+P2+P3)+0(P4+P5+P6)+28(P7+P8+P9);

PointsMod43={{P0+2*P9,P1,P2+32*P9,P3+23*P9,P4+3*P9,P5+32*P9,P6+21*P9,P7+22*P9,P8+27*P9},{P0+33*P9,P1+6*P9,P2+32*P9,P3+27*P9,P4+22*P9,P5+30*P9,P6+20*P9,P7+7*P9,P8+40*P9},{P0+25*P9,P1+10*P9,P2+38*P9,P3+4*P9,P4+28*P9,P5+19*P9,P6+24*P9,P7+10*P9,P8+30*P9},{P0+27*P9,P1+12*P9,P2+32*P9,P3+37*P9,P4+39*P9,P5+34*P9,P6+24*P9,P7+15*P9,P8+40*P9},{P0+25*P9,P1+35*P9,P2+15*P9,P3+16*P9,P4+11*P9,P5+P9,P6+18*P9,P7+33*P9,P8+38*P9},{P0+26*P9,P1+20*P9,P2+33*P9,P3+36*P9,P4+8*P9,P5+34*P9,P6+3*P9,P7+25*P9,P8+37*P9},{P0+26*P8,P1+29*P8,P2+7*P8,P3+14*P8,P4+5*P8,P5+9*P8,P6+4*P8,P7+34*P8,P9},{P0+39*P9,P1+22*P9,P2+40*P9,P3,P4+22*P9,P5+P9,P6+37*P9,P7+32*P9,P8+2*P9},{P0+38*P9,P1+6*P9,P2+9*P9,P3+15*P9,P4+36*P9,P5+12*P9,P6+14*P9,P7+35*P9,P8+28*P9}};
Pns={Pn0,Pn1,Pn2,Pn3,Pn4,Pn5,Pn6,Pn7,Pn8,Pn9};
(*Keeping track of coefficients *)
Coeffs=Table[ToExpression[StringJoin["c",ToString[m]]],{m,1,2000}];
maxdepth=40;
slsroot=Solve[x^2+7==0,Modulus->prime^maxdepth];
Print[Mod[x/.slsroot,prime]];

(*Starting to build tangent spaces and keeping them the same *)
(*pointscut -> PointsMod29*)
GoodCuts = {};
ptssub = {}; (*base points of tangent space*)
tanvecs={}; (*tangent vectors corresponding to points in ptssub*)
For[k=1,k<=Length[PointsMod43],k++,
pt=Solve[Join[Table[PointsMod43[[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[((EqsFPPwithrrKeum/.rootvalue)/.P0-> 1)/.tanvector, t],{Pn0}];
(*from Borisov SevenfoldCover: 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 -> cut, currentptsub -> ptssub, currenttanvecs -> tanvecs*)
currentcut=cut;
currentptssub=ptssub;
currenttanvecs=tanvecs;

For[depth=1,depth<=30,depth++,
Print["Depth: ", depth];

rootvalue=rr->Mod[x/.slsroot[[2]],prime^(depth+1)];
EqsFPPnewdepthKeum=PolynomialMod[(EqsFPPwithrrKeum/.P0->1)/.rootvalue,prime^(depth+1)];
count=1; (*Equations mod prime^{depth+1}*)

(*make extensions to p^(depth+1), e.g. fit the cuts to the depth+1*)
newcut=currentcut + prime^depth (a1 (P1+P2+P3)+a2 (P4+P5+P6)+a3(P7+P8+P9));
newptssub={};
newtanvecs={};

(*calculating new base points and tangent vectors*)
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=(EqsFPPnewdepthKeum[[k]]/.newtanvecs[[l]][[1]])/prime^depth;
AppendTo[Eqs,PolynomialMod[Factor[tbz/.t->0],prime]];
AppendTo[Eqs,PolynomialMod[Factor[Coefficient[tbz,t]],prime]];

tbz=(EqsFPPnewdepthKeum[[k]]/.newtanvecs[[1]][[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]]; (*This command is causing problems*)
currentcut=newcut/.sls;
currentptssub=newptssub/.sls;
currenttanvecs=newtanvecs/.sls;
];
AppendTo[GoodCuts,currentcut];
GoodCuts

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

r7=Mod[rr/.rootvalue, prime^31];

whatis=Coefficient[GoodCuts[[1]],P1];
thebasis=Append[Join[Table[Mod[whatis^k,prime^31],{k,0,2}],Table[Mod[whatis^k r7,prime^31],{k,0,2}]],prime^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[[1]],P4];
thebasis=Append[Join[Table[Mod[whatis^k,prime^31],{k,0,2}],Table[Mod[whatis^k r7,prime^31],{k,0,2}]],prime^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[[1]],P7];
thebasis=Append[Join[Table[Mod[whatis^k,prime^31],{k,0,2}],Table[Mod[whatis^k r7,prime^31],{k,0,2}]],prime^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 use this to solve for the coefficients of the cut*)
whatis=Coefficient[GoodCuts[[1]],P1];
thebasis=Append[Join[Table[Mod[whatis^k,prime^31],{k,0,1}],Table[Mod[whatis^k r7,prime^31],{k,0,1}]],prime^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[[1]],P4];
thebasis=Append[Join[Table[Mod[whatis^k,prime^31],{k,0,1}],Table[Mod[whatis^k r7,prime^31],{k,0,1}]],prime^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[[1]],P7];
thebasis=Append[Join[Table[Mod[whatis^k,prime^31],{k,0,1}],Table[Mod[whatis^k r7,prime^31],{k,0,1}]],prime^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]];

CorrectCut= RootReduce[P0+a1(P1+P2+P3)+a2(P4+P5+P6)+ a3(P7+P8+P9)] (*The correct cut we were aiming for*)

(*Now we want to verify that the cut is nonreduced. We will do it numerically.*)
sols=NSolve[Join[Table[0==EqsFPPwithrrKeum[[k]]/.rr->I (Sqrt[7]),{k,1,84}],{0==CorrectCut,P0==1,P1==RandomInteger[10]+1}],50];
Length[DeleteDuplicates[sols]]

fname="TenFunctionsKeum"; 
OpenRead[fname];
TenFunctions=Read[fname];
Close[fname];

QuadricsW=Flatten[ToExpression[Import["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;

UWGpts={};

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

prec=1000;
Do[If[Mod[n,20]==0, Print[n]];
sols=NSolve[Join[{u0==1,w1==RandomInteger[100]+ I RandomInteger[100],u1==RandomInteger[100]+ I RandomInteger[100],w2==RandomInteger[100]+ I RandomInteger[100]},Table[0==QuadricsW[[k]],{k,1,9}]],prec];
gvalue=g->(g/.Solve[g^7== g7/.sols[[1]]][[1]]);
AppendTo[UWGpts,Append[sols[[1]],gvalue]];
,{n,40}];

Ppoints=Table[Table[Ps[[n]]-> (TenFunctions[[n]]/.UWGpts[[k]]),{n,1,10}],{k,1,Length[UWGpts]}];

(*Below is the 2-torsion cut, computed previously. We write it here for future convenience*)
CorrectCut=P0+(1/2) (-1+I (Sqrt[7])) (P1+P2+P3)+(1/7) (272-848 I (Sqrt[7])) (P4+P5+P6)+(1/7) (832-192 I (Sqrt[7])) (P7+P8+P9);
(*Below is code to make a coordinate change from P to Q*)
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};
(*Coefficient c2 is chosen by inspection from the resulting equations*)
solcs=RootReduce[Solve[{c2==((15/1856-(I Sqrt[7])/1856)),Coefficient[CorrectCut/.PtoQ,Q1]==1 ,Coefficient[CorrectCut/.PtoQ,Q4]==0,Coefficient[CorrectCut/.PtoQ,Q9]==1}]][[1]]

PtoQ = Simplify[PtoQ/.solcs]
QtoP =Simplify[Solve[Table[Ps[[k]]==(Ps[[k]]/.PtoQ),{k,1,10}],Qs][[1]]]

FullSimplify[CorrectCut/.PtoQ]

(*Now, we will substitute into the original equations for the FPP*)

(*Change P-points to Q-points*)
Qpoints={};
For[k=1,k≤Length[Ppoints],k++,
AppendTo[Qpoints,Table[Qs[[m]]->(( Qs[[m]]/.QtoP)/.Ppoints[[k]]),{m,1,Length[Qs]}]];
];
Accuracy[Ppoints]
Accuracy[Qpoints]
Length[Qpoints]
(*Now, we compute the new equations*)
deg=3;
EqsFPPqKeum={};
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];
(*Weights mod 7*)
Mod7Weights={0,-1,-4,-2,-3,-5,-6,-3,-5,-6};

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

(*Generate general form of the degree 3 polynomials*)
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]]];
];
];
(*Generate table of values of the polynomials, with Qpoints filled in*)
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[Qpoints]},{k,1,Length[Monomials]}];
(*Approximate the zero set of polynomials with RootApproximant*)

ns=Block[{$MaxExtraPrecision=1000},RootApproximant[Chop[RowReduce[Chop[NullSpace[TableOfValues],10^(-900)]],10^(-900)]]];
(*The final equations*)
EqsFPPqKeum=Join[EqsFPPqKeum,Table[Sum[Monomials[[k]]ns[[i]][[k]],{k,1,Length[Monomials]}],{i,1,Length[ns]}]];
];

EqsFPPqKeum=DeleteDuplicates[EqsFPPqKeum];
Length[EqsFPPqKeum]
Table[ByteCount[EqsFPPqKeum[[k]]],{k,1,Length[EqsFPPqKeum]}]

(*Rationalize the denominators*)
EqsFPPqnKeum=Numerator[Factor[EqsFPPqKeum]];
ByteCount[EqsFPPqnKeum]
Table[ByteCount[EqsFPPqnKeum[[k]]],{k,1,Length[EqsFPPqnKeum]}]
Print[EqsFPPqnKeum];

Success! Now, we will run the beautification code

Let’s start with 0-graded C_7 part i.e. indexes 1-12.

GradedEqs=Take[EqsFPPKeum,{1,12}];
cs=Table[ToExpression[StringJoin["c",ToString[i]]],{i,1,Length[GradedEqs]}];
GradedMonoms=MonomialList[Sum[GradedEqs[[i]],{i,1,Length[GradedEqs]}],Qs]/(MonomialList[Sum[GradedEqs[[i]],{i,1,Length[GradedEqs]}],Qs]/.Table[Qs[[i]]->1,{i,1,Length[Qs]}])

NumGradedMonoms=Length[GradedMonoms]

(*Each column represents a polynomial.*)
CoeffMatrix=Table[Coefficient[GradedEqs,GradedMonoms[[i]]],{i,1,Length[GradedMonoms]}];

GenLinComb=GradedEqs.cs;

AllSolutions={};
DesiredMaxMonoms=18;

For[iteration=1,iteration<=3000,iteration++,

If[Mod[iteration,200]==0, Print[iteration]];

RandomToBeOne=RandomInteger[{1,Length[GradedMonoms]}];

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[MatrixCoeffsInUse]<Length[GradedEqs], i++,
RandomToBeZero=RandomChoice[Complement[Range[Length[GradedMonoms]],IndicesInUse]];
AppendTo[IndicesInUse,RandomToBeZero];

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

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,GradedMonoms[[i]]]] GradedMonoms[[i]],{i,1,Length[GradedMonoms]}];
AppendTo[AllSolutions,NewEquation];
(*Print["Iteration: ", ToString[iteration]];
Print[NumMonoms];
Print[NewEquation];*)
];
];

C3Aut={Q0->Q0,Q1->Q2, Q2->Q3,Q3->Q1,Q4->Q5,Q5->Q6,Q6->Q4,Q7->Q8,Q8->Q9,Q9->Q7};

(*Sort by ByteCount*)
Length[AllSolutions]
SortedSolutionsIndex=Sort[Table[{ByteCount[AllSolutions[[i]]],i},{i,1,Length[AllSolutions]}]];
BestEqs={};
currentrank=0;
For[i=1, i<=Length[SortedSolutionsIndex],i++,
index=SortedSolutionsIndex[[i]][[2]];

(*Try the C3 action*)
BestEqsTry=Append[BestEqs,AllSolutions[[index]]];
MatrixCoeffsBest=Table[Coefficient[BestEqsTry,GradedMonoms[[j]]],{j,1,Length[GradedMonoms]}];
If[MatrixRank[MatrixCoeffsBest]>currentrank,
BestEqs=BestEqsTry;
currentrank++;
];

(*Order 1*)
BestEqsTry=Append[BestEqs,AllSolutions[[index]]/.C3Aut];
MatrixCoeffsBest=Table[Coefficient[BestEqsTry,GradedMonoms[[j]]],{j,1,Length[GradedMonoms]}];
If[MatrixRank[MatrixCoeffsBest]>currentrank,
BestEqs=BestEqsTry;
currentrank++;
];

(*Order 2*)
BestEqsTry=Append[BestEqs,(AllSolutions[[index]]/.C3Aut)/.C3Aut];
MatrixCoeffsBest=Table[Coefficient[BestEqsTry,GradedMonoms[[j]]],{j,1,Length[GradedMonoms]}];
If[MatrixRank[MatrixCoeffsBest]>currentrank,
BestEqs=BestEqsTry;
currentrank++;
];

];
currentrank
BestEqs
Sort[Table[Length[MonomialList[BestEqs[[i]]]],{i,1,Length[BestEqs]}]]

Let’s start with 1-graded C_7 part i.e. indexes 13-24.

Grade1Eqs=Take[EqsFPPKeum,{13,24}];
cs=Table[ToExpression[StringJoin["c",ToString[i]]],{i,1,Length[Grade1Eqs]}];
Grade1Monoms=MonomialList[Sum[Grade1Eqs[[i]],{i,1,Length[Grade1Eqs]}],Qs]/(MonomialList[Sum[Grade1Eqs[[i]],{i,1,Length[Grade1Eqs]}],Qs]/.Table[Qs[[i]]->1,{i,1,Length[Qs]}])

NumGrade1Monoms=Length[Grade1Monoms]

(*Each column represents a polynomial. Rows have same monomial*)
CoeffMatrix=Table[Coefficient[Grade1Eqs,Grade1Monoms[[i]]],{i,1,Length[Grade1Monoms]}];

GenLinComb=Grade1Eqs.cs;

AllSolutions={};
DesiredMaxMonoms=18;

For[iteration=1,iteration<=7000,iteration++,

If[Mod[iteration,200]==0, Print["Iteration: ",iteration]];

RandomToBeOne=RandomInteger[{1,Length[Grade1Monoms]}];

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[MatrixCoeffsInUse]<Length[Grade1Eqs], i++,
RandomToBeZero=RandomChoice[Complement[Range[Length[Grade1Monoms]],IndicesInUse]];
AppendTo[IndicesInUse,RandomToBeZero];

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

If[MatrixRank[TestMatrixCoeffsInUse]>MatrixRank[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,Grade1Monoms[[i]]]] Grade1Monoms[[i]],{i,1,Length[Grade1Monoms]}];
AppendTo[AllSolutions,NewEquation];
(*Print["Iteration: ", ToString[iteration], " ", NumMonoms];*)
(*Print[NumMonoms];*)
(*Print[NewEquation];*)
];
];

(*Sort by ByteCount*)
Length[AllSolutions]
SortedSolutionsIndex=Sort[Table[{ByteCount[AllSolutions[[i]]],i},{i,1,Length[AllSolutions]}]];
BestEqs={};
currentrank=0;
For[i=1, i<=Length[SortedSolutionsIndex],i++,
index=SortedSolutionsIndex[[i]][[2]];
BestEqsTry=Append[BestEqs,AllSolutions[[index]]];
MatrixCoeffsBest=Table[Coefficient[BestEqsTry,Grade1Monoms[[j]]],{j,1,Length[Grade1Monoms]}];
If[MatrixRank[MatrixCoeffsBest]>currentrank,
BestEqs=BestEqsTry;
currentrank++;
];
];
currentrank


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

Let’s start with 3-graded C_7 part i.e. indexes 37-48.

GradedEqs=Take[EqsFPPKeum,{37,48}];
cs=Table[ToExpression[StringJoin["c",ToString[i]]],{i,1,Length[GradedEqs]}];
GradedMonoms=MonomialList[Sum[GradedEqs[[i]],{i,1,Length[GradedEqs]}],Qs]/(MonomialList[Sum[GradedEqs[[i]],{i,1,Length[GradedEqs]}],Qs]/.Table[Qs[[i]]->1,{i,1,Length[Qs]}])

NumGradedMonoms=Length[GradedMonoms]

(*Each column represents a polynomial. Rows have same monomial*)
CoeffMatrix=Table[Coefficient[GradedEqs,GradedMonoms[[i]]],{i,1,Length[GradedMonoms]}];

GenLinComb=GradedEqs.cs;
Length[GradedEqs]

AllSolutions={};
DesiredMaxMonoms=20;

For[iteration=1,iteration<=1000,iteration++,

If[Mod[iteration,200]==0, Print[iteration]];

RandomToBeOne=RandomInteger[{1,Length[GradedMonoms]}];

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[MatrixCoeffsInUse]<Length[GradedEqs], i++,
RandomToBeZero=RandomChoice[Complement[Range[Length[GradedMonoms]],IndicesInUse]];
AppendTo[IndicesInUse,RandomToBeZero];

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

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,GradedMonoms[[i]]]] GradedMonoms[[i]],{i,1,Length[GradedMonoms]}];
AppendTo[AllSolutions,NewEquation];
(*Print["Iteration: ", ToString[iteration]];
Print[NumMonoms];
Print[NewEquation];*)
];
];

(*Sort by ByteCount*)
Length[AllSolutions]
SortedSolutionsIndex=Sort[Table[{ByteCount[AllSolutions[[i]]],i},{i,1,Length[AllSolutions]}]];
BestEqs={};
currentrank=0;
For[i=1, i<=Length[SortedSolutionsIndex],i++,
index=SortedSolutionsIndex[[i]][[2]];
BestEqsTry=Append[BestEqs,AllSolutions[[index]]];
MatrixCoeffsBest=Table[Coefficient[BestEqsTry,GradedMonoms[[j]]],{j,1,Length[GradedMonoms]}];
If[MatrixRank[MatrixCoeffsBest]>currentrank,
BestEqs=BestEqsTry;
currentrank++;
];
];
currentrank
BestEqs
Sort[Table[Length[MonomialList[BestEqs[[i]]]],{i,1,Length[BestEqs]}]]

(*The code above generates new, random polynomials that generate the FPP each time it is run. We have collected the equations below in the files "Beautification 0 ByteCount.txt" "Beautification 1 ByteCount.txt" and "Beautification 3 ByteCount.txt". Rerunning the program will require writing to these files or editing the code above to save the equations to Graded0Eqs, Graded1Eqs, Graded3Eqs. *)

(*Load Beautified 0-graded equations*)
fname="Beautification 0 ByteCount.txt";
OpenRead[fname];
Graded0Eqs=Read[fname];
Close[fname];

(*Load beautified 1-graded equations*)
fname="Beautification 1 ByteCount.txt";
OpenRead[fname];
Graded1Eqs=Read[fname];
Close[fname];

(*Load beautified 3-graded equations*)
fname="Beautification 3 ByteCount.txt";
OpenRead[fname];
Graded3Eqs=Read[fname];
Close[fname];

Length[Graded0Eqs]
Length[Graded1Eqs]
Length[Graded3Eqs]

Graded4Eqs=Graded1Eqs/.C3Aut;
Graded2Eqs=Graded4Eqs/.C3Aut;

Graded5Eqs=Graded3Eqs/.C3Aut;
Graded6Eqs=Graded5Eqs/.C3Aut;

Length[Graded2Eqs]
Length[Graded4Eqs]
Length[Graded5Eqs]
Length[Graded6Eqs]

BeautifiedEqsFPP=Join[Graded0Eqs,Graded1Eqs, Graded2Eqs,Graded3Eqs,Graded4Eqs,Graded5Eqs,Graded6Eqs];
Length[BeautifiedEqsFPP]
ByteCount[BeautifiedEqsFPP]

(*Let's also store it in magma format with rr to check the Hilbert polynomial*)
BeautifiedEqsFPPwithrr={};
AllMonomials=MonomialList[Sum[BeautifiedEqsFPP[[i]],{i,1,Length[BeautifiedEqsFPP]}]]/(MonomialList[Sum[BeautifiedEqsFPP[[i]],{i,1,Length[BeautifiedEqsFPP]}]]/.Table[Qs[[i]]->1,{i,1,Length[Qs]}]);
For[numeq=1,numeq≤Length[BeautifiedEqsFPP],numeq++,
eq=Expand[BeautifiedEqsFPP[[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[BeautifiedEqsFPPwithrr,eqwithrr];
];
(* just checking *)
Simplify[(BeautifiedEqsFPPwithrr/.rr->((Sqrt[-7])))-BeautifiedEqsFPP]

(*Save the output*)
fname="BeautifiedFPPwithrr.txt";
OpenWrite[fname];
Write[fname, BeautifiedEqsFPPwithrr];
Close[fname];

(*Set numerators*)
fname="BeautifiedFPPwithrr.txt";
OpenRead[fname];
BeautifiedEqsFPPwithrr=Read[fname];
Close[fname];

BeautifiedEqsFPPnwithrr=Numerator[Factor[BeautifiedEqsFPPwithrr]];
ByteCount[BeautifiedEqsFPPnwithrr]

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