
// Weights of the minimal canonical Fano simplices up to dim 3
weights:=[
	[
		[1,1]
	],
	[
		[1,1,1],[1,1,2]
	],
	[
		[1,1,1,1],[1,1,1,2],[1,1,1,3],[1,1,2,2],[1,1,2,3],
		[1,1,2,4],[1,1,3,4],[1,1,3,5],[1,1,4,6],[1,2,3,5],
		[1,3,4,5],[2,3,5,7],[3,4,5,7]
	]
];

// Number and dimensions of the simplices in the decomposition
t:=3;
di:=[2,2,2];

// Create the lattices
Ni:=[ToricLattice(d+1) : d in di];
N,embs,projs:=DirectSum(Ni);

// All the permutations of a sequence
function SequencePermutations(s)
	S:=[[s[i],i] : i in [1..#s]];
	return {[p[j,1] : j in [1..#p]] : p in Permutations(SequenceToSet(S))};
end function;

// All possible identifications between the vertices of the (three, in this
// case) simplices appearing in the decomposition
all_Vij:={};
for aa in SequencePermutations([0..di[1]]) do
	for bb in SequencePermutations([0..di[2]]) do
		Vij:=[[],[]];
		Vij[1,2]:=&cat[[[aa[i],bb[i]]] : i in [1..di[1]-1]];
		for cc in SequencePermutations([0..di[2]]) do
			for dd in SequencePermutations([0..di[1]+di[2]+1]) do
				Vij[1,3]:=&cat[[[dd[i],cc[i]]] : i in [1..di[1]-1] | dd[i] le di[1]];
				Vij[2,3]:=&cat[[[dd[i]-di[1]-1,cc[i]]] : i in [1..di[1]-1] | dd[i] gt di[1]];
				Include(~all_Vij,Vij);
			end for;
		end for;
	end for;
end for;

// Return true if the polytope P is minimal
function is_minimal(P)
    pts:=Points(P);
    for v in Vertices(P) do
        Q:=Polytope(SetToSequence(Exclude(pts,v)));
        if Dimension(Q) eq Dimension(P) and IsFano(Q) then
            return false;
        end if;
    end for;
    return true;
end function;

// Check that the built polytope P has the parameters we need
function is_good(P)
    if not Dimension(P) eq di[1]+t-1 or not IsIntegral(P) or
            not NumberOfInteriorPoints(P) eq 1 or not IsFano(P) then
        return false;
    end if;
    if not is_minimal(P) then
        return false;
    end if;
    return true;
end function;

// Return the volume of the dual of a WPS given its weights wts
function vol_from_wts(wts)
	return (&+wts)^(#wts-1)/&*wts;
end function;

// Returns the i-th value of the Sylvester sequence
function sylv(i)
	if i eq 1  then
		return 2;
	else
		return &*[sylv(j) : j in [1..i-1]]+1;
	end if;
end function;

// Build of the possible P and store them in the set "list"
tot:=#all_Vij;
ii:=0;
list:={};
for Vij in all_Vij do
	V:=&cat[&cat[[embs[i](Basis(Ni[i])[v[1]+1])-embs[j](Basis(Ni[j])[v[2]+1]) : v in Vij[i,j]] : j in [i+1..t]] : i in [1..t]];
	for w1 in weights[di[1]] do
		for w2 in [w : w in weights[di[2]] | vol_from_wts(w) le vol_from_wts(w1)] do
			for w3 in [w : w in weights[di[3]] | vol_from_wts(w) le vol_from_wts(w2)] do
				// define the quotient
				w:=[w1,w2,w3];
				W:=[embs[i](Ni[i]!w[i]) : i in [1..t]];
				_,phi:=Quotient(W cat V);
				// build P
				PP:=Polytope(&cat[embs[i](Basis(Ni[i])) : i in [1..t]]);
				P:=Image(phi,PP);
				if is_good(P) then 
					Include(~list,[Eltseq(v) : v in NormalForm(P)]);
				end if;
			end for;
		end for;
	end for;
	ii+:=1;
	printf "%o/%o - found: %o \n",ii,tot,#list;
end for;

// It returns the following list of 4 4-dimensional polytopes:
//[[1,0,0,0],[0,1,0,0],[0,0,1,0],[-1,0,-1,0],[-1,-1,0,0],[0,0,0,1],[-2,0,0,-1]]
//[[1,0,0,0],[0,1,0,0],[0,0,1,0],[-1,0,-1,0],[-1,-1,0,0],[0,0,0,1],[-1,0,0,-1]]
//[[1,0,0,0],[0,1,0,0],[-1,-1,0,0],[0,0,1,0],[-2,0,-1,0],[0,0,0,1],[-2,0,0,-1]]
//[[1,0,0,0],[0,1,0,0],[0,0,1,0],[-2,0,-1,0],[-2,-1,0,0],[0,0,0,1],[-2,0,0,-1]]
//For each of them we verify their dual volume satisfies Theorem 1.1
&and[Volume(Dual(Polytope(s))) le 2*(sylv(4)-1)^2 : s in list];
// It returns true