In this section, we will try to construct the FPP from points of Z

Us = {U0, U1, U2, U3, U4, U5, U6, U7, U8, U9, U10, U11, U12, U13, U14, U15, U16, U17, U18, U19};
Ps = {P0, P1, P2, P3, P4, P5, P6, P7, P8, P9, P10, P11, P12, P13, P14, P15, P16, P17, P18, P19};
Zs = {Z0, Z1, Z2, Z3, Z4, Z5, Z6, Z7, Z8, Z9, Z10, Z11, Z12};

PsToUs = {};
Do [
AppendTo[PsToUs, Ps[[i]]->Us[[i]]];
, {i, 1, 10}];

Do [
AppendTo[PsToUs, Ps[[i]]->Us[[i]]/(U10)^(1/2)];
, {i, 11, 20}];

PsToUs

Mod7Weights={0,6,5,3,1 ,2,4,1,2,4, 0,6,5,3,1 ,2,4,1,2,4}; 

Order3Ps={P1->P2,P2->P3,P3->P1,P4->P5,P5->P6,P6->P4,P7->P8,P8->P9,P9->P7, 
P11->P12, P12->P13, P13->P11, P14->P15, P15->P16, P16->P14, P17->P18, P18->P19, P19->P17};

Order3Zs = {Z1 -> Z2, Z2->Z3, Z3->Z1, Z4->Z5, Z5->Z6, Z6->Z4, Z7->Z8, Z8->Z9, Z9->Z7, Z10->Z11, Z11->Z12, Z12->Z10};

C2AutZs = {Z0->-Z0,
Z1->Z4, Z2->Z5, Z3->Z6, 
Z4->Z1, Z5->Z2, Z6->Z3,
Z7->Z10, Z8->Z11, Z9->Z12,
Z10->Z7, Z11->Z8, Z12->Z9};

ZWeights = {{0, 0}, 
{1, 1}, {4, 2}, {2, 4},
{-1, 1}, {-4, 2}, {-2, 4},
{1, 3}, {4, 6}, {2, 5},
{-1, 3}, {-4, 6}, {-2, 5}
};

C2AutPs = {P10->-P10, P11->-P11, P12->-P12,P13->-P13,P14->-P14,P15->-P15,P16->-P16,P17->-P17,P18->-P18, P19->-P19};

(*Import files*)
SetDirectory[NotebookDirectory[]];
FileNames[]

filename = "Equations/Z_singular_equations.txt";
ZSingularEquations=ToExpression[Import[filename]]/.rr->(Sqrt[-7]);

filename = "Dependency/points_on_Z.txt";
ZPoints=ToExpression[Import[filename]][[320;;480]];

ToRR[x_]:=Ret=Rationalize[Chop[N[Re[x],1900],10^(-1800)]]+rr Rationalize[Chop[N[Im[x]/(Sqrt[7]),1900],10^(-1800)]];


(* this module finds simple linear compbinations of rows of matrices that are defined over quadratic imaginary field *)
NiceLinComb[CoeffMatrix_,root_,prec_]:=
Module[
{BackToCxReduced,CoeffMatrixReIm,i,j,k,u,r,v,SmallLinComb,rk,NewRow1,NewRow2,BackToCxMatrix},
(* Building the rational matrix of double the size by looking at real and imaginary parts of rows and  (I root) times rows. *)
CoeffMatrixReIm={};
For[i=1,i≤ Length[CoeffMatrix],i++,
NewRow1={};
NewRow2={};
For[j=1,j≤ Length[CoeffMatrix[[i]]],j++,
AppendTo[NewRow1,Re[CoeffMatrix[[i]][[j]]]];
AppendTo[NewRow1,Im[CoeffMatrix[[i]][[j]]]/root];
AppendTo[NewRow2,Im[CoeffMatrix[[i]][[j]]](-root)];
AppendTo[NewRow2,Re[CoeffMatrix[[i]][[j]]]];
];
AppendTo[CoeffMatrixReIm,NewRow1];
AppendTo[CoeffMatrixReIm,NewRow2];
];

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

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

(* 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]][[2k-1]]+ I root SmallLinComb[[i]][[2k]],{k,1,Length[SmallLinComb[[1]]]/2}],{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]]],prec]]>rk,
rk++;
BackToCxReduced=Append[BackToCxReduced,BackToCxMatrix[[i]]];
(* Print[i," ",rk];*)
If[rk==Length[CoeffMatrix],i=Length[BackToCxMatrix]+1];
];
];
BackToCxReduced
];

We first find the 7 C_14 invariant quadratics in terms of Zi’s

deg=2;

Powers={b1,b2,b3,b4,b5,b6,b7,b8,b9,b10, b11, b12, b13};
AllPowers=Solve[Join[Table[Powers[[k]]≥0,{k,1,Length[Powers]}],{Sum[Powers[[k]],{k,1,Length[Powers]}]==deg}],Integers];

count = 0;
Equations = {};
Equationsrr = {};

RankTable = ConstantArray[0, {7, 7}];

For[mod7A=1,mod7A≤6,mod7A++,
mod7B=0;
(*Print["Weight=",{mod7A, mod7B}];*)

CorrectPowers={};
Do[
If[Mod[Sum[ZWeights[[k]] Powers[[k]],{k,1,Length[ZWeights]}]/.AllPowers[[k]],7]=={mod7A, mod7B},
AppendTo[CorrectPowers,AllPowers[[k]]];
];
,{k,1,Length[AllPowers]}];

(* These are monomials of weight (mod7A, 0) *)
Monomials = Table[Product[Zs[[n]]^(Powers[[n]]/.CorrectPowers[[k]]),{n,1,Length[Zs]}],{k,1,Length[CorrectPowers]}];
(*Print[Monomials];*)

Do [
AppendTo[Equations,  mon + (mon/.C2AutZs)];
, {mon, Monomials}];
];

AppendTo[Equations, Z0^2];
C14Quadratics = DeleteDuplicates[Equations]

We want to find 10 linearly independent on Z cubics  F satisfying:
1. F has second weight 0
2. C_2(F) = -F
3. F vanishes on z0 = 0

7 of them should be Z0*(C_14 quadratics computed earlier).

(* We will look for a C3- and C2- invariant set of cubic monomials which are {C7} x {1} invariant. 
We also don't want to use Z0. *)

ListCoeffs=Table[ToExpression[StringJoin["q",ToString[i]]],{i,1,400}];

(* Useful cubic monomials are chosen to be C3 invariant and give C14 quadratics  - TBD *)
AllMonomials = Flatten[{{ (Z6^2) Z8,(Z5^2) Z7,(Z1^2) Z12,Z1 Z2 Z3,Z11 (Z3^2),Z10 (Z2^2),(Z4^2) Z9,Z4 Z5 Z6},
{Z4 (Z7^2),(Z3^2) Z8,Z2 Z5 Z7,Z2 Z4 Z6,(Z12^2) Z6,Z1 Z10 Z7},
{(Z2^2) Z7,Z1 Z5 Z6,Z11 Z12 Z7,(Z11^2) Z5,Z1 Z4 Z9,Z12 Z3 Z9},
{Z1 (Z7^2),Z1 Z2 Z6,Z12 Z7 Z8,Z12 (Z4^2),Z11 Z5 Z8,Z11 Z3 Z6},
{(Z10^2) Z4,Z3 Z4 Z5,Z10 Z11 Z9,(Z1^2) Z9,Z11 Z2 Z8,Z3 Z6 Z8},
{Z10 (Z5^2),Z2 Z3 Z4,Z10 Z8 Z9,Z2 (Z8^2),Z1 Z12 Z4,Z12 Z6 Z9},
{Z1 (Z10^2),Z11 (Z6^2),Z10 Z2 Z5,Z1 Z3 Z5,Z3 (Z9^2),Z10 Z4 Z7}}];

Length[AllMonomials];
cubic = AllMonomials.ListCoeffs[[1;;Length[AllMonomials]]];

(* We first add the relation that the cubic is anti-invariant with respect to C2 *)
cubicC2 = cubic/.C2AutZs;
C2AntiInvariance = Solve[Coefficient[-1*cubic - cubicC2, AllMonomials] == 0,{q27,q28,q29,q30,q31,q32,q33,q34,q35,q36,q37,q38,q40,q41,q42,q43,q44,q5,q6,q7,q8,q9}][[1]];

antiInvariantCubic = cubic/.C2AntiInvariance
coeffQs = Variables[Coefficient[antiInvariantCubic, AllMonomials]];

(* Here we impose the relation that this vanishes on z0. Our precomputed points were in Ps, so we will instead recompute them in Zs. We also precomputed them. *)

filename = "Dependency/points_on_z0_zs.txt";
Z0PointsZs=ToExpression[Import[filename]];

{Dimensions[Z0PointsZs],Accuracy[Z0PointsZs]}

(*

iter =10;
PointsOnZ0InZs={};
For[k=1,k≤iter,k++,
Print["Iteration ", k];

current = NSolve[Join[{Z0==0,Z1==RandomChoice[Join[-Range[100],Range[100]]],Z2==RandomChoice[Join[-Range[100],Range[100]]]},Table[0==ZSingularEquations[[i]],{i,1,Length[ZSingularEquations]}]],1100];

Print[Accuracy[current[[1;;10]]]];

If [Length[current] ≠ 0,
PointsOnZ = Join[PointsOnZ0InZs,current[[1;;10]]];
];
];
*)

AllRelations = {};
Do [
AppendTo[AllRelations, antiInvariantCubic/.pt];
, {pt, Z0PointsZs}];
Length[AllRelations];

coeffMatrix = {};
Do [
AppendTo[coeffMatrix, Coefficient[rel, coeffQs]];
, {rel,AllRelations}];

{Dimensions[coeffMatrix],MatrixRank[coeffMatrix],Accuracy[coeffMatrix]}

nullCoeff = NullSpace[coeffMatrix];
Accuracy[nullCoeff];
reducedCoeff =RowReduce[nullCoeff,Tolerance-> 10^(-500)];
Accuracy[reducedCoeff];
exactCoeff=RootApproximant[Chop[reducedCoeff,10^(-500)], 2];
exactCoeffinRR = ToRR[exactCoeff];


(* Here we translate the solution back to assignments of coefficients *)
AllAssignments = {};

Do [
assignment = {};
Do [
AppendTo[assignment, coeffQs[[i]]->row[[i]]];
, {i, 1, Length[coeffQs]}];

AppendTo[AllAssignments, assignment];
, {row, exactCoeffinRR}];

AllAssignments;
Length[AllAssignments];

resultCubics = {};

Do [
AppendTo[resultCubics, antiInvariantCubic/.assign];
, {assign, AllAssignments}];

Length[resultCubics]

(* We only need one new cubic, provided that its C3 translates complete the basis of sections. Trial and error shows that taking the second cubic works. *)

new=resultCubics[[2]]

TenCubicsC3 = {
(Z0^3), Z0 Z12 Z5+Z0 Z2 Z9, Z0 Z10 Z6+Z0 Z3 Z7, Z0 Z11 Z4+Z0 Z1 Z8, Z0 Z10 Z3+Z0 Z6 Z7, Z0 Z1 Z11+Z0 Z4 Z8, Z0 Z12 Z2+Z0 Z5 Z9,
new,
new/.Order3Zs,
new/.Order3Zs/.Order3Zs 
};

MatrixRank[TenCubicsC3/.rr->(Sqrt[-7])/.ZPoints] 
(* must be 10 *)

fname=OpenWrite["Equations/ten_cubics.txt"];
Write[fname,TenCubicsC3];
Close[fname];

(* Now we change rr back to square root of -7 *)
TenCubicsC3 = TenCubicsC3/.rr->(Sqrt[-7]);

Finally, we compute the equations of the new FPP

Qs ={Q0, Q1, Q2, Q3, Q4, Q5, Q6, Q7, Q8, Q9};
QsToZs = {};
Do [
AppendTo[QsToZs, Qs[[i]]->TenCubicsC3[[i]]];
, {i, 1, 10}];

(*Do [
Print[elt];
,{elt, QsToZs}];*)

Monomials=Flatten[Table[Qs[[i]]Qs[[j]]Qs[[k]],{i,1,10},{j,1,i},{k,1,j}]];
Length[Monomials]

TableOfValues=Table[Monomials[[k]]/.QsToZs/.ZPoints[[l]],{l,1,Length[ZPoints]},{k,1,Length[Monomials]}];

{Dimensions[TableOfValues],Accuracy[TableOfValues]}

result = NullSpace[TableOfValues];
Print["nullspace done"];
Print[Dimensions[result]];
Print[Accuracy[result]];

rowResult = RowReduce[result, Tolerance->10^(-400)];
Print["Rowreduce done"];
Print[Accuracy[rowResult]];

exactResult = RootApproximant[Chop[rowResult, 10^(-390)], 2];
rowResultrr =Table[ToRR[exactResult[[i]]],{i,1,Length[exactResult]}];
newEquationsrr = rowResultrr.Monomials;
newEquations = exactResult.Monomials;

ByteCount[newEquationsrr]

(* We use the simplicaton module *)

exactResultNice=NiceLinComb[exactResult,7^(1/2),100];
rowResultrr =Table[ToRR[exactResultNice[[i]]],{i,1,Length[exactResultNice]}];
newEquationsrr = rowResultrr.Monomials;
newEquations = exactResult.Monomials;

fname=OpenWrite["Equations/Target_FPP_simple.txt"];
Write[fname,newEquationsrr ];
Close[fname];
Length[newEquationsrr ]

Now we reduce the equations modulo p for verification purposes. 
We will choose a prime p such that square root of -7 and cubic roots of 1 both exist (for example p = 37) for some further simplifications.

(* We convert Qs to Ts as follows. Here w is thought of as a non-trivial cubic root of 1. *)
Ts={T0,T1,T2,T3,T4,T5,T6,T7,T8,T9};
QstoTs = {
Q0->T0,
Q1->T1 + T2 + T3,
Q2-> T1 + w T2 + w^2 T3,
Q3 -> T1 + w^2 T2 + w T3,
Q4->T4 + T5 + T6,
Q5-> T4 + w T5 + w^2 T6,
Q6-> T4 + w^2 T5 + w T6,
Q7->T7 + T8 + T9,
Q8-> T7 + w T8 + w^2 T9,
Q9-> T7 + w^2 T8 + w T9
};

Weight3Ts = {0, 0, 1, 2, 0, 1, 2, 0, 1, 2};
Order3Ts = {T2-> w T2, T3->w^2 T3, T5-> w T5, T6->w^2 T6,T8-> w T8, T9->w^2 T9};

FPPinTs =newEquationsrr/.QstoTs/.w^3->1;

(* Now we find the Degree 3 monomials in terms of Ts *)
deg=3;

Powers={b1,b2,b3,b4,b5,b6,b7,b8,b9, b10};
AllPowers=Solve[Join[Table[Powers[[k]]≥0,{k,1,Length[Powers]}],{Sum[Powers[[k]],{k,1,Length[Powers]}]==deg}],Integers];

Mod3Monomials = {};
Do [
Print[mod3];

CorrectPowers={};
Do[
If[Mod[Sum[Weight3Ts[[k]] Powers[[k]],{k,1,Length[Weight3Ts]}]/.AllPowers[[k]],3]==mod3,
AppendTo[CorrectPowers,AllPowers[[k]]];
];
,{k,1,Length[AllPowers]}];

Monomials=Table[Product[Ts[[n]]^(Powers[[n]]/.CorrectPowers[[k]]),{n,1,10}],{k,1,Length[CorrectPowers]}];
(*Print[Monomials];
Print[Length[Monomials]];*)

AppendTo[Mod3Monomials, Monomials];

, {mod3, 0, 2}];

(* We split the 84 equations in Ts to 3*84 equations separated by its C3 weight *)
ExpandedEquations = {};
Do [
weight0Eq = 0;
weight1Eq = 0;
weight2Eq = 0;

Do [
weight0Eq += Coefficient[eq, term]*term;
, {term, Mod3Monomials[[1]]}];

Do [
weight1Eq += Coefficient[eq, term]*term;
, {term, Mod3Monomials[[2]]}];

Do [
weight2Eq += Coefficient[eq, term]*term;
, {term, Mod3Monomials[[3]]}];

AppendTo[ExpandedEquations, {weight0Eq, weight1Eq, weight2Eq}];
, {eq, FPPinTs}];

(* We collect the coefficient matrix of the 252 equations. This matrix will have rank 84, and we will find the linearly independent set of 84 equations from this *)
CoeffMatrix = {};
Do [
row = Coefficient[eq, Flatten[Mod3Monomials]];
AppendTo[CoeffMatrix, row];
, {eq, Flatten[ExpandedEquations]}];


(* First we check that the matrix has rank 84. We will use a numerical approximation of this matrix to save computation time *)
approx = N[CoeffMatrix/.rr->(Sqrt[-7])/.w->Exp[2*π*I/3], 500];
MatrixRank[approx]

(* We precomputed where the 84 linearly independent vectors should be, we verify it here *)
input = Join[approx[[1;;80]], approx[[81;;82]], approx[[85;;85]], approx[[250;;250]]];
MatrixRank[input]

(* Here we reduce the equations modulo 37. *)
exactInput = Join[CoeffMatrix[[1;;80]], CoeffMatrix[[81;;82]], CoeffMatrix[[85;;85]], CoeffMatrix[[250;;250]]];

p = 37
reducedMatrix = {};
Do [
reducedRow = row/.{rr->17,w->10};
newRow = {};
Do [
item = Mod[Numerator[elt], p]/Mod[Denominator[elt], p];
AppendTo[newRow, item];
, {elt, reducedRow}];

AppendTo[reducedMatrix, newRow];
, {row, exactInput}];

Dimensions[reducedMatrix]

reducedEquations = reducedMatrix.Flatten[Mod3Monomials];
Length[reducedEquations]

(* We also compute 3 minors from the equations. *)
JacobianT =PolynomialMod[Factor[Table[D[reducedEquations[[i]],Ts[[k]]],{i,1,Length[reducedEquations]},{k,1,Length[Ts]}]], p]; 
Dimensions[JacobianT]

(* Minor 1 *)
Minor1 = JacobianT[[{12,24,35,45,59,72,83},{1,2,3,4,6,7,8}]];
Dimensions[Minor1]
M1 = Det[Minor1];

(* Minor 2 *)
Minor2 = JacobianT[[{11,23,36,48,52,71,81},{1,2,3,4,5,6,10}]];
Dimensions[Minor2]
M2= Det[Minor2];

(* Minor 3 *)
Minor3 = JacobianT[[{10,16,35,47,60,69,84},{1,2,3,4,7,8,9}]];
Dimensions[Minor3]
M3= Det[Minor3];

fname=OpenWrite["Equations/reduced_equations_37.txt"];
Write[fname, reducedEquations];
Close[fname];

fname=OpenWrite["Equations/minors_37.txt"];
Write[fname, {M1, M2, M3}];
Close[fname];