(* ::Package:: *)

(* ::Title:: *)
(*Simplifying Keum's FPP*)


(* ::Subsection:: *)
(*In this notebook, we will present the code to simplify the Keum fake projective plane. *)


(* ::Input:: *)
(*(*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[]*)


(* ::Input:: *)
(*(*Loading equations for Keum's FPP from Borisov21*)*)
(*Ps={P0,P1,P2,P3,P4,P5,P6,P7,P8,P9};*)


(* ::Input:: *)
(*fname="EqsFPPwithrrKeum.txt";*)
(*OpenRead[fname];*)
(*EqsFPPwithrrKeum=Read[fname];*)
(*Close[fname];*)


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


(* ::Input:: *)
(*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]];*)


(* ::Input:: *)
(*(*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}];*)
(*];*)
(**)


(* ::Input:: *)
(*(*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*)
(**)


(* ::Input:: *)
(*(*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}];*)


(* ::Input:: *)
(*(*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*)*)


(* ::Input:: *)
(*(*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]]*)


(* ::Input:: *)
(**)
(*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;*)


(* ::Input:: *)
(*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]}];*)
(**)


(* ::Input:: *)
(*(*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]*)


(* ::Input:: *)
(*(*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];*)


(* ::Subsection:: *)
(*Success! Now, we will run the beautification code*)
(**)


(* ::Subsection:: *)
(*Let's start with 0-graded C_7 part i.e. indexes 1-12. *)
(**)


(* ::Input:: *)
(*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;*)
(**)


(* ::Input:: *)
(*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];*)*)
(*];*)
(*];*)


(* ::Input:: *)
(*C3Aut={Q0->Q0,Q1->Q2, Q2->Q3,Q3->Q1,Q4->Q5,Q5->Q6,Q6->Q4,Q7->Q8,Q8->Q9,Q9->Q7};*)


(* ::Input:: *)
(*(*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]}]]*)
(**)


(* ::Subsection:: *)
(*Let's start with 1-graded C_7 part i.e. indexes 13-24. *)
(**)


(* ::Input:: *)
(*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;*)
(**)


(* ::Input:: *)
(*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];*)*)
(*];*)
(*];*)


(* ::Input:: *)
(*(*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]*)


(* ::Input:: *)
(* *)


(* ::Subsection:: *)
(*Let's start with 3-graded C_7 part i.e. indexes 37-48. *)
(**)


(* ::Input:: *)
(*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]*)
(**)


(* ::Input:: *)
(*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];*)*)
(*];*)
(*];*)


(* ::Input:: *)
(**)


(* ::Input:: *)
(*(*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]}]]*)
(**)


(* ::Input:: *)
(**)
(*(*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. *)*)


(* ::Input:: *)
(*(*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]*)
(**)


(* ::Input:: *)
(*Graded4Eqs=Graded1Eqs/.C3Aut;*)
(*Graded2Eqs=Graded4Eqs/.C3Aut;*)
(**)
(*Graded5Eqs=Graded3Eqs/.C3Aut;*)
(*Graded6Eqs=Graded5Eqs/.C3Aut;*)
(**)
(*Length[Graded2Eqs]*)
(*Length[Graded4Eqs]*)
(*Length[Graded5Eqs]*)
(*Length[Graded6Eqs]*)


(* ::Input:: *)
(*BeautifiedEqsFPP=Join[Graded0Eqs,Graded1Eqs, Graded2Eqs,Graded3Eqs,Graded4Eqs,Graded5Eqs,Graded6Eqs];*)
(*Length[BeautifiedEqsFPP]*)
(*ByteCount[BeautifiedEqsFPP]*)


(* ::Input:: *)
(*(*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]*)


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


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


(* ::Input:: *)
(*BeautifiedEqsFPPnwithrr=Numerator[Factor[BeautifiedEqsFPPwithrr]];*)
(*ByteCount[BeautifiedEqsFPPnwithrr]*)
(**)
(*fname="BeautifiedFPPwithrrn.txt";*)
(*OpenWrite[fname];*)
(*Write[fname, BeautifiedEqsFPPnwithrr];*)
(*Close[fname];*)


(* ::Input:: *)
(**)
