/* c9_example.m

   Compute cubic surface with prescribed Galois group.

   This file constructs an example for the cyclic Group of order 9.
   For other groups, just change the input data. */


/* Solution of the equation problem -- generic I_40 <> 0 case. */
ring_z4<x,y,z,w> := PolynomialRing(IntegerRing(),4);

/* Given a separable polynomial pol5 of degree 5, this function returns the cubic surface having as pentahedral coefficients the zeroes of pol5. */
penta_pol_2_gl := function(pol5)
 alg := quo<Parent(pol5) | pol5>;

 nf_r4 := PolynomialRing(alg,4);
 /* Basis of trace-zero-subspace */
 bm :=  BasisMatrix(Kernel(Transpose(Matrix(
                         [[Trace(alg.1^i) : i in [0..4]]] ))));
 bas := [&+ [alg.1^(i-1) * bm[j,i] : i in [1..5]] : j in [1..4]];
 lf := &+[bas[i] * nf_r4.i : i in [1..4]];
 kf := lf^3 * alg.1; /* The pentahedral coefficient */
 gl := &+[Monomials(a)[1] * Trace(Coefficients(a)[1]) : a in Terms(kf)];
 gl := PolynomialRing(RationalField(),4)!gl;
 gl2 := ring_z4!(gl * Denominator(gl));
 gl2 := gl2 div GCD(Coefficients(gl2));
 return gl2;
end function;

/* Given an invariant vector, this function returns the corresponding cubic surface described by 20 coefficients. */
inv_2_gl := function(inv2)
 inv3 := [inv2[i] * (inv2[5]^3)^i : i in [1..5]];
 suc,s5 := IsPower(inv3[5],8);
 s2 := inv3[4] / s5^6;  s4 := inv3[3] / s5^4;
 s1 := inv3[2] / s5^3;  s3 := (s4^2 - inv3[1])/( 4* s5);

 pol5 := Polynomial([-s5,s4,-s3,s2,-s1,1]);
 return penta_pol_2_gl(pol5), pol5;
end function;


/* Coble's 40 irrational invariants. */
minor_ijk := function(mat,i,j,k)
  return Determinant(Submatrix(mat,[i,j,k],[1,2,3]));
end function;

gamma_3 := function(mat,i,j,k,l,m,n)
 return [minor_ijk(mat,i,k,l),  minor_ijk(mat,j,k,l),
         minor_ijk(mat,k,m,n),  minor_ijk(mat,l,m,n),
         minor_ijk(mat,m,i,j),  minor_ijk(mat,n,i,j)];
end function;

gamma_2 := function(mat, i,j,k,l,m,n)
  return [minor_ijk(mat,i,j,k), minor_ijk(mat,l,m,n),
          minor_ijk(mat,3,4,1)* minor_ijk(mat,5,6,1)*minor_ijk(mat,5,3,2)*minor_ijk(mat,4,6,2)
         -minor_ijk(mat,3,4,2)* minor_ijk(mat,5,6,2)*minor_ijk(mat,5,3,1)*minor_ijk(mat,4,6,1)];
end function;

all_gamma := function(mat)
 /* mat must be 3x6-matrix. */
 res1 := [&*gamma_3 (mat, 1, i,j,k,l,m) : i,j,k,l,m in [2..6] |
            (#{i,j,k,l,m} eq 5) and (j lt k) and (l lt m)];
 res2 := [&*gamma_2(mat,1,i,j,k,l,m) : i,j,k,l,m in [2..6] |
            (#{i,j,k,l,m} eq 5) and (i lt j) and (k lt l) and (l lt m)];
 return res1 cat res2;
end function;

/* End of function definitions. Code to be immediately executed starts here. */


/* Construct linear and cubic relations of the gamma variety by interpolation.
Start with a list of blow-up points. 300 items are hopefully sufficient. */
bu_ll := [[[1,0,0],[0,1,0],[0,0,1],[1,1,1],
           [Random([-100..100]),Random([-100..100]),Random([-100..100])],
           [Random([-100..100]),Random([-100..100]),Random([-100..100])]] : k in [1..300]];

/* Points in P^39 to start the interpolation. */
printf"Evaluating minors\n";
gamma_ll := [all_gamma(Matrix(a)) : a in bu_ll];

/* The linear subspace. */
printf"Reconstructing relations\n";
bm := BasisMatrix(RowSpace(Matrix(RationalField(),gamma_ll)));
assert NumberOfRows(bm) eq 10;
r10 := PolynomialRing(RationalField(),10);
se := ElementToSequence(Vector([r10.i : i in [1..10]]) * ChangeRing(bm,r10));
/* se consists of linear forms defining a morphism from P^9 into the P^39 with the 40 gammas as coordinates.
bm is the corresponding 10x40-matrix. Both will be needed below. */

/* Computing the pre-image in P^9 of the cubic relations. */
ul := [ElementToSequence(Solution(bm, Vector(RationalField(),a))) : a in gamma_ll];
mon3 := Monomials(&+[r10.i : i in [1..10]]^3);
m3_mat := [[Evaluate(f,p) : f in mon3] : p in ul];
ker := Kernel( Transpose(Matrix(m3_mat)));
rel3_co := BasisMatrix(ker);
/* List of cubic relations */
rel3 := [&+[mon3[i] * rel3_co[k,i] : i in [1..#mon3]] : k in [1..NumberOfRows(rel3_co)]];
assert #rel3 eq 30;


/* Constructing generators of W(E_6) by interpolation.
   Applying a list of generators to the point bo. */

printf"Constructing representation of W(E_6)\n";
/* Points in general position. Fixed for the calculation to follow. */
bo :=  [[1,0,0],[0,1,0],[0,0,1],[1,1,1],[13,19,-31],[5,81,-31]];
/* Applying a Cremona transformation */
bo2 := [bo[i] : i in [1..4]] cat [[1 / a : a in bo[i]] : i in [5,6]];
/* Applying two generators of S_6. */
bo3 := [bo[i] : i in [2,3,4,5,6,1]];
bo4 := [bo[i] : i in [2,1,3,4,5,6]];

/* The four corresponding gamma vectors */
g_ll_gen := [all_gamma(Matrix(a)) : a in [bo, bo2, bo3, bo4]];
/* Normalization */
mul_l := [ LCM([Denominator(b) : b in a]) : a in g_ll_gen];
g_ll_gen := [[IntegerRing()!(mul_l[i] * a) : a in g_ll_gen[i]] : i in [1..#g_ll_gen]];
ggt_l := [GCD(a) : a in g_ll_gen];
g_ll_gen := [[a div ggt_l[i] : a in g_ll_gen[i]] : i in [1..#g_ll_gen]];

/* Recovering the signed permutation on the 40 gamma-coordinates. */
mat_l := [];
for i := 2 to #g_ll_gen do
 mat := ZeroMatrix(IntegerRing(),40,40);
 for j := 1 to 40 do
  k := Index(g_ll_gen[i], g_ll_gen[1][j]);
  if k eq 0 then k := Index(g_ll_gen[i], -g_ll_gen[1][j]); end if;
  mat[k,j] := g_ll_gen[1][j] div g_ll_gen[i][k];
 end for;
 assert IsRegular(mat);
 Append(~mat_l,mat);
 /* mat_l consists of the three matrices for bo->bo2, bo->bo3, and bo->bo4. */
end for;

we6 := sub<GL(40,IntegerRing()) | mat_l>;
assert Order(we6) eq 51840;
/* 40-dimensional linear representation of W(E_6) that acts on the gamma variety. */


/* Construct a cyclic group of order 9. There is only one subgroup of this type. */
cl := SubgroupClasses(we6);
c9 := [a`subgroup : a in cl | IsCyclic(a`subgroup) and Order(a`subgroup) eq 9][1];

/* A number field with Galois group c9. */
printf"Initializing Numberfield\n";
q<t> := PolynomialRing(RationalField());
pol9 := t^9 - 19*t^8 + 152*t^7 - 665*t^6 + 1729*t^5 - 2717*t^4 + 2508*t^3 - 1254*t^2 + 285*t - 19;
nf := NumberField(pol9);
nf_aut, aut_set, aut_action := AutomorphismGroup(nf);

suc, phi := IsIsomorphic(nf_aut, c9);
/* Explicit automorphism from the Galois group to the permutation group. */
assert suc;

gen_gal := GeneratorsSequence(nf_aut);
/* For c9, there is only one generator. */
gen_c9 := [ChangeRing(phi(a),RationalField()) : a in gen_gal];
/* The generator(s) of c9 as a 40x40-matrix. */


/* Constructing the lattice in \\cal{O}^10. */
std_bas := RowSequence(IdentityMatrix(IntegerRing(),10));
/* The standard base of ZZ^10. */

printf"Constructing lattice in \\cal{O}^10\n";
ord := MaximalOrder(nf);

/* The two operations of c9 on \\cal{O}^10, described by (10*#c9)x(10*#c9)-matrices. */
gal_l := []; we6_mat_l := [];
for i := 1 to #gen_gal do
 gal_ord_mat := Matrix([ElementToSequence(ord!(aut_action(gen_gal[i])(ord.k)))
                        : k in [1..Degree(ord)]]);
 vec := Vector([Evaluate(f,std_bas[1]) : f in se]);
 we6_dim_10_mat := [Solution(bm, Vector([Evaluate(f,std_bas[k]) : f in se]) * gen_c9[i])
                    : k in [1..10]];
 Append(~gal_l, DiagonalJoin([gal_ord_mat : i in [1..10]]));
 m1 := Matrix(we6_dim_10_mat);

 Append(~we6_mat_l,
 VerticalJoin([HorizontalJoin([IdentityMatrix(RationalField(),Degree(ord)) *
                       m1[k,j] : j in [1..10]]) : k in [1..10]]) );
end for;

/* The rank-10 ZZ-lattice in \\cal{O}^10. */
z_ur := &meet [Kernel(ChangeRing(we6_mat_l[i] - gal_l[i], IntegerRing())) : i in [1..#gen_gal]];
lat := Lattice(z_ur);
lat := BasisReduction(lat);
z_bm := BasisMatrix(lat); /* ZZ-basis of the lattice, 10 vectors in ZZ^(10*#c9). */

/* ZZ-basis of the lattice, 10 vectors in \\cal{O}^10. */
z_bm_o_10 := [[ord![z_bm[i,j + Degree(ord) * k] : j in [1..Degree(ord)]] : k in [0..9]] : i in [1..10]];


/* LLL-reduction */
/* The Minkowski basis */
mink_bas := [&cat [&cat [[Re(c), Im(c)] : c in Conjugates(a: Precision := 300)] : a in b] : b in z_bm_o_10];
e1,e2,e3 := LLL(Matrix(mink_bas));
/* The LLL-minimal basis for the lattice in \\cal{O}^10. */
z_bm_o_10 := RowSequence( e2 * Matrix(z_bm_o_10));

/* The coordinate transformation on P^9 to LLL basis. */
nf_r10 := PolynomialRing(nf, 10);
para := [&+[nf_r10.i * z_bm_o_10[i][j] : i in [1..10]] : j in [1..10] ];

printf"Transforming cubic relations\n";
/* The naive transformation */
time kub_tran := [Evaluate(f, para) : f in rel3];
/* Decomposing into components defined over Q. */
kub_des := [];
for i := 1 to #kub_tran do
 co_seq := [ElementToSequence(MonomialCoefficient(kub_tran[i],nf_r10!m3 )) : m3 in mon3];
 for j := 1 to Degree(nf) do
  Append(~kub_des, &+[mon3[i] * co_seq[i][j] : i in [1..#mon3]]);
 end for;
end for;

/* Eliminating redundancy in cubic relations. Determining 30 linearly independent relations. */
co_mat := [[MonomialCoefficient(f,m3) : m3 in mon3 ] : f in kub_des];
co_mat_2 := BasisMatrix(RowSpace(Matrix(co_mat)));
kub_des := [&+[co_mat_2[i][j] * mon3[j] : j in [1..#mon3]] : i in [1..30] ];


/* Point search. Very naive O(N^10)-algorithm. Search bound 1. */
printf"Searching for rational points on the twist\n";
sol1 := [];
r4<x1,x2,x3,x4> := PolynomialRing(RationalField(),4);
time
for i1 in [0..1] do
 for i2 in [-1..1] do
  for i3 in [-1..1] do
   for i4 in [-1..1] do
    for i5 in [-1..1] do
     for i6 in [-1..1] do
      spez1 := Evaluate(kub_des[1],[i1,i2,i3,i4,i5,i6,x1,x2,x3,x4]);
      sol1 := sol1 cat [ [i1,i2,i3,i4,i5,i6] cat [j1,j2,j3,j4] : j1,j2,j3,j4 in [-1..1] |
                          Evaluate(spez1,[j1,j2,j3,j4]) eq 0 ];
     end for;
    end for;
   end for;
  end for;
 end for;
end for;

sol2 := [a : a in sol1 | &and [Evaluate(f,a) eq 0 : f in kub_des]];
sol3 := [a  : a in sol2 | #[c : c in a | c eq 0] lt 9]; /* Excluding the trivial points. */

printf"Using solution %o\n",sol3[1];
printf"Mapping back to the original gamma-variety\n";
gz_10 := [Evaluate(f,sol3[1]) : f in para];  /* Back to the original coordinates in P^9. */
gz_40 := [Evaluate(f,gz_10) : f in se];      /* Back to P^39. */


/* Calculating the Clebsch-Salmon-invariants from the values of gamma. */
printf"Computing Clebsch-Salmon-invariants\n";
P_2 := RationalField()!&+[a^2 : a in gz_40];
P_4 := RationalField()!&+[a^4 : a in gz_40];
P_6 := RationalField()!&+[a^6 : a in gz_40];
P_8 := RationalField()!&+[a^8 : a in gz_40];
P_10 := RationalField()!&+[a^10 : a in gz_40];

i8 := (-6) * P_2;
i16 := -24 * P_4  + (41/16) * P_2^2;
i24 := 576/13 * P_6 -396/13 * P_4 * P_2 +  29/13 * P_2^3;
i32 := -62208/1171 *P_8 + 54864/1171*P_6 * P_2
          + 203616/1171 * P_4^2 -61287/1171 * P_4 * P_2^2
          + 13393/4684 * P_2^4;

i40 := 41472/155 * P_10 -4605984/36301 *P_8 * P_2
          -106272/403 * P_6 * P_4 +19990440/471913 * P_6 * P_2^2
       + 47719206/471913 * P_4^2 * P_2 -7468023/471913 * P_4 * P_2^3
       + 10108327/18876520 * P_2^5;

printf"Invariants %o found. Computing Equation\n",[i8,i16,i24,i32,i40];

/* Determining an equation. */
gl := inv_2_gl([i8,i16,i24,i32,i40]);
printf"Minimizing and reducing of result\n";
time gl2 := MinimizeReduceCubicSurface(gl);
gl2;
