(*This Mathematica notebook forms part of the proof of Proposition 4.6 in the paper "On the arithmetic of simple singularities of type E" by Beth Romano and Jack A. Thorne*)
(*The following program begins with the root datum of the group E8. It implements the procedure outlined in Proposition 4.6 to generate a list of cusp data, and then verifies for each cusp datum (M_0, M_1) the existence of a function f satisfying the given conditions*)

(*Positive roots of E8 in a given basis; the first 8 roots are the simple roots in the Bourbaki numbering.*)
PRH={{1,0,0,0,0,0,0,0},{0,1,0,0,0,0,0,0},{0,0,1,0,0,0,0,0},{0,0,0,1,0,0,0,0},{0,0,0,0,1,0,0,0},{0,0,0,0,0,1,0,0},{0,0,0,0,0,0,1,0},{0,0,0,0,0,0,0,1},{1,0,1,0,0,0,0,0},{0,1,0,1,0,0,0,0},{0,0,1,1,0,0,0,0},{0,0,0,1,1,0,0,0},{0,0,0,0,1,1,0,0},{0,0,0,0,0,1,1,0},{0,0,0,0,0,0,1,1},{1,0,1,1,0,0,0,0},{0,1,1,1,0,0,0,0},{0,1,0,1,1,0,0,0},{0,0,1,1,1,0,0,0},{0,0,0,1,1,1,0,0},{0,0,0,0,1,1,1,0},{0,0,0,0,0,1,1,1},{1,1,1,1,0,0,0,0},{1,0,1,1,1,0,0,0},{0,1,1,1,1,0,0,0},{0,1,0,1,1,1,0,0},{0,0,1,1,1,1,0,0},{0,0,0,1,1,1,1,0},{0,0,0,0,1,1,1,1},{1,1,1,1,1,0,0,0},{1,0,1,1,1,1,0,0},{0,1,1,2,1,0,0,0},{0,1,1,1,1,1,0,0},{0,1,0,1,1,1,1,0},{0,0,1,1,1,1,1,0},{0,0,0,1,1,1,1,1},{1,1,1,2,1,0,0,0},{1,1,1,1,1,1,0,0},{1,0,1,1,1,1,1,0},{0,1,1,2,1,1,0,0},{0,1,1,1,1,1,1,0},{0,1,0,1,1,1,1,1},{0,0,1,1,1,1,1,1},{1,1,2,2,1,0,0,0},{1,1,1,2,1,1,0,0},{1,1,1,1,1,1,1,0},{1,0,1,1,1,1,1,1},{0,1,1,2,2,1,0,0},{0,1,1,2,1,1,1,0},{0,1,1,1,1,1,1,1},{1,1,2,2,1,1,0,0},{1,1,1,2,2,1,0,0},{1,1,1,2,1,1,1,0},{1,1,1,1,1,1,1,1},{0,1,1,2,2,1,1,0},{0,1,1,2,1,1,1,1},{1,1,2,2,2,1,0,0},{1,1,2,2,1,1,1,0},{1,1,1,2,2,1,1,0},{1,1,1,2,1,1,1,1},{0,1,1,2,2,2,1,0},{0,1,1,2,2,1,1,1},{1,1,2,3,2,1,0,0},{1,1,2,2,2,1,1,0},{1,1,2,2,1,1,1,1},{1,1,1,2,2,2,1,0},{1,1,1,2,2,1,1,1},{0,1,1,2,2,2,1,1},{1,2,2,3,2,1,0,0},{1,1,2,3,2,1,1,0},{1,1,2,2,2,2,1,0},{1,1,2,2,2,1,1,1},{1,1,1,2,2,2,1,1},{0,1,1,2,2,2,2,1},{1,2,2,3,2,1,1,0},{1,1,2,3,2,2,1,0},{1,1,2,3,2,1,1,1},{1,1,2,2,2,2,1,1},{1,1,1,2,2,2,2,1},{1,2,2,3,2,2,1,0},{1,2,2,3,2,1,1,1},{1,1,2,3,3,2,1,0},{1,1,2,3,2,2,1,1},{1,1,2,2,2,2,2,1},{1,2,2,3,3,2,1,0},{1,2,2,3,2,2,1,1},{1,1,2,3,3,2,1,1},{1,1,2,3,2,2,2,1},{1,2,2,4,3,2,1,0},{1,2,2,3,3,2,1,1},{1,2,2,3,2,2,2,1},{1,1,2,3,3,2,2,1},{1,2,3,4,3,2,1,0},{1,2,2,4,3,2,1,1},{1,2,2,3,3,2,2,1},{1,1,2,3,3,3,2,1},{2,2,3,4,3,2,1,0},{1,2,3,4,3,2,1,1},{1,2,2,4,3,2,2,1},{1,2,2,3,3,3,2,1},{2,2,3,4,3,2,1,1},{1,2,3,4,3,2,2,1},{1,2,2,4,3,3,2,1},{2,2,3,4,3,2,2,1},{1,2,3,4,3,3,2,1},{1,2,2,4,4,3,2,1},{2,2,3,4,3,3,2,1},{1,2,3,4,4,3,2,1},{2,2,3,4,4,3,2,1},{1,2,3,5,4,3,2,1},{2,2,3,5,4,3,2,1},{1,3,3,5,4,3,2,1},{2,3,3,5,4,3,2,1},{2,2,4,5,4,3,2,1},{2,3,4,5,4,3,2,1},{2,3,4,6,4,3,2,1},{2,3,4,6,5,3,2,1},{2,3,4,6,5,4,2,1},{2,3,4,6,5,4,3,1},{2,3,4,6,5,4,3,2}};

(*The Cartan matrix of E8*)
CartanH={{2,0,-1,0,0,0,0,0},{0,2,0,-1,0,0,0,0},{-1,0,2,-1,0,0,0,0},{0,-1,-1,2,-1,0,0,0},{0,0,0,-1,2,-1,0,0},{0,0,0,0,-1,2,-1,0},{0,0,0,0,0,-1,2,-1},{0,0,0,0,0,0,-1,2}};
MatrixForm[CartanH]

(*The root of E8*)
RH= Join[PRH,-PRH];
(* Odd and even parts of grading *)
ROdd = {}; REven={};
For[i=1,i≤Length[RH],i++,
If[OddQ[{1,1,1,1,1,1,1,1}.RH[[i]]],AppendTo[ROdd,RH[[i]]],AppendTo[REven,RH[[i]]]
];
];

(*A basis of simple roots for G, as linear combinations of the simple roots of H = E7*)
SG={{0,1,1,1,1,0,0,0},{0,0,0,0,0,1,1,0},{0,0,0,1,1,0,0,0},{1,0,1,0,0,0,0,0},{0,1,0,1,0,0,0,0},{0,0,0,0,1,1,0,0},{0,0,0,0,0,0,1,1},{0,0,1,1,0,0,0,0}};
(*The Cartan matrix of G - type D8 in Bourbaki numbering*)
CartanG = SG.CartanH.Transpose[SG];
MatrixForm[CartanG]

(*The weights of V (i.e. elements of ROdd) in the basis .5*SG*)
vwts=2*ROdd.Inverse[SG];
(*Positive roots of G in the basis SG*)
PRG={{1,0,0,0,0,0,0,0},{0,1,0,0,0,0,0,0},{0,0,1,0,0,0,0,0},{0,0,0,1,0,0,0,0},{0,0,0,0,1,0,0,0},{0,0,0,0,0,1,0,0},{0,0,0,0,0,0,1,0},{0,0,0,0,0,0,0,1},{1,1,0,0,0,0,0,0},{0,1,1,0,0,0,0,0},{0,0,1,1,0,0,0,0},{0,0,0,1,1,0,0,0},{0,0,0,0,1,1,0,0},{0,0,0,0,0,1,1,0},{0,0,0,0,0,1,0,1},{1,1,1,0,0,0,0,0},{0,1,1,1,0,0,0,0},{0,0,1,1,1,0,0,0},{0,0,0,1,1,1,0,0},{0,0,0,0,1,1,1,0},{0,0,0,0,1,1,0,1},{0,0,0,0,0,1,1,1},{1,1,1,1,0,0,0,0},{0,1,1,1,1,0,0,0},{0,0,1,1,1,1,0,0},{0,0,0,1,1,1,1,0},{0,0,0,1,1,1,0,1},{0,0,0,0,1,1,1,1},{1,1,1,1,1,0,0,0},{0,1,1,1,1,1,0,0},{0,0,1,1,1,1,1,0},{0,0,1,1,1,1,0,1},{0,0,0,1,1,1,1,1},{0,0,0,0,1,2,1,1},{1,1,1,1,1,1,0,0},{0,1,1,1,1,1,1,0},{0,1,1,1,1,1,0,1},{0,0,1,1,1,1,1,1},{0,0,0,1,1,2,1,1},{1,1,1,1,1,1,1,0},{1,1,1,1,1,1,0,1},{0,1,1,1,1,1,1,1},{0,0,1,1,1,2,1,1},{0,0,0,1,2,2,1,1},{1,1,1,1,1,1,1,1},{0,1,1,1,1,2,1,1},{0,0,1,1,2,2,1,1},{1,1,1,1,1,2,1,1},{0,1,1,1,2,2,1,1},{0,0,1,2,2,2,1,1},{1,1,1,1,2,2,1,1},{0,1,1,2,2,2,1,1},{1,1,1,2,2,2,1,1},{0,1,2,2,2,2,1,1},{1,1,2,2,2,2,1,1},{1,2,2,2,2,2,1,1}};

(*The partial order on weights, with respect to the root basis SG*)
GeQ[x_,y_]:=(
For[q=1,q<=8,q++,
If[vwts[[x]][[q]]<vwts[[y]][[q]],
Return [False]
]
];
Return[True]
); 
LLeQ[x_,y_]:=(
For[q=1,q<=8,q++,
If[vwts[[x]][[q]]>vwts[[y]][[q]],
Return [False]
]
]; 
Return[True]
); 
(*For a subset M_0 of {1..128} (identified with the rows of vwts), the function min returns the set of upper bounds of {1..128} - M_0 with respect to the above partial order*)
min[x_]:=( 
l = Complement[Range[128],x];
z={}; 
For[k=1,k<Length[l]+1,k++,
For[p=1,p<Length[l]+1,p++,
If[GeQ[l[[k]],l[[p]]] && (k ≠p), 
AppendTo[z, l[[p]]]
]
]
]; 
For[k=1,k<Length[z]+1,k++, 
l=DeleteCases[l, z[[k]]]
]; 
Return[l];
);
(*We verify that the highest weight of V equals the highest weight of H*)
Print[ROdd[[min[{}][[1]]]]=={2,3,4,6,5,4,3,2}]

(*The following function returns True if a set M_0 satisfies 2. of Proposition 2.15*)
UnstableQ[x_]:=(
y=Complement[Range[128],x];
Return[Reduce[
Exists[{b1,b2,b3,b4,b5,b6,b7,b8},{b1,b2,b3,b4,b5,b6,b7,b8}≠{0,0,0,0,0,0,0,0}&&AllTrue[Range[Length[y]],{b1,b2,b3,b4,b5,b6,b7,b8}.vwts[[y[[#]]]]≤0&]
]
]
];
);

(*The following function returns True if a set M_0 satisfies 1. of Proposition 2.15*)
KostantQ[x_]:= (
If[ContainsAll[x, Range[9,64]],Return[True],Return[False];];
);
(*The following function of (M0, M1) returns True if a set M_0 satisfies 3. of Proposition 2.15 for some alpha in M1*)
sl2UnstableQ[x_]:=(
M0=x[[1]];
M1=x[[2]];
M2=Complement[Range[128], M0];
rootlist={};
For[i=1,i≤8,i++,
counter=0;
For[j=1,j≤Length[M0],j++,If[MemberQ[vwts[[M2]],vwts[[M0[[j]]]]+Normal[SparseArray[{i->2},8]]]∨MemberQ[vwts[[M2]],vwts[[M0[[j]]]]-Normal[SparseArray[{i->2},8]]],counter=1]];
If[counter==0,AppendTo[rootlist,i]];];
goodrootlist={};
For[i=1,i≤Length[rootlist],i++,counter=0;
For[j=1,j≤Length[M1],j++,If[MemberQ[vwts[[M2]],vwts[[M1[[j]]]]-Normal[SparseArray[{rootlist[[i]]->2},8]]]&&UnstableQ[Join[M0,{M1[[j]]}]],counter=1];];
If[counter==1,AppendTo[goodrootlist,rootlist[[i]]]];];
If[Length[goodrootlist]>0,Return[True],Return[False]];
);

(* 2 * the character \delta_G which appears in the description of the Haar measure on G in terms of the Iwasawa decomposition*)
delta=2*ConstantArray[1,56].PRG;

(*We now generate a list of cusp data, using UnstableQ to remove those (M_0, M_1) for which a cocharacter witnesses the non-stability of the elements of V(M_0) (observe the existence of a cocharacter is hereditary under inclusions M_0 \subset M_0')*)
allcusps={{{64},min[{64}],delta-vwts[[64]]}};(*The final list is stored in allcusps. allcusps is initialised with the cusp datum corresponding to the highest weight of V = highest root of H*)
case1={{{64},min[{64}]}};
marker=1; (*Set to 0 once we create no new cusp data in a given step*)
timesaround=1; (*The size of the sets M_0 created during the current step*)
counter=0; (*The number of cusp data so far added to the list*)
While[marker==1,
(marker=0;
case={};
For[i=1,i<Length[case1]+1,i++, 
For[j=1,j<Length[case1[[i]][[2]]]+1,j++,
temp=Sort[Append[case1[[i]][[1]],case1[[i]][[2]][[j]]]]; 
If[¬ UnstableQ[temp],(AppendTo[case,{temp,min[temp]}];marker=1;counter=counter+1; )
];
];
];
 case1=DeleteDuplicates[case];
For[i=1,i<Length[case1]+1,i++,temp=delta-Sum[vwts[[case1[[i]][[1]][[k]]]],{k,1,Length[case1[[i]][[1]]]}];
AppendTo[allcusps,Join[case1[[i]],{temp}]]
];
);
timesaround = timesaround+1;
];

(*200 minutes*)
Length[allcusps]

(*We get 9437 cusp data. This is the number which appears in the introduction to our paper*)

(*The next function tests a cusp datum (M_0, M_1)) for the existence of a function f(a) (a in M_1) that witnesses a power-saving on the number of points in S(M_0, M_1) of bounded height*)
H[c_]:=(
m=Length[c[[2]]];
var=Take[{a1,a2,a3,a4,a5,a6,a7,a8},m];
l=ConstantArray[1,m];
Return [Reduce[Exists[{a1,a2,a3,a4,a5,a6,a7,a8},If[Sum[If[(var.vwts[[c[[2]]]]+c[[3]])[[i]]>0,1,0],{i,1,8}]==8,True,False]&&(var.l<Length[c[[1]]])]]]
);

(*We test all the elements of allcusps to see whether either they satisfy the condition for existence of a function f; or if they satisfy the remaining conditions 1. or 3. of Proposition 2.15*)
For[w=1,w≤Length[allcusps],w++,
If[¬(H[allcusps[[w]]]∨ sl2UnstableQ[allcusps[[w]]]∨KostantQ[allcusps[[w]][[1]]]),Print[w, " ", allcusps[[w]]]];
];

(*We find that all cusp data have now been eliminated!*)