(* ::Package:: *)

(* ::Title:: *)
(*Map of Keum's FPP to P5 Q*)


(* ::Subsection:: *)
(*In this notebook, we will present the code to map Keum's FPP to P5 via |5H|*)


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


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


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


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


(* ::Subsection:: *)
(*Now we compute the order 2 neighborhood for p2*)


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


(* ::Subsection:: *)
(*We can now generate S3*)
(**)


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


(* ::Subsection:: *)
(* We will now compute random points on the FPP*)


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


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


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


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


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


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


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


(* ::Input:: *)
(*fname=OpenWrite["Coeffs Vanishing.txt"];*)
(*Write[fname, tbzrr];*)
(*Close[fname];*)


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


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


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


(* ::Subsection:: *)
(*We will now try to lift the solutions to the p-adics, and then to the complex numbers*)
(* *)


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


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


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


(* ::Input:: *)
(*(*The second point Point3mod37 worked! *)*)


(* ::Subsection:: *)
(*Let's plug in the solutions to S3, S5, S6, D3.*)
(**)


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


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


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


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


(* ::Input:: *)
(*(*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=\!\(\**)
(*TagBox[*)
(*RowBox[{"{", *)
(*RowBox[{*)
(*RowBox[{"Q0", "->", "P0"}], ",", *)
(*RowBox[{"Q1", "->", *)
(*FractionBox[*)
(*RowBox[{"4", " ", "I", " ", "P1"}], *)
(*RowBox[{*)
(*RowBox[{"-", "I"}], "+", *)
(*SqrtBox["7"]}]]}], ",", *)
(*RowBox[{"Q2", "->", *)
(*FractionBox[*)
(*RowBox[{"4", " ", "I", " ", "P2"}], *)
(*RowBox[{*)
(*RowBox[{"-", "I"}], "+", *)
(*SqrtBox["7"]}]]}], ",", *)
(*RowBox[{"Q3", "->", *)
(*FractionBox[*)
(*RowBox[{"4", " ", "I", " ", "P3"}], *)
(*RowBox[{*)
(*RowBox[{"-", "I"}], "+", *)
(*SqrtBox["7"]}]]}], ",", *)
(*RowBox[{"Q4", "->", *)
(*FractionBox[*)
(*RowBox[{"1856", " ", "I", " ", "P4"}], *)
(*RowBox[{*)
(*RowBox[{"15", " ", "I"}], "+", *)
(*SqrtBox["7"]}]]}], ",", *)
(*RowBox[{"Q5", "->", *)
(*FractionBox[*)
(*RowBox[{"1856", " ", "I", " ", "P5"}], *)
(*RowBox[{*)
(*RowBox[{"15", " ", "I"}], "+", *)
(*SqrtBox["7"]}]]}], ",", *)
(*RowBox[{"Q6", "->", *)
(*FractionBox[*)
(*RowBox[{"1856", " ", "I", " ", "P6"}], *)
(*RowBox[{*)
(*RowBox[{"15", " ", "I"}], "+", *)
(*SqrtBox["7"]}]]}], ",", *)
(*RowBox[{"Q7", "->", *)
(*RowBox[{"-", *)
(*FractionBox[*)
(*RowBox[{"464", " ", *)
(*RowBox[{"(", *)
(*RowBox[{*)
(*RowBox[{*)
(*RowBox[{"(", *)
(*RowBox[{*)
(*RowBox[{"67", " ", "I"}], "+", *)
(*RowBox[{"47", " ", *)
(*SqrtBox["7"]}]}], ")"}], " ", "P4"}], "+", *)
(*RowBox[{"4", " ", *)
(*RowBox[{"(", *)
(*RowBox[{*)
(*RowBox[{"15", " ", "I"}], "+", *)
(*SqrtBox["7"]}], ")"}], " ", "P7"}]}], ")"}]}], *)
(*RowBox[{"7", " ", *)
(*RowBox[{"(", *)
(*RowBox[{*)
(*RowBox[{*)
(*RowBox[{"-", "27"}], " ", "I"}], "+", *)
(*RowBox[{"4", " ", *)
(*SqrtBox["7"]}]}], ")"}]}]]}]}], ",", *)
(*RowBox[{"Q8", "->", *)
(*RowBox[{"-", *)
(*FractionBox[*)
(*RowBox[{"464", " ", *)
(*RowBox[{"(", *)
(*RowBox[{*)
(*RowBox[{*)
(*RowBox[{"(", *)
(*RowBox[{*)
(*RowBox[{"67", " ", "I"}], "+", *)
(*RowBox[{"47", " ", *)
(*SqrtBox["7"]}]}], ")"}], " ", "P5"}], "+", *)
(*RowBox[{"4", " ", *)
(*RowBox[{"(", *)
(*RowBox[{*)
(*RowBox[{"15", " ", "I"}], "+", *)
(*SqrtBox["7"]}], ")"}], " ", "P8"}]}], ")"}]}], *)
(*RowBox[{"7", " ", *)
(*RowBox[{"(", *)
(*RowBox[{*)
(*RowBox[{*)
(*RowBox[{"-", "27"}], " ", "I"}], "+", *)
(*RowBox[{"4", " ", *)
(*SqrtBox["7"]}]}], ")"}]}]]}]}], ",", *)
(*RowBox[{"Q9", "->", *)
(*RowBox[{"-", *)
(*FractionBox[*)
(*RowBox[{"464", " ", *)
(*RowBox[{"(", *)
(*RowBox[{*)
(*RowBox[{*)
(*RowBox[{"(", *)
(*RowBox[{*)
(*RowBox[{"67", " ", "I"}], "+", *)
(*RowBox[{"47", " ", *)
(*SqrtBox["7"]}]}], ")"}], " ", "P6"}], "+", *)
(*RowBox[{"4", " ", *)
(*RowBox[{"(", *)
(*RowBox[{*)
(*RowBox[{"15", " ", "I"}], "+", *)
(*SqrtBox["7"]}], ")"}], " ", "P9"}]}], ")"}]}], *)
(*RowBox[{"7", " ", *)
(*RowBox[{"(", *)
(*RowBox[{*)
(*RowBox[{*)
(*RowBox[{"-", "27"}], " ", "I"}], "+", *)
(*RowBox[{"4", " ", *)
(*SqrtBox["7"]}]}], ")"}]}]]}]}]}], "}"}],*)
(*Short[#, 5]& ]\);*)
(*QPointsOnDoubleFiber=QtoP/.PPointsOnDoubleFiber;*)
(**)


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


(* ::Input:: *)
(*fname="Deg3EqsDoubleFiberQ.txt";*)
(*OpenWrite[fname];*)
(*Write[fname, Deg3EqsDoubleFiber];*)
(*Close[fname];*)
(**)


(* ::Subsection:: *)
(**)
(*We now have all the ingredients necessary to compute |5H|*)


(* ::Input:: *)
(*DoubleFiber=Deg3EqsDoubleFiber;*)


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


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


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


(* ::Input:: *)
(*fname="ZsIn5H.txt";*)
(*OpenWrite[fname];*)
(*Write[fname, ZFunctions];*)
(*Close[fname];*)


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


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


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


(* ::Input:: *)
(*ZFunctionNew ={};*)
(*AppendTo[ZFunctionNew,ZFunctions[[1]]];*)


(* ::Subsubsection:: *)
(*Fixing second section*)


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


(* ::Input:: *)
(*FactorInteger[1927561216 ]*)


(* ::Input:: *)
(*denom=2*7;*)
(*MinimalPolynomial[FixedVal2*denom]*)


(* ::Input:: *)
(*FactorInteger[37256083456]*)


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


(* ::Input:: *)
(*minpolyarr2[[1]]*)
(*sols2[[1]]*)


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


(* ::Subsubsection:: *)
(*Fixing third section*)


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


(* ::Input:: *)
(*FactorInteger[729 ]*)


(* ::Input:: *)
(*denom=3;*)
(*MinimalPolynomial[FixedVal3*denom]*)


(* ::Input:: *)
(*FactorInteger[780448]*)


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


(* ::Input:: *)
(*minpolyarr3[[1]]*)
(*sols3[[1]]*)


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


(* ::Subsubsection:: *)
(*Fixing fourth section*)


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


(* ::Input:: *)
(*FactorInteger[194672 ]*)


(* ::Input:: *)
(*denom=2 23;*)
(*MinimalPolynomial[FixedVal4*denom]*)


(* ::Input:: *)
(*FactorInteger[34187517952]*)


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


(* ::Input:: *)
(*minpolyarr4[[1]]*)
(*sols4[[1]]*)


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


(* ::Subsubsection:: *)
(*Fixing fifth section*)


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


(* ::Input:: *)
(*FactorInteger[5488 ]*)


(* ::Input:: *)
(*denom=2 7;*)
(*MinimalPolynomial[FixedVal5*denom]*)


(* ::Input:: *)
(*FactorInteger[524386566144]*)


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


(* ::Input:: *)
(*minpolyarr5[[1]]*)
(*sols5[[1]]*)


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


(* ::Subsubsection:: *)
(*Fixing sixth section*)


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


(* ::Input:: *)
(*FactorInteger[2916 ]*)


(* ::Input:: *)
(*denom=2 3;*)
(*MinimalPolynomial[FixedVal6*denom]*)


(* ::Input:: *)
(*FactorInteger[232202043392]*)


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


(* ::Input:: *)
(*minpolyarr6[[1]]*)
(*sols6[[1]]*)


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


(* ::Subsubsection:: *)
(*We're done fixing the group action.  Let's see what the new values look like. *)


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


(* ::Subsection:: *)
(*Let's compute the 56 degree 6 equations*)
(**)


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


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


(* ::Subsubsection:: *)
(*Let's check the group action*)


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


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


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


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


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


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


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


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


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


(* ::Input:: *)
(*fname=OpenWrite["SixSections.txt"]*)
(*Write[fname,SixSections];*)
(*Close[fname];*)


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


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


(* ::Input:: *)
(*fname="EqsFPPInP5WithZCoords.txt";*)
(*OpenWrite[fname];*)
(*Write[fname, EqsFPPInP5];*)
(*Close[fname];*)
(**)


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


(* ::Input:: *)
(*fname="EqsFPPInP5WithZCoordsrr.txt";*)
(*OpenWrite[fname];*)
(*Write[fname, EqsFPPInP5withrr];*)
(*Close[fname];*)


(* ::Input:: *)
(*EqsFPPInP5withrrn=Numerator[Factor[EqsFPPInP5withrr]];*)
(*fname="EqsFPPInP5WithZCoordsrrn.txt";*)
(*OpenWrite[fname];*)
(*Write[fname, EqsFPPInP5withrrn];*)
(*Close[fname];*)


(* ::Input:: *)
(*(*Recomputing the action*)*)
(*C3Action={Z1->Z2,Z2->(1+3Sqrt[-7])/8 Z4,Z4->Z1,Z3->Z6,Z6->(1+3Sqrt[-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]*)


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


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


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


(* ::Input:: *)
(*BeautifiedEqsFPPInP5=Flatten[BeautifiedEqsFPPInP5];*)


(* ::Input:: *)
(*fname="BeautifiedEqsFPPInP5WithZCoords.txt";*)
(*OpenWrite[fname];*)
(*Write[fname, BeautifiedEqsFPPInP5];*)
(*Close[fname];*)
(**)


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


(* ::Input:: *)
(*fname="BeautifiedEqsFPPInP5WithZCoordsrr.txt";*)
(*OpenWrite[fname];*)
(*Write[fname, BeautifiedEqsFPPInP5withrr];*)
(*Close[fname];*)


(* ::Input:: *)
(*BeautifiedEqsFPPInP5withrrn=Numerator[Factor[BeautifiedEqsFPPInP5withrr]];*)
(*fname="BeautifiedEqsFPPInP5WithZCoordsrrn.txt";*)
(*OpenWrite[fname];*)
(*Write[fname, BeautifiedEqsFPPInP5withrrn];*)
(*Close[fname];*)
