// copy and paste the first 210 lines as one block into a Magma, then
// play with the worked examples from line 212 ffoll.

function Pts(r,A)
  return [[a*k mod r : a in A] : k in [1..r-1]];
end function;

function JunSuff(r,A)
  Points := Pts(r,A);
  Jun := [P : P in Points | &+P eq r];
  Mid := [P : P in Points | &+P eq 2*r];
  Sums := [[P1[i]+P2[i] : i in [1..#A]] : P1 in Jun, P2 in Jun];
    if &and[P in Sums : P in Mid] then return true;
    else return false, [P : P in Mid | P notin Sums][1];
    end if;
end function;

function Solutions(r);
  TempSolutions := [[0,0,0,0]]; // initiate list with nonsense
  for a1 in [i : i in [0..r div 4]],
    a2 in [i : i in [a1..(r-a1) div 3]],
      a3 in [i : i in [a2..(r-a1-a2) div 2]] do
     a4 := r-a1-a2-a3;
     A := [a1,a2,a3,a4];
      if
    (a4 ne r) and (GreatestCommonDivisor([r] cat A) eq 1) and
    (&and[A ne Sort(B) : B in Pts(r,C), C in TempSolutions])
    and (JunSuff(r,A)) then
      Append(~TempSolutions,A);
     end if;
  end for;
  return Remove(TempSolutions,1); // take out initial nonsense
end function;

/* r := 23; A := [1,2,4,16]; //try it for 1/23(1,2,4,16)
  Jun := [P : P in Pts(r,A) | &+P eq r]; Jun;
  [B : B in Jun | &and [JunSuff(B[i], [-r]
                           cat Remove(B,i)) : i in [1..4]]];
*/

KK := FiniteField(1009);
RR<x,y,z,t> := PolynomialRing(KK,4);
function InvariantMonomials(r,A)
  Inv := [x^r,y^r,z^r,t^r,x*y*z*t];
    for i in [0..r-1], j in [0..r-i-1], k in [0..r-i-j-1] do
      if i*j*k eq 0 then
    Append(~Inv,x^((-A[2]*i-A[3]*j-A[4]*k) mod r)*y^i*z^j*t^k);
      elif (A[2]*i+A[3]*j+A[4]*k) mod r eq 0 then
        Append(~Inv,y^i*z^j*t^k);
      end if;
    end for;
  Exclude(~Inv,1); // Omit 1 before testing divisibility.
S := [[RR!1]]; // S is a list of lists Si of mons of deg i.
  for i in [1..r] do
    Si := [m : m in Inv | TotalDegree(m) eq i];
    Append(~S,Si); // Omit those divisible by mons of deg i.
    Inv := [m : m in Inv |
             not &or[IsDivisibleBy(m,n) : n in Si]];
  end for;
    if #Inv ne 0 then
        error("Invariant monomial has degree > r");
    end if;
  return S;
end function;

function EigenSp(r,A)
  Inv := InvariantMonomials(r,A); 
  Eig := [&cat[S : S in Inv]]; // Collate the list of lists.
  Rbar := quo<RR | Exclude(Eig[1],1)>;
    for s in [1..r-1] do
  temp := [x^((s-A[2]*i-A[3]*j-A[4]*k) mod r)*y^i*z^j*t^k :
       k in [0..r-i-j-1], j in [0..r-i-1], i in [0..r-1] |
            i*j*k eq 0] cat [y^i*z^j*t^k :
       k in [0..r-i-j-1], j in [0..r-i-1], i in [0..r-1] |
       (i*j*k ne 0) and ((A[2]*i+A[3]*j+A[4]*k) mod r eq s)];
  temp := [m : m in temp | Rbar!m ne 0];
  Append(~Eig,temp);
  end for;
return Rotate(Eig,-1); // Put invariants in final rth place.
end function;

function NonZeroMonomials(Rbar,S)
 return [m : m in S | Rbar!m ne 0];
end function;

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 constructions
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.
  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;
  // I,M; // Uncomment this line for a diagnostic display.
  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;
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;

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;

 // given matrix N, ask for index of positive entry in a row
 // with other entries negative, or vice versa. Return 0 if none.
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;

// Hirzebruch-Jung continued fraction expansion. e.g. 20/3 = [7,3] = 7 - 1/3
function HJ(x);
CF := [];
while x ne Ceiling(x) do
Append(~CF,Ceiling(x));
x := 1/(Ceiling(x)-x);
end while;
Append(~CF,Ceiling(x));
return CF;
end function;


//
// typical way of using
//

Sol20:= [A : A in Solutions(20) | A[1] eq 1];
ASets20 := [ASets(20,A) : A in Sol20];
for i in [1..#Sol20] do A := Sol20[i]; Eig := EigenSp(20,A); i, #ASets20[i], &and[IsNonSingularPiece(20,A,ASets20[i][j],Eig) : j in [1..#ASets20[i]]]; end for;
i:=8; r:=20; A:= Sol20[i]; Eig:= EigenSp(r,A); AS := ASets20[i];
for j in [1..5] do IsNonSingularPiece(r,A,AS[j],Eig); end for;
for j in [15,17] do IsNonSingularPiece(r,A,AS[j],Eig); end for;

Sol30:= [A : A in Solutions(30) | A[1] eq 1];
ASets30 := [ASets(30,A) : A in Sol30];
i:=27; A := Sol30[i]; Eig := EigenSp(30,A);
for S in ASets30[i] do IsNonSingularPiece(30,A,S,Eig); end for;
// 27th case 1/30(1,5,8,16), 56 affine pieces, only No.11 and No.21 singular

r:=15; Solr:= [A : A in Solutions(r) | A[1] eq 1]; ASetsr := [ASets(r,A) : A in Solr];
for i in [1..#Solr] do A := Solr[i]; Eig := EigenSp(r,A); i, #ASetsr[i], &and[IsNonSingularPiece(r,A,ASetsr[i][j],Eig) : j in [1..#ASetsr[i]]]; end for;

First interesting case is r = 12, with 7 cases, of which 2 have typical mild discrepancy.
First singular case is 1/15(1,3,5,6), 28 affines
r = 20: 3 singular case: 1/20(1,2,5,12) with 42 affines. 1/20(1,3,4,12) with 33. 1/20(1,4,5,10) with 30
r = 24: 23 case, of which 15 are nonsingular and 8 singular

Do r = 24, Sol24[9] = 1/24(1,2,5,16). The 4 singular affine pieces are given by
> [j : j in [1..#ASetsr[i]] | not IsNonSingularPiece(r,A,ASetsr[i][j],Eig)];   
[ 19, 27, 28, 42 ]

Case 19:
> M0;
[ 0 -1  2  1] m*c
[ 0  2 -4  1] m*b
[ 2 -1  0  0] m*l
[ 1  2 -1  0] m*b*l
[ 0  2  4  0] m*b*c*l
[ 1  1  1  1] m^2*b*c*l    disc m = 0  (presumably (m = d = 0) + (m = l = 0))
[-1 -1 -1  2] d t^2 = d*x*y*z
[ 1 -1 -3  1] m x*t=m*y*z^3;
[-1  3 -1  0] b y^3=b*x*z
[-1  0  5  0] c z^5 = c*x
[ 1  0  3 -1] l x*z^3 = l*t
> KernelMatrix(M0);
[ 1  0  0  0  0  0  0 -1  0 -1  0]  z^2*t = m*c*y because y*z^3
[ 0  1  0  0  0  0  0 -1 -1  0  0]  y^2*t = m*b*z^4 because y*z^4
[ 0  0  1  0  0  0  0 -1  0  0 -1]  x^2 = m*l*y   because y*z^3
[ 0  0  0  1  0  0  0 -1 -1  0 -1]  x*y^2 = m*b*l*z because z*t
[ 0  0  0  0  1  0  0 -1 -1 -1 -1]  y^2*z^4 = m*b*c*l because t
[ 0  0  0  0  0  1  0 -2 -1 -1 -1]  x*y*z*t = m^2*b*c*l because y
[ 0  0  0  0  0  0  1 -1  0 -1  1]  d*l = m*c because x*y*z
m*c*x*y*z = m*y*z^6 = x*z^3*t = l*t^2 = l*d*x*y*z

Case 27:
> M0;
[ 0  5 -2  0]
[ 1 -1 -3  1]
[ 1  2 -1  0]
[ 2 -1  0  0]
[ 0 -1  2  1]
[ 0  2  4  0]
[ 1  1  1  1]
[ 0 -4  0  2] d  t^2 = d*y^4
[ 0  3  2 -1] mu  y^3*z^2 = mu*t
[ 0  2 -4  1] la    y^2*t = la*z^4
[ 1 -3  1  0] l  x*z = l*y^3
[-1  0  5  0] c  z^5 = c*x
> KernelMatrix(M0);
[ 1  0  0  0  0  0  0  0 -1 -1  0  0]  y^5 = mu*la*z^2 because z^4
[ 0  1  0  0  0  0  0  0  0 -1 -1  0]  x*t = la*l*y*z^3 because y*z^4
[ 0  0  1  0  0  0  0  0 -1 -1 -1  0]  x*y^2 = mu*la*l*z because z*t
[ 0  0  0  1  0  0  0  0 -1 -1 -2  0]  x^2 =  mu*la*l^2*y because y^3
[ 0  0  0  0  1  0  0  0  0 -1 -1 -1]  z^2*t = la*l*c*y because x*y
[ 0  0  0  0  0  1  0  0 -1 -1 -1 -1]  y^2*z^4 = mu*la*l*c*y because y
[ 0  0  0  0  0  0  1  0 -1 -2 -2 -1]  x*y*z*t = mu*la^2*l^2*c*y not basic
[ 0  0  0  0  0  0  0  1  1 -1 -1 -1]  d*mu = la*l*c because x*y
z^2*t = d*mu*y because y^4: d*mu*y^4 = y^3*z^2*t

Case 28
> M0;
[ 0  3  2 -1]
[ 1 -3  1  0]
[ 1  2 -1  0]
[ 2 -1  0  0]
[ 0 -1  2  1]
[ 0  4  0  1]
[ 1  1  1  1]
[ 0 -4  0  2]
[ 0 -2  4 -1]
[ 0  5 -2  0]
[ 1 -1 -3  1]
[-1  2  1  1]
> KernelMatrix(M0);
[ 1  0  0  0  0  0  0  0 -1 -1  0  0]
[ 0  1  0  0  0  0  0  0 -1  0 -1  0]
[ 0  0  1  0  0  0  0  0 -1 -1 -1  0]
[ 0  0  0  1  0  0  0  0 -2 -1 -2  0]
[ 0  0  0  0  1  0  0  0 -1  0 -1 -1]
[ 0  0  0  0  0  1  0  0 -1 -1 -1 -1]
[ 0  0  0  0  0  0  1  0 -2 -1 -2 -1]
[ 0  0  0  0  0  0  0  1  0  1 -1 -1]

Case 42;
> M0;
[ 0 -1  2  1]
[ 0  5 -2  0]
[ 0  2  4  0]
[ 1  1  1  1]
[ 0 -4  0  2]  t^2 = d*y^4
[ 0  3  2 -1]  y^3*z^2 = mu*t
[ 1  0 -5  0]  x = a*z^5
[ 0  2 -4  1]  y^2*t = la*z^4
[ 0 -3  6  0]  z^6 = c*y^3
> KernelMatrix(M0);
[ 1  0  0  0  0  0  0 -1 -1]  z^2*t = la*c*y  because y^3
[ 0  1  0  0  0 -1  0 -1  0]  y^5 = mu*la*z^2 because z^4
[ 0  0  1  0  0 -1  0 -1 -1]  y^2*z^4 = mu*la*c because z^2
[ 0  0  0  1  0 -1 -1 -2 -2]  x*y*z*t = a*mu*la^2*c^2 because 1
[ 0  0  0  0  1  1  0 -1 -1]  d*mu = la*c because y^4

y^3*z^2*t = mu*t^2 = d*mu*y^4

Sol30:= [A : A in Solutions(30) | A[1] eq 1];
ASets30 := [ASets(30,A) : A in Sol30];
for i in [26..30] do A := Sol30[i]; Eig := EigenSp(30,A); i, #ASets30[i], &and[IsNonSingularPiece(30,A,ASets30[i][j],Eig) : j in [1..#ASets30[i]]]; end for;

i:=27; A := Sol30[i]; Eig := EigenSp(30,A);
for S in ASets30[i] do IsNonSingularPiece(30,A,S,Eig); end for;
27th case 1/30(1,5,8,16), 56 affine pieces of which No.11 and No.21 are singular

//

// Application 1: AHilb(1/23(1,2,4,16) is nonsingular with
// a single exceptional divisor of discrepancy 1
// corresponding to P9 = 1/23(9,18,13,6)
// appearing on affine pieces numbers 4,5,9,10,11,12

r := 23; A := [1,2,4,16];
Eig := EigenSp(r,A);
AH23 := ASets(r,A);
[#S : S in AH23];
  // 26
[IsNonSingularPiece(r,A,Si,Eig) : Si in AH23];
  // all true
for i in [1..5] do IsNonSingularPiece(r,A,AH23[i],Eig); end for;

// Application 2:

r := 20;
Sol20:= [A : A in Solutions(20) | A[1] eq 1];
ASets20 := [ASets(r,A) : A in Sol20];
[#S : S in ASets20];
  // [ 20, 24, 20, 20, 20, 20, 20, 42, 27, 33, 20, 30 ]

for A in Sol20 do
   Eig:=EigenSp(20,A);
   &and[IsNonSingularPiece(20,A,S,Eig) : S in ASets(r,A)];
end for;

The three singular guys are:
> Sol20[8];
[ 1, 2, 5, 12 ]
> Sol20[10];
[ 1, 3, 4, 12 ]
> Sol20[12];
[ 1, 4, 5, 10 ]

//////////////////
A := Sol20[12];
AH12 := ASets20[12];
Eig := EigenSp(20,A);
[i : i in [1..#AH12] | not IsNonSingularPiece(r,A,AH12[i],Eig)];
 // [18,21] // i.e. only two affine pieces out of 30 are sing

> IsNonSingularPiece(r,A,AH12[18],Eig);
false
[ 3 -2 -3  0]
[ 2 -3  2  0]
[-1  4 -3  0]
[-2  3  2  0]
[ 0  0 -2  1]

[ 1 -1 -1  1  0]

>   IsNonSingularPiece(r,A,AH12[21],Eig);
false
[ 3 -2 -1 -1]
[ 2 -3  0  1]
[-1  4 -1 -1]
[-2  3  0  1]
[ 0  0  2 -1]

[ 1 -1 -1  1  0]

//////////////////
A := Sol20[10];
AH10 := ASets20[10];
Eig := EigenSp(20,A);
[i : i in [1..#AH10] | not IsNonSingularPiece(r,A,AH10[i],Eig)];
 // [24,28] // i.e. only two affine pieces out of 33 are sing

> IsNonSingularPiece(r,A,AH10[24],Eig);
false
[ 2 -2 -4  0]
[ 1 -3  2  0]
[ 0  4 -3  0]
[-1  3  3  0]
[ 0  0 -3  1]

[ 1 -1 -1  1  0]


//////////////////
AH8 := ASets20[8];
[#S : S in AH8];
  /* [ 4, 5, 7, 7, 4, 4, 9, 11, 9, 10, 8, 7, 5, 4, 12,
  10, 11, 9, 10, 11, 10, 7, 9, 10, 9, 5, 12, 11, 11, 7,
  9, 9, 7, 5, 4, 4, 7, 4, 4, 4, 7, 5 ] */
> [IsNonSingularPiece(r,A,S,Eig) : S in AH8];
  // [ true, true, .. true, false, true, false, true, .. true ]
> [i : i in [1..#AH8] | not IsNonSingularPiece(r,A,AH8[i],Eig)]; 
  // [ 15, 17 ]

//////////////////
// All the calc for r = 20
r := 20;
Sol20:= [A : A in Solutions(20) | A[1] eq 1];
ASets20 := [ASets(r,A) : A in Sol20];
[#S : S in ASets20];
for A in Sol20 do
   Eig:=EigenSp(20,A);
   &and[IsNonSingularPiece(20,A,S,Eig) : S in ASets(r,A)];
end for;
A := Sol20[12];
AH12 := ASets20[12];
Eig := EigenSp(20,A);
[i : i in [1..#AH12] | not IsNonSingularPiece(r,A,AH12[i],Eig)];
  IsNonSingularPiece(r,A,AH12[18],Eig);
  IsNonSingularPiece(r,A,AH12[21],Eig);
A := Sol20[10];
AH10 := ASets20[10];
Eig := EigenSp(20,A);
[i : i in [1..#AH10] | not IsNonSingularPiece(r,A,AH10[i],Eig)];
  IsNonSingularPiece(r,A,AH10[24],Eig);
  IsNonSingularPiece(r,A,AH10[28],Eig);
A := Sol20[8];
AH8 := ASets20[8];
Eig := EigenSp(20,A);
[i : i in [1..#AH10] | not IsNonSingularPiece(r,A,AH8[i],Eig)];
  IsNonSingularPiece(r,A,AH8[15],Eig);
  IsNonSingularPiece(r,A,AH8[17],Eig);

Case 1/20(1,2,5,12). Has 42 affines, of which 2 are singular,
making PP^1 x (3-fold node).
for i in [1..5] do IsNonSingularPiece(r,A,AH8[i],Eig); end for;

discrepant divisors correspond to the vectors:
  P6 = 1/20(6,12,10,12)   disc = 1  on 3,4,9,10,11,18,19,23,24,25
  P9 = 1/20(9,18,5,8)   disc = 1  on 8,10,11,12,13,18
  P14 = 1/20(14,8,10,8)   disc = 1  on 16,19,20,21,24,25,27,28
  P1 + e1 = 1/20(21,2,5,12) disc = 1 on 29,30,32,33,39,40,41,42
  P2 + e1 = 1/20(22,4,10,4) disc = 1 on 20,21,22,27,28,29,31,32,37,38,41,42
  P11 + e2 = 1/20(11,22,15,12)   disc = 2 on 7,8,9,10,16,19


Disc := [];
for i in [1..#AH8] do 
  nonsing, M, De := IsNonSingularPiece(r,A,AH8[i],Eig);
  if nonsing then
    Disc := Disc cat [M[i] : i in [1..4] | De[i] ne 0];
  end if;
end for;

> DiscDiv := SetToSequence(SequenceToSet(Disc));
> [&+[A[i] : i in [1..4]] div 20 : A in DiscDiv];
  // [ 2, 2, 3, 2, 2, 2 ]

> DiscDiv;
[
    ( 6 12 10 12),
    ( 9 18  5  8),
    (11 22 15 12),
    (14  8 10  8),
    (22  4 10  4),
    (21  2  5 12)
]

i.e. The discrepant divisors on the nonsing pieces are P6, P9,
P14, P2+e1, P1 + e1 and  P11+e2 (disc 2).

Disc := [];
for i in [1..#AH8] do 
  nonsing, M, De := IsNonSingularPiece(r,A,AH8[i],Eig);
  if nonsing then
    Disc := Disc cat [M[i] : i in [1..4] | De[i] ne 0];
  end if;
end for;

P6 := DiscDiv[1];
P9 := DiscDiv[2];
Q11 := DiscDiv[3];
P14 := DiscDiv[4];
Q2 := DiscDiv[5];
Q1 := DiscDiv[6];
IB1302525034

for i in [17..20] do IsNonSingularPiece(r,A,AH8[i],Eig); end for;


The affine pieces involving P6 are
3,4,5,10,11,18,19,23,24,25 and 17 (singular).

The star of P6 in the fan of A-Hilb is
e2 over the P6-centred pentagon P4 P1 P2 Q11 P9
i.e. the 5 simplexes
4. P6 e2 P4 P1
3. P6 e2 P1 P2
9. P6 e2 P2 Q11
10. P6 e2 Q11 P9
11. P6 e2 P9 P4

and

23. P6 P1 P4 P12
25. P6 P14 P1 P2
19. P6 P14 P2 Q11
24. P6 P14 P12 P1
17. P6 P14 Q11 P12 P9 (not simplex)
18. P6 P12 P9 P4