d=3; (* ambient dimension *)
k=3; (* size of the considered hypercube *)
p=2; (* number of vertices of P *)
q=2; (* number of vertices of Q *)

(* computes the list of difference vectors *)

ListDiff={};

For[n=1,n<(k+1)^(p+q),n++,

(* computes a pair (p,q) *)

m=n;
pq={};
cte=0;
For[i=1,i≤p+q,i++,
mm=Mod[m,k+1];
m=(m-mm)/(k+1);
pq=Append[pq,mm];
];

(* computes the difference vector from pq *)
If[pq[[p]]≤k/2,
diff={};
For[i=1,i<=p+q,i++,
If[i<p,
diff=Append[diff,pq[[i]]-pq[[p]]];
];
If[i>p && i<p+q,
diff=Append[diff,-pq[[i]]+pq[[p+q]]];
];
If[i==p+q,
diff=Append[diff,pq[[p]]-pq[[p+q]]];
];
];

(* only consider this row if the p+q-2 first coordinates are not all zero *)

f=0;
For[i=1,i<=p+q-2,i++,
If[diff[[i]]!=0,f=1];
];

If[f==1,

(* negates the row if the first non-zero coefficient is negative *)

i=1;
While[diff[[i]]==0,
i++;
];
If[diff[[i]]<0,
diff=-diff;
];

ListDiff=Append[ListDiff,diff];
];
];
];

(* discards the difference vectors that are not primitive or whose first non-zero coordinate is negative *)

For[i=1,i≤Length[ListDiff],i++,
ListDiff[[i]]=ListDiff[[i]]/Apply[GCD,ListDiff[[i]]];
];
ListDiff=DeleteDuplicates[ListDiff];
Print["Number of difference vectors: ",Length[ListDiff]];

(* extracts the left-hand side from the difference vectors *)
ListDiffRHS={};
For[i=1,i≤Length[ListDiff],i++,
ListDiffRHS=Append[ListDiffRHS,ListDiff[[i]][[1;;p+q-2]]];
];
ListDiffRHS=DeleteDuplicates[ListDiffRHS];
Print["Number of left-hand sides of difference vectors: ",Length[ListDiffRHS]];

(* Main computation, see article *)
(* The function to optimize has been explicitely computed when d=3 and 4 for effciency *)

ListOpt=Subsets[ListDiff,{d}];
Print["Number of linear systems: ",Length[ListOpt]];

SingularMatrices=0;

MinGlob=1000;
imin=-1;
For[i=1,i≤Length[ListOpt],i++,
If[Mod[i,10000]==0,
Print["Done: ",i,"/", Length[ListOpt], " Minimal Distance: ",MinGlob];
];

(* Here we treat separately the cases when we have a formula (see BuildFormula.nb) *)

If[d==3,
a11=ListOpt[[i]][[1]][[1]];
a12=ListOpt[[i]][[1]][[2]];
a21=ListOpt[[i]][[2]][[1]];
a22=ListOpt[[i]][[2]][[2]];
a31=ListOpt[[i]][[3]][[1]];
a32=ListOpt[[i]][[3]][[2]];
det=(a12^2) ((a21^2)+(a31^2))+(((a22 a31-a21 a32))^2)-2 a11 a12 (a21 a22+a31 a32)+(a11^2) ((a22^2)+(a32^2));
If[det≠0,
b1=ListOpt[[i]][[1]][[3]];
b2=ListOpt[[i]][[2]][[3]];
b3=ListOpt[[i]][[3]][[3]];
LL=Abs[a22 a31 b1-a21 a32 b1-a12 a31 b2+a11 a32 b2+a12 a21 b3-a11 a22 b3]/Sqrt[det];
If[LL≠0&&LL<MinGlob,
MinGlob=LL;
imin=i;
];
,
SingularMatrices=SingularMatrices+1;
];
,
If[d==4,
a11=ListOpt[[i]][[1]][[1]];
a12=ListOpt[[i]][[1]][[2]];
a13=ListOpt[[i]][[1]][[3]];
a21=ListOpt[[i]][[2]][[1]];
a22=ListOpt[[i]][[2]][[2]];
a23=ListOpt[[i]][[2]][[3]];
a31=ListOpt[[i]][[3]][[1]];
a32=ListOpt[[i]][[3]][[2]];
a33=ListOpt[[i]][[3]][[3]];
a41=ListOpt[[i]][[4]][[1]];
a42=ListOpt[[i]][[4]][[2]];
a43=ListOpt[[i]][[4]][[3]];
det=(a11^2) (a23^2) (a32^2)-2 (a11^2) a22 a23 a32 a33+(a11^2) (a22^2) (a33^2)+(a23^2) (a32^2) (a41^2)-2 a22 a23 a32 a33 (a41^2)+(a22^2) (a33^2) (a41^2)-2 (a23^2) a31 a32 a41 a42+2 a22 a23 a31 a33 a41 a42+2 a21 a23 a32 a33 a41 a42-2 a21 a22 (a33^2) a41 a42+(a11^2) (a23^2) (a42^2)+(a23^2) (a31^2) (a42^2)-2 a21 a23 a31 a33 (a42^2)+(a11^2) (a33^2) (a42^2)+(a21^2) (a33^2) (a42^2)+(a13^2) ((a22^2) ((a31^2)+(a41^2))+(((a32 a41-a31 a42))^2)-2 a21 a22 (a31 a32+a41 a42)+(a21^2) ((a32^2)+(a42^2)))-2 ((a22 a31-a21 a32) (-a23 a32+a22 a33) a41+(a22 a31-a21 a32) (a23 a31-a21 a33) a42+(a11^2) (a22 a23+a32 a33) a42) a43+((((a22 a31-a21 a32))^2)+(a11^2) ((a22^2)+(a32^2))) (a43^2)+(a12^2) ((a23^2) ((a31^2)+(a41^2))+(((a33 a41-a31 a43))^2)-2 a21 a23 (a31 a33+a41 a43)+(a21^2) ((a33^2)+(a43^2)))+2 a11 a12 (-(a23^2) (a31 a32+a41 a42)-(a33 a41-a31 a43) (a33 a42-a32 a43)+a23 (a22 a31 a33+a21 a32 a33+a22 a41 a43+a21 a42 a43)-a21 a22 ((a33^2)+(a43^2)))+2 a13 (a12 (-a22 a23 ((a31^2)+(a41^2))+a21 a23 (a31 a32+a41 a42)-(a32 a41-a31 a42) (a33 a41-a31 a43)+a21 a22 (a31 a33+a41 a43)-(a21^2) (a32 a33+a42 a43))+a11 (-a21 a23 ((a32^2)+(a42^2))-(a32 a41-a31 a42) (-a33 a42+a32 a43)-(a22^2) (a31 a33+a41 a43)+a22 (a23 a31 a32+a21 a32 a33+a23 a41 a42+a21 a42 a43)));
If[det≠0,
b1=ListOpt[[i]][[1]][[4]];
b2=ListOpt[[i]][[2]][[4]];
b3=ListOpt[[i]][[3]][[4]];
b4=ListOpt[[i]][[4]][[4]];
LL=Abs[-a21 a33 a42 b1+a21 a32 a43 b1+a13 a32 a41 b2-a12 a33 a41 b2-a13 a31 a42 b2+a11 a33 a42 b2+a12 a31 a43 b2-a11 a32 a43 b2+a13 a21 a42 b3-a12 a21 a43 b3-a13 a21 a32 b4+a12 a21 a33 b4+a23 (-a32 a41 b1+a31 a42 b1+a12 a41 b3-a11 a42 b3-a12 a31 b4+a11 a32 b4)+a22 (a33 a41 b1-a31 a43 b1-a13 a41 b3+a11 a43 b3+a13 a31 b4-a11 a33 b4)]/Sqrt[det];
If[LL≠0&&LL<MinGlob,
MinGlob=LL;
imin=i;
];
,
SingularMatrices=SingularMatrices+1;
];
,
Mat=ListOpt[[i]][[1;;d,1;;p+q-2]];
b=-ListOpt[[i]][[1;;d,p+q-1]];
MM=Transpose[Mat].Mat;
If[MatrixRank[MM]==d-1,
sol=(Inverse[MM].Transpose[Mat]).b;
f=1;
For[j=1,j≤Length[sol],j++,
If[sol[[j]]<0||sol[[j]]>1,
f=0;
];
];
If[f==1,
LL=Norm[Mat.sol-b];
If[LL≠0&&LL<MinGlob,
MinGlob=LL;
imin=i;
];
];
,
SingularMatrices=SingularMatrices+1;
];
];
];
];
Print["epsilon(",d,",",k,")=",MinGlob];
Print["Optimal configuration: ",ListOpt[[imin]]];
Opt=ListOpt[[imin]];
Print["Number of singular matrices: ",SingularMatrices];