This notebook is a supporting file for Section 4 of the paper.

1. Find dimension 19 space of sections of H^0(X,4K). Find the dimension 13 subspace that has vanishing conditions 2A_1+B_1 and A_2+2B_2 on the exceptional curves of the blowup of singular points. Find the dimension 6 space of C_3 invariant sections of H^0(X,4K) that pass through the singular points (easy).

2. Recheck the 2A+B conditions on the 13-dimensional space. Store the restriction maps to the curves, which are used later.

3. Find/verify the basis of the C3 invariant sections of 12K_FPP/C3. Set up 67 cubic equations. Reduce to 58 equations.

4. Make a choice and set up three linear equations. Simplify 58 equations to reduce to 55 cubic equations.

5. Find 12 quadratic equations, simplify them. Write out the equations in Magma format.

6. Solve the equations to find f and d.

SetDirectory[$InitialDirectory];
SetDirectory["DeskTop/EnricoFPP/paperfiles"];

fname="EquationsOfXintermsofW";
OpenRead[fname];
AllEqs23Wgood=Read[fname];
Close[fname];

1. Find dimension 19 space of sections of H^0(X,4K). Find the dimension 13 subspace that has vanishing conditions 2A_1+B_1 and A_2+2B_2 on the exceptional curves of the blowup of singular points. Find the dimension 6 space of of invariant sections of H^0(X,4K) that pass through the singular points.

(* we start by finding some random points on X with high accuracy *)
VarsW={W0,W1,W2,W3,W4,W5,W6,W7};
prec=6000;
SetOfPointsW = {};
numpts=200;
For[i=1,i≤numpts,i++,
If[Mod[i,numpts/10]==0,Print[i]];
rpt1=(RandomInteger[{1,100000}]+I RandomInteger[{1,100000}])/100000;
rpt2=(RandomInteger[{1,100000}]+I RandomInteger[{1,100000}])/100000;
Eqs=Join[{W0==1,W1+W2+W3==rpt1,W4+W5+W6==rpt2},Table[0==AllEqs23Wgood[[j]],{j,1,9}]];
sols=NSolve[Eqs,200];
fr=FindRoot[Table[Eqs[[k]],{k,1,8}],Table[{VarsW[[j]],VarsW[[j]]/.sols[[1]]},{j,1,8}],MaxIterations->100,WorkingPrecision->prec,AccuracyGoal->prec];
AppendTo[SetOfPointsW,fr];
];

(* quality control *)
For[i=1,i≤numpts,i++, If[Chop[Sum[Abs[AllEqs23Wgood[[j]]]^2/.SetOfPointsW[[i]],{j,1,9}]]≠0,Print["warning, point ",i," looks fishy"];];];

(* I will now find points on the curve W4+W5+W6=1 *)
SetOfPointsWwithW456={};
For[numpt=1,numpt≤100,numpt++,
If[Mod[numpt,10]==0,Print[numpt]];
rpt1= numpt + I RandomInteger[1000]/1000;
RelsToSolve=Join[{W2==rpt1,W0==1,W4+W5+W6==0},Table[AllEqs23Wgood[[i]]==0,{i,1,9}]];
sols=NSolve[RelsToSolve,300];
fr=FindRoot[Table[RelsToSolve[[k]],{k,1,8}],Table[{VarsW[[j]],VarsW[[j]]/.sols[[1]]},{j,1,8}],MaxIterations->100,WorkingPrecision->prec,AccuracyGoal->prec];
AppendTo[SetOfPointsWwithW456,N[fr,prec]];
];

(* Looking for sections of 6K that are zero on W4+W5+W6=0 *)
maxlength=500;
CoeffList=Table[ToExpression[StringJoin["A",ToString[i]]],{i,1,maxlength}];

AllMons=MonomialList[(W1+W2+W3+W0)^2(W4+W5+W6+W7) + (W4+W5+W6)^3];
count=Length[AllMons];
AllMons=Table[AllMons[[i]]/(AllMons[[i]]/.Table[VarsW[[k]]->1,{k,1,8}]),{i,1,count}];

UnivPoly=Sum[AllMons[[k]] CoeffList[[k]],{k,1,count}];
Eqs=Table[0==UnivPoly/.SetOfPointsWwithW456[[k]],{k,1,count+20}];
sols=Chop[ Solve[Eqs][[1]]];

AllEqs={};
For[j=1,j≤count,j++,
UEq=Coefficient[UnivPoly/.sols,CoeffList[[j]]];
If[ByteCount[UEq]>100,
Eq=Sum[RootApproximant[Coefficient[UEq,AllMons[[n]]],2] AllMons[[n]],{n,1,Length[AllMons]}];
AppendTo[AllEqs,Eq];
];
];
Length[AllEqs]

(* Checking to see which of these could be added to the degree two polynomials *)
For[j=1,j≤Length[AllEqs],j++,
D4P=Append[{(W0^2),W0 W1,W0 W2,W0 W3,(W1^2),(W2^2),(W3^2),W1 W2,W2 W3,W1 W3,W4 W7,W5 W7,W6 W7,W4 W5,W5 W6,W4 W6,W7^2},AllEqs[[j]]/(W4+W5+W6)];
AppendTo[D4P,D4P[[18]]/.{W1->W2,W2->W3,W3->W1,W4->W5,W5->W6,W6->W4}];
Print[j," ",ByteCount[AllEqs[[j]]]," ",MatrixRank[Table[D4P/.SetOfPointsW[[k]],{k,1,50}]]];
];

(* writing into the file *)
fname="Deg4inPluckerW";
OpenWrite[fname];
Write[fname,D4P];
Close[fname];

(* reading from the file *)
fname="Deg4inPluckerW";
OpenRead[fname];
Deg4inPluckerW=Read[fname];
Close[fname];

(* finding the conditions of passing through the first two singular points *)
VarsC={c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,c14,c15,c16,c17,c18,c19};
UnivDeg4=VarsC.Deg4inPluckerW;
pt1={W0->0,W1->1,W2->0,W3->0,W4->1,W5->0,W6->0,W7->0};
pt2={W0->0,W1->0,W2->1,W3->0,W4->0,W5->1,W6->0,W7->0};
pt3={W0->0,W1->0,W2->0,W3->1,W4->0,W5->0,W6->1,W7->0};
EqsOnC={0==UnivDeg4/.pt1,0==UnivDeg4/.pt2};
solsc=Simplify[Solve[EqsOnC,{c5,c6}]][[1]];
(* tell us which coefficients are solved for *)
Print[Table[solsc[[i]][[1]],{i,1,Length[solsc]}]];
UnivDeg4throughpt1and2=UnivDeg4/.solsc;

We now find the exceptional locus of the blowup at the first point and one of the lines there.

VarT={T0,T2,T3,T4,T5,T6,T7};
(* T0 = W0/W1, T2 = W2/W0, ..., T4 = (W4-W1)/W0,... , T7 = W7/W0 are 
affine coordinates on one of the open sets on the cover of blowup of P^7 at (0:1:0:0:1:0:0:0) *)
TtoW ={W0->T0,W1->1,W2->T0 T2 ,W3->T0 T3,W4->1+T0 T4 ,W5->T0 T5 ,W6->T0 T6,W7->T0 T7 };
EqsW5=Table[AllEqs23Wgood[[k]],{k,1,5}];(* the first five equatons cut out the surface near pt1 *)
(* all equations are divisible by T0, which we remove *)
RelsT=Table[Sum[SeriesCoefficient[(AllEqs23Wgood[[j]]/.TtoW),{T0,0,k}] T0^(k-1),{k,1,5}],{j,1,5}];
(* coordinates on the exceptional divisor of blowup of the point in P^7 *)
VarT2={T2,T3,T4,T5,T6,T7};
solsingeq=Simplify[Solve[Table[0==Coefficient[{1,a2,a3,a4,a5}.RelsT/.T0->0,VarT2[[k]]],{k,1,6}]]];
(* This is a cut that's singular at pt1 *)
singeq =Simplify[({1,a2,a3,a4,a5}.EqsW5/.solsingeq[[1]])];
singeqT=Simplify[singeq/.TtoW];

(* finding a plane in which the exceptional divisor lies *)
solsexclin=Simplify[Solve[Table[0==RelsT[[j]]/.T0->0,{j,1,5}]][[1]]]

(* finding the two exceptional lines *)
sollines=Simplify[Solve[SeriesCoefficient[(singeqT/.solsexclin),{T0,0,2}]==0,T2]];
(**************************************************************************************)
(* There is a choice to be made here. A wrong choice results in no solutions later on. *)
(* This choice also causes some modification of the later code, since Mathematica solves *)
(* for different variables and we need to have the lines correctly chosen later *)

Line2=sollines[[1]]  (*{T2->((21-61 I (Sqrt[15])+252 I (Sqrt[15]) T3)/1512)} *)
Line1=sollines[[2]]  (*{T2->(1/504) (-349-203 I (Sqrt[15]))+(-(9/2)-I (Sqrt[15])) T3}*)

(* finding a coordinate chart on the blowup so that the complement contains the second line
by replacing W0 by W0n *)
W0n=(((T2/.Line2)/.T3->0)) W0- W2 + (Coefficient[T2/.Line2,T3]) W3;
(* recall the original chance {W0->T0,W1->1,W2->T0 T2 ,W3->T0 T3,W4->1+T0 T4 ,W5->T0 T5 ,W6->T0 T6,W7->T0 T7 }*)
TtoW=Solve[{W1==1,W0n==T0, W2==T0  T2, W3==T0 T3, W4==T0 T4+1, W5== T0 T5, W6==T0 T6, W7==T0 T7},VarsW][[1]];
RelsT=Table[Sum[SeriesCoefficient[(AllEqs23Wgood[[j]]/.TtoW),{T0,0,k}] T0^(k-1),{k,1,5}],{j,1,5}];
solsexc=Simplify[Solve[Table[0==RelsT[[j]]/.T0->0,{j,1,5}]]];

(* coordinates on the exceptional divisor of blowup of the point in P^7 *)
VarT2={T2,T3,T4,T5,T6,T7};
(* finding singular equation and its restriction to the exceptional divisor of the blowup of the point in P^7 *)
solsingeq=Solve[Table[0==Coefficient[{1,a2,a3,a4,a5}.RelsT/.T0->0,VarT2[[k]]],{k,1,6}]];
singeq =({1,a2,a3,a4,a5}.Table[AllEqs23Wgood[[k]],{k,1,5}]/.solsingeq[[1]]);
singeqT=Sum[SeriesCoefficient[singeq/.TtoW,{T0,0,k}] T0^{k-2},{k,2,6}];
solsexclin=Simplify[Solve[0==Table[RelsT[[j]]/.T0->0,{j,1,5}]]][[1]]
singeqTline=FullSimplify[singeqT/.solsexclin/.T0->0]

(* We think of the neighborhood of the line as a map from Spec C[y_1][[y_2]] to the blowup. *)
Line1=Simplify[Solve[0==singeqTline,T2][[1]][[1]]]
UptoZeroiny2=Simplify[({T0->0,T2->(T2/.Line1/.T3->y1),T3->y1,T4->(T4/.solsexclin/.Line1/.T3->y1),T5->(T5/.solsexclin/.Line1/.T3->y1),T6->(T6/.solsexclin/.Line1/.T3->y1),T7->(T7/.solsexclin/.Line1/.T3->y1)})];

(* Look at the restriction of the section of O(4H) with condition of passing through the two points at the exceptional line. *)
tbz=FullSimplify[(Sum[SeriesCoefficient[UnivDeg4throughpt1and2/.TtoW,{T0,0,k}]T0^(k-1),{k,1,5}])/.UptoZeroiny2];
Eqs2AB={0==tbz/.y1->0,0==tbz/.y1->1};

We now find the exceptional locus of the blowup at the second point and an appropriate line there.
The procedure is similar.

EqsW5=Table[AllEqs23Wgood[[k]],{k,1,5}];
VarT={T0,T1,T3,T4,T5,T6,T7};
TtoW ={W0->T0,W1->T0 T1,W2->1 ,W3->T0 T3,W4->T0 T4 ,W5->1+T0 T5 ,W6->T0 T6,W7->T0 T7 };
RelsT=Table[Sum[SeriesCoefficient[(AllEqs23Wgood[[j]]/.TtoW),{T0,0,k}] T0^(k-1),{k,1,5}],{j,1,5}];
VarT2={T1,T3,T4,T5,T6,T7};
solsingeq=Simplify[Solve[Table[0==Coefficient[{1,a2,a3,a4,a5}.RelsT/.T0->0,VarT2[[k]]],{k,1,6}]]];
singeq =Simplify[({1,a2,a3,a4,a5}.EqsW5/.solsingeq[[1]])];
singeqT=Simplify[singeq/.TtoW];
solsexclin=Simplify[Solve[Table[0==RelsT[[j]]/.T0->0,{j,1,5}]][[1]]];
sollines=Simplify[Solve[SeriesCoefficient[(singeqT/.solsexclin),{T0,0,2}]==0,T1]];

(* it is a bit tricky to know which line is the C3 translate of the previous first line *)
(* we want to use the other line *)
Line1=sollines[[1]] (*{T1->((I (9231 I-1129 (Sqrt[15])+504 (9 I+2 (Sqrt[15])) T3))/35532)} *)
Line2=sollines[[2]] (*{T1->(61/252)+(I/(12 (Sqrt[15])))-2 I (Sqrt[(3/5)]) T3}*)
Expand[Simplify[Solve[T1==(61/252)+(I/(12 (Sqrt[15])))-2 I (Sqrt[(3/5)]) T3,T3]]] (* this confirms that our choice is correct: the new Line 2 is the translate of old line 1 *)


W0n=(((T1/.Line2)/.T3->0)) W0- W1 + (Coefficient[T1/.Line2,T3]) W3;TtoW=Solve[{W2==1,W0n==T0, W1==T0  T1, W3==T0 T3, W4==T0 T4, W5== T0 T5+1, W6==T0 T6, W7==T0 T7},VarsW][[1]];

TtoW2nd=TtoW;

RelsT=Table[Sum[SeriesCoefficient[(AllEqs23Wgood[[j]]/.TtoW),{T0,0,k}] T0^(k-1),{k,1,5}],{j,1,5}];
solsexc=Simplify[Solve[Table[0==RelsT[[j]]/.T0->0,{j,1,5}]]];

solsingeq=Solve[Table[0==Coefficient[{1,a2,a3,a4,a5}.RelsT/.T0->0,VarT2[[k]]],{k,1,6}]];
singeq =({1,a2,a3,a4,a5}.Table[AllEqs23Wgood[[k]],{k,1,5}]/.solsingeq[[1]]);
singeqT=Sum[SeriesCoefficient[singeq/.TtoW,{T0,0,k}] T0^{k-2},{k,2,6}];
solsexclin=Solve[0==Table[RelsT[[j]]/.T0->0,{j,1,5}]][[1]];
singeqTquad=Simplify[singeqT/.solsexclin/.T0->0];
singeqTquad=singeqTquad/Coefficient[singeqTquad,T1];
singeqTquad=RootReduce[singeqTquad/.{T1->0,T3->0} ]+ RootReduce[Coefficient[singeqTquad,T1]]T1+ RootReduce[Coefficient[singeqTquad,T3]]T3;
Line1=Simplify[Solve[0==singeqTquad,T1][[1]][[1]]];
UptoZeroiny2=Simplify[({T0->0,T1->(T1/.Line1/.T3->y1),T3->y1,T4->(T4/.solsexclin/.Line1/.T3->y1),T5->(T5/.solsexclin/.Line1/.T3->y1),T6->(T6/.solsexclin/.Line1/.T3->y1),T7->(T7/.solsexclin/.Line1/.T3->y1)})];

tbz=(Sum[SeriesCoefficient[UnivDeg4throughpt1and2/.TtoW,{T0,0,k}]T0^(k-1),{k,1,5}])/.UptoZeroiny2;
Eqs2ABboth=Join[Eqs2AB,{0==tbz/.y1->0,0==tbz/.y1->1}];
sols2AB=FullSimplify[Solve[Eqs2ABboth,{c10,c11,c8,c9}]][[1]];
Print["solved for both 2A+B conditions"];

(* tell us which coefficients are solved for *)
Print[Table[sols2AB[[i]][[1]],{i,1,Length[sols2AB]}]];

(* find dimension 13 space *)
UnivDeg4with2ABboth=UnivDeg4throughpt1and2/.sols2AB;

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

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

(* find C3 invariant sections of 4K that vanish at the singular points *)
InvDeg2through={W0^2, W0(W1+W2+W3),W1^2+W2^2+W3^2-(W4^2+W5^2+W6^2), W1 W2 +W1 W3 + W2 W3, W7^2, W7( W4 + W5 + W6)}

2. Recheck the 2A+B conditions on the 13-dimensional space. Store the restriction maps to the curves, which are used later.

Varcd19={c1,c2,c3,c4,c7,c12,c13,c14,c15,c16,c17,c18,c19,d1,d2,d3,d4,d5,d6};
MatrixRank[Table[N[Coefficient[UnivDeg4with2ABboth,Varcd19[[k]]]/.SetOfPointsW[[l]],400],{k,1,13},{l,1,40}]]

I am switching from c7 to c0 that’s the value at the third point. Then we would be eventually looking for c0 nonzero.

singpt3={W0->0,W1->0,W2->0,W3->1,W4->0,W5->0,W6->1,W7->0};
Varcd19={c1,c2,c3,c4,c7,c12,c13,c14,c15,c16,c17,c18,c19,d1,d2,d3,d4,d5,d6};Solve[c0==UnivDeg4with2ABboth/.singpt3,c7][[1]]
UnivDeg4with2ABbothc0 = UnivDeg4with2ABboth/.(Solve[c0==UnivDeg4with2ABboth/.singpt3,c7][[1]]);
(* these will be the homogenous coordinates on the space of f *)
Varallc0={c0,c1,c2,c3,c4,c12,c13,c14,c15,c16,c17,c18,c19,d1,d2,d3,d4,d5,d6};
MatrixRank[Table[N[Coefficient[UnivDeg4with2ABbothc0,Varallc0[[k]]]/.SetOfPointsW[[l]],300],{k,1,13},{l,1,30}]]

This code rechecks that the 13 equations are vanishing to the 2A+B at the singular points in an appropriate way.
It also calculates their restrictions to the lines, and the restriction of the C3 invariant equation that vanishes at singular points.
We will not try to have our parameters on the exceptional P^1 set up so as to have the intersection point at y1=infty,
i.e. we are a Mobius transformation away from the previous approach.

UnivSymmc0=UnivDeg4with2ABbothc0 (UnivDeg4with2ABbothc0/.{W1->W2,W2->W3,W3->W1,W4->W5,W5->W6,W6->W4})(UnivDeg4with2ABbothc0/.{W1->W3,W2->W1,W3->W2,W4->W6,W5->W4,W6->W5});
InvDeg2through={W0^2, W0(W1+W2+W3),W1^2+W2^2+W3^2-(W4^2+W5^2+W6^2), W1 W2 +W1 W3 + W2 W3, W7^2, W7( W4 + W5 + W6)};
UnivInvDeg2through={d1,d2,d3,d4,d5,d6}.InvDeg2through;

(* Now getting restrictions *)
EqsW5=Table[AllEqs23Wgood[[k]],{k,1,5}];
VarT={T0,T2,T3,T4,T5,T6,T7};
VarT1st = VarT;
TtoW ={W0->T0,W1->1,W2->T0 T2 ,W3->T0 T3,W4->1+T0 T4 ,W5->T0 T5 ,W6->T0 T6,W7->T0 T7 };
RelsT=Table[Sum[SeriesCoefficient[(AllEqs23Wgood[[j]]/.TtoW),{T0,0,k}] T0^(k-1),{k,1,5}],{j,1,5}];
VarT2={T2,T3,T4,T5,T6,T7};

solsingeq=FullSimplify[Solve[Table[0==Coefficient[{1,a2,a3,a4,a5}.RelsT/.T0->0,VarT2[[k]]],{k,1,6}]]]
singeq =Simplify[({1,a2,a3,a4,a5}.EqsW5/.solsingeq[[1]])];
singeqT=Simplify[singeq/.TtoW];
solsexclin=Simplify[Solve[Table[0==RelsT[[j]]/.T0->0,{j,1,5}]][[1]]];
sollines=Simplify[Solve[SeriesCoefficient[(singeqT/.solsexclin),{T0,0,2}]==0,T2]];
Line2=sollines[[1]]
Line1=sollines[[2]]

(* writing a parametric equation of the first line with T3 as the parameter *)
param=FullSimplify[Table[VarT[[k]]->(VarT[[k]]/.solsexclin/.Line1/.T0->0/.T3->y1),{k,1,7}]];
(* the singular equation has double vanishing *)
RelsTplus=Append[RelsT,Simplify[( singeqT/T0^2)]];
Simplify[RelsTplus/.param] (* get zeros as expected *)

(* first infinitesimal neighborhood of the first exceptional line *)
param1 = {T0->y2, T2->(T2/.param) + t2 y2, T3-> y1, T4->(T4/.param) + t4 y2, T5->(T5/.param) + t5 y2, T6->(T6/.param) + t6 y2, T7->(T7/.param) + t7 y2};
solvenext=FullSimplify[Solve[Table[0==SeriesCoefficient[RelsTplus[[k]]/.param1, {y2,0,1}],{k,1,6}],{t2,t4,t5,t6,t7}]][[1]];
param1solved=param1/.solvenext;

online1=Simplify[SeriesCoefficient[(UnivDeg4with2ABbothc0/.TtoW)/.param1solved,{y2,0,2}]]; (* degree three numerator and degree one denominator in y1 *)
Num=Numerator[Factor[online1]];
SeriesCoefficient[Num,{y1,0,4}] (* should be 0 *)
online1good=Sum[RootReduce[SeriesCoefficient[Num,{Varallc0[[k]],0,1},{ y1,0,l}]/Coefficient[Num,c0 y1^3]] Varallc0[[k]]y1^l, {l,0,3},{k,1,19}];
UnivInvDeg2throughLine1=Simplify[SeriesCoefficient[(UnivInvDeg2through/.TtoW)/.param1solved,{y2,0,1}]];(* is exact, degree one in y1, no denominator *)

(* now the second line *)
param=Table[VarT[[k]]->(VarT[[k]]/.solsexclin/.Line2/.T0->0/.T3->y1),{k,1,7}];
RelsTplus=Append[RelsT,Simplify[( singeqT/T0^2)]]; 
Simplify[RelsTplus/.param](* get zeros as expected *)

(* first infinitesimal neighborhood of the second exceptional line *)
param1 = {T0->y2, T2->(T2/.param) + t2 y2, T3-> y1, T4->(T4/.param) + t4 y2, T5->(T5/.param) + t5 y2, T6->(T6/.param) + t6 y2, T7->(T7/.param) + t7 y2};
solvenext=FullSimplify[Solve[Table[0==SeriesCoefficient[RelsTplus[[k]]/.param1, {y2,0,1}],{k,1,6}],{t2,t4,t5,t6,t7}]][[1]];
param1solved=param1/.solvenext;

(* to get the opposite line at the other point, we simply cycle the variables; no need to do it for the C3 invariant ones *)
online2=Simplify[SeriesCoefficient[((UnivDeg4with2ABbothc0/.{W1->W3,W2->W1,W3->W2,W4->W6,W5->W4,W6->W5})/.TtoW)/.param1solved,{y2,0,2}]];
Num=Numerator[Factor[online2]];
SeriesCoefficient[Num,{y1,0,4}] (* should be 0 *)
online2good=Sum[RootReduce[SeriesCoefficient[Num,{Varallc0[[k]],0,1},{ y1,0,l}]/Coefficient[Num,c0 y1^3]] Varallc0[[k]]y1^l, {l,0,3},{k,1,19}];
UnivInvDeg2throughLine2=Simplify[SeriesCoefficient[(UnivInvDeg2through/.TtoW)/.param1solved,{y2,0,1}]];

3. Find/verify the basis of the C3 invariant sections of 12K_FPP/C3. Set up 67 cubic equations. Reduce to 58 equations.

(* this is a day long calculation which is best not to repeat *)
deg=6;
VarsW={W0,W1,W2,W3,W4,W5,W6,W7};
AllMons=MonomialList[Sum[(W0+W1+W2+W3)^(2k) (W4+W5+W6+W7)^(deg-2k),{k,0,deg/2}]];
count=Length[AllMons]
AllMons=Table[AllMons[[k]]/(AllMons[[k]]/.Table[VarsW[[l]]->1,{l,1,8}]),{k,1,count}];

Basis12H={};
MatrixValues={};
rk=0;
For[i=1,i≤count,i++,
newbasiselt=AllMons[[i]]+(AllMons[[i]]/.{W1->W2,W2->W3,W3->W1,W4->W5,W5->W6,W6->W4})+(AllMons[[i]]/.{W1->W3,W2->W1,W3->W2,W4->W6,W5->W4,W6->W5});
newrow=Table[newbasiselt/.SetOfPointsW[[k]],{k,1,Length[SetOfPointsW]}];
If[MatrixRank[Table[Append[MatrixValues,newrow][[j]][[k]],{k,1,100},{j,1,rk+1}]]==rk+1,
rk++;
Print[rk," ",i];
AppendTo[Basis12H,newbasiselt];
AppendTo[MatrixValues,newrow];
];
];

Basis12H

(* This is a basis of the space of even C3-symmetric sections of degree 6 on X. The dimension of the space of even degree 6 sections is 12K 11K/2 + 1 = 66 K^2 + 1 = 66*3 + 1. The trace of the generator of C3 is 1, so we get 67. *)
Basis12H={3 (W0^6),(W0^5) W1+(W0^5) W2+(W0^5) W3,(W0^4) (W1^2)+(W0^4) (W2^2)+(W0^4) (W3^2),(W0^4) W1 W2+(W0^4) W1 W3+(W0^4) W2 W3,(W0^4) (W4^2)+(W0^4) (W5^2)+(W0^4) (W6^2),(W0^4) W4 W5+(W0^4) W4 W6+(W0^4) W5 W6,(W0^4) W4 W7+(W0^4) W5 W7+(W0^4) W6 W7,(W0^3) (W1^3)+(W0^3) (W2^3)+(W0^3) (W3^3),(W0^3) (W1^2) W2+(W0^3) (W2^2) W3+(W0^3) W1 (W3^2),(W0^3) W1 (W2^2)+(W0^3) (W1^2) W3+(W0^3) W2 (W3^2),3 (W0^3) W1 W2 W3,(W0^3) W1 (W4^2)+(W0^3) W2 (W5^2)+(W0^3) W3 (W6^2),(W0^3) W1 W4 W5+(W0^3) W3 W4 W6+(W0^3) W2 W5 W6,(W0^3) W2 W4 W5+(W0^3) W1 W4 W6+(W0^3) W3 W5 W6,(W0^3) W1 W4 W7+(W0^3) W2 W5 W7+(W0^3) W3 W6 W7,(W0^3) W3 (W4^2)+(W0^3) W1 (W5^2)+(W0^3) W2 (W6^2),(W0^2) (W1^4)+(W0^2) (W2^4)+(W0^2) (W3^4),(W0^2) (W1^3) W2+(W0^2) (W2^3) W3+(W0^2) W1 (W3^3),(W0^2) W1 (W2^3)+(W0^2) (W1^3) W3+(W0^2) W2 (W3^3),(W0^2) (W1^2) (W2^2)+(W0^2) (W1^2) (W3^2)+(W0^2) (W2^2) (W3^2),(W0^2) (W1^2) W2 W3+(W0^2) W1 (W2^2) W3+(W0^2) W1 W2 (W3^2),(W0^2) (W1^2) (W4^2)+(W0^2) (W2^2) (W5^2)+(W0^2) (W3^2) (W6^2),(W0^2) (W1^2) W4 W5+(W0^2) (W3^2) W4 W6+(W0^2) (W2^2) W5 W6,(W0^2) (W2^2) W4 W5+(W0^2) (W1^2) W4 W6+(W0^2) (W3^2) W5 W6,(W0^2) (W1^2) W4 W7+(W0^2) (W2^2) W5 W7+(W0^2) (W3^2) W6 W7,(W0^2) (W3^2) (W4^2)+(W0^2) (W1^2) (W5^2)+(W0^2) (W2^2) (W6^2),(W0^2) (W3^2) W4 W5+(W0^2) (W2^2) W4 W6+(W0^2) (W1^2) W5 W6,(W0^2) (W3^2) W4 W7+(W0^2) (W1^2) W5 W7+(W0^2) (W2^2) W6 W7,(W0^2) W1 W2 (W4^2)+(W0^2) W2 W3 (W5^2)+(W0^2) W1 W3 (W6^2),W0 (W1^5)+W0 (W2^5)+W0 (W3^5),W0 (W1^4) W2+W0 (W2^4) W3+W0 W1 (W3^4),W0 W1 (W2^4)+W0 (W1^4) W3+W0 W2 (W3^4),W0 (W1^3) (W2^2)+W0 (W2^3) (W3^2)+W0 (W1^2) (W3^3),W0 (W1^3) W2 W3+W0 W1 (W2^3) W3+W0 W1 W2 (W3^3),W0 (W1^2) (W2^3)+W0 (W1^3) (W3^2)+W0 (W2^2) (W3^3),W0 (W1^3) (W4^2)+W0 (W2^3) (W5^2)+W0 (W3^3) (W6^2),W0 (W1^3) W4 W5+W0 (W3^3) W4 W6+W0 (W2^3) W5 W6,W0 (W2^3) W4 W5+W0 (W1^3) W4 W6+W0 (W3^3) W5 W6,W0 (W1^3) W4 W7+W0 (W2^3) W5 W7+W0 (W3^3) W6 W7,W0 (W3^3) (W4^2)+W0 (W1^3) (W5^2)+W0 (W2^3) (W6^2),W0 (W3^3) W4 W5+W0 (W2^3) W4 W6+W0 (W1^3) W5 W6,W0 (W3^3) W4 W7+W0 (W1^3) W5 W7+W0 (W2^3) W6 W7,W0 (W1^2) (W2^2) W3+W0 (W1^2) W2 (W3^2)+W0 W1 (W2^2) (W3^2),W0 (W1^2) W2 (W4^2)+W0 (W2^2) W3 (W5^2)+W0 W1 (W3^2) (W6^2),W0 (W1^2) W2 W4 W5+W0 W1 (W3^2) W4 W6+W0 (W2^2) W3 W5 W6,W0 (W2^2) W3 W4 W5+W0 (W1^2) W2 W4 W6+W0 W1 (W3^2) W5 W6,(W1^6)+(W2^6)+(W3^6),(W1^5) W2+(W2^5) W3+W1 (W3^5),W1 (W2^5)+(W1^5) W3+W2 (W3^5),(W1^4) (W2^2)+(W2^4) (W3^2)+(W1^2) (W3^4),(W1^4) W2 W3+W1 (W2^4) W3+W1 W2 (W3^4),(W1^2) (W2^4)+(W1^4) (W3^2)+(W2^2) (W3^4),(W1^4) (W4^2)+(W2^4) (W5^2)+(W3^4) (W6^2),(W1^4) W4 W5+(W3^4) W4 W6+(W2^4) W5 W6,(W2^4) W4 W5+(W1^4) W4 W6+(W3^4) W5 W6,(W1^4) W4 W7+(W2^4) W5 W7+(W3^4) W6 W7,(W3^4) (W4^2)+(W1^4) (W5^2)+(W2^4) (W6^2),(W3^4) W4 W5+(W2^4) W4 W6+(W1^4) W5 W6,(W3^4) W4 W7+(W1^4) W5 W7+(W2^4) W6 W7,(W1^3) (W2^3)+(W1^3) (W3^3)+(W2^3) (W3^3),(W1^3) (W2^2) W3+W1 (W2^3) (W3^2)+(W1^2) W2 (W3^3),(W1^2) (W2^3) W3+(W1^3) W2 (W3^2)+W1 (W2^2) (W3^3),(W1^3) W2 (W4^2)+(W2^3) W3 (W5^2)+W1 (W3^3) (W6^2),(W1^3) W2 W4 W5+W1 (W3^3) W4 W6+(W2^3) W3 W5 W6,(W2^3) W3 W4 W5+(W1^3) W2 W4 W6+W1 (W3^3) W5 W6,(W1^3) W2 W4 W7+(W2^3) W3 W5 W7+W1 (W3^3) W6 W7,W1 (W3^3) (W4^2)+(W1^3) W2 (W5^2)+(W2^3) W3 (W6^2)};

MatrixRank[Table[Basis12H[[k]]/.SetOfPointsW[[l]],{k,1,67},{l,1,100}]] (* get 67 as expected *)

Accuracy[SetOfPointsW]
Beautify[x_]:=
Ret=Expand[RootApproximant[N[Re[x],4000],2]+RootApproximant[N[Im[x]/(Sqrt[15]),4000],2] I(Sqrt[15])];

(* calculating the coefficients of f sigma(f) sigma^2(f) in the above basis; this takes a few hours *)
CoeffList=Table[ToExpression[StringJoin["A",ToString[j]]],{j,1,67}];
UnivPoly=Sum[CoeffList[[j]] Basis12H[[j]],{j,1,67}];
Deg=3;
Deg3contrib=0;
AllMons=MonomialList[Sum[Varallc0[[k]],{k,1,13}]^3];
AllMons=Table[AllMons[[i]]/(AllMons[[i]]/.Table[Varallc0[[k]]->1,{k,1,Length[Varallc0]}]),{i,1,Length[AllMons]}];
Print[Length[AllMons]];
For[num1=1,num1≤Length[AllMons],num1++,
 If[Mod[num1,50]==0,Print[num1]];
tosym=Coefficient[UnivSymmc0,AllMons[[num1]]];
SymMon=tosym+(tosym/.{W1->W2,W2->W3,W3->W1,W4->W5,W5->W6,W6->W4})+(tosym/.{W1->W3,W2->W1,W3->W2,W4->W6,W5->W4,W6->W5});
Eqs=Table[0==(SymMon-UnivPoly)/.SetOfPointsW[[j]],{j,1,67}]; 
sols=Chop[Solve[Eqs][[1]]];
newproj=Table[Beautify[CoeffList[[j]]/.sols],{j,1,67}];
Deg3contrib+=newproj AllMons[[num1]];
];

UnivSymInBasis=Deg3contrib; (* the coefficients are large *)
Table[ByteCount[UnivSymInBasis[[i]]],{i,1,67}]

(* the same thing for UnivInvDeg2through^3 *)
Deg=3;
InvCubeInBasis=0;
AllMons=MonomialList[Sum[Varallc0[[k]],{k,14,19}]^3];
AllMons=Table[AllMons[[i]]/(AllMons[[i]]/.Table[Varallc0[[k]]->1,{k,1,Length[Varallc0]}]),{i,1,Length[AllMons]}];
Print[Length[AllMons]];
For[num1=1,num1≤Length[AllMons],num1++,
If[Mod[num1,50]==0, Print[num1];];
tosym=Coefficient[UnivInvDeg2through^3,AllMons[[num1]]];
SymMon=tosym+(tosym/.{W1->W2,W2->W3,W3->W1,W4->W5,W5->W6,W6->W4})+(tosym/.{W1->W3,W2->W1,W3->W2,W4->W6,W5->W4,W6->W5});
Eqs=Table[0==(SymMon-UnivPoly)/.SetOfPointsW[[j]],{j,1,67}];
sols=Chop[Solve[Eqs][[1]]];
newproj=Table[Beautify[CoeffList[[j]]/.sols],{j,1,67}];
InvCubeInBasis+=newproj AllMons[[num1]];
];

(* sanity check *)
Chop[Simplify[N[(Sum[InvCubeInBasis[[i]] Basis12H[[i]],{i,1,67}] -3  UnivInvDeg2through^3)/.SetOfPointsW[[170]],2000]]]
Chop[Simplify[N[(Sum[UnivSymInBasis[[i]] Basis12H[[i]],{i,1,67}] -3  UnivSymmc0)/.SetOfPointsW[[170]],2000]]]

We can now set up the equations, in equation and just polynomial format. The “aton” refers to having too many of them.

EqsAton=Table[InvCubeInBasis[[i]]==UnivSymInBasis[[i]],{i,1,67}];
RelsAton=Table[InvCubeInBasis[[i]]-UnivSymInBasis[[i]],{i,1,67}];

(* we will find out what monomials are present in these equations *)
ExpList={};
For[i=1,i≤67,i++,
If[Mod[i,20]==0,Print[i," ",Length[ExpList]]];
ML=MonomialList[RelsAton[[i]]];
For[j=1,j≤Length[ML],j++,
newmon=ML[[j]];
newexp=Table[Exponent[newmon,Varallc0[[k]]],{k,1,19}];
ExpList=Union[ExpList,{newexp}];
];
];
MonsList=Table[Product[Varallc0[[k]]^(ExpList[[i]][[k]]),{k,1,19}],{i,1,Length[ExpList]}];
Length[MonsList]

(* we put the coefficients in a matrix *)
BigMatrix={};
For[i=1,i≤Length[RelsAton],i++,
If[Mod[i,15]==0,Print[i]];
AppendTo[BigMatrix,Table[RootReduce[Coefficient[RelsAton[[i]],MonsList[[k]]]],{k,1,Length[MonsList]}]];
];
Print["The rank of the matrix is ",MatrixRank[N[BigMatrix,200]]];

Now we will clean up the matrix.

ByteCount[BigMatrix]
For[k=1,k≤67,k++,
firsteq=BigMatrix[[k]];
lcm=Denominator[firsteq[[1]]];
For[i=2,i≤Length[firsteq],i++,
lcm=LCM[lcm,Denominator[firsteq[[i]]]];
];
BigMatrix[[k]]=Expand[firsteq lcm];
];
ByteCount[BigMatrix]

(* We write down the matrix in terms of real and imaginary parts *)
BigMatrixReIm={};
For[i=1,i≤67,i++,
NewRow1={};
NewRow2={};
For[j=1,j≤511,j++,
AppendTo[NewRow1,Re[BigMatrix[[i]][[j]]]];
AppendTo[NewRow1,Im[BigMatrix[[i]][[j]]]/Sqrt[15]];
AppendTo[NewRow2,Im[BigMatrix[[i]][[j]]](-Sqrt[15])];
AppendTo[NewRow2,Re[BigMatrix[[i]][[j]]]];
];
AppendTo[BigMatrixReIm,NewRow1];
AppendTo[BigMatrixReIm,NewRow2];
];
ByteCount[BigMatrixReIm]
BestRels=LatticeReduce[BigMatrixReIm];  (* takes awhile *)
ByteCount[BestRels]

Length[BestRels]
For[k=1,k≤116,k++,
firsteq=BestRels[[k]];
gcd=firsteq[[1]];
For[i=2,i≤Length[firsteq],i++,
gcd=GCD[gcd,Denominator[firsteq[[i]]]];
];
BestRels[[k]]=Expand[firsteq/gcd];
];
ByteCount[BestRels]

{u,r, v} = SmithDecomposition[BestRels]; (* this takes awhile, ten hours or so *)
vt=Transpose[v];
SmallBestRels={};
For[i=1,i≤116,i++,
If[Mod[i,10]==0,Print[i];];
basisvec=Table[0,{k,1,1022}];
basisvec[[i]]=1;
AppendTo[SmallBestRels,LinearSolve[vt,basisvec]];
];

ByteCount[SmallBestRels] 
SmallerBestRels=LatticeReduce[SmallBestRels];
ByteCount[SmallerBestRels]

BackToCxMatrix={};
For[i=1,i≤116,i++,
newrow=Table[SmallerBestRels[[i]][[2k-1]]+ I Sqrt[15] SmallerBestRels[[i]][[2k]],{k,1,511}];
AppendTo[BackToCxMatrix,newrow];
];
MatrixRank[N[BackToCxMatrix,400]]

BackToCxReduced={};
rk=0;
For[i=1,i≤116,i++,
newred=Append[BackToCxReduced,BackToCxMatrix[[i]]];
If[MatrixRank[N[newred,1000]]>rk, rk++;Print[i," ",rk];BackToCxReduced=newred;];
If[rk==58,i=120];
];

ByteCount[BackToCxReduced]
Length[BackToCxReduced]
Table[ByteCount[BackToCxReduced[[i]]],{i,1,58}]

EqsAton={};
For[i=1,i≤58,i++,
AppendTo[EqsAton,0==Sum[BackToCxReduced[[i]][[k]] MonsList[[k]],{k,1,Length[MonsList]}]];
];
RelsAton=Table[EqsAton[[i]][[2]],{i,1,58}];

fname="58CubicRels";
OpenWrite[fname];
Write[fname, RelsAton];
Close[fname];

fname="58CubicRels";
OpenRead[fname];
RelsAton=Read[fname];
Close[fname];

4. Make a choice and set up three linear equations. Simplify 58 equations to reduce to 55 cubic equations.

w=Exp[2 Pi I /3]
solfixed=NSolve[Join[Table[0==AllEqs23Wgood[[k]],{k,1,9}],{W0==1,W1==W2,W1==W3,W4==W5,W5==W6}],300]; (* find fixed points of C3 action *)
N[W1/.solfixed]
solfixedexact=Table[VarsW[[k]]->Beautify[VarsW[[k]]/.solfixed[[i]]],{i,1,Length[solfixed]},{k,1,8}] (* exact values *)

(* we will be solving the linear equations for c1, c2, c3 *)
Varallred={c0,c4,c12,c13,c14,c15,c16,c17,c18,c19,d1,d2,d3,d4,d5,d6};
Varsolved={c1,c2,c3};


UnivSymmSum=UnivDeg4with2ABbothc0 +(UnivDeg4with2ABbothc0/.{W1->W2,W2->W3,W3->W1,W4->W5,W5->W6,W6->W4})+(UnivDeg4with2ABbothc0/.{W1->W3,W2->W1,W3->W2,W4->W6,W5->W4,W6->W5}); (* f+ sigma(f) + sigma^2(f) *)

(* we will be making equations that are defined over the Q(sqrt[-15]) *)
(* this takes longer than expected *)
LinRel1=Sum[RootReduce[Coefficient[(3  UnivInvDeg2through - UnivSymmSum)/.solfixedexact[[6]], Varallc0[[k]]]]Varallc0[[k]],{k,1,19}];
LinRel2=Sum[RootReduce[Coefficient[((3 w UnivInvDeg2through - UnivSymmSum)/.solfixedexact[[3]]) +((3 w^2 UnivInvDeg2through - UnivSymmSum)/.solfixedexact[[4]]), Varallc0[[k]]]] Varallc0[[k]],{k,1,19}];
LinRel3=Sum[RootReduce[Coefficient[((3 w UnivInvDeg2through - UnivSymmSum)/.solfixedexact[[3]])  w + w^2((3 w^2 UnivInvDeg2through - UnivSymmSum)/.solfixedexact[[4]]),Varallc0[[k]]]] Varallc0[[k]],{k,1,19}];
LinRel4=Sum[RootReduce[Coefficient[((3 w^2 UnivInvDeg2through - UnivSymmSum)/.solfixedexact[[3]]) +((3 w UnivInvDeg2through - UnivSymmSum)/.solfixedexact[[4]]), Varallc0[[k]]]] Varallc0[[k]],{k,1,19}];
LinRel5=Sum[RootReduce[Coefficient[((3 w^2 UnivInvDeg2through - UnivSymmSum)/.solfixedexact[[3]])  w ^2+ w((3 w UnivInvDeg2through - UnivSymmSum)/.solfixedexact[[4]]), Varallc0[[k]]]] Varallc0[[k]],{k,1,19}];

There are two possible choices of linear equations.

LinEqOne={0==LinRel1,0==LinRel2,0==LinRel3};
LinEqTwo={0==LinRel1,0==LinRel4,0==LinRel5};

We were successful when using LinEqTwo.

(* these can be reordered as needed *)
sols=Solve[LinEqOne,Varsolved][[1]];
sols=Solve[LinEqTwo,Varsolved][[1]];



solsclean={};
For[k=1,k≤3,k++,
clean=Sum[RootReduce[Coefficient[Varsolved[[k]]/.sols, Varallred[[i]]]]Varallred[[i]],{i,1,16}];
AppendTo[solsclean,Varsolved[[k]]->clean];
];
N[solsclean]
RelsAtonred={};
For[i=1,i≤58,i++,
AppendTo[RelsAtonred,RelsAton[[i]]/.solsclean];
];

(* repeat the procedure for finding 58 equations *)
ExpListRed={};
For[i=1,i≤58,i++,
If[Mod[i,10]==0,Print[i," ",Length[ExpListRed]]];
ML=MonomialList[RelsAtonred[[i]]];
For[j=1,j≤Length[ML],j++,
newmon=ML[[j]];
newexp=Table[Exponent[newmon,Varallred[[k]]],{k,1,16}];
ExpListRed=Union[ExpListRed,{newexp}];
];
];
Length[ExpListRed]
MonsListRed=Table[Product[Varallred[[k]]^(ExpListRed[[i]][[k]]),{k,1,16}],{i,1,Length[ExpListRed]}];

BigMatrix={};
For[i=1,i≤58,i++,
If[Mod[i,5]==0,Print[i]];
AppendTo[BigMatrix,Table[RootReduce[Coefficient[RelsAtonred[[i]]  ,MonsListRed[[j]] ]],{j,1,816}]];
];
BM=N[BigMatrix,400];
ByteCount[BM]
MatrixRank[BM]

ByteCount[BigMatrix]
For[k=1,k≤58,k++,
firsteq=BigMatrix[[k]];
lcm=Denominator[firsteq[[1]]];
For[i=2,i≤Length[firsteq],i++,
lcm=LCM[lcm,Denominator[firsteq[[i]]]];
];
BigMatrix[[k]]=Expand[firsteq lcm];
];

ByteCount[BigMatrix]
BigMatrixReIm={};
For[i=1,i≤58,i++,
NewRow1={};
NewRow2={};
For[j=1,j≤816,j++,
AppendTo[NewRow1,Re[BigMatrix[[i]][[j]]]];
AppendTo[NewRow1,Im[BigMatrix[[i]][[j]]]/Sqrt[15]];
AppendTo[NewRow2,Im[BigMatrix[[i]][[j]]](-Sqrt[15])];
AppendTo[NewRow2,Re[BigMatrix[[i]][[j]]]];
];
AppendTo[BigMatrixReIm,NewRow1];
AppendTo[BigMatrixReIm,NewRow2];
];

ByteCount[BigMatrixReIm]
BestRels=LatticeReduce[BigMatrixReIm];
ByteCount[BestRels]

{u,r, v} = SmithDecomposition[BestRels];
vt=Transpose[v];
SmallBestRels={};
For[i=1,i≤110,i++,
If[Mod[i,10]==0,Print[i];];
basisvec=Table[0,{k,1,1632}];
basisvec[[i]]=1;
AppendTo[SmallBestRels,LinearSolve[vt,basisvec]];
];

SmallerBestRels=LatticeReduce[SmallBestRels];
Length[SmallerBestRels]
ByteCount[SmallerBestRels]
SmallerBestRels=LatticeReduce[SmallerBestRels];
Length[SmallerBestRels]
ByteCount[SmallerBestRels]
Table[ByteCount[SmallerBestRels[[i]]],{i,1,110}]


BackToCxMatrix={};
For[i=1,i≤110,i++,
newrow=Table[SmallerBestRels[[i]][[2k-1]]+ I Sqrt[15] SmallerBestRels[[i]][[2k]],{k,1,816}];
AppendTo[BackToCxMatrix,newrow];
];

BackToCxReduced={};
rk=0;
For[i=1,i≤110,i++,
newred=Append[BackToCxReduced,BackToCxMatrix[[i]]];
If[MatrixRank[N[newred,1000]]>rk, rk++;Print[i," ",rk];BackToCxReduced=newred;];
If[rk==55,i=120];
];

EqsAtonRed={};
RelsAtonRed={};
For[i=1,i≤55,i++,
AppendTo[EqsAtonRed,0==Sum[BackToCxReduced[[i]][[k]] MonsListRed[[k]],{k,1,Length[MonsListRed]}]];
AppendTo[RelsAtonRed,Sum[BackToCxReduced[[i]][[k]] MonsListRed[[k]],{k,1,Length[MonsListRed]}]];
];
fname="55CubicEqs";
OpenWrite[fname];
Write[fname, RelsAtonRed];
Close[fname];

5. Find 12 quadratic equations, simplify them. Write out the equations in Magma format.

a0=SeriesCoefficient[online1good, {y1,0,0}];
a1=SeriesCoefficient[online1good,{y1,0,1}];
a2=SeriesCoefficient[online1good, {y1,0,2}];
a3=SeriesCoefficient[online1good,{y1,0,3}];
b0=SeriesCoefficient[UnivInvDeg2throughLine1, {y1,0,0}];
b1=SeriesCoefficient[UnivInvDeg2throughLine1, {y1,0,1}];
EqsOnLine1={3a3 b0 - a2 b1, a2 b0 - a1 b1, a1 b0 - 3 a0 b1, 9 a0 a3 - a1 a2, 3 a0 a2 - a1^2, 3 a1 a3 - a2^2};

a0=SeriesCoefficient[online2good, {y1,0,0}];
a1=SeriesCoefficient[online2good,{y1,0,1}];
a2=SeriesCoefficient[online2good, {y1,0,2}];
a3=SeriesCoefficient[online2good,{y1,0,3}];
b0=SeriesCoefficient[UnivInvDeg2throughLine2, {y1,0,0}];
b1=SeriesCoefficient[UnivInvDeg2throughLine2, {y1,0,1}];
EqsOnLine2={3a3 b0 - a2 b1, a2 b0 - a1 b1, a1 b0 - 3 a0 b1, 9 a0 a3 - a1 a2, 3 a0 a2 - a1^2, 3 a1 a3 - a2^2};

EqsOnLines=Join[EqsOnLine1,EqsOnLine2]/.solsclean;

ExpListDeg2={};
For[i=1,i≤12,i++,
ML=MonomialList[EqsOnLines[[i]]];
For[j=1,j≤Length[ML],j++,
newmon=ML[[j]];
newexp=Table[Exponent[newmon,Varallred[[k]]],{k,1,16}];
ExpListDeg2=Union[ExpListDeg2,{newexp}];
];
];
Length[ExpListDeg2]
MonsListDeg2=Table[Product[Varallred[[k]]^(ExpListDeg2[[i]][[k]]),{k,1,16}],{i,1,Length[ExpListDeg2]}];

BigMatrix={};
For[i=1,i≤12,i++,
AppendTo[BigMatrix,Table[RootReduce[Coefficient[EqsOnLines[[i]]  ,MonsListDeg2[[j]] ]],{j,1,136}]];
];

ByteCount[BigMatrix]
For[k=1,k≤12,k++,
firsteq=BigMatrix[[k]];
lcm=Denominator[firsteq[[1]]];
For[i=2,i≤Length[firsteq],i++,
lcm=LCM[lcm,Denominator[firsteq[[i]]]];
];
BigMatrix[[k]]=Expand[firsteq lcm];
];

ByteCount[BigMatrix]
BigMatrixReIm={};
For[i=1,i≤12,i++,
NewRow1={};
NewRow2={};
For[j=1,j≤136,j++,
AppendTo[NewRow1,Re[BigMatrix[[i]][[j]]]];
AppendTo[NewRow1,Im[BigMatrix[[i]][[j]]]/Sqrt[15]];
AppendTo[NewRow2,Im[BigMatrix[[i]][[j]]](-Sqrt[15])];
AppendTo[NewRow2,Re[BigMatrix[[i]][[j]]]];
];
AppendTo[BigMatrixReIm,NewRow1];
AppendTo[BigMatrixReIm,NewRow2];
];

ByteCount[BigMatrixReIm]
BestRels=LatticeReduce[BigMatrixReIm];
ByteCount[BestRels]

{u,r, v} = SmithDecomposition[BestRels];
vt=Transpose[v];
SmallBestRels={};
For[i=1,i≤24,i++,
Print[i];
basisvec=Table[0,{k,1, 272}];
basisvec[[i]]=1;
AppendTo[SmallBestRels,LinearSolve[vt,basisvec]];
];

SmallerBestRels=LatticeReduce[SmallBestRels];
Length[SmallerBestRels]
ByteCount[SmallerBestRels]

BackToCxMatrix={};
For[i=1,i≤24,i++,
newrow=Table[SmallerBestRels[[i]][[2k-1]]+ I Sqrt[15] SmallerBestRels[[i]][[2k]],{k,1,136}];
AppendTo[BackToCxMatrix,newrow];
];

BackToCxReduced={};
rk=0;
For[i=1,i≤24,i++,
newred=Append[BackToCxReduced,BackToCxMatrix[[i]]];
If[MatrixRank[N[newred,1000]]>rk, rk++;Print[i," ",rk];BackToCxReduced=newred;];
];

EqsDeg2Red={};
RelsDeg2Red={};
For[i=1,i≤12,i++,
AppendTo[EqsDeg2Red,0==Sum[BackToCxReduced[[i]][[k]] MonsListDeg2[[k]],{k,1,Length[MonsListDeg2]}]];
AppendTo[RelsDeg2Red,Sum[BackToCxReduced[[i]][[k]] MonsListDeg2[[k]],{k,1,Length[MonsListDeg2]}]];
];

ByteCount[EqsDeg2Red]
Table[ByteCount[EqsDeg2Red[[i]]],{i,1,12}]

fname="12QuadEqs";
OpenWrite[fname];
Write[fname,EqsDeg2Red];
Close[fname];

Rels23Withrr={};
For[i=1,i≤12,i++,
Print[i];
newrel=0;
For[j=1,j≤Length[MonsListDeg2],j++,
cf=Coefficient[RelsDeg2Red[[i]] ,MonsListDeg2[[j]]];
newrel+=(Re[cf]+Im[cf] (rr/Sqrt[15])) MonsListDeg2[[j]] ;
];
AppendTo[Rels23Withrr,newrel];
];

For[i=1,i≤55,i++,
Print[i];
newrel=0;
For[j=1,j≤Length[MonsListRed],j++,
cf=Coefficient[RelsAtonRed[[i]],MonsListRed[[j]]];
newrel+=(Re[cf]+Im[cf] (rr/Sqrt[15])) MonsListRed[[j]];
];
AppendTo[Rels23Withrr,newrel];
];

fname="Rels23withrr";
OpenWrite[fname];
Write[fname,Rels23Withrr];
Close[fname];

6. Solve the equations to find f and d.

fname="Rels23withrr";
OpenRead[fname];
Rels23Withrr=Read[fname];
Close[fname];

Varallred={c0,c4,c12,c13,c14,c15,c16,c17,c18,c19,d1,d2,d3,d4,d5,d6};

Point19=Solve[{0==18*c0+c4,0==14*c0+c12,0==14*c0+c13,0==13*c0+c14,0==14*c0+c15,0==13*c0+c16,0==3*c0+c17,0==16*c0+c18,0==11*c0+c19,0==15*c0+d1,0==6*c0+d2,0==3*c0+d3,0==16*c0+d4,0==10*c0+d5,0==13*c0+d6, 1==d1},Modulus->19][[1]]
Mod[Rels23Withrr/.rr->2/.Point19,19]

(* finding solutions modulo powers of 19 *)

NewPt={v1,v2,v3,v4,v5,v6,v7,v8,v9,v10,v11,v12,v13,v14,v15,v16};
Pt=Varallred/.Point19
rvalue=2;
maxdep=200;
For[dep=1,dep≤maxdep,dep++,
If[Mod[dep,20]==0,Print[dep," "];];
(* computing the square root of (-15) *)
For[i=0,i≤18,i++,
If[Mod[(rvalue+19^dep i)^2+15,19^(dep+1)]==0,rvalue+=19^dep i;i=20;]; 
];
(* trying on new point *)
Ptplus=Table[Varallred[[i]]->(Pt[[i]]+19^dep  NewPt[[i]]),{i,1,16}];
Eqs={0==v11}; (* we want to normalize d1=1 in the solution *)
For[i=1,i≤30,i++,
AppendTo[Eqs,0==Simplify[(Rels23Withrr[[i]]/.rr->rvalue/.Ptplus)/19^dep]]
];
sols=Solve[Eqs,Modulus->19];
Pt=Varallred/.Ptplus/.sols[[1]];
];

(* we will now try to translate from approximate p-adic to complex *)
smallsols={};
rootscx={};
For[i=1,i≤16,i++,
FourNumbers={Pt[[i]],1, Mod[rvalue ,19^(dep-1)], 19^(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[-15] smallsols[[i]][[3]])/(smallsols[[i]][[1]])];
];

rootscx
Rels23=Rels23Withrr/.rr->Sqrt[-15];
Simplify[Rels23/.Table[Varallred[[i]]->rootscx[[i]],{i,1,16}]] (* we got it! *)

This is simply to present the formula for d in the paper.

SectionD=UnivInvDeg2through/.Table[Varallred[[i]]->rootscx[[i]],{i,1,16}]

42318123032 (W0^2)+42318123032(((-23709-33355 I (Sqrt[15])) W0 (W1+W2+W3))/18758)+42318123032(((-369999-115101 I (Sqrt[15])) (W1 W2+W1 W3+W2 W3))/9379)+42318123032(((-134064-81144 I (Sqrt[15])) ((W1^2)+(W2^2)+(W3^2)-(W4^2)-(W5^2)-(W6^2)))/9379)+42318123032(((-60345558-90294750 I (Sqrt[15])) (W4+W5+W6) W7)/7043629)+42318123032(((457763819877-572077298835 I (Sqrt[15])) (W7^2))/42318123032)

This is to write the equations to a file, to be used in the next section.

SectionF=Factor[(UnivDeg4with2ABbothc0/.Solve[LinEqTwo,Varsolved][[1]])/.Table[Varallred[[i]]->rootscx[[i]],{i,1,16}]];
fname="Sectionsfandd";
OpenWrite[fname];
Write[fname,{SectionF,SectionD}];
Close[fname];