/* A Magma routine calculates A-Hilb AA^4 for a finite diagonal
subgroup A in SL(4, CC).  The function ASets has INPUT an expression
1/r(a1,a2,a3,a4) with r a positive integer, and a1,..a4 integers mod r
representing characters of the group mu_r of rth roots of unity. (The
SL(4) condition is that sum ai == 0 mod r. Not all ai should have a
common divisors with r.) The output consists of a list of sets of
monomials that generate an ideal IZ of CC[x,y,z,t] with quotient ring
A_Z = CC[x,y,z,t]/IZ the regular representation of mu_r. This means
that A_Z has a monomial basis MB with exactly one monomial m in each
eigenspace (character space) a mod r.

The output AS = ASets(r,A) consists of a list of lists of monomials
that generate the ideal IZ of an A-cluster. Further routines are
required to translate this to the toric fan of A-Hilb. The first
step is the routine RelMatrix(r,A,Si,Eig) that takes input a
particular A-Set Si, together with the specification of the current
group 1/r(A) and the data Eig specifying the generating monomials
for each character space as a module over the invariants (this data
is computationally expensive). The output RelMatrix has rows the
relations that express the monomials n not in MB as a multiple of m
in MB. (For example, [a,0,-b,-c] says that x^a must be a scalar
multiple of the basis monomial z^b*t^c -- the scalar coefficient is
a function on the affine piece of AHilb.) These rows are points of
the dual lattice L generating the dual cone that make up the fan of
AHilb. */

RR<x,y,z,t> := PolynomialRing(Rationals(),4);
function InvariantMonomials(r,A)
  assert IsDivisibleBy(&+A,r);
  S := [ [RR!1] ];
/* == S is a list of lists of A-invariant polynomials of
  degree 0 .. r, not products of smaller degree invariants == */
 for d in [1..r] do
   Invd := [m : m in MonomialsOfDegree(RR,d) |
   IsDivisibleBy(&+[A[j]*Degree(m,RR.j) : j in [1..4]],r) ];
   Sd := [ m : m in Invd | &and[ not IsDivisibleBy(m, n)
      : n in &cat[S[e] : e in [2..d]] ]];
   Append(~S,Sd);
end for;
return S;
end function;

/* Eig[s] is the basis of the sth eigenspace (character space)
of the A-action. This calculation is expensive. */
function EigenSp(r,A)
Inv := InvariantMonomials(r,A); 
Eig := [&cat[S : S in Inv]]; // Concatenate the list of lists.
Rbar := quo< RR | Exclude(Eig[1],1) >;
MB := [RR!m : m in MonomialBasis(Rbar) ];
   for s in [1..r-1] do
     Append(~Eig, [m : m in MB |
     &+[A[j]*Degree(m,RR.j) : j in [1..4]] mod r eq s] );
   end for;
return Rotate(Eig,-1); // Put invariants in final rth place.
end function;
 
/* The main function is the tree search for A-regular lists
that (1) are monomial basis for a quotient ring R/I by a monomial
ideal, and (2) have exactly one monomial in each eigenspace. E is
a list under construction of the filled eigenspaces so far. */
function NonZeroMonomials(Rbar,S)
  return [m : m in S | Rbar!m ne 0];
end function;
/* Determines the branch points in the tree with a choice to make. */
function NotOne(R,E)
  return [[i] cat L : i in [1..#E] | #L ne 1 where
     L is NonZeroMonomials(R,E[i])];
end function;

function ASets(r,A)
Eig := EigenSp(r,A);
 // initialise dynamic construction
Bag := []; // To store the fruit.
finished := false;
I := []; // List of indexes i1, i2, ... where I make choices.
C := []; // List of lists of candidates mi in Eig[i], i in I.
       // Cycle through monomials in last list then delete it.
M := []; // List of number of current choices from C.
   // In other words, current choices are C[i,M[i]] for i in I.
 
   while not finished do
   // Make the current quotient ring by killing all monomials
   // in all Eig[i] except the current choices mi for i in I.
   Rbar := quo< RR | Exclude(Eig[r],1) cat
         &cat[Exclude(Eig[I[i]],C[i,M[i]]) : i in [1..#I]] >;
   // Analyse whether it is an A-cluster (A-regular list).
   tlist := NotOne(Rbar,Eig);
     if tlist eq [] then
     Include(~Bag,Basis(DivisorIdeal(Rbar))); end if;
     // If so, bag the fruit, intending to backtrack.
     // If each tlist[i] is [i, plus >= 1 monomial], so no
     // contradiction or result, go forward: add new index to I, 
     // new list of candidates to C, and first choice to M.
        if (tlist ne []) and (&and[#l ne 1 : l in tlist]) then 
        Append(~I,Integers()!tlist[1,1]);
        Append(~C,Remove(tlist[1],1)); Append(~M,1);
        end if;
/* Uncomment followimg line for a diagnostic display of where
we are in the tree -- the branches and choices made */
  // I,M;
   if (tlist eq [])
       or ((tlist ne []) and (&or[#l eq 1 : l in tlist])) then
  // Backtrack, include setting finished to true if can't.
      ss := #I;
        if &and[M[i] eq #(C[i]) : i in [1..ss]] then
        finished := true;
        else tt := Max([i : i in [1..ss] | M[i] ne #(C[i])]);
          I:=I[1..tt]; M:=M[1..tt]; M[tt]:=M[tt]+1; C:=C[1..tt];
        end if;
   end if;
   end while;
return Bag, Eig;
 // Include the expensive Eig as second output of the function.
end function;

function EigenBasis(S,E)
  Rbar := quo< RR | S >;
  MonB := MonomialBasis(Rbar);
  return [[n : n in MonB | n in E[i]][1] : i in [1..#E]];
 end function;

/* A basis of ideal IZ: for monomial n not in MonB, find the
unique m in MonB in the same eigenspace. The rows of RelMatrix
have positive exponents of n and negative of m */

function RelMatrix(rr,AA,SS,Eig)
  Rbar := quo<RR|SS>;
  EB := EigenBasis(SS,Eig);
  M := Matrix(4,[Integers()|1,1,1,1]); // add invariant x*y*z*t
 for m in SS do
   if m eq x*y*z*t then; // don't add invariant x*y*z*t twice
   else
   B := [Degree(m,RR.i) : i in [1..4]];
     c := &+[AA[i]*B[i] : i in [1..4]] mod rr;
        if c eq 0 then c := rr; end if;
        n := EB[c];
        C := [Degree(n,RR.i) : i in [1..4]];
        M := VerticalJoin(M,Matrix(4,[B[i]-C[i] : i in [1..4]])); 
        end if;
    end for;
return M;
end function;