// Tue 5th Aug repaired the dependence on a1 = 1. Still to test.
// The ASets routine goes to line 155
// It includes IsNonSingularPiece that works when it works..

RR<x,y,z,t> := PolynomialRing(Rationals(),4);
function InvariantMonomials(r,A)
   assert #A eq 4 and IsDivisibleBy(&+A,r);
   Inv0 := [x^i*y^j*z^k*t^l : i in [0..r div GCD(r,A[1])],
j in [0..r div GCD(r,A[2])], k in [0..r div GCD(r,A[3])],
l in [0..r div GCD(r,A[4])] |
IsDivisibleBy(i*A[1]+j*A[2]+k*A[3]+l*A[4],r)];
   return [RR!1] cat MinimalBasis(I) where I is Ideal(Remove(Inv0,1));
end function;

function EigenSp(r,A)
   cha := func< m | &+[A[j]*Degree(m,RR.j) : j in [1..4]] mod r >;
   Inv := InvariantMonomials(r,A); 
   Eig := [Inv]; // Concatenate the list of lists.
   Rbar := quo< RR | Remove(Inv,1) >;  // don't set 1 = 0.
   MB := [RR!m : m in MonomialBasis(Rbar) ];
// MB;
      for s in [1..r-1] do
         Append(~Eig, [m : m in MB | cha(m) 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;

function convex_row(N)
  n := NumberOfRows(N);
  J := 0;
  for i in [1..n] do
      if J eq 0 then
         pos := [j : j in [1..n+4] | N[i,j] gt 0];
         neg := [j : j in [1..n+4] | N[i,j] lt 0];
           if (#neg eq 1) and (#pos gt 1) then
              found := true; J := neg[1];
           elif (#pos eq 1) and (#neg gt 1) then
              found := true; J := pos[1];
           end if;
      end if;
  end for;
  return J;
end function;

function IsNonSingularPiece(r,A,Si,Eig)
M := RelMatrix(r,A,Si,Eig);
  M0 := M;
  good := true;
  while good and (NumberOfRows(M0) gt 4) do
    N := KernelMatrix(M0);
    j := convex_row(N);
       if j ne 0 then RemoveRow(~M0,j);
       else good := false;
       end if;
  end while;
if good then
  N0 := Transpose(Adjoint(M0)); // dual cone
    if &and[N0[i,j] le 0 : j in [1..4], i in [1..4]]
      then N0 := -1*N0;
    end if;
  Discr := [1/r*&+[N0[i,j] : j in [1..4]]-1 : i in [1..4]];
  // return true, M0, N0, Discr;
  return true, N0, Discr;
else
  return false, M0, KernelMatrix(M0);
end if;
end function;

// for the 1-2 symmetric case
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]] >;

// Example
for r in [25..37] do
A := [1,1,5,r-7];
printf "=== Group 1/%o(%o,%o,%o,%o)===\n", r, A[1], A[2], A[3], A[4]; 
AS, Eig := ASets(r,A);
#AS; [#S : S in AS]; [i : i in [1..#AS] | not IsNonSingularPiece(r,A,AS[i],Eig)];
IsNonSingularPiece(r,A,AS[j],Eig) where j is $1[1];
end for;