/* Explicit example for the Brauer-Manin obstruction.
   Surface with Galois-Group of Order 27.
   To compute other examples just change the input data. */

Z_4<x,y,z,w> := PolynomialRing(IntegerRing(),4);
Q_4<k1,k2,k3,k4> := PolynomialRing(RationalField(), 4);
Q_u<T> := PolynomialRing(RationalField());

/* Equation of the surface */
surf := -3*x^3 - 6*x^2*y - 3*x^2*z + 3*x^2*w - 3*x*y^2 + 3*x*y*w + 3*x*z^2 + 6*x*w^2 + 2*y^3
- 4*y^2*z - y^2*w + 10*y*z^2 - 4*y*z*w - 9*y*w^2 + 6*z^3 - 8*z^2*w
- 8*z*w^2 + 4*w^3;

/* Splitting field of the Brauer class, leading to 3 pairwise skew lines */
spp := T^3 + T^2 - 10*T - 8;
nf_sp := NumberField(spp);

/* Explicit description of the lines: */
ideal_lines := ideal<Q_4 | [Evaluate(surf,[1,i,k1 + k2 * i, k3 + k4 * i]) : i in [-2..2]] >;
printf"Searching for lines:\n";
gb_lines := GroebnerBasis(ideal_lines);
p_27 := Evaluate(gb_lines[4],[0,0,0,T]);

/* search for 3 pairwise skew lines to blow them down */
printf"Searching for orbit of 3 pairwise skew line:\n";
fac := Factorization(p_27,nf_sp);
bdf := false;
bdl := false;
i := 1;
while (i le #fac) and (bdf cmpeq false) do
 if Degree(fac[i][1]) eq 3 then
  nf_2 := ext<nf_sp | fac[i][1]>;
  nf_2_r<u> := PolynomialRing(nf_2);
  rt := Roots(fac[i][1],nf_2);
  if #rt eq Degree(fac[i][1]) then /* All lines of the orbit are defined over nf_2 */
   lines_3 := [ [ Roots(Evaluate(gb_lines[l],[u,u,u,rt[j][1]]))[1][1] : l in [1..3]]
          cat [ rt[j][1]]: j in [1..#rt]];
   p_l := [ [[ 1, 0, lines_3[j][1], lines_3[j][3]],
             [ 0, 1, lines_3[j][2], lines_3[j][4]]] : j in [1..#lines_3]];
   if {Rank(Matrix(p_l[j] cat p_l[k])) : j,k in [1..#lines_3] | j lt k} eq {4} then
    bdf := nf_2;
    bdl := lines_3;
   end if;
  end if;
 end if;
 i := i + 1;
end while;

/* Test against unexpected input */
if bdf cmpeq false then
 print"No number field for blow-down found\n";
 assert false; /* Terminate script, unexpected input-data */
end if;

/* Construct triliear equation of surface by using the lines in lines_3 */

/* search for points to interpolate the equation */
printf"Searching for rational points:\n";
time
points := PointSearch(Scheme(ProjectiveSpace(RationalField(),3),surf),50);
points := [ElementToSequence(a) : a in points];

/* Sort rational points by height */
mul_l := [LCM([Denominator(a) : a in pt]) : pt in points];
points := [[mul_l[i] * points[i][j]  : j in [1..4]] : i in [1..#points]];
Sort(~points,func<a,b | Max([AbsoluteValue(c) : c in a]) - Max([AbsoluteValue(c) : c in b])>);

assert #points ge 20; /* To less points for further computations */

/* Construct a trilinear equation in P^1 x P^1 x P^1 */
/* Point in (A^1)^3 on the surface */
p_3A1_l := [];

p_l := [ [[ 1, 0, lines_3[j][1], lines_3[j][3]],
          [ 0, 1, lines_3[j][2], lines_3[j][4]]] : j in [1..#lines_3]];

nf_2_r<u> := PolynomialRing(bdf);
printf"Mapping points to trilinear surface:\n";

for j := 1 to #points do
 Append(~p_3A1_l,[Roots(Determinant(Matrix(nf_2_r,[[1,0,0,u],points[j]] cat p_l[i])))[1][1] : i in [1..#bdl]]);
end for;

/* Linear system for the coefficients of the trilinear equation */
mat := Matrix(bdf,[[a[1]*a[2]*a[3],a[1]*a[2],a[1]*a[3],a[2]*a[3],a[1],a[2],a[3],1] : a in p_3A1_l]);
printf"Computing trilinear form:\n";

ker := BasisMatrix(Kernel(Transpose(mat))); /* Coefficients of the trilinear form. */

r3<u1,u2,u3> := PolynomialRing(bdf, 3); /* Ring for affine computations */
r6<v1,v2,v3,v4,v5,v6> := PolynomialRing(bdf, 6); /* Ring for projective computations  */

/* Affine and projective trilinear equation */
tril_eq_pro := ker[1,1]*v1*v2*v3 + ker[1,2]*v1*v2*v6 + ker[1,3]*v1*v3*v5 + ker[1,4]*v2*v3*v4
             + ker[1,5]*v1*v5*v6 + ker[1,6]*v2*v4*v6 + ker[1,7]*v3*v4*v5 + ker[1,8]*v4*v5*v6;

tril_eq := Evaluate(tril_eq_pro,[u1,u2,u3,1,1,1] );

/* Exceptional curves on the trilinear surface */
printf"Searching for exceptional curves on degree 6 del Pezzo:\n";

ex_c_tril := [ideal<r3 | [ Evaluate(tril_eq,j,i) : i in [-2..2]]> : j in [1..3]];

gb_c_tril := [GroebnerBasis(a) : a in ex_c_tril];

ex_l := [];

for i := 1 to 3 do
 e1,e2,e3 := IsUnivariate(gb_c_tril[i][2]);
 ra := Roots(e2);
 for j := 1 to #ra do
  f1,f2,f3 := IsUnivariate(Evaluate(gb_c_tril[i][1],e3,ra[j][1]));
  rb := Roots(f2);
  if i eq 1 then    akt := [[0, rb[1][1], ra[j][1]], [1, rb[1][1], ra[j][1]] ];   end if;
  if i eq 2 then    akt := [[rb[1][1], 0, ra[j][1]] , [rb[1][1], 1, ra[j][1]]];   end if;
  if i eq 3 then    akt := [[rb[1][1], ra[j][1], 0], [rb[1][1], ra[j][1], 1] ];   end if;
  Append(~ex_l,akt);
 end for;
end for;

/* Compute the intersection matrix of the 6 exceptional curves */
printf"Computing intersection matrix:\n";
sm := ZeroMatrix(IntegerRing(), 6,6);

for i := 1 to 6 do for j := 1 to 6 do
 if i ne j then
  tm := [[ex_l[i][1][k] - ex_l[i][2][k] : k in [1..3]],
         [- ex_l[j][1][k] + ex_l[j][2][k] : k in [1..3]]  ];
  mat1 := Matrix(tm);
  mat2 := Matrix(tm cat [[ex_l[j][2][k] - ex_l[i][2][k] : k in [1..3]] ]);
  if Rank(mat1) eq Rank(mat2)  then sm[i,j] := 1; end if;
 end if;
end for; end for;

printf"Computing image of line in P^2:\n";
/* indices of triples of pairwise skew lines */
ind := [[i,j,k] : i in [1,2], j in [3,4], k in [5,6] | sm[i,j] + sm[i,k] + sm[j,k] eq 0][1];

/* Projective linear transformation to normalize the trilinear form */
m1 := Matrix(bdf,[[ ex_l[ind[2]][1][1], 0,0,ex_l[ind[3]][1][1],0,0],
               [ 0,ex_l[ind[3]][1][2],0,0, ex_l[ind[1]][1][2],0],
               [ 0,0,ex_l[ind[1]][1][3], 0,0,ex_l[ind[2]][1][3] ],
               [1,0,0,1,0,0],[0,1,0, 0,1,0],[0,0,1,0,0, 1]]);

/* normalizes trilinear equation */
tril_eq_pro_nf := tril_eq_pro^m1;

/* Matrix to normalize the points on the trilinear surface */
m1ti := Transpose(m1^(-1));

/* Points on the normalized trilinear surface */
pts_tri_pro_nf := [ElementToSequence(Vector(pt cat [1,1,1]) * m1ti) : pt in p_3A1_l];

c1 := MonomialCoefficient(tril_eq_pro_nf,v1*v2*v3);
c2 := MonomialCoefficient(tril_eq_pro_nf,v4*v5*v6);

m2 := DiagonalMatrix([1/c1,1,1,-1/c2,1,1]);
m2ti := Transpose(m2^(-1));

/* Normalize leading coefficients -- trilinear equation gets form v1*v2*v3 - v4*v5*v6, this surface is birational to P^2 by
   [ x  :  y :         z     ] --> ([  y :  z ], [  z :  x ], [  x :  y ])

   [ v3 : v6 : (v4 / v1) *v6 ] <-- ([ v1 : v4 ], [ v2 : v5 ], [ v3 : v6 ])  */

tril_eq_pro_nf2 := tril_eq_pro_nf^m2;
pts_tri_pro_nf2 := [ElementToSequence(Vector(pt) * m2ti) : pt in pts_tri_pro_nf];

/* Image of two rational points on the trilinear surface: */
pp1 := pts_tri_pro_nf2[1];
pp2 := pts_tri_pro_nf2[2];

/* Image of the points in P^2: */
pts_p2_start  := [[pp1[3],pp1[6],  pp1[4] /pp1[1] * pp1[6]],
                  [pp2[3],pp2[6],  pp2[4] /pp2[1] * pp2[6]] ];

/* Points in P^2 on the line joining pts_p2_start: */
l_pts_p2 := [ [pts_p2_start[1][j] + i * pts_p2_start[2][j] : j in [1..3]]   : i in [-5..5]];

/* Image of the points on the trilinear surface -- normalized, initial and inital affine */
l_pts_tril_nf := [[a[2], a[3], a[1], a[3], a[1], a[2]] : a in l_pts_p2 ];
l_pts_tril := [ Vector(a) * Transpose(m2) * Transpose(m1) : a in l_pts_tril_nf  ];
l_pts_tril_aff := [[a[1] / a[4] , a[2]/a[5], a[3] / a[6]] : a in l_pts_tril | 0 ne a[4]*a[5]*a[6]];

/* Image of the points on the line in P^3 */
l_pts_p3 := [ElementToSequence(BasisMatrix(
               &meet [Rowspace(Matrix(nf_2_r,[[1,0,0,a[i]]] cat p_l[i])) : i in [1..3]]
              )) : a in l_pts_tril_aff];


printf"Interpolating cubic form\n";

/* Compute the space of cubic forms (with rational coefficients) vanishing on the image of the line */
mon_3 := Monomials((x+y+z+w)^3); /* Monomials als a basis of the space of cubic forms. */

/* Expressing the trace correctly requires explicit declaration of all fields... */
tower_up := RelativeField(nf_sp, bdf);
tower_lo := RelativeField(RationalField(), nf_sp);

/* Each point on the line gives us a linear equation for the coefficients of the cubic forms */
mat_cubic := [];
printf"Building linear system for cubic form:\n";
time
for i := 1 to #l_pts_p3 do
 akt := [Evaluate(f, l_pts_p3[i]) : f in mon_3];
 for j := 0 to 17 do
  zq := [Trace(tower_lo!Trace(tower_up!(a * (bdf.1^j)))) : a in akt];
  mul := LCM([Denominator(a) : a in zq]);
  Append(~mat_cubic,[IntegerRing()!(a * mul) : a in zq]);
 end for;
end for;

/* Two-dimensional coefficient lattice */
cf_sp :=  BasisMatrix(Lattice(Kernel(Transpose(Matrix(IntegerRing(),mat_cubic)))));
printf"Linear system solved\n";
if NumberOfRows(cf_sp) ne 2 then
 printf"Space of cubic forms must be two-dimensional.\n";
 assert false;
end if;
cf1 := &+[cf_sp[1,i] * mon_3[i] : i in [1..20]];
cf2 := &+[cf_sp[2,i] * mon_3[i] : i in [1..20]];

/* The Brauer-Manin Form */
if {Evaluate(cf1,a) : a in points} eq {0} then
 bm_form := cf2;
else
 bm_form := cf1;
end if;
printf"Computation of Brauer-Manin-Form finished.\n";

/* Now we check by computing the local evaluation map on rational points */
printf"Checking result:\n";
p1 := 3;   /* For this example only 3 and 7 are BM-primes. Both are inert. */
p2 := 7;


bm_form_values := [ Evaluate(bm_form,a) : a  in points | 0 ne Evaluate(bm_form,a) ];
{[Valuation(a,p1) mod 3, Valuation(a,p2) mod 3] : a in bm_form_values};

{&+ [Valuation(a,p1), Valuation(a,p2)] mod 3 : a in bm_form_values};