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;

// 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]
	]
];

// Make them into barycentric coordinates
weights:=[[[w[i]/&+[w[i] : i in [1..d+1]] : i in [1..d+1]] : w in weights[d]] : d in [1..3]];

// barycentric coordinates of 4-simplices
file:="minimal4.txt";
alls4:=[eval(s) : s in FileProcess(file)];
dweights:={};
for s in alls4 do
  w:=Eltseq(Kernel(Matrix(s)).1);
  sum:=&+[w[i] : i in [1..5]];
  w:=[w[i]/sum : i in [1..5]];
  Include(~dweights,w);
end for;
weights[4]:=SetToSequence(dweights);

// 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;

// 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;

// Given two (ordered) baricentric coordinates w1 and w2 it calculates the
// right hand side of 5.7
function bound(w1,w2)
  d1:=#w1-1;
  d2:=#w2-1;
  d:=d1+1;
  q:=d2-1;
  // The baricentric coordinates of the common vertices
  V1:=w1[1..q];
  V2:=w2[1..q];
  // We can exclude the lowest of the remaining barycentric coordinates.
  // In the paper this is done when we coarsen the lattice
  W1:=Remove(Sort(w1[q+1..d1+1]),1);
  W2:=Remove(Sort(w2[q+1..d2+1]),1);
  return Factorial(d)*(1/2*(Factorial(2*(d1-q))/(Factorial(q+2*(d1-q))*(Factorial(d1-q))^2))*&*[1/w : w in V1]*&*[1/(w^2) : w in W1] +  1/2*(Factorial(2*(d2-q))/(Factorial(q+2*(d2-q))*(Factorial(d2-q))^2))*&*[1/w : w in W2]*&*[1/(w^2) : w in W2]);
end function;

// case d1=3 d2=2
d:=4;
errors:=[];
i:=0;
for w3 in weights[3] do
	i+:=1;
	i;
	for w2 in weights[2] do
		for ww3 in SequencePermutations(w3) do
			for ww2 in SequencePermutations(w2) do
				vol:=bound(ww3,ww2);
				if vol ge 2*(sylv(d)-1)^2 then
					Include(~errors,[ww3,ww2]);
				end if;
			end for;
		end for;
	end for;
end for;

#errors;
// Returns 0

// case d1=4 d2=3
d:=5;
errors:=[];
i:=0;
for w4 in weights[4] do
	i+:=1;
	printf "%o - %o\n",i,#errors;
	for w3 in weights[3] do
		for ww4 in SequencePermutations(w4) do
			for ww3 in SequencePermutations(w3) do
				vol:=bound(ww4,ww3);
				if vol ge 2*(sylv(d)-1)^2 then
					Include(~errors,[ww4,ww3]);
				end if;
			end for;
		end for;
	end for;
end for;

#errors;
// Returns 0



