/* This applies the Magma functions from my McKay webpage
  www.warwick.ac.uk/~masda/McKay/ASets_v3.txt
in particular ASets(r,A) and Relmatrix(r,A,Si,Eig) to analyse the
A-Hilb of a 12-sym group 1/r(a,a,b,c) with a | r and 2*a+b+c = r,
and GCD([r,a,b,c]) = 1. The outcome is a list of toric cones that
make up the fan of A-Hilb. (This assume a current conjecture, but
it works so far in most cases with r up to about 150. If it
crashes with the error message "assert fails", it may indicate
a counter-example.)

The formatted output lists the cones starting with the link of C,
then the cones of ASets for which x is a basic monomial. The cones
with y a basic monomials repeat the same information. The routine
omits these by default. As an experimental result, the output seems
to follow 2 overall patterns: either AHilb is smooth, and its
exceptional divisors have discrepancy 0 and 1. Or there are
exceptional divisors of large discrepancy on affine pieces that
have a conifold singularity a*d = b*c.

To run the script, copy and past the whole of the file
  www.warwick.ac.uk/~masda/McKay/ASets_v3.txt
into a Magma, then copy the file below with your own choice
of r and A replacing the top line of code. */

// Let Magma choose a (1,2)-symmetric group for you:
b := Random([1..12]); c := Random([5..16]);
a := Random([al : al in Divisors(b+c) | GCD([al,b,c]) eq 1]); // require GCD([r,a,b,c]) eq 1;
A := [a,a,b,c]; r := &+A; assert GCD([r,a,b,c]) eq 1;

// or more usefully input 12-sym group of your choice, for example
// r := 31; a := 7; A := [a,a,1,r-2*a-1];
// r := 51; A := [17,17,12,5];

if IsEven(r) and IsEven(A[3]) and IsEven(A[4]) then
printf "Warning: Group 1/%o(%o,%o,%o,%o). Centroid C is in lattice L, so \
A-Hilb is a crepant resolution, and there are better ways of computing it.\n",
r, A[1], A[2], A[3], A[4];
end if;
AS0, Eig := ASets(r,A);
AS := [Si : Si in AS0 | x in Si and y in Si]
      cat [Si : Si in AS0 | x notin Si and y in Si]
      cat [Si : Si in AS0 | x in Si and y notin Si];
kk := #AS; printf "A-Hilb has %o affine pieces.\n", kk;
if kk eq r then printf("A-Hilb is a crepant resolution.\n"); end if;

L := ToricLattice(4);
  // cross product of two 3-vectors, but first entry doubled
cross := func< A,B | [A[2]*B[3]-A[3]*B[2],
A[2]*B[3]-A[3]*B[2],
A[3]*B[1]-A[1]*B[3],
A[1]*B[2]-A[2]*B[1]] >;

/* Specifically for (1,2)-symmetric case, analyse each
A-set as cone of fan. The observed result, that I hope
to prove, is that every cone is simplicial or cone over
square tile. The line R := Rays(C); assert #R le 5;
verifies this or would causes a crash. */

for i in [1..kk] do
   M0 := Matrix(Rationals(),4, Eltseq(RelMatrix(r,A,AS[i],Eig)));
      if x in AS[i] and y in AS[i] then
      k := NumberOfRows(M0);
      C := Cone([L | M0[i] : i in [1..k]]); // C;
      R := Rays(C); assert #R eq 4;
      M1 := Submatrix(M0, [i : i in [1..NumberOfRows(M0)] |
                       &or [L!M0[i] in Cone(r): r in R]], [1..4]);
      N0 := r*Transpose(M1^-1);
printf "==== Group 1/%o(%o,%o,%o,%o), Affine piece No %o out of %o ====\n",
   r, A[1], A[2], A[3], A[4], i, kk;
printf "In the link of C, Basic Cone with vertices:\n%o,\nDiscrepancies:\n%o.\n\n",
    N0, [1/r*&+[N0[i,j] : j in [1..4]]-1 : i in [1..4]];
    elif x notin AS[i] and y in AS[i] then
    k := NumberOfRows(M0);
    C := Cone([L | M0[i] : i in [1..k]]); // C;
    R := Rays(C); assert #R le 5;
    M1 := Submatrix(M0, [i : i in [1..NumberOfRows(M0)] |
                     &or [L!M0[i] in Cone(r): r in R]], [1..4]);
if #R eq 4 then N0 := r*Transpose(M1^-1);
printf "==== Group 1/%o(%o,%o,%o,%o), Affine piece No %o out of %o ====\n",
   r, A[1], A[2], A[3], A[4], i, kk;
printf "Monomial x is Basic, Basic Cone with vertices:\n%o,\nDiscrepancies:\n%o.\n\n",
   N0, [1/r*&+[N0[i,j] : j in [1..4]]-1 : i in [1..4]];
elif #R eq 5 then N := Submatrix(M1,[1,2,4,5],[1,3,4]);
printf "==== Group 1/%o(%o,%o,%o,%o), Affine piece No %o out of %o ====\n",
   r, A[1], A[2], A[3], A[4], i, kk;
printf "Monomial x is Basic, Cone over Square Tile with vertices\n%o\n%o\n%o\n%o,\nDiscrepancies:\n[%o, %\o, %o, %o].\n\n",
 Matrix(4, Eltseq(cross(N[1],N[3]))),
 Matrix(4, Eltseq(cross(N[3],N[4]))),
 Matrix(4, Eltseq(cross(N[4],N[2]))),
 Matrix(4, Eltseq(cross(N[2],N[1]))),
 (&+cross(N[1],N[3])/r)-1,
 (&+cross(N[3],N[4])/r)-1,
 (&+cross(N[4],N[2])/r)-1,
 (&+cross(N[2],N[1])/r)-1
         where N is Submatrix(M1,[1,2,4,5],[1,3,4]);
      end if;
   end if;
end for;