My kid’s math teacher recently posted a puzzle for the parents that got me down a rabbit hole of sorts and made me reflect on a number of things. This post is about that … and I mix the math and programming bits. I describe the puzzle in the first section and then the rest reveal aspects of the solution. So if you want to avoid the spoilers, you should stop reading after the first section.
The Puzzle
Consider an \(N \times N\) square grid, where each cell may contain a letterdigit combination  like \(A1, A2, B4, D5\) and so on. The letters are to be drawn from \(ABCDEFGHI\) and the digits from \(123456789\). For the purpose of this discussion, we’ll limit ourselves to \(N \leq 9\), but you can consider larger grids on your own by extension.
The task is to fill the grid with combinations such that 
 Digits that feature in each column or row do not appear more than once.
 Letters that feature in each column or row do not appear more than once.
 Each letterdigit combination does not appear more than once in the whole grid.
WARNING: Spoilers ahead, so if you want to try and solve this on your own, you might want to stop reading right now and get back after you’ve attempted it.
The obvious and the easy
It is obvious that \(N = 2\) has no solution. \(N = 3\) is rather easy 
$$
\begin{array}{ccc}
A1 & B2 & C3 \\
B3 & C1 & A2 \\
C2 & A3 & B1
\end{array}
$$
With that out of the way, now try \(N = 4\). Really, go and try it out now. If it feels strange that the tactic for \(N = 3\) doesn’t seem to work for \(N = 4\), try alternatives. Yes, there are solutions for \(N = 4\) so your attempts won’t be fruitless.
One thing to notice is that if you have one solution, you can generate other solutions by 
 Exchanging any two columns
 Exchanging any two rows
 Permuting the letters across the whole board. Ex: substitute \(ABC\to BCA\).
 Permuting the digits across the whole board. Ex: substitute \(123\to 231\).
Due to the above, it is easy to see that if a solution exists, you can manipulate it such that the first row will read \(A1\ B2\ C3\ldots\) and the first coumn’s letters are also in the order \(A\ B\ C\ldots\) when read from top to bottom. This is because you can permute the letters in such a way that the first row letters read in order, likewise for the digits and then swap rows in the appropriate order so that the left column letters read in order.
Modulo arithmetic
At first, I tried to produce the required permutation series using modulo arithmetic … using the Julia function 
grid(N,a,b,c,d) = 1 .+ [(mod(a*i+b*j,N), mod(c*i+d*j,N)) for i in 1:N, j in 1:N]
In the above program, mod(m,n)
is the remainder when you divide m
by n
,
even considering negative values for m
. For example, mod(5,3)
is 2
because \(5=3+2\), but mod(5,3)
is 1
since \(5 = 6 + 1\). Note that
we’ll get a matrix whether we range over 1:N
or 0:(N1)
for i
and j
.
If you try that out, you’ll find that grid(N, 1, 1, 1, 2)
produces a solution
for all odd values of \(N\). However, no values of a,b,c,d
work to
solve for even \(N\). You might’ve noticed that when you tried to solve for
\(N = 4\) by hand.
In the above program, I’m representing the entry for each cell as a pair of
numbers (i,j)
where i
and j
are in the range \([1,N]\). It is a trivial
matter to translate that to a matrix of letterdigit combinations with a
function like below 
function printgrid(grid)
n = size(grid)[1]
letters = "ABCDEFGHI"
digits = "123456789"
ld(i,j) = letters[grid[i,j][1]] * digits[grid[i,j][2]]
for i in 1:n
println(join([ld(i,j) for j in 1:n], " "))
end
end
Do we have a solution?
It’s good to have a test function on a grid that we can use to check whether the
grid solves the puzzle or not. To do that, we code a function that produces a
boolean value of true
or false
depending on whether the grid solves the
puzzle or not.
# We assume that the given grid is NxN in structure and
# has tuples with values in the range 0:(N1). We won't test
# for those conditions for simplicity.
function solves(N, grid)
colsokL = all(allunique(grid[i,j][1] for j in 1:N) for i in 1:N)
colsokD = all(allunique(grid[i,j][2] for j in 1:N) for i in 1:N)
rowsokL = all(allunique(grid[j,i][1] for j in 1:N) for i in 1:N)
rowsokD = all(allunique(grid[j,i][2] for j in 1:N) for i in 1:N)
colsokL && colsokD && rowsokL && rowsokD && allunique(grid)
end
We can now use the solves
function to test whether the previous grid
function produces solutions for various values of N,a,b,c,d
.
julia> solves(5, grid(5,1,1,1,2))
true
julia> solves(4, grid(4,1,1,1,2))
false
julia> solves(5, grid(5,1,1,1,1))
false
julia> printgrid(grid(5,1,1,1,2))
C4 D1 E3 A5 B2
D5 E2 A4 B1 C3
E1 A3 B5 C2 D4
A2 B4 C1 D3 E5
B3 C5 D2 E4 A1
Play around with the values for a,b,c,d
to see what works and what doesn’t.
Finding solutions by searching
We saw that the 4x4 case wasn’t solvable using our modulo arithmetic grid
function. So how do we search for solutions in such cases?
First off, it is useful to narrow down possibilities since otherwise we’ll have to look through \(16! = 20922789888000\) cases. We’ll use the “standard form” we came up with earlier for that  i.e. where 
 The first row is assumed to read \([A1,B2,C3,\ldots]\), and
 the first column’s letters have to be in order \([A,B,C,\ldots]\) when read topdown.
We’ll try and fill our the grid from top left corner and move right and down. So we’ll need some way to enumerate the possibilities at any given position assuming that all previous positions are filled in a consistent manner.
# We assume that all the elements up to (i,j) are consistent
# with the constraints.
function satisfies(N, grid, i, j, a, b)
colokL = a ∉ (grid[r,j][1] for r in 1:i1)
colokD = b ∉ (grid[r,j][2] for r in 1:i1)
rowokL = a ∉ (grid[i,c][1] for c in 1:j1)
rowokD = b ∉ (grid[i,c][2] for c in 1:j1)
ix = i*N+j
lduniq = (a,b) ∉ (grid[r,c] for r in 1:i for c in 1:N if r*N+c < i*N+j)
return colokL && colokD && rowokL && rowokD && lduniq
end
function possibilities(N, grid, i, j)
if i == 1
return [(j,j)]
end
[(r,c) for r in 1:N for c in 1:N
if (j > 1  r == i) && satisfies(N, grid, i, j, r, c)]
end
Now we just need to explore the possibilities one by one and produce a solution if one exists.
function explore(N, grid, i, j)
if i > N return true end
for rc in possibilities(N, grid, i, j)
grid[i,j] = rc
jnext = mod(j,N)+1
inext = i + (if (jnext == 1) 1 else 0 end)
if explore(N, grid, inext, jnext)
return true
end
end
return false
end
If the above explore
function returns true
, then the grid
matrix contains
a solution. If it returns false
, then there is no solution. Let’s try them
in a few known situations 
julia> g = [(0,0) for i in 1:3, j in 1:3]
3×3 Matrix{Tuple{Int64, Int64}}:
(0, 0) (0, 0) (0, 0)
(0, 0) (0, 0) (0, 0)
(0, 0) (0, 0) (0, 0)
julia> explore(3, g, 1, 1)
true
julia> printgrid(g)
A1 B2 C3
B3 C1 A2
C2 A3 B1
julia> g = [(0,0) for i in 1:4, j in 1:4]
4×4 Matrix{Tuple{Int64, Int64}}:
(0, 0) (0, 0) (0, 0) (0, 0)
(0, 0) (0, 0) (0, 0) (0, 0)
(0, 0) (0, 0) (0, 0) (0, 0)
(0, 0) (0, 0) (0, 0) (0, 0)
julia> explore(4, g, 1, 1)
true
julia> printgrid(g)
A1 B2 C3 D4
B3 A4 D1 C2
C4 D3 A2 B1
D2 C1 B4 A3
So yeah! It was able to find the solution for the 4x4 grid.
The biggie!
Try to solve the 6x6 grid by hand. Go ahead, really give it a try. When I did it, it felt like this problem was one of those life consuming problems like the Collatz conjecture. I later learnt (from my kid’s math teacher) that Euler struggled with the 6x6 too, so if you’re struggling, don’t be disheartened.
… or if you’re like me, you can run the explore
function on it … just
that you’ll have to wait a bit to get the solution.
julia> g = [(0,0) for i in 1:6, j in 1:6]
6×6 Matrix{Tuple{Int64, Int64}}:
(0, 0) (0, 0) (0, 0) (0, 0) (0, 0) (0, 0)
(0, 0) (0, 0) (0, 0) (0, 0) (0, 0) (0, 0)
(0, 0) (0, 0) (0, 0) (0, 0) (0, 0) (0, 0)
(0, 0) (0, 0) (0, 0) (0, 0) (0, 0) (0, 0)
(0, 0) (0, 0) (0, 0) (0, 0) (0, 0) (0, 0)
(0, 0) (0, 0) (0, 0) (0, 0) (0, 0) (0, 0)
julia> @time explore(6, g, 1, 1)
42.486919 seconds (62.89 M allocations: 5.005 GiB, 0.31% gc time)
false
What? Oh! The 6x6 has no solutions! … and we have a computational proof for that too? Apparently this took Euler ages to conjecture, which is how it started for me initially – i.e. trying to get \(6\times 6\) to work by hand – but then I increasingly suspected there is no solution based on the patterns I was seeing, then decided to write a program … etc.
Note that we weren’t particularly aiming for efficient code, just obvious code,
so I expect we can make this faster as well. Maybe that’s the topic of another
post, but you can get a quick win by just resorting to shortcircuit evaluation
for the satisfies
function 
function satisfies(N, grid, i, j, a, b)
ix = i*N+j
(a ∉ (grid[r,j][1] for r in 1:i1)
&& b ∉ (grid[r,j][2] for r in 1:i1)
&& a ∉ (grid[i,c][1] for c in 1:j1)
&& b ∉ (grid[i,c][2] for c in 1:j1)
&& (a,b) ∉ (grid[r,c] for r in 1:i for c in 1:N if r*N+c < i*N+j))
end
The idea is that if the first condition is false, then the subsequent ones don’t need to be checked, thereby saving time. This version, on my machine runs in < 10 seconds  i.e. gives about a 4x speed up already. A larger thought worth reflecting on is that what would’ve taken Euler perhaps a lifetime to check can now be checked in under 10 seconds! We have a computer proof on our hands … as long as we’re confident enough that our program is correct.
The 4x4 pattern
So now we have an easy solution for odd \(N\), and a weird looking solution for \(4\times 4\) and we also know that there is no solution for \(6\times 6\). Is there a solution for \(8\times 8\)?
When considering this, my thinking took a drastic turn. I tried to extract a pattern out of the 4x4 solution  a single instance  and generalize it for the 8x8. What was odd was that it worked!
Let’s look at the 4x4 solution again 
$$
\begin{array}{cccc}
A1 & B2 & C3 & D4 \\
B3 & A4 & D1 & C2 \\
C4 & D3 & A2 & B1 \\
D2 & C1 & B4 & A3
\end{array}
$$
If we focus on the letter patterns, we see the permutations  \(ABCD\), \(BADC\), \(CDBA\) and \(DCBA\). These looked very familiar to me. They reminded me of … the WalshHadamard transform matrix –
$$
\begin{array}{rrrr}
1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1
\end{array}
$$
… but that matrix didn’t seem to have anything to do with permutations. Ok let’s look at the digit patterns, which read – \(1234\), \(3412\), \(4321\) and \(2143\). This is essentially the same permutation sequence as with the letters, but in a different order. If we number the permutation sequence of the letters as \(0123\), this order is \(0231\).
So, what does the WalshHadamard matrix for an 8x8 transform look like?
$$
\begin{array}{rrrrrrrr}
1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1 & 1 & 1 & 1 & 1
\end{array}
$$
How to make sense of this as a set of permutations? One thing that can help is to note that the rows of the matrix can be generated from the following fundamental patterns 
$$
\begin{array}{lcrrrrrrrr}
P1 & = & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\
P2 & = & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\
P3 & = & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\
P4 & = & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1
\end{array}
$$
Given that we’d take the first row there  \(P1\)  to mean “no permutations”, we can set that aside for the moment.
The rows of the 8x8 matrix can be expressed in terms of these as 
$$
\begin{array}{l}
P1 \\
P4 \\
P3 \\
P3 \times P4\ \text{(elementwise multiplication)} \\
P2 \\
P2 \times P4 \\
P2 \times P3 \\
P2 \times P3 \times P4
\end{array}
$$
Note that the second half is the same as the first half, except with \(P2 \times\) as an added permutation. If we take the four permutations \(P1,P2,P3,P4\) to mean \(ABCDEFGH\), \(EFGHABCD\), \(CDABGHEF\) and \(BADCFEHG\) respectively, we can use the multiplication to mean application of permutations (which in this case are also commutative) to get the following letter sequence –
$$
\begin{array}{lcccccccc}
0) & A & B & C & D & E & F & G & H \\
1) & B & A & D & C & F & E & H & G \\
2) & C & D & A & B & G & H & E & F \\
3) & D & C & B & A & H & G & F & E \\
4) & E & F & G & H & A & B & C & D \\
5) & F & E & H & G & B & A & D & C \\
6) & G & H & E & F & C & D & A & B \\
7) & H & G & F & E & D & C & B & A
\end{array}
$$
The above pattern meets the “no repeated letters in columns and rows” criterion. We need to pick a different ordering of this pattern for the digits. With some fooling around, I tried \(04627315\), which gave me –
$$
\begin{array}{cccccccc}
A1 & B2 & C3 & D4 & E5 & F6 & G7 & H8 \\
B5 & A6 & D7 & C8 & F1 & E2 & H3 & G4 \\
C7 & D8 & A5 & B6 & G3 & H4 & E1 & F2 \\
D3 & C4 & B1 & A2 & H7 & G8 & F5 & E6 \\
E8 & F7 & G6 & H5 & A4 & B3 & C2 & D1 \\
F4 & E3 & H2 & G1 & B8 & A7 & D6 & C5 \\
G2 & H1 & E4 & F3 & C6 & D5 & A8 & B7 \\
H6 & G5 & F8 & E7 & D2 & C1 & B4 & A3
\end{array}
$$
.. which solves the 8x8 case. Now what exactly is the \(04627315\) pattern? If we write those numbers out in binary (3 bits), we get –
$$
\begin{array}{c}
000 \\
100 \\
110 \\
010 \\
\\
111 \\
011 \\
001 \\
101
\end{array}
$$
We see that the first four and the last four in the above sequence are bitwise inverses. Within each set of four as well, if you drop the last bit, you can again divide it into two where the first half and the second half are bitwise flipped – i.e. \(00, 10\) and \(11, 01\) are bitwise flipped. So we have an algorithm to construct the sequence for the digits given the sequence for the letters for \(N = 2^k\) (for which the WalshHadamard matrix is defined). Let’s write it out for \(16\times 16\), for example –
$$
\begin{array}{ccc}
0000 & = & 0 \\
1000 & = & 8 \\
1100 & = & 12 \\
0100 & = & 4 \\
1110 & = & 14 \\
0110 & = & 6 \\
0010 & = & 2 \\
1010 & = & 10 \\
1111 & = & 15 \\
0111 & = & 7 \\
0011 & = & 3 \\
1011 & = & 11 \\
0001 & = & 1 \\
1001 & = & 9 \\
1101 & = & 13 \\
0101 & = & 5
\end{array}
$$
Just to be sure, we need to check whether this indeed solves the \(16\times 16\) case. So let’s work that out –
A1 B2 C3 D4 E5 F6 G7 H8 I9 J10 K11 L12 M13 N14 O15 P16
B9 A10 D11 C12 F13 E14 H15 G16 J1 I2 L3 K4 N5 M6 P7 O8
C13 D14 A15 B16 G7 H8 E9 F10 K5 L6 I7 J8 O1 P2 M3 N4
D5 C6 B7 A8 H1 G2 F3 E4 L13 K14 J15 I16 P7 O8 N9 M10
E15 F16 G13 H14 A9 B10 C11 D12 M7 N8 O5 P6 K3 L4 I1 J2
F7 E8 H5 G6 B3 A4 D1 C2 N15 M16 P13 O14 L9 K10 J11 I12
G3 H4 E1 F2 C7 D8 A5 B6 O11 P12 M9 N10 I15 J16 K13 L14
H11 G12 F9 E10 D15 C16 B13 A14 P3 O4 N1 M2 J7 I8 L5 K6
I16 J15 K14 L13 M10 N9 O12 P11 A8 B7 C6 D5 E4 F3 G2 H1
J8 I7 L6 K5 N4 M3 P2 O1 B16 A15 D14 C13 F12 E11 H10 G9
K4 L3 I2 J1 O8 P7 M6 N5 C12 D11 A10 B9 G16 H15 E14 F13
L12 K11 J9 I10 P16 O15 N14 M13 D4 C3 B2 A1 H8 G7 F6 E5
M2 N1 O4 P3 K8 L7 I6 J5 E10 F9 G12 H11 A14 B13 C16 D15
N10 M9 P12 O11 L14 K13 J16 I15 F2 E1 H4 G3 B6 A5 D8 C7
O14 P13 M16 N15 I12 J11 K10 L9 G6 H5 E8 F7 C2 D1 A4 B3
P6 O5 N8 M7 J2 I1 L4 K3 H13 G13 F16 E15 D12 C11 B10 A9
Well, that looks fine to me, but I admit I may miss repetitions there and so we need to really check it with code … which I’ll set aside for now.^{1} ^{2}
So now, how do we compute that permutation sequence for the digits – i.e. produce the index table above.
function dmap(N,i)
if N == 2 return i end
N2 = div(N,2)
if i < N2
dmap(N2, i) * 2
else
xor(N21, dmap(N2,iN2)) * 2 + 1
end
end
Now we can test whether we get the same sequence for 16 –
julia> dmap.(16,0:15)
16element Vector{Int64}:
0
8
12
4
14
6
2
10
15
7
3
11
1
9
13
5
So there you go. It looks like we have a way to get a solution for the case where \(N = 2^k\).
Proofs, Prooofs, Proooofs!
What we have here is a bunch of conjectures. How do we know that this scheme for \(N = 2^k\) will work?
If we consider how we built up the permutations from the Hadamard matrix, we can rewrite it rather systematically as below (for \(8\times 8\)) –
$$
\begin{array}{lccc}
\text{permutation} & \text{bit pattern} & \text{decimal} & \text{dmap}(8,i) \\
P2^0 \times{} P3^0 \times{} P4^0 & 000 & 0 & 0 \\
P2^0 \times{} P3^0 \times{} P4^1 & 001 & 1 & 4 \\
P2^0 \times{} P3^1 \times{} P4^0 & 010 & 2 & 6 \\
P2^0 \times{} P3^1 \times{} P4^1 & 011 & 3 & 2 \\
P2^1 \times{} P3^0 \times{} P4^0 & 100 & 4 & 7 \\
P2^1 \times{} P3^0 \times{} P4^1 & 101 & 5 & 3 \\
P2^1 \times{} P3^1 \times{} P4^0 & 110 & 6 & 1 \\
P2^1 \times{} P3^1 \times{} P4^1 & 111 & 7 & 5
\end{array}
$$
So the permutation set for the letters is straight forward, and since we know that the permutations \(P2\), \(P3\) and \(P4\) commute and are also orthogonal, we can see that the above set will produce no overlaps for the letters when used to make a matrix.
So, given a sequence of length \(N = 2^k\), how do we produce the “WalshHadamard permutation sequence” for it? To start with, we need to compute the “basis permutations” –
# k ranges over 2^k < N and N = 2^m
function pk(N, k)
p = 2^k
p2 = 2*p
if p2 == N
return vcat(p+1:N, 1:p)
end
[i  mod(i,p2) + pk(p2,k)[mod(i,p2)+1] for i in 0:N1]
end
# N must be a power of 2
function walsh_hadamard_permutations(N)
perms = []
seq = vcat(1:N)
for i in 0:(N1)
s = copy(seq)
(j,n) = (i,0)
while 2^n < N
s = if mod(j,2) == 1 s[pk(N,n)] else s end
(j,n) = (div(j,2), n+1)
end
push!(perms, s)
end
perms
end
Julia hint: Each entry of the result of
walsh_hadamard_permutations
will be an index array that you can use to index a sequence to permute it.
That was a bit circuitous actually, as the permutations can be generated more easily recursively like this, where we’re producing a matrix instead of an array of arrays 
function walsh_hadamard_permutations(N)
if N == 1
return [1]
end
N2 = div(N,2)
p = walsh_hadamard_permutations(N2)
vcat(hcat(p, p .+ N2), hcat(p .+ N2, p))
end
Now, why would using the dmap
approach produce a digit pattern such that we
get no repeats for the letterdigit combinations? I’m not sure I know how to
prove that for the general \(N = 2^k\) case yet. Maybe this is indeed a life
spoiler problem and not knowing the reason will continue to haunt me for a
while. For now, I’m settling for “it looks reasonable and within reach” and am
satisfied with the fact that I found a reordering sequence for \([0,2^k)\)
that I didn’t know about before this, purely by playing around with
patterns.
(21 Feb 2022) Continuing the journey
In ^{1} and ^{2}, I noted after writing the material above that the 16x16 I had posted there was wrong and wasn’t a solution and also conjectured that there are no solutions for \(N = 2^{3k+1}\). Well, I explored the WalshHadamard based technique some more and after one more generalization step was able to produce (and verify) a solution for the 16x16 which I give below without much fanfare (so much for that conjecture!) –
A1 B2 C3 D4 E5 F6 G7 H8 I9 J10 K11 L12 M13 N14 O15 P16
B5 C6 D7 A8 F9 G10 H11 E12 J13 K14 L15 I16 N1 O2 P3 M4
C9 D10 A11 B12 G13 H14 E15 F16 K1 L2 I3 J4 O5 P6 M7 N8
D13 A14 B15 C16 H1 E2 F3 G4 L5 I6 J7 K8 P9 M10 N11 O12
E6 F7 G8 H5 I10 J11 K12 L9 M14 N15 O16 P13 A2 B3 C4 D1
F10 G11 H12 E9 J14 K15 L16 I13 N2 O3 P4 M1 B6 C7 D8 A5
G14 H15 E16 F13 K2 L3 I4 J1 O6 P7 M8 N5 C10 D11 A12 B9
H2 E3 F4 G1 L6 I7 J8 K5 P10 M11 N12 O9 D14 A15 B16 C13
I11 J12 K9 L10 M15 N16 O13 P14 A3 B4 C1 D2 E7 F8 G5 H6
J15 K16 L13 I14 N3 O4 P1 M2 B7 C8 D5 A6 F11 G12 H9 E10
K3 L4 I1 J2 O7 P8 M5 N6 C11 D12 A9 B10 G15 H16 E13 F14
L7 I8 J5 K6 P11 M12 N9 O10 D15 A16 B13 C14 H3 E4 F1 G2
M16 N13 O14 P15 A4 B1 C2 D3 E8 F5 G6 H7 I12 J9 K10 L11
N4 O1 P2 M3 B8 C5 D6 A7 F12 G9 H10 E11 J16 K13 L14 I15
O8 P5 M6 N7 C12 D9 A10 B11 G16 H13 E14 F15 K4 L1 I2 J3
P12 M9 N10 O11 D16 A13 B14 C15 H4 E1 F2 G3 L8 I5 J6 K7
I’ll take a moment there to share the combined thrill of finding a solution and
the disappointment of having made a failed conjecture. Well, apparently it’s
been proved that solutions exist for all N except for \(N = 2\ \text{and}
6\)!! Isn’t that weird? For the moment I’m going to ignore that and pretend
that there is no such proof available .. at least not one that proves by
providing a construction.
So how did that solution come about??
A generalization of the “WalshHadamard permutations”
I wanted to see whether the “WalshHadamard permutations” can be extended to nonpowers of 2. To get there, let’s look at the powerof2 one I’d given earlier 
function walsh_hadamard_permutations(N)
if N == 1
return [1]
end
N2 = div(N,2)
p = walsh_hadamard_permutations(N2)
vcat(hcat(p, p .+ N2), hcat(p .+ N2, p))
end
If we look at what’s happening in the recursive step, we see a division by 2 and a replication of the “N/2” permutations in the last line in a 2x2 pattern. What if we parameterize this “2” and control what is used in every recursive descent step? To get at that, we can change the definition of the parameter \(N\) to be a “vector of bases”, each entry indicating what value to use at each successive recursive descent step. Also, it would be nice to be able to calculate single entries of the resultant permutation matrix, or single rows, without having to calculate the entire matrix. So I remodelled the function to produce a function of two indices (a “row index” and a “column index”) which when called with the appropriate indices yields the matrix’s value.
For non2 values, I used a generalization of the bit swap operation “1,2;2,1” that you get with a conventional WalshHadamard, to be a successive rotation  like “1,2,3;2,3,1;3,1,2” for example, when we produce a 3x3 permutation matrix. So putting all of that together, I generalized the above approach to yield 
# "gwhp = generalized WalshHadamard permutations"
function gwhp(N)
if length(N) == 0
return (i,j) > 0
# Returning 1 in this function will produce a matrix whose entries
# for a 1based index instead of a 0based index.
end
bn = N[end]
N2 = N[1:end1]
bn2 = prod(N2)
wh2 = gwhp(N2) # Calcuating it here instead of inside the function
# will avoid the recursive step when calculating the
# value of each matrix entry.
(r,c) > wh2(mod(r,bn2), mod(c,bn2)) + bn2 * mod(div(c,bn2) + div(r,bn2), bn)
# The first part of the expression produces the same pattern for
# each block of `bn2 x bn2`. The second part varies each block according
# to the rotation pattern. Note that the values of r and c in the second
# part will take on constant values for each `bn2 x bn2` block. Together,
# they achieve the same effect as the `vcat` and `hcat` in the earlier
# presented WalshHadamard permutations function.
end
# We need to take n as a separate parameter since when given
# function representation we can't tell what its intended size is.
function printwhp(whp, n)
for i in 0:n1
for j in 0:n1
print(whp(i,j), " ")
end
println()
end
end
One additional thing to note is that the indices are zerobased and not onebased … since we’re using modulo arithmetic here and zerobased is easier to work with when using modulo arithmetic. The entries of the permutation matrix are also therefore zerobased. Let’s see how to get the 8x8 WalshHadamard permutations matrix 
julia> printwhp(gwhp([2,2,2]), 8)
0 1 2 3 4 5 6 7
1 0 3 2 5 4 7 6
2 3 0 1 6 7 4 5
3 2 1 0 7 6 5 4
4 5 6 7 0 1 2 3
5 4 7 6 1 0 3 2
6 7 4 5 2 3 0 1
7 6 5 4 3 2 1 0
To see what it does for odd values, here are a few examples –
julia> printwhp(gwhp([5]), 5)
0 1 2 3 4
1 2 3 4 0
2 3 4 0 1
3 4 0 1 2
4 0 1 2 3
julia> printwhp(gwhp([2,3]), 6)
0 1 2 3 4 5
1 0 3 2 5 4
2 3 4 5 0 1
3 2 5 4 1 0
4 5 0 1 2 3
5 4 1 0 3 2
julia> printwhp(gwhp([3,2]), 6)
0 1 2 3 4 5
1 2 0 4 5 3
2 0 1 5 3 4
3 4 5 0 1 2
4 5 3 1 2 0
5 3 4 2 0 1
Note: I’ve added extra spaces to make the patterns stand out in the above examples.
printwhp
itself won’t print those spaces. You can try to modify it to do that if you’re interested.
So having generalized the WHPs thus, we also need to generalize some of the
other functions before we can get a letterdigit matrix computed for
nonpowerof2 values of N – dmap
, and a function that puts together gdmap
(the generalization of dmap
) and gwhp
to produce the result matrix.
The generalization of dmap
is interesting. We need to find a generalization
for the “bit flip” operation that is recursively applied. To do this, we treat
the vector \(N\) as the “bases” of a “generalized place value system” \(N =
[b_1,b_2,b_3,\ldots]\) where we express any number \(n\) as a unique series
of “generalized digits” \([d_1,d_2,d_3,\ldots]\) such that \(0 \leq d_k \lt
b_k\) and \(n = \sum_{k}{d_k \prod_{0}^{k1}{b_m}}\) where we take \(b_0 =
1\). If all our \(b_k = 2\), we can see that the “bit flip” operation is a
digitwise increment \(\mod b_k\). So that’s what we’ll take as the
generalization of the “bit flip”. This gives
us the following “generalized dmap” –
"""
brot(N, m, n, f=1)
N is a "generalized number base" vector.
m is the number of rotation steps to applu to each digit
(modulo the corresponding base).
n is the number whose digits are to be "rotated".
If `N` is a power of 2 and $fm != 0 \mod 2$, then this is
effectively just `xor(N1, n)` which flips all the bits
in the binary representation of `n`.
"""
function brot(N, m, n, f=1)
if m == 0 return n end
rn = 0
bs = 1
for bk in N
dk = mod(n, bk)
# rn rebuilds the number n but with the digits
# "rotated" by m steps
rn += bs * mod(dk+f*m,bk)
n = div(n, bk)
bs *= bk
end
rn
end
function gdmap(N, f=1)
if length(N) == 1 return (i > mod(f*i,N[1])) end
N2 = N[1:end1]
bn2 = prod(N2)
g2 = gdmap(N2,f)
function mi(i)
m = div(i, bn2)
brot(N2, m, g2(i  m * bn2)) * N[end] + m
end
end
We’ve also used the same approach of producing a function instead of a vector, so we can compute individual offsets separately without having to compute potentially large vectors and matrices.
Correspondingly, we can define the generalized version of the function that
combines the gwhp
and gdmap
to produce the letterdigit matrix as –
function gmat(N,f=1)
w = gwhp(N)
d = gdmap(N,f)
# The digits are the same pattern as the letters except for a reordering
# of the rows.
(r,c) > (w(r,c), w(d(r),c))
end
Some examples of what gmat
produces –
julia> printgrid(gmat([2,2]),4)
A1 B2 C3 D4
B3 A4 D1 C2
C4 D3 A2 B1
D2 C1 B4 A3
julia> printgrid(gmat([2,3]),6)
A1 B2 C3 D4 E5 F6
B4 A3 D6 C5 F2 E1
C5 D6 E1 F2 A3 B4
D2 C1 F4 E3 B6 A5
E3 F4 A5 B6 C1 D2
F6 E5 B2 A1 D4 C3
julia> printgrid(gmat([3,2]),6)
A1 B2 C3 D4 E5 F6
B3 C1 A2 E6 F4 D5
C5 A6 B4 F2 D3 E1
D6 E4 F5 A3 B1 C2
E2 F3 D1 B5 C6 A4
F4 D5 E6 C1 A2 B3
… where we adapt the previous printgrid
for the new function representation
like this 
function printgrid(m, n)
letters = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
digits = string.(1:26)
ld(i,j) = letters[1+m(i,j)[1]] * digits[1+m(i,j)[2]]
for i in 0:n1
println(join([ld(i,j) for j in 0:n1], " "))
end
end
So the gmat
function produces the following for gmat([4,4])
–
julia> printgrid(gmat([4,4], 1),16)
# .... what I showed at the start of this section ...
julia> printgrid(gmat([4,4]), 16)
A1 B2 C3 D4 E5 F6 G7 H8 I9 J10 K11 L12 M13 N14 O15 P16
B5 C6 D7 A8 F9 G10 H11 E12 J13 K14 L15 I16 N1 O2 P3 M4
C9 D10 A11 B12 G13 H14 E15 F16 K1 L2 I3 J4 O5 P6 M7 N8
D13 A14 B15 C16 H1 E2 F3 G4 L5 I6 J7 K8 P9 M10 N11 O12
E14 F15 G16 H13 I2 J3 K4 L1 M6 N7 O8 P5 A10 B11 C12 D9
F2 G3 H4 E1 J6 K7 L8 I5 N10 O11 P12 M9 B14 C15 D16 A13
G6 H7 E8 F5 K10 L11 I12 J9 O14 P15 M16 N13 C2 D3 A4 B1
H10 E11 F12 G9 L14 I15 J16 K13 P2 M3 N4 O1 D6 A7 B8 C5
I11 J12 K9 L10 M15 N16 O13 P14 A3 B4 C1 D2 E7 F8 G5 H6
J15 K16 L13 I14 N3 O4 P1 M2 B7 C8 D5 A6 F11 G12 H9 E10
K3 L4 I1 J2 O7 P8 M5 N6 C11 D12 A9 B10 G15 H16 E13 F14
L7 I8 J5 K6 P11 M12 N9 O10 D15 A16 B13 C14 H3 E4 F1 G2
M8 N5 O6 P7 A12 B9 C10 D11 E16 F13 G14 H15 I4 J1 K2 L3
N12 O9 P10 M11 B16 C13 D14 A15 F4 G1 H2 E3 J8 K5 L6 I7
O16 P13 M14 N15 C4 D1 A2 B3 G8 H5 E6 F7 K12 L9 I10 J11
P4 M1 N2 O3 D8 A5 B6 C7 H12 E9 F10 G11 L16 I13 J14 K15
The second 16x16 also satisfies the constraints and is a solution.
We can also adapt the solves
function for the new representation to use less
memory and rely more on compute instead, as follows –
function solves(m, N)
ixs = 0:N1
row(r) = [m(r,c) for c in ixs]
(all(allunique(m(i,j)[1] for j in ixs) for i in ixs)
&& all(allunique(m(i,j)[2] for j in ixs) for i in ixs)
&& all(allunique(m(j,i)[1] for j in ixs) for i in ixs)
&& all(allunique(m(j,i)[2] for j in ixs) for i in ixs)
&& all(allunique(vcat(row(r1),row(r2))) for r1 in ixs for r2 in ixs if r2 < r1))
end
Using the above solves
, we can check whether we have a solution –
julia> solves(gmat([4,4]), 16)
true
julia> solves(gmat([2,3]), 6)
false
julia> solves(gmat([2,2,2,2]), 16)
false
julia> solves(gmat([2,2,2]), 8)
true
julia> solves(gmat([3,5]), 15)
true
julia> solves(gmat([5]), 5)
false
The last result is rather disappointing. the matrix produced reads 
A1 B2 C3 D4 E5
B2 C3 D4 E5 A1
C3 D4 E5 A1 B2
D4 E5 A1 B2 C3
E5 A1 B2 C3 D4
Aha  while our brot
corresponds to “bit flip” when the bases are 2,
it doesn’t quite flip anything for other bases. All we have to do to
get a solution is to pass in an f
value of 1
… which produces –
julia> printgrid(gmat([5],1),5)
A1 B2 C3 D4 E5
B5 C1 D2 E3 A4
C4 D5 E1 A2 B3
D3 E4 A5 B1 C2
E2 A3 B4 C5 D1
julia> solves(gmat([5],1),5)
true
So we now have a generalized approach that solves more than just the odd sized matrix.
I’m getting more interested in proofs after this much exploration and explanation … so the journey into patterns shall continue I suppose!
(22022022) Proving the “odd N” case
Let’s start with the simple case  of odd \(N\), i.e. \(N \equiv 1 \mod 2\). We constructed the matrix using the function 
$$ \text{grid}(i,j) = ((i+j\mod N), (i+2j\mod N)) $$
The first and second components generate the whole range \([0,N)\) only for odd \(N\). So for even \(N\), the rows and/or columns will have repeated values due to the stepping by 2. Or to put it in modulo arithmetic terms,
$$
\begin{array}{lcccr}
& 2i & \equiv & 2j & \mod N \\
\implies & i & \equiv & j & \mod N \\
\text{iff} & \text{gcd}(2,N) & = & 1
\end{array}
$$
So what we need to find out is whether \(\text{grid}(i,j) = \text{grid}(m,n)\) implies \(i = m\) and \(j = n\).
The equation gives us –
$$
\begin{array}{llcr}
& i+j & \equiv & m+n\mod N \\
\text{and} & i+2j & \equiv & m+2n\mod N \\
\implies\text{ (by subtraction)} & j & \equiv & n\mod N
\end{array}
$$
And since \(0 \leq j \lt N\) and \(0 \leq n \lt N\), we have $$j = n$$. Similarly, we can show that $$i = m$$. So we now know why \(\text{grid}(i,j)\) produces a solution.
You might’ve noticed that we skipped the row and column uniqueness criteria for letters and digits. Those are trivially ensured by the function since for a given \(i\), the columns all take on different values and likewise for a given \(j\), the rows all take on different values – for both components.
Why can’t we apply that for the “even N” case?
From the previous argument for “odd N”, it looks like we need to find appropriate \(a,b,c,d\) that meet the criteria for even \(N\). For example, if we take \(N = 8\), we have two numbers \(3\) and \(5\) which are both coprime to \(8\). So what do we get if we construct a \(\text{grid}(i,j)\) function based on these two numbers like shown below –
$$ \text{grid}_8(i,j) = (i + 3j\ \mod 8, i + 5j\ \mod 8) $$
With the above function, it is easy to see that we’ll be able to meet all of the row and column uniqueness criteria, since \(\text{gcd}(3,8) = 1\) and \(\text{gcd}(5,8) = 1\).
However, when we set out to show that \(\text{grid}_8(i,j) = \text{grid}_8(m,n) \implies i = m\ \text{and}\ j = n\), we will need to show that from the following equations 
$$
\begin{array}{rclr}
i + 3j & \equiv & m + 3n & \mod 8\\
i + 5j & \equiv & m + 5n & \mod 8
\end{array}
$$
When we subtract the two equations, we get \(2j \equiv 2n\ \mod 8\). That equation does not imply \(j = n\) because \(\text{gcd}(2,8) \neq 1\). This means there exist two different values for \(j\) and \(n\) which can be used to satisfy \(2j \equiv 2n\ \mod 8\), even with the constraints \(0 \leq j \lt 8\) and \(0 \leq n \lt 8\). For example, \(2\times 3 \equiv 2\times 7\ \mod 8\) because \(2\times 7 = 14\) which leaves a remainder \(6\) when divided by \(8\). Similarly, \(2\times 2 \equiv 2\times 6\ \mod 8\) since \(12\) leaves a remainder of \(4\) when divided by \(8\).
Given that \(\text{gcd}(k,N) = 1\) where \(N\) is even, implies that \(k\) is odd (otherwise the GCD will be at least 2), it implies that using two odd coefficients can never yield a uniqueelement matrix because the difference between two odd numbers is always even which will have a GCD of at least 2 with an even number.
So what about using values other than 1 for \(a\) and \(c\)?
Our equations then become 
$$
\begin{array}{rclr}
ai + bj & \equiv & am + bn & \mod N \\
ci + dj & \equiv & cm + dn & \mod N
\end{array}
$$
We know that \(i \equiv j\ \mod N\) implies \(ai \equiv aj\ \mod N\) for all \(a\). So we can multiply the equations by opposite coefficients to get –
$$
\begin{array}{rclr}
cai + cbj & \equiv & cam + cbn & \mod N \\
aci + adj & \equiv & acm + adn & \mod N
\end{array}
$$
If we then subtract the equations, we get –
$$ (adbc)j \equiv (adbc)n\ \mod N $$
For that to imply \(j \equiv n\ \mod N\), we require \(\text{gcd}(adbc,N) = 1\). That implies that both \(ad\) and \(bc\) cannot be odd, since the difference will not be coprime to \(N\). One of them will have to be even. If one of them, say \(ad\) is even, then that implies either \(a\) is even or \(d\) is even. In either case we will fail the column/row uniqueness critertion for even \(N\).
Therefore it is not possible to use the modulo arithmetic approach to construct a solution for even \(N\).
So the initial exploration we did amounts to treating computation as “experimental mathematics” :)

PS: There are errors there that my eye missed (ex: H8 in the first and 3rd rows). When I checked using the
solves
function, it looks like the 16 case isn’t solved by this approach, though the 8 case is. But all is not lost. After all, we know that it doesn’t solve the case of \(2^1\) either, for which there is no solution. When I checked for a few powers of 2 up to \(2^{12}\), I found that the approach doesn’t work for all powers of the form \(2^{4+3k}\) for integer \(k\) and (to my relief) works for other powers of 2, at least up to \(2^{12}\). This gets more and more interesting! ↩︎ 
PS2: Given that \(2^1\) is also of the same pattern and the approach certainly does not work for it and there is no solution for it, I’m at this point conjecturing that there are no solutions for \(N = 2^{3k+1}\). Why those? ↩︎