In this notebook we begin to construct a 8:1 cover of the FPP, trying to keep the basis reasonably nice. We will save all of the relevant results so that subsequent steps do not need to use them.

SetDirectory["/Users/borisov/Desktop/2025/FPPOct2025/"] (* This needs to be replaced with whatever directory you store the original equations and the results in.*)

(* reading in various data computed by Step23.nb *)

fname=OpenRead["EqsRsmallWithRR"];
EqsRsmall=Read[fname]/.{rr1->Sqrt[-2],rr2->Sqrt[-3]};
Close[fname];

fname=OpenRead["UtoR"];
UtoR=Read[fname];
Close[fname];

fname=OpenRead["C3actionU"];
C3actionU=Read[fname];
Close[fname];

fname=OpenRead["C3actionU2"];
C3actionU2=Read[fname];
Close[fname];

fname=OpenRead["pointsC2C2cover"];
pointsC2C2cover=Read[fname];
Close[fname];

fname=OpenRead["pointsFPP"];
pointsFPP=Read[fname];
Close[fname];

Repeating the same standard simplification commands as in Steps23.nb.

NormalizedTrace[x_]:=AlgebraicNumberTrace[x]/Exponent[MinimalPolynomial[x,t],t];
To4[x_]:=NormalizedTrace[{x,x/Sqrt[-2],x/Sqrt[-3],x/Sqrt[6]}];
ToRR[x_]:=To4[x].{1,rr1,rr2,-rr1 rr2};

SimplifyEquationsRootM2M3[eqs_]:=Module[{SmallLinComb,coeffmatrix,eqs4,NiceEquations,mons,ml,vars,varsto1,k,i,t,rk,rkfull,MatrixLarge,newrow,BackToCxMatrix,BackToCxReduced,m},

vars=Variables[eqs];
varsto1=Table[vars[[k]]->1,{k,1,Length[vars]}];
(* find monomials *)
mons={};
For[k=1,k<=Length[eqs],k++,
ml=MonomialList[eqs[[k]]];
ml=Table[ml[[i]]/(ml[[i]]/.varsto1),{i,1,Length[ml]}];
mons=Union[mons,ml];
];

(* getting the matrix *)
coeffmatrix=Table[Coefficient[eqs[[k]],mons],{k,1,Length[eqs]}];

eqs4=Join[Join[coeffmatrix,Sqrt[-2] coeffmatrix],Join[Sqrt[-3] coeffmatrix,Sqrt[6] coeffmatrix]];
MatrixLarge={};
For[k=1,k<=Length[eqs4],k++,
newrow={};
For[m=1,m<=Length[eqs4[[k]]],m++,
newrow=Join[newrow,To4[eqs4[[k]][[m]]]];
];
AppendTo[MatrixLarge,newrow];
];

Print["computing rank"];
rkfull =MatrixRank[MatrixLarge]/4;

(* Getting rid of denominators.*)
Print["getting rid of denominators"];
MatrixLarge=Table[(Denominator[MatrixLarge[[k]]]/.
List->LCM) MatrixLarge[[k]],{k,1,Length[MatrixLarge]}];

(*Saturating the lattice using Smith decomposition.*)
Print["starting Smith decomposition"]; 
{u,r,v}=SmithDecomposition[MatrixLarge];
Print["Smith decomposition is done"];
SmallLinComb=Table[LinearSolve[Transpose[v],Table[KroneckerDelta[k,i],{k,1,Length[MatrixLarge[[1]]]}]],{i,1,MatrixRank[MatrixLarge]}];

(*Key step:using the lattice reduction algorithm.*)
Print["starting lattice reduction"]; 
SmallLinComb=LatticeReduce[SmallLinComb];
Print["lattice reduction is done"];

(*Getting back to the complex matrix.*)
BackToCxMatrix=Table[Table[SmallLinComb[[i]][[4 k-3]]+Sqrt[-2]  SmallLinComb[[i]][[4k-2]]+Sqrt[-3] SmallLinComb[[i]][[4 k-1]]+Sqrt[6]  SmallLinComb[[i]][[4k]],{k,1,Length[SmallLinComb[[1]]]/4}],{i,1,Length[SmallLinComb]}];

(*Picking the first rows that give the correct rank.*)
BackToCxReduced={};
rk=0;
For[i=1,i<=Length[BackToCxMatrix],i++,
If[MatrixRank[N[Append[BackToCxReduced,BackToCxMatrix[[i]]],800]]>rk,rk++;
BackToCxReduced=Append[BackToCxReduced,BackToCxMatrix[[i]]];
(*Print[i," ",rk];*)
If[rk==rkfull,i=Length[BackToCxMatrix]+1];
];
];

NiceEquations=BackToCxReduced.mons;
NiceEquations
];

EquationsTorrM2M3[eqs_,rr1_,rr2_]:=Module[{k,t,i,vars,varsto1,ml,eqsrr},

vars=Variables[eqs];
varsto1=Table[vars[[k]]->1,{k,1,Length[vars]}];

eqsrr={};
For[k=1,k<=Length[eqs],k++,
ml=MonomialList[eqs[[k]]];
ml=Table[ml[[i]]/(ml[[i]]/.varsto1),{i,1,Length[ml]}];
AppendTo[eqsrr,Sum[ToRR[Coefficient[eqs[[k]],ml[[i]]]]ml[[i]],{i,1,Length[ml]}]]
];

(* output *)
eqsrr
];

There should be one, up to scaling, section of 3H on the Q8 cover of FPP which is odd for the involution to the C2C2 cover, weights (1,0) for C3xC3; and similar for (2,0). Their product f0 and squares f1, f2 will be linear combinations of Us.

UnknownCoeffs=Table[ToExpression[StringJoin["c",ToString[i]]],{i,1,100}];
Us=Table[ToExpression[StringJoin["U",ToString[i]]],{i,0,39}];
w=Exp[2Pi I/3];

(* f0 has weight 0 for both C3 actions. It is therefore a linear combination of the following.*)
f0 = {c1,c2,c3,c4,c5,c6}.{U0, U1, U11+U21+U31, U12+U22+U32, U13+U23+U33,U10+U20+U30};
f1={c7,c8,c9,c10,c11}.{U7,U11+w U21+w^2U31, U12+w U22+w^2U32, U13+w U23+w^2U33,U10+w U20+w^2U30};
f2={c12,c13,c14,c15,c16}.{U4,U11+w^2 U21+w U31, U12+w^2 U22+w U32, U13+w^2 U23+ w U33,U10+w^2 U20+w U30};

tbz=f0^2 - f1 f2;
Mon2C=DeleteDuplicates[Flatten[Table[UnknownCoeffs[[i]]UnknownCoeffs[[j]],{i,1,16},{j,1,16}]]];
tbzvalues=Table[Coefficient[tbz,Mon2C[[k]]]/.pointsC2C2cover,{k,1,Length[Mon2C]}];
rels=RowReduce[Transpose[tbzvalues]];

Accuracy[rels]
relsexact=Complement[RootApproximant[Chop[rels,10^(-700)],4].Mon2C,{0}];
simplerelsexact = SimplifyEquationsRootM2M3[relsexact];

One needs to be careful so as not to pick up trivial solutions that come from the even sections.

extra={c8==1,c13==1};
sols=NSolve[Join[extra,Table[0==simplerelsexact[[k]],{k,1,Length[simplerelsexact]}]],800];
{Length[sols],Accuracy[sols]}
solsexact=Table[RootApproximant[Chop[UnknownCoeffs[[k]]/.sols[[1]],10^(-400)]],{k,1,16}];
F1=EquationsTorrM2M3[(-1573 I (Sqrt[3])+11 (Sqrt[6])){f1}/.Table[UnknownCoeffs[[k]]->solsexact[[k]],{k,1,16}],rr1,rr2][[1]]/.{rr1->Sqrt[-2],rr2->Sqrt[-3]} (* scaled for future needs *)
F0=EquationsTorrM2M3[{RootReduce[Collect[(f0/solsexact[[1]])/.Table[UnknownCoeffs[[k]]->solsexact[[k]],{k,1,16}],Us]]},rr1,rr2][[1]]/.{rr1->Sqrt[-2],rr2->Sqrt[-3]}

We now want to compute two more sections of H^0(3H) on the 8:1 cover. Specifically, these will be odd for the covering involution and C3-fixed. If we have one such v, then f0 = v Sqrt(f1) and f2 = v^2 satisfy  f0^2 = f1 f2. We will have f0 of the weight (2,0) and f2 of the weight (0,0).

f0 = {c1,c2,c3,c4,c5}.{U4,U11+w^2 U21+w U31, U12+w^2 U22+w U32, U13+w^2 U23+ w U33,U10+w^2 U20+w U30};
f1=c12 F1;
f2={c6,c7,c8,c9,c10,c11}.{U0, U1, U11+U21+U31, U12+U22+U32, U13+U23+U33,U10+U20+U30};

tbz=f0^2 - f1 f2;
Mon2C=DeleteDuplicates[Flatten[Table[UnknownCoeffs[[i]]UnknownCoeffs[[j]],{i,1,12},{j,1,12}]]];
tbzvalues=Table[Coefficient[tbz,Mon2C[[k]]]/.pointsC2C2cover,{k,1,Length[Mon2C]}];
rels=RowReduce[Transpose[tbzvalues]];

Accuracy[rels]
relsexact=Complement[RootApproximant[Chop[rels,10^(-700)],4].Mon2C,{0}];
simplerelsexact = SimplifyEquationsRootM2M3[relsexact];

extra={c12==1,c3==1,c2==0};
sols=NSolve[Join[extra,Table[0==simplerelsexact[[k]],{k,1,Length[simplerelsexact]}]],800];
{Length[sols],Accuracy[sols]}
solsexact=Table[RootApproximant[Chop[UnknownCoeffs[[k]]/.sols[[1]],10^(-400)]],{k,1,12}];
F2=EquationsTorrM2M3[{f0}/.Table[UnknownCoeffs[[k]]->solsexact[[k]],{k,1,12}],rr1,rr2][[1]]/.{rr1->Sqrt[-2],rr2->Sqrt[-3]}

extra={c12==1,c3==0,c2==1};
sols=NSolve[Join[extra,Table[0==simplerelsexact[[k]],{k,1,Length[simplerelsexact]}]],800];
{Length[sols],Accuracy[sols]}
solsexact=Table[RootApproximant[Chop[UnknownCoeffs[[k]]/.sols[[1]],10^(-400)]],{k,1,12}];
F3=EquationsTorrM2M3[{f0}/.Table[UnknownCoeffs[[k]]->solsexact[[k]],{k,1,12}],rr1,rr2][[1]]/.{rr1->Sqrt[-2],rr2->Sqrt[-3]}

We now have seven sections of H^0(8:1, 3H). Specifically, these are Sqrt(CutR), Sqrt(Cut1), Sqrt(Cut2), F_i/Sqrt(F_1), for i=0,1,2,3. We find equations on these sections.

Rs={R0,R1,R2,R3,R4,R5,R6,R7,R8,R9};
CutR =R1+R4+R7;
wtsC3C3={{0,0},{0,0},{0,1},{0,2},{1,0},{1,1},{1,2},{2,0},{2,1},{2,2}};
C3action=Table[Rs[[k]]->w^wtsC3C3[[k]][[1]]Rs[[k]],{k,1,10}]
C3action2=Table[Rs[[k]]->w^wtsC3C3[[k]][[2]]Rs[[k]],{k,1,10}]
CutR=R1+R4+R7; 
Cut1= CutR/.C3action;
Cut2=Cut1/.C3action;

VtoR={Sqrt[CutR],Sqrt[Cut1],Sqrt[Cut2],F0/Sqrt[F1],F1/Sqrt[F1],F2/Sqrt[F1],F3/Sqrt[F1]};

valuespoints8={};
For[k=1,k<=Length[pointsC2C2cover],k++,
AppendTo[valuespoints8, {Sqrt[CutR],Sqrt[Cut1],Sqrt[Cut2],F0/Sqrt[F1],F1/Sqrt[F1],F2/Sqrt[F1],F3/Sqrt[F1]}/.Join[pointsFPP[[k]],pointsC2C2cover[[k]]]];
];
Accuracy[valuespoints8]
Vs={V0,V1,V2,V3,V4,V5,V6};

eqs8={};
eqs8rr={};
For[deg=3,deg<=4,deg++,
Print["deg=",deg];
ml=MonomialList[Sum[Vs[[k]],{k,1,7}]^deg];
ml=Table[ml[[k]]/(ml[[k]]/.Table[Vs[[m]]->1,{m,1,7}]),{k,1,Length[ml]}];
mleven={};mlodd={};
For[m=1,m<=Length[ml],m++,
If[0==Mod[Sum[Exponent[ml[[m]],Vs[[i]]],{i,4,7}],2],
AppendTo[mleven,ml[[m]]],
AppendTo[mlodd,ml[[m]]]
]
];

ns=NullSpace[Table[mleven[[i]]/.Table[Vs[[k]]->valuespoints8[[m]][[k]],{k,1,7}],{m,1,Length[valuespoints8]},{i,1,Length[mleven]}]];
If[Length[ns]==0,Print["not at this degree, even"],
ns=RowReduce[Chop[ns,10^(-700)]];
Print["row reduction done"];
Print[{Accuracy[ns],Length[ns]}];
eqs8=Join[eqs8,SimplifyEquationsRootM2M3[RootApproximant[Chop[ns,10^(-600)],4].mleven]];
eqs8rr=Join[eqs8rr,EquationsTorrM2M3[eqs8,rr1,rr2]];
];

ns=NullSpace[Table[mlodd[[i]]/.Table[Vs[[k]]->valuespoints8[[m]][[k]],{k,1,7}],{m,1,Length[valuespoints8]},{i,1,Length[mlodd]}]];
If[Length[ns]==0,Print["not at this degree, odd"],
ns=RowReduce[Chop[ns,10^(-700)]];
Print["row reduction done"];
Print[{Accuracy[ns],Length[ns]}];
eqs8=Join[eqs8,SimplifyEquationsRootM2M3[RootApproximant[Chop[ns,10^(-600)],4].mlodd]];
eqs8rr=Join[eqs8rr,EquationsTorrM2M3[eqs8,rr1,rr2]];
];

];

(* verifying: should get 24 *)
sols=NSolve[Join[Table[0==eqs8[[k]],{k,1,Length[eqs8]}],{1==V0,0==Sum[Vs[[k]] RandomInteger[1000],{k,1,7}],0==Sum[Vs[[k]] RandomInteger[1000],{k,1,7}]}],500];
Length[sols] (* get 24, as expected, this is the quotient of 8:1 cover of FPP by a C3 *)

It is possible to reduce the number of degree 4 equations, because some of them are just degree 3 times linear, but there is no clear advantage to doing it at this point, as we would still have at least 30 degree 4 equations left.

C3actionV={V0->V1,V1->V2,V2->V0,V3-> w^2 V3, V4-> w V4, V5->V5, V6->V6};
Length[eqs8]
eqs8deg3=Table[eqs8[[k]],{k,1,4}];
eqs8deg4=Table[eqs8[[k]],{k,5,62}];

MonV3=Flatten[Table[Vs[[i]]Vs[[j]]Vs[[k]],{i,1,7},{j,1,7},{k,1,7}]];
MonV4=Flatten[Table[Vs[[i]]Vs[[j]]Vs[[k]]Vs[[l]],{i,1,7},{j,1,7},{k,1,7},{l,1,7}]];

eqs=Join[eqs8deg3,eqs8deg3/.C3actionV];
{Length[eqs8deg3],MatrixRank[Table[Coefficient[eqs[[i]],MonV3[[l]]],{i,1,Length[eqs]},{l,1,Length[MonV3]}]]}
eqs=Join[eqs8deg4,eqs8deg4/.C3actionV];
{Length[eqs8deg4],MatrixRank[N[Table[Coefficient[eqs[[i]],MonV4[[l]]],{i,1,Length[eqs]},{l,1,Length[MonV4]}],300]]}

We would like to understand the action of the group on the space above. Specifically, we would like to know the action of the order four element of the octonion  group that negates V1, V2 and keeps V0. The action on V3,...,V6 should have two eigenvalues i and two eigenvalues -i.

C3actionV={V0->V1,V1->V2,V2->V0,V3-> w^2 V3, V4-> w V4, V5->V5, V6->V6};
NewV3=c1 V3 + c2 V4 + c5 V5 + c6 V6;
EigenBasisOdd={NewV3,NewV3/.C3actionV,(NewV3/.C3actionV)/.C3actionV,c7 V5 + c8 V6};
MatrixBasisOdd=Table[Coefficient[EigenBasisOdd,Vs[[k]]],{k,4,7}];
MatrixBasisOdd//MatrixForm

C4actionVansatz={V0->V0,V1->-V1,V2->-V2,
V3->c1 V3 + c2 V4 + c3 V5 + c4 V6,
V4->c5 V3 + c6 V4 + c7 V5 + c8 V6,
V5->c9 V3 + c10 V4 + c11 V5 + c12 V6,
V6->c13 V3 + c14 V4 + c15 V5 + c16 V6
};
C4actionEqs={c1+c6+c11+c16==0};
For[k=4,k<=7,k++,
tbz=((Vs[[k]]/.C4actionVansatz)/.C4actionVansatz)+Vs[[k]];
For[l=4,l<=7,l++,AppendTo[C4actionEqs,0==Coefficient[tbz,Vs[[l]]]]]
]

(* from the character table it looks like the image of V3 should have zero coefficient by V4 and vice versa *)
AppendTo[C4actionEqs,0==c2];
AppendTo[C4actionEqs,0==c5];

sols=Solve[C4actionEqs]; (* we will just do these first *)
sols1=FullSimplify[sols[[1]]] (* there are other solutions, but this is the main branch *)
C4actionVansatz=C4actionVansatz/.sols1;

The idea now is to find points with coordinates (1:0:0:...). Then the action should be permuting these points, and we would be able to find the action.

sols=NSolve[Join[Table[0==eqs8[[k]],{k,1,4}],{1==V0,0==V1,0==V2}],500];

sols1={};
For[k=1,k<=Length[sols],k++,
For[l=1,l<=62,l++,
If[Chop[eqs[[l]]/.sols[[k]],10^(-300)]!=0,
l=100;]
];
If[l<100,AppendTo[sols1,sols[[k]]]];
];
Length[sols1]

For[k=2,k<=Length[sols1],k++,
For[l=3,l<=Length[sols1],l++,
tbz=Join[ ((Vs/.C4actionVansatz)/.sols1[[1]]) - (Vs/.sols1[[k]]),((Vs/.C4actionVansatz)/.sols1[[3]]) - (Vs/.sols1[[l]])];
maybe=
NSolve[Table[tbz[[m]]==0,{m,1,Length[tbz]}]];
If[Length[maybe]!=0,Print[{k,l}]];
];];

Let us pick one of these and see if we get an automorphism. We just run through all of the above pairs until one clicks.

tbz=Join[ ((Vs/.C4actionVansatz)/.sols1[[1]]) - (Vs/.sols1[[5]]),((Vs/.C4actionVansatz)/.sols1[[3]]) - (Vs/.sols1[[22]])];
maybe=
NSolve[Table[tbz[[m]]==0,{m,1,Length[tbz]}]];
valuesc=RootApproximant[{c1,c11,c16,c3,c4,c7}/.maybe[[1]]];
C4actionV=FullSimplify[ C4actionVansatz/.Table[{c1,c11,c16,c3,c4,c7}[[k]]->valuesc[[k]],{k,1,6}]]
N[Table[eqs8[[l]]/.C4actionV/.Table[Vs[[k]]->valuespoints8[[33]][[k]],{k,1,7}],{l,1,1}]] 
(* this one is wrong! *)

tbz=Join[ ((Vs/.C4actionVansatz)/.sols1[[1]]) - (Vs/.sols1[[24]]),((Vs/.C4actionVansatz)/.sols1[[3]]) - (Vs/.sols1[[13]])];
maybe=
NSolve[Table[tbz[[m]]==0,{m,1,Length[tbz]}]];
valuesc=Table[ToRR[RootApproximant[{c1,c11,c16,c3,c4,c7}[[k]]/.maybe[[1]]]]/.{rr1->Sqrt[-2],rr2->Sqrt[-3]},{k,1,6}];
C4actionV=Collect[ C4actionVansatz/.Table[{c1,c11,c16,c3,c4,c7}[[k]]->valuesc[[k]],{k,1,6}],Vs];
C4actionV=Table[Vs[[l]]->Sum[Vs[[k]]ToRR[Coefficient[Vs[[l]]/.C4actionV,Vs[[k]]]]/.{rr1->Sqrt[-2],rr2->Sqrt[-3]},{k,1,7}],{l,1,7}]
DeleteDuplicates[Flatten[N[Table[eqs8[[l]]/.C4actionV/.Table[Vs[[k]]->valuespoints8[[m]][[k]],{k,1,7}],{l,1,62},{m,10,20}]]]] (* this one is good! *)

We will now compute the remaining 40 basis elements of H^0(8FPP, 6H). They will be split according to the action of the  C3 x C3 group. Specifically, they will be given as quadratic equations in U1,...,U40, which are zero on the V0 + V1 + V2 = 0 and V3 = 0. To make this work, we need to construct multiple points  on these curves. 

Prior to this, we want to express all 40 Us and 10 Rs in terms of Vs. Of course, we need to take one cubic root to do that.

(* get a database of points *)
PointsV={};
For[num=1,num<=20,num++,
blah=Join[{1==V0,0==Sum[Vs[[k]]RandomInteger[1000],{k,1,7}],0==Sum[Vs[[k]]RandomInteger[1000],{k,1,7}]},Table[0==eqs8[[k]],{k,1,6}]];
(* don't need to use all 62 equations *)
sols=NSolve[blah,800,WorkingPrecision->100];
If[Length[sols]==24,
For[m=1,m<=5,m++,
AppendTo[PointsV,sols[[m]]];
]
];
];
{Length[PointsV],Accuracy[PointsV]}

UtoRV=Table[Us[[k]]->Rs[[k]],{k,1,10}];
For[k=11,k<=20,k++,
AppendTo[UtoRV,Us[[k]]->(Us[[k]]/.UtoR)((Sqrt[R1+(E^((2 I π)/3)) R4+(E^(-((2 I π)/3))) R7]) (Sqrt[R1+(E^(-((2 I π)/3))) R4+(E^((2 I π)/3)) R7]))/V1/V2]
]
For[k=21,k<=30,k++,
AppendTo[UtoRV,Us[[k]]->(Us[[k]]/.UtoR)((Sqrt[R1+R4+R7]) (Sqrt[R1+(E^(-((2 I π)/3))) R4+(E^((2 I π)/3)) R7]))/V0/V2]
]
For[k=31,k<=40,k++,
AppendTo[UtoRV,Us[[k]]->(Us[[k]]/.UtoR)((Sqrt[R1+R4+R7]) (Sqrt[R1+(E^((2 I π)/3)) R4+(E^(-((2 I π)/3))) R7]))/V1/V0]
]

PointsVRU={};
For[numpt=1,numpt<=Length[PointsV], numpt++,
If[Mod[numpt,10]==1,Print[{numpt,Length[PointsVRU],Accuracy[PointsVRU]}]];
solsR1R4R7=Solve[{(V0/.PointsV[[numpt]])^2 == CutR,(V1/.PointsV[[numpt]])^2 == Cut1,(V2/.PointsV[[numpt]])^2 == Cut2}][[1]];
sols=NSolve[Append[Table[0==EqsRsmall[[k]],{k,1,20}],0==(V3 V4-(F0/.UtoRV))]/.Join[PointsV[[numpt]],solsR1R4R7]]; (* it is faster to not use all of the equations *)
If[Length[sols]==3, (* sanity check *)
For[k=1,k<=3,k++,
ptVR = Join[PointsV[[numpt]],sols[[k]],solsR1R4R7];
AppendTo[PointsVRU, Join[ptVR,UtoRV/.ptVR]];
]
]
];

Print[{numpt,Length[PointsVRU],Accuracy[PointsVRU]}];

It will be better for us to have the U_10, ..., U_39 in a C3-invariant basis. We will call them Utildes.

Utildes= Table[ToExpression[StringJoin["Utilde",ToString[k]]],{k,0,39}];
UtildesToU={};
For[k=1,k<=10,k++,
AppendTo[UtildesToU,Utildes[[k]]->Us[[k]]];
AppendTo[UtildesToU,Utildes[[10+k]]->Us[[10+k]]+Us[[20+k]]+Us[[30+k]]];
AppendTo[UtildesToU,Utildes[[20+k]]->Us[[10+k]]+w Us[[20+k]]+w^2 Us[[30+k]]];
AppendTo[UtildesToU,Utildes[[30+k]]->Us[[10+k]]+w^2 Us[[20+k]]+w Us[[30+k]]];
]

I will now write R0, R1, R4, R7 in terms of V-s. Then I will do it for R2^3, and express other R-s in terms of R2.

(* R1, R4, R7 are easy *)
R147toV=RootReduce[FullSimplify[Solve[{CutR==V0^2, Cut1==V1^2, Cut2==V2^2},{R1,R4,R7}]]][[1]];

(* R0 also turns out to be a degree 2 polynomial in V-s *)
Vto1=Table[Vs[[l]]->1,{l,1,7}];
Mons2V=MonomialList[Sum[Vs[[k]],{k,1,7}]^2];
Mons2V=Table[Mons2V[[k]]/(Mons2V[[k]]/.Vto1),{k,1,Length[Mons2V]}];
goodMons2V={};
For[k=1,k<=Length[Mons2V],k++,
If[Mod[Sum[Exponent[Mons2V[[k]],Vs[[l]]],{l,4,7}],2]==0,
AppendTo[goodMons2V,Mons2V[[k]]];
];
];
Mons2V=goodMons2V;


blah=Join[{R0}, Mons2V];
Length[blah]
values=blah/.PointsVRU;
ns=NullSpace[values];
mxns=ReverseSortBy[ns[[1]],Abs][[1]];
ns=ns/mxns;
acc=Round[Accuracy[ns]];
sols=Solve[0==RootApproximant[Chop[ns[[1]],10^(-acc+50)],4].blah,R0];
R0147toV=Append[R147toV,sols[[1]]];
R0147toV[[4]]=R0->(R0/.RootReduce[Collect[R0147toV[[4]],Vs]]);

I also want to write R2^3 as a rational function in V-s. The approach is to write it as a rational function in R0, R1, R4, R7. I will also do it for R2 R3, R2 R6, R2 R9, R2^2 R5, R2^2 R8 so that we are then able to write these in terms of R2.

PointsR={};
For[k=1,k<=30,k++,
sols=NSolve[Join[{1==R0,0==Sum[RandomInteger[200]Rs[[l]],{l,1,10}],0==Sum[RandomInteger[200]Rs[[l]],{l,1,10}]},
Table[0==EqsRsmall[[l]],{l,1,84}]],600];
For[l=1,l<=10,l++,
If[Accuracy[sols[[l]]]>350, AppendTo[PointsR,sols[[l]]]]];
];
{Length[PointsR],Accuracy[PointsR]}

(* finding a formula for R2^3 *)
deg=9;
sols=Solve[{a0>=0,a1>=0,a4>=0,a7>=0,a0+a1+a4+a7==deg,Mod[a4-a7,3]==0},Integers];
Mons1=R0^a0 R1^a1 R4^a4 R7^a7/.sols;
deg+=3;
sols=Solve[{a0>=0,a1>=0,a4>=0,a7>=0,a0+a1+a4+a7==deg,Mod[a4-a7,3]==0},Integers];
Mons2=R0^a0 R1^a1 R4^a4 R7^a7/.sols;
maybe=Join[R2^3 Mons1, Mons2];
Length[maybe]
values=Table[maybe[[k]]/.PointsR,{k,1,Length[maybe]}];
ns=NullSpace[Transpose[values]];
Length[ns]
rr=RowReduce[ns,Tolerance->10^(-400)];
Accuracy[rr]
eqs=RootApproximant[Chop[rr,10^(-400)],4].maybe;
eqs=SimplifyEquationsRootM2M3[RootReduce[eqs]]; 
R2cubeToR=-(eqs[[1]]/.R2->0)/Coefficient[eqs[[1]],R2^3];

(* I will now do it for R2 R3 *)
deg=9;
sols=Solve[{a0>=0,a1>=0,a4>=0,a7>=0,a0+a1+a4+a7==deg,Mod[a4-a7,3]==0},Integers];
Mons1=R0^a0 R1^a1 R4^a4 R7^a7/.sols;
deg+=2;
sols=Solve[{a0>=0,a1>=0,a4>=0,a7>=0,a0+a1+a4+a7==deg,Mod[a4-a7,3]==0},Integers];
Mons2=R0^a0 R1^a1 R4^a4 R7^a7/.sols;
maybe=Join[R2 R3  Mons1, Mons2];
Length[maybe]
values=Table[maybe[[k]]/.PointsR,{k,1,Length[maybe]}];
ns=NullSpace[Transpose[values]];
Length[ns]
rr=RowReduce[ns,Tolerance->10^(-400)];
Accuracy[rr]
eqs=RootApproximant[Chop[rr,10^(-400)],4].maybe;
eqs=SimplifyEquationsRootM2M3[RootReduce[eqs]]; 
R2R3toR=-(eqs[[1]]/.R2->0)/Coefficient[eqs[[1]],R2 R3];

(* I will now do it for R2 R6 *)
deg=9;
sols=Solve[{a0>=0,a1>=0,a4>=0,a7>=0,a0+a1+a4+a7==deg,Mod[a4-a7,3]==0},Integers];
Mons1=R0^a0 R1^a1 R4^a4 R7^a7/.sols;
deg+=2;
sols=Solve[{a0>=0,a1>=0,a4>=0,a7>=0,a0+a1+a4+a7==deg,Mod[a4-a7,3]==1},Integers]; (* needs to accomodate the second grading *)
Mons2=R0^a0 R1^a1 R4^a4 R7^a7/.sols;
maybe=Join[R2 R6  Mons1, Mons2];
Length[maybe]
values=Table[maybe[[k]]/.PointsR,{k,1,Length[maybe]}];
ns=NullSpace[Transpose[values]];
Length[ns]
rr=RowReduce[ns,Tolerance->10^(-400)];
Accuracy[rr]
eqs=RootApproximant[Chop[rr,10^(-400)],4].maybe;
eqs=SimplifyEquationsRootM2M3[RootReduce[eqs]]; 
R2R6toR=-(eqs[[1]]/.R2->0)/Coefficient[eqs[[1]],R2 R6];

(* I will now do it for R2 R9 *)
deg=9;
sols=Solve[{a0>=0,a1>=0,a4>=0,a7>=0,a0+a1+a4+a7==deg,Mod[a4-a7,3]==0},Integers];
Mons1=R0^a0 R1^a1 R4^a4 R7^a7/.sols;
deg+=2;
sols=Solve[{a0>=0,a1>=0,a4>=0,a7>=0,a0+a1+a4+a7==deg,Mod[a4-a7,3]==2},Integers]; (* needs to accomodate the second grading *)
Mons2=R0^a0 R1^a1 R4^a4 R7^a7/.sols;
Mons2=R0^a0 R1^a1 R4^a4 R7^a7/.sols;
maybe=Join[R2 R9  Mons1, Mons2];
Length[maybe]
values=Table[maybe[[k]]/.PointsR,{k,1,Length[maybe]}];
ns=NullSpace[Transpose[values]];
Length[ns]
rr=RowReduce[ns,Tolerance->10^(-400)];
Accuracy[rr]
eqs=RootApproximant[Chop[rr,10^(-400)],4].maybe;
eqs=SimplifyEquationsRootM2M3[RootReduce[eqs]]; 
R2R9toR=-(eqs[[1]]/.R2->0)/Coefficient[eqs[[1]],R2 R9];

(* I will now do it for R2^2 R5 *)
deg=9;
sols=Solve[{a0>=0,a1>=0,a4>=0,a7>=0,a0+a1+a4+a7==deg,Mod[a4-a7,3]==0},Integers];
Mons1=R0^a0 R1^a1 R4^a4 R7^a7/.sols;
deg+=3;
sols=Solve[{a0>=0,a1>=0,a4>=0,a7>=0,a0+a1+a4+a7==deg,Mod[a4-a7,3]==1},Integers]; (* needs to accomodate the second grading *)
Mons2=R0^a0 R1^a1 R4^a4 R7^a7/.sols;
maybe=Join[R2^2 R5  Mons1, Mons2];
Length[maybe]
values=Table[maybe[[k]]/.PointsR,{k,1,Length[maybe]}];
ns=NullSpace[Transpose[values]];
Length[ns]
rr=RowReduce[ns,Tolerance->10^(-400)];
Accuracy[rr]
eqs=RootApproximant[Chop[rr,10^(-400)],4].maybe;
eqs=SimplifyEquationsRootM2M3[RootReduce[eqs]]; 
R22R5toR=-(eqs[[1]]/.R2->0)/Coefficient[eqs[[1]],R2^2 R5];

(* I will now do it for R2^2 R8 *)
deg=9;
sols=Solve[{a0>=0,a1>=0,a4>=0,a7>=0,a0+a1+a4+a7==deg,Mod[a4-a7,3]==0},Integers];
Mons1=R0^a0 R1^a1 R4^a4 R7^a7/.sols;
deg+=3;
sols=Solve[{a0>=0,a1>=0,a4>=0,a7>=0,a0+a1+a4+a7==deg,Mod[a4-a7,3]==2},Integers]; (* needs to accomodate the second grading *)
Mons2=R0^a0 R1^a1 R4^a4 R7^a7/.sols;
maybe=Join[R2^2 R8  Mons1, Mons2];
Length[maybe]
values=Table[maybe[[k]]/.PointsR,{k,1,Length[maybe]}];
ns=NullSpace[Transpose[values]];
Length[ns]
rr=RowReduce[ns,Tolerance->10^(-400)];
Accuracy[rr]
eqs=RootApproximant[Chop[rr,10^(-400)],4].maybe;
eqs=SimplifyEquationsRootM2M3[RootReduce[eqs]]; 
R22R8toR=-(eqs[[1]]/.R2->0)/Coefficient[eqs[[1]],R2^2 R8];

We will now do a test of our results. We will compute some points for the V-s, then compute R-s and then check that we get a point on the original FPP. We will also time this process.

Timing[
sols=NSolve[Join[{1==V0,0==Sum[Vs[[k]]RandomInteger[1000],{k,1,7}],0==Sum[Vs[[k]]RandomInteger[1000],{k,1,7}]},Table[0==eqs8[[k]],{k,1,6}]],800,WorkingPrecision->200];
]

Timing[
pt=sols[[1]];
pt=Join[pt, R0147toV/.pt];
AppendTo[pt,R2->(R2cubeToR/.pt)^(1/3)];
AppendTo[pt,R3-> ((R2R3toR/R2)/.pt)];
AppendTo[pt,R6-> ((R2R6toR/R2)/.pt)];
AppendTo[pt,R9-> ((R2R9toR/R2)/.pt)];
AppendTo[pt,R5-> ((R22R5toR/R2^2)/.pt)];
AppendTo[pt,R8-> ((R22R8toR/R2^2)/.pt)];
pt=Join[pt, UtildesToU/.UtoRV/.pt];
]
DeleteDuplicates[N[EqsRsmall/.pt]]

(* Now we will build a database of points *)
PointsVRUtilde={};

For[num=1,num<=30,num++,
If[Mod[num,7]==3,Print[num]];
sols=NSolve[Join[{1==V0,0==Sum[Vs[[k]]RandomInteger[1000],{k,1,7}],0==Sum[Vs[[k]]RandomInteger[1000],{k,1,7}]},Table[0==eqs8[[k]],{k,1,6}]],1200,WorkingPrecision->300];
For[l=1,l<=7,l++,
pt=sols[[l]];
pt=Join[pt, R0147toV/.pt];
AppendTo[pt,R2->(R2cubeToR/.pt)^(1/3)];
AppendTo[pt,R3-> ((R2R3toR/R2)/.pt)];
AppendTo[pt,R6-> ((R2R6toR/R2)/.pt)];
AppendTo[pt,R9-> ((R2R9toR/R2)/.pt)];
AppendTo[pt,R5-> ((R22R5toR/R2^2)/.pt)];
AppendTo[pt,R8-> ((R22R8toR/R2^2)/.pt)];
pt=Join[pt,UtildesToU/.UtoRV/.pt];
If[Accuracy[pt]>800,
AppendTo[PointsVRUtilde,pt];
]
];
];

{Length[PointsVRUtilde],Accuracy[PointsVRUtilde]}

We will now compute the remaining 40 sections of H^0(8FPP, 6H). We will write them as rational functions in V-s and R2,  of total V-degree 2, which are zero on the V0+V1+V2=0 union V3=0. We first get a database of points on these two curves.

PointsV012={};
For[num=1,num<=20,num++,
If[Mod[num,7]==3,Print[num]];
sols=NSolve[Join[{1==V0,0==V0+V1+V2,0==Sum[Vs[[k]]RandomInteger[1000],{k,1,7}]},Table[0==eqs8[[k]],{k,1,6}]],800,WorkingPrecision->300];
For[l=1,l<=5,l++,
pt=sols[[l]];
pt=Join[pt, R0147toV/.pt];
AppendTo[pt,R2->(R2cubeToR/.pt)^(1/3)];
AppendTo[pt,R3-> ((R2R3toR/R2)/.pt)];
AppendTo[pt,R6-> ((R2R6toR/R2)/.pt)];
AppendTo[pt,R9-> ((R2R9toR/R2)/.pt)];
AppendTo[pt,R5-> ((R22R5toR/R2^2)/.pt)];
AppendTo[pt,R8-> ((R22R8toR/R2^2)/.pt)];
pt=Join[pt,Table[Utildes[[k]]->(((Utildes[[k]]/. UtildesToU)/.UtoRV)/.pt),{k,1,40}]];
AppendTo[PointsV012,pt];
];
];

PointsV3={};
For[num=1,num<=20,num++,
If[Mod[num,7]==3,Print[num]];
sols=NSolve[Join[{1==V0,0==V3,0==Sum[Vs[[k]]RandomInteger[1000],{k,1,7}]},Table[0==eqs8[[k]],{k,1,6}]],800,WorkingPrecision->300];
For[l=1,l<=5,l++,
pt=sols[[l]];
pt=Join[pt, R0147toV/.pt];
AppendTo[pt,R2->(R2cubeToR/.pt)^(1/3)];
AppendTo[pt,R3-> ((R2R3toR/R2)/.pt)];
AppendTo[pt,R6-> ((R2R6toR/R2)/.pt)];
AppendTo[pt,R9-> ((R2R9toR/R2)/.pt)];
AppendTo[pt,R5-> ((R22R5toR/R2^2)/.pt)];
AppendTo[pt,R8-> ((R22R8toR/R2^2)/.pt)];
pt=Join[pt,Table[Utildes[[k]]->(((Utildes[[k]]/. UtildesToU)/.UtoRV)/.pt),{k,1,40}]];
AppendTo[PointsV3,pt];
];
];

{Length[PointsV012],Accuracy[PointsV012]}
{Length[PointsV3],Accuracy[PointsV3]}

We will now compute the weight 0 mod 3 new sections.

Vto1=Table[Vs[[k]]->1,{k,1,7}];
Mons4V=MonomialList[Sum[Vs[[k]],{k,1,7}]^4];
Mons4V=Table[Mons4V[[k]]/(Mons4V[[k]]/.Vto1),{k,1,Length[Mons4V]}];
(* we want the odd ones *)
goodMons4V={};
For[k=1,k<=Length[Mons4V],k++,
If[Mod[Sum[Exponent[Mons4V[[k]],Vs[[l]]],{l,4,7}],2]==0, (* we will be dividing by (V0+V1+V2)V3, so we want this even to get odd sections *)
AppendTo[goodMons4V,Mons4V[[k]]];
]
]
Mons4V=goodMons4V;

Length[Mons4V]
MatrixRank[Mons4V/.PointsVRUtilde]
MatrixRank[Mons4V/.Join[PointsV3,PointsV012]]

(* finding the space of linear combinations that are zero on both curves *)
ns=NullSpace[Mons4V/.Join[PointsV3,PointsV012]];
Length[ns]
rr=RowReduce[ns,Tolerance->10^(-500)];
rrexact=RootApproximant[Chop[rr,10^(-500)],4];
polys=SimplifyEquationsRootM2M3[RootReduce[rrexact.Mons4V]];

goodpolys={};
rk=0;
values={};
For[k=1,k<=Length[polys],k++,
newrow=polys[[k]]/.PointsVRUtilde;
If[MatrixRank[Append[values,newrow]]>rk,
AppendTo[values,newrow];
rk++;
Print["rk=",rk];
AppendTo[goodpolys,polys[[k]]];
]
]

NewUs=goodpolys/(V0+V1+V2)/V3;

Now we want to get the odd ones of C3 weight 1. We will try to write them as linear combinations of R2,R5,R8,\tilde U17, \tilde U18, \tilde U19, times degree two polys in V-s, divided by (V0+V1+V2)/V3 where P-s are of degree two and even.

Mons2V=MonomialList[Sum[Vs[[k]],{k,1,7}]^2];
Mons2V=Table[Mons2V[[k]]/(Mons2V[[k]]/.Vto1),{k,1,Length[Mons2V]}];
(* we want the odd ones *)
goodMons2V={};
For[k=1,k<=Length[Mons2V],k++,
If[Mod[Sum[Exponent[Mons2V[[k]],Vs[[l]]],{l,4,7}],2]==0, (* we will be dividing by (V0+V1+V2)V3, so we want this even to get odd sections *)
AppendTo[goodMons2V,Mons2V[[k]]];
]
]
Mons2V=goodMons2V;

Mons=Join[R2 Mons2V, R5 Mons2V, R8 Mons2V, Utildes[[18]]Mons2V, Utildes[[19]] Mons2V,Utildes[[20]]Mons2V];
Length[Mons]
MatrixRank[Mons/.PointsVRUtilde]
MatrixRank[Mons/.Join[PointsV012,PointsV3]]


(* finding the space of linear combinations that are zero on both curves *)

ns=NullSpace[Mons/.Join[PointsV012,PointsV3]];
Length[ns]
rr=RowReduce[ns,Tolerance->10^(-500)];
rrexact=RootApproximant[Chop[rr,10^(-500)],4];
polys=SimplifyEquationsRootM2M3[RootReduce[rrexact.Mons]]; 

goodpolys={};
rk=0;
values={};
For[k=1,k<=Length[polys],k++,
newrow=polys[[k]]/.PointsVRUtilde;
If[MatrixRank[Append[values,newrow]]>rk,
AppendTo[values,newrow];
rk++;
Print["rk=",rk];
AppendTo[goodpolys,polys[[k]]];
]
]

MatrixRank[Join[Utildes,goodpolys/V3/(V0+V1+V2)]/.PointsVRUtilde]

NewUswt1=goodpolys/(V0+V1+V2)/V3;

Now we want to get the odd ones of C3 weight 2. We will try to write them as linear combinations of R3,R6,R9, \tilde U14, \tilde U15, \tilde U16, times degree two polys in V-s, divided by (V0+V1+V2)/V3 where P-s are of degree two and even.

Mons2V=MonomialList[Sum[Vs[[k]],{k,1,7}]^2];
Mons2V=Table[Mons2V[[k]]/(Mons2V[[k]]/.Vto1),{k,1,Length[Mons2V]}];
(* we want the odd ones *)
goodMons2V={};
For[k=1,k<=Length[Mons2V],k++,
If[Mod[Sum[Exponent[Mons2V[[k]],Vs[[l]]],{l,4,7}],2]==0, (* we will be dividing by (V0+V1+V2)V3, so we want this even to get odd sections *)
AppendTo[goodMons2V,Mons2V[[k]]];
]
]
Mons2V=goodMons2V;

Mons=Join[R3 Mons2V, R6 Mons2V, R9 Mons2V, Utildes[[15]]Mons2V, Utildes[[16]] Mons2V,Utildes[[17]]Mons2V];
Length[Mons]
MatrixRank[Mons/.PointsVRUtilde]
MatrixRank[Mons/.Join[PointsV012,PointsV3]]

(* finding the space of linear combinations that are zero on both curves *)

ns=NullSpace[Mons/.Join[PointsV012,PointsV3]];
Length[ns]
rr=RowReduce[ns,Tolerance->10^(-500)];
rrexact=RootApproximant[Chop[rr,10^(-500)],4];
polys=SimplifyEquationsRootM2M3[RootReduce[rrexact.Mons]]; 

goodpolys={};
rk=0;
values={};
For[k=1,k<=Length[polys],k++,
newrow=polys[[k]]/.PointsVRUtilde;
If[MatrixRank[Append[values,newrow]]>rk,
AppendTo[values,newrow];
rk++;
Print["rk=",rk];
AppendTo[goodpolys,polys[[k]]];
]
]


NewUswt2=goodpolys/(V0+V1+V2)/V3;
MatrixRank[Join[Utildes,NewUs,NewUswt1,NewUswt2]/.PointsVRUtilde]

We now want to compute the action of two C3-s on the 80-dimensional space H^0(8FPP,6H).

w=Exp[2Pi I /3];

C3actionU
C3action=Table[Rs[[k]]->w^wtsC3C3[[k]][[1]]Rs[[k]],{k,1,10}]
C3actionV={V0->V1,V1->V2,V2->V0,V3-> w^2 V3, V4-> w V4, V5->V5, V6->V6};

Table[wtsC3C3[[i]][[1]],{i,1,10}]
Table[wtsC3C3[[i]][[2]],{i,1,10}]

Table[FullSimplify[((Utildes[[k]]/.UtildesToU)/.C3actionU)/(Utildes[[k]]/.UtildesToU)],{k,1,40}]

Table[Simplify[((Us[[k]]/.C3actionU)/.UtoR) -(( Us[[k]]/.UtoR)/.C3action)],{k,1,40}]

(* getting two actions on V, R, U, Utilde up to 40 *)

C3action1UVR=Join[C3actionV,C3actionU,C3action];
wtsC3C3Utilde1={};
For[k=1,k<=40,k++,
For[l=0,l<=2,l++,
If[w^l == FullSimplify[((Utildes[[k]]/.UtildesToU)/.C3action1UVR)/(Utildes[[k]]/.UtildesToU)],
AppendTo[wtsC3C3Utilde1,l]
]
];
];
wtsC3C3Utilde1
C3action1UVRUtilde=Join[C3action1UVR,Table[Utildes[[k]]->Utildes[[k]]w^wtsC3C3Utilde1[[k]],{k,1,40}]];


C3action2UVR=Join[C3actionU2,Table[Rs[[k]]->Rs[[k]] w^wtsC3C3[[k]][[2]],{k,1,10}]]; (* action on V is trivial *)
wtsC3C3Utilde2={};
For[k=1,k<=40,k++,
For[l=0,l<=2,l++,
If[w^l == FullSimplify[((Utildes[[k]]/.UtildesToU)/.C3action2UVR)/(Utildes[[k]]/.UtildesToU)],
AppendTo[wtsC3C3Utilde2,l]
]
];
];
wtsC3C3Utilde2
C3action2UVRUtilde=Join[C3action2UVR,Table[Utildes[[k]]->Utildes[[k]]w^wtsC3C3Utilde2[[k]],{k,1,40}]];

(* now constructing the action on new U-s *)
(* we see that the second action is easy *)
Simplify[(NewUs/.C3action2UVRUtilde)-NewUs]
Simplify[(NewUswt1/.C3action2UVRUtilde)- w  NewUswt1]
Simplify[(NewUswt2/.C3action2UVRUtilde)- w^2  NewUswt2]

(* for the first action *)
Ts={T1,T2,T3,T4,T5,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15,T16};
actionNewUs={};
For[l=1,l<=16,l++,
tbz=((NewUs[[l]]/.C3action1UVRUtilde)==Sum[UnknownCoeffs[[k]]NewUs[[k]],{k,1,16}])/.Table[PointsVRUtilde[[k]],{k,1,20}];
AppendTo[actionNewUs,Ts[[l]]->Sum[Ts[[k]]RootApproximant[Chop[UnknownCoeffs[[k]]/.NSolve[tbz,550][[1]],10^(-500)],4],{k,1,16}]]
(* it was somehow sensitive to solving with the correct accuracy, with 550 working *)
];

actionNewUswt1={};
For[l=1,l<=12,l++,
tbz=((NewUswt1[[l]]/.C3action1UVRUtilde)==Sum[UnknownCoeffs[[k]]NewUswt1[[k]],{k,1,12}])/.Table[PointsVRUtilde[[k]],{k,1,20}];
AppendTo[actionNewUswt1,Ts[[l]]->Sum[Ts[[k]]RootApproximant[Chop[UnknownCoeffs[[k]]/.NSolve[tbz][[1]],10^(-500)],4],{k,1,12}]]
];

actionNewUswt2={};
For[l=1,l<=12,l++,
tbz=((NewUswt2[[l]]/.C3action1UVRUtilde)==Sum[UnknownCoeffs[[k]]NewUswt2[[k]],{k,1,12}])/.Table[PointsVRUtilde[[k]],{k,1,20}];
AppendTo[actionNewUswt2,Ts[[l]]->Sum[Ts[[k]]RootApproximant[Chop[UnknownCoeffs[[k]]/.NSolve[tbz][[1]],10^(-500)],4],{k,1,12}]]
];

(* this verifies that we have results in the correct field *)
blah=Table[Coefficient[(Ts[[k]]/.actionNewUs),Ts[[l]]],{k,16},{l,16}];
DeleteDuplicates[Flatten[Table[RootReduce[To4[blah[[k]][[l]]].{1,Sqrt[-2],Sqrt[-3],Sqrt[6]}]==blah[[k]][[l]],{k,1,16},{l,1,16}]]]
actionNewUssimple=Table[Ts[[k]]->Sum[To4[blah[[k]][[l]]].{1,Sqrt[-2],Sqrt[-3],Sqrt[6]} Ts[[l]],{l,1,16}],{k,1,16}];
blah=Table[Coefficient[(Ts[[k]]/.actionNewUswt1),Ts[[l]]],{k,12},{l,12}];
DeleteDuplicates[Flatten[Table[RootReduce[To4[blah[[k]][[l]]].{1,Sqrt[-2],Sqrt[-3],Sqrt[6]}]==blah[[k]][[l]],{k,1,12},{l,1,12}]]]
actionNewUswt1simple=Table[Ts[[k]]->Sum[To4[blah[[k]][[l]]].{1,Sqrt[-2],Sqrt[-3],Sqrt[6]} Ts[[l]],{l,1,12}],{k,1,12}];
blah=Table[Coefficient[(Ts[[k]]/.actionNewUswt2),Ts[[l]]],{k,12},{l,12}];
DeleteDuplicates[Flatten[Table[RootReduce[To4[blah[[k]][[l]]].{1,Sqrt[-2],Sqrt[-3],Sqrt[6]}]==blah[[k]][[l]],{k,1,12},{l,1,12}]]]
actionNewUswt2simple=Table[Ts[[k]]->Sum[To4[blah[[k]][[l]]].{1,Sqrt[-2],Sqrt[-3],Sqrt[6]} Ts[[l]],{l,1,12}],{k,1,12}];

We now want to build the eigenspaces in each of the new subspaces. These will be given names Utilde40,...,Utilde79.

Utildes= Table[ToExpression[StringJoin["Utilde",ToString[k]]],{k,0,79}];
WtC3C3all80Utildes2=Join[wtsC3C3Utilde2,Table[0,{i,1,16}],Table[1,{i,1,12}],Table[2,{i,1,12}]]
WtC3C3all80Utildes1=Join[wtsC3C3Utilde1,Table[0,{i,1,6}],Table[1,{i,1,5}],Table[2,{i,1,5}],Table[0,{i,1,4}],Table[1,{i,1,4}],Table[2,{i,1,4}],Table[0,{i,1,4}],Table[1,{i,1,4}],Table[2,{i,1,4}]]

NewUtildesToAllother={};
count=41; (* this is the first new Utilde *)

(* first the NewUS, weight 0 *)
rk=0;
values={};
TtoPoly=Table[Ts[[k]]->NewUs[[k]],{k,1,16}];

For[k=1,k<=16,k++,
new=Together[(Ts[[k]]+ (Ts[[k]]/.actionNewUssimple)+((Ts[[k]]/.actionNewUssimple)/.actionNewUssimple))/.TtoPoly];
newvalues=new/.PointsVRUtilde;
If[MatrixRank[Append[values,newvalues]]>rk,
rk++;
AppendTo[values,newvalues];
AppendTo[NewUtildesToAllother,Utildes[[count]]->new];
count++;
]
]
Print[rk];

(* wt 1*)
rk=0;
values={};
For[k=1,k<=16,k++,
new=Together[(Ts[[k]]+ w^2(Ts[[k]]/.actionNewUssimple)+w((Ts[[k]]/.actionNewUssimple)/.actionNewUssimple))/.TtoPoly];
newvalues=new/.PointsVRUtilde;
If[MatrixRank[Append[values,newvalues]]>rk,
rk++;
AppendTo[values,newvalues];
AppendTo[NewUtildesToAllother,Utildes[[count]]->new];
count++;
]
]
Print[rk];

(* wt 2*)
rk=0;
values={};
For[k=1,k<=16,k++,
new=Together[(Ts[[k]]+ w(Ts[[k]]/.actionNewUssimple)+w^2((Ts[[k]]/.actionNewUssimple)/.actionNewUssimple))/.TtoPoly];
newvalues=new/.PointsVRUtilde;
If[MatrixRank[Append[values,newvalues]]>rk,
rk++;
AppendTo[values,newvalues];
AppendTo[NewUtildesToAllother,Utildes[[count]]->new];
count++;
]
]
Print[rk];

(* now the NewUSwt1 *)
rk=0;
values={};
TtoPoly=Table[Ts[[k]]->NewUswt1[[k]],{k,1,12}];

For[k=1,k<=12,k++,
new=Together[(Ts[[k]]+ (Ts[[k]]/.actionNewUswt1simple)+((Ts[[k]]/.actionNewUswt1simple)/.actionNewUswt1simple))/.TtoPoly];
newvalues=new/.PointsVRUtilde;
If[MatrixRank[Append[values,newvalues]]>rk,
rk++;
AppendTo[values,newvalues];
AppendTo[NewUtildesToAllother,Utildes[[count]]->new];
count++;
]
]
Print[rk];

rk=0;
values={};
For[k=1,k<=12,k++,
new=Together[(Ts[[k]]+ w^2 (Ts[[k]]/.actionNewUswt1simple)+ w((Ts[[k]]/.actionNewUswt1simple)/.actionNewUswt1simple))/.TtoPoly];
newvalues=new/.PointsVRUtilde;
If[MatrixRank[Append[values,newvalues]]>rk,
rk++;
AppendTo[values,newvalues];
AppendTo[NewUtildesToAllother,Utildes[[count]]->new];
count++;
]
]
Print[rk];

rk=0;
values={};
For[k=1,k<=12,k++,
new=Together[(Ts[[k]]+ w (Ts[[k]]/.actionNewUswt1simple)+ w^2((Ts[[k]]/.actionNewUswt1simple)/.actionNewUswt1simple))/.TtoPoly];
newvalues=new/.PointsVRUtilde;
If[MatrixRank[Append[values,newvalues]]>rk,
rk++;
AppendTo[values,newvalues];
AppendTo[NewUtildesToAllother,Utildes[[count]]->new];
count++;
]
]
Print[rk];

(* now the NewUSwt2 *)
rk=0;
values={};
TtoPoly=Table[Ts[[k]]->NewUswt2[[k]],{k,1,12}];

For[k=1,k<=12,k++,
new=Together[(Ts[[k]]+ (Ts[[k]]/.actionNewUswt2simple)+((Ts[[k]]/.actionNewUswt2simple)/.actionNewUswt2simple))/.TtoPoly];
newvalues=new/.PointsVRUtilde;
If[MatrixRank[Append[values,newvalues]]>rk,
rk++;
AppendTo[values,newvalues];
AppendTo[NewUtildesToAllother,Utildes[[count]]->new];
count++;
]
]
Print[rk];

rk=0;
values={};
For[k=1,k<=12,k++,
new=Together[(Ts[[k]]+ w^2 (Ts[[k]]/.actionNewUswt2simple)+ w((Ts[[k]]/.actionNewUswt2simple)/.actionNewUswt2simple))/.TtoPoly];
newvalues=new/.PointsVRUtilde;
If[MatrixRank[Append[values,newvalues]]>rk,
rk++;
AppendTo[values,newvalues];
AppendTo[NewUtildesToAllother,Utildes[[count]]->new];
count++;
]
]
Print[rk];

rk=0;
values={};
For[k=1,k<=12,k++,
new=Together[(Ts[[k]]+ w (Ts[[k]]/.actionNewUswt2simple)+ w^2((Ts[[k]]/.actionNewUswt2simple)/.actionNewUswt2simple))/.TtoPoly];
newvalues=new/.PointsVRUtilde;
If[MatrixRank[Append[values,newvalues]]>rk,
rk++;
AppendTo[values,newvalues];
AppendTo[NewUtildesToAllother,Utildes[[count]]->new];
count++;
]
]
Print[rk];

We will now construct a lot of points with high accuracy for all 80 Utilde (and V and R).

(* Now we will build a database of points *)
PointsAllUtilde={};

numpts=300;
While[Length[PointsAllUtilde]<numpts,
If[Mod[Length[PointsAllUtilde],36]==6,
Print[{Length[PointsAllUtilde],Accuracy[PointsAllUtilde]}]
];
sols=NSolve[Join[{1==V0,0==Sum[Vs[[k]]RandomInteger[100],{k,1,7}],0==Sum[Vs[[k]]RandomInteger[100],{k,1,7}]},Table[0==eqs8[[k]],{k,1,6}]],4200,WorkingPrecision->800];
For[l=1,l<=7,l++,
pt=sols[[l]];
pt=Join[pt, R0147toV/.pt];
AppendTo[pt,R2->(R2cubeToR/.pt)^(1/3)];
AppendTo[pt,R3-> ((R2R3toR/R2)/.pt)];
AppendTo[pt,R6-> ((R2R6toR/R2)/.pt)];
AppendTo[pt,R9-> ((R2R9toR/R2)/.pt)];
AppendTo[pt,R5-> ((R22R5toR/R2^2)/.pt)];
AppendTo[pt,R8-> ((R22R8toR/R2^2)/.pt)];
pt=Join[pt,UtildesToU/.UtoRV/.pt];
pt=Join[pt,NewUtildesToAllother/.pt];
If[Accuracy[pt]>2500,
AppendTo[PointsAllUtilde,pt];
]
];
];

{Length[PointsAllUtilde],Accuracy[PointsAllUtilde]}

We will now save a lot of this information, so that we can start with a fresh file.

fname=OpenWrite["VariablesAndRelationsStep4"];
Write[fname,{Rs,Us,Utildes,UtildesToU,NewUtildesToAllother,UtoRV,R0147toV,R2cubeToR,R22R5toR,R22R8toR,R2R3toR,R2R6toR,R2R9toR}];
Close[fname];

fname=OpenWrite["PointsAllUtilde"];
Write[fname,PointsAllUtilde];
Close[fname];

fname=OpenWrite["ActionsStep4"];
Write[fname,{WtC3C3all80Utildes1,WtC3C3all80Utildes2,C4actionV}];
Close[fname];

fname=OpenWrite["PointsVRU"];
Write[fname,PointsVRU];
Close[fname];

fname=OpenWrite["eqs8"];
Write[fname,eqs8];
Close[fname];

fname=OpenWrite["PointsVRUtilde"];
Write[fname,PointsVRUtilde];
Close[fname];