interface(quiet=true):
# graph number 3 in Table 1 on page 8 of [Chua-Kummer-St],

with(linalg): with(simplex):
die := rand(-100..100):

#choose a cycle basis for the totally cyclic graph, define the matrix
B := array([
[1,1, 1, 0, 0, 0, 0, 0],
[1, 0, 1, 1, 1, 0, 0, 0],
[ 0, 0, 1, 0, 1, 1, 1, 0],
[0,0, 0, 0, 0, 0, 1, 1]]):

#print Riemann matrix and cycle basis matrix
print(multiply(B,transpose(B)), Riemann);
print(B);

#define vector for defining linear functionals later
X := array([[x,y,z,w]]):

#this finds the vertices in the Voronoi polytope (stored in points)
nonunique := {}:
points := {}:
for i1 in {-1,1} do
for i2 in {-1,1} do
for i3 in {-1,1} do
for i4 in {-1,1} do
for i5 in {-1,1} do
for i6 in {-1,1} do
for i7 in {-1,1} do
for i8 in {-1,1} do

M := multiply(X,multiply(B,array([[i1],[i2],[i3],[i4],[i5],[i6], [i7], [i8]])))[1,1]:
if (member(M,points)) then nonunique := nonunique union {M}:
    else points := points union {M}: fi:
od:od:od:od:od:od:od:od:

points := points minus nonunique:
numelems(points);

#polyts will store the different resulting Delaunay polytopes (as lists of vertices in Polymake-readable format)
#nums[i] will store the number of polytopes with i vertices
polyts := {}:
nums := Vector[row](16):

#this finds the Delaunay polytopes corresponding to each vertex of the Voronoi polytope.
for i1 in {-1,1} do
for i2 in {-1,1} do
for i3 in {-1,1} do
for i4 in {-1,1} do
for i5 in {-1,1} do
for i6 in {-1,1} do
for i7 in {-1,1} do
for i8 in {-1,1} do

M := multiply(X,multiply(B,array([[i1],[i2],[i3],[i4],[i5],[i6], [i7], [i8]])))[1,1]:
if (member(M,points)) then
lprint(i1,i2,i3,i4,i5,i6, i7, i8); #which vertex the Delaunay polytope corresponds to
sys := {}: #list of inequalities satisfied by our polytope
j:=1:ii:=i1: sys := sys union {ii*(B[1,j]*x+B[2,j]*y+B[3,j]*z+B[4,j]*w) >= 0, ii*(B[1,j]*x+B[2,j]*y+B[3,j]*z+B[4,j]*w) <= 1}:
j:=2:ii:=i2: sys := sys union {ii*(B[1,j]*x+B[2,j]*y+B[3,j]*z+B[4,j]*w) >= 0, ii*(B[1,j]*x+B[2,j]*y+B[3,j]*z+B[4,j]*w) <= 1}:
j:=3:ii:=i3: sys := sys union {ii*(B[1,j]*x+B[2,j]*y+B[3,j]*z+B[4,j]*w) >= 0, ii*(B[1,j]*x+B[2,j]*y+B[3,j]*z+B[4,j]*w) <= 1}:
j:=4:ii:=i4: sys := sys union {ii*(B[1,j]*x+B[2,j]*y+B[3,j]*z+B[4,j]*w) >= 0, ii*(B[1,j]*x+B[2,j]*y+B[3,j]*z+B[4,j]*w) <= 1}:
j:=5:ii:=i5: sys := sys union {ii*(B[1,j]*x+B[2,j]*y+B[3,j]*z+B[4,j]*w) >= 0, ii*(B[1,j]*x+B[2,j]*y+B[3,j]*z+B[4,j]*w) <= 1}:
j:=6:ii:=i6: sys := sys union {ii*(B[1,j]*x+B[2,j]*y+B[3,j]*z+B[4,j]*w) >= 0, ii*(B[1,j]*x+B[2,j]*y+B[3,j]*z+B[4,j]*w) <= 1}:
j:=7:ii:=i7: sys := sys union {ii*(B[1,j]*x+B[2,j]*y+B[3,j]*z+B[4,j]*w) >= 0, ii*(B[1,j]*x+B[2,j]*y+B[3,j]*z+B[4,j]*w) <= 1}:
j:=8:ii:=i8: sys := sys union {ii*(B[1,j]*x+B[2,j]*y+B[3,j]*z+B[4,j]*w) >= 0, ii*(B[1,j]*x+B[2,j]*y+B[3,j]*z+B[4,j]*w) <= 1}:
vertices := {}:
lprint(sort(sys));
for run from 1 to 200 do
vertices := vertices union {subs(minimize(die()*x+die()*y+die()*z+die()*w,sys),[1,x,y,z,w])}:
od:
verts := nops(vertices):
nums[verts] := nums[verts] + 1:
polyt := [nops(vertices), sort(vertices)]: #a polytope is stored as the number of vertices, followed by a list of vertices in polymake readable format
polyts := polyts union {polyt}:
print(polyt);
print(` `);
fi:
od:od:od:od:od:od:od:od:

#print all polytopes -- these can then be plugged into polymake to find out which ones they correspond to
for polyt in polyts do
print(polyt);
print(`   `);
od:

#print nums, which tells us how many polytopes of each type we have. In these cases, the Delaunay polytopes with the same number of vertices which arise from the same voronoi polytope are the same.
print(nums);