Using FindInstance for Identifying Feasible Planar Solutions in a 3D Latin Hypercube

Question

Consider a $3d$ lattice latin hypercube with $n$ steps in each dimension, so it has $n^3$ positions. Coordinates $X, Y, Z \in \{1,2,...n\}$. I want to find all of the permutations of them where they add to $S=3(n+1)/2$.

The number of permutations is only non-zero when $n$ is odd, and is A002047 from the Sloane sequence library. I've brute forced my way to generating the permutations via Permutations and Subsets and it works alright, but wanted to see if FindInstance could do better. It's easier to define the problem.

Here is what I have:

nn = 3;
ss = 3 (num + 1)/2;
yvars = Table[y[i], {i, nn}]; (* instantiate the y[] and z[] variables *)
zvars = Table[z[i], {i, nn}];
xvars = Table[x[i] = i, {i, nn}]; (* just set these to {1,...,nn} *)

(* Force Z to be a permutation *)
betweenZ = 1 <= zvars <= nn; 
unequalZ = Unequal @@ zvars;

(* same for Y *)
betweenY = 1 <= yvars <= nn;
unequalY = Unequal @@ yvars;

(* force the sums to equal S *)
sumeqns = xvars + yvars + zvars == Table[ss, {nn}]; (* sum properly *)

rules = FindInstance[{betweenZ,unequalZ, betweenY, unequalY, sumeqns}, 
                      Join @@ {yvars, zvars}, 
                      Integers, 
                      20]
 Length@rules
 (* 2 *)

This works fine, except it is slow. With $n=7$ it takes 22 seconds on my machine to get a single permutation. Using my brutish force approach, I can get all 244 permutations for when $n=9$ in about 3 seconds.

I was wondering if there were other ways to pass the constraints into FindInstance that might be more suitable for processing in its innards. Thoughts?

Edit

Turns out this is the same problem as explored in 1966, where the authors used CDC3200 computer to explore the number of permutations.

ARRAYS AND BROOKS B. T. BENNETT and R. B. POTTS

https://www.cambridge.org/core/services/aop-cambridge-core/content/view/4DF66D255AB9ACA1E9CE9A0034490105/S144678870000505Xa.pdf/arrays_and_brooks.pdf

Answer 1

One way to get a single solution faster is to set it up as a feasibility problem and use FindMinimum. This in turn uses fast machine precision linear programming under the hood.

I recast as a 0-1 ILP. So each y variable in the original formulation becomes a vector, with component j equal to 1 if the original had that var set to j, and other components are 0.

nn = 9;
ss = 3 (nn + 1)/2;
yvars = Array[y, {nn, nn}];
zvars = Array[z, {nn, nn}];
allvars = Flatten[{yvars, zvars}];
cons = Flatten[{Table[
     j + yvars[[j]].Range[nn] + zvars[[j]].Range[nn] == ss, {j, nn}], 
    Map[0 <= # <= 1 &, allvars], Thread[Total[yvars] == 1], 
    Thread[Total[zvars] == 1], Thread[Total[Transpose@yvars] == 1], 
    Thread[Total[Transpose@zvars] == 1]}];

Timing[{min, vals} = 
   FindMinimum[{1, Join[cons, {Element[allvars, Integers]}]}, 
    allvars];]

(* Out[87]= {0.12, Null} *)

Now one recovers the permutations as below.

Transpose[Position[yvars /. vals, 1]]
Transpose[Position[zvars /. vals, 1]]

(* Out[94]= {{1, 2, 3, 4, 5, 6, 7, 8, 9}, {7, 5, 9, 2, 6, 8, 3, 1, 4}}

Out[95]= {{1, 2, 3, 4, 5, 6, 7, 8, 9}, {7, 8, 3, 9, 4, 1, 5, 6, 2}} *)

In the original formulation, this would mean: y[1]==7,y[2]==5,... and z[1]==7,z[2]==8,... (no need to transpose if it is easier to see that way, I just did so because the above is a common way to represent a permutation).