Patterns

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 letter-digit 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 -

  1. Digits that feature in each column or row do not appear more than once.
  2. Letters that feature in each column or row do not appear more than once.
  3. Each letter-digit 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 -

  1. Exchanging any two columns
  2. Exchanging any two rows
  3. Permuting the letters across the whole board. Ex: substitute \(ABC\to BCA\).
  4. 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:(N-1) 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 letter-digit 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:(N-1). 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 -

  1. The first row is assumed to read \([A1,B2,C3,\ldots]\), and
  2. the first column’s letters have to be in order \([A,B,C,\ldots]\) when read top-down.

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:i-1)
    colokD = b ∉ (grid[r,j][2] for r in 1:i-1)
    rowokL = a ∉ (grid[i,c][1] for c in 1:j-1)
    rowokD = b ∉ (grid[i,c][2] for c in 1:j-1)
    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 short-circuit 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:i-1)        
     && b ∉ (grid[r,j][2] for r in 1:i-1)     
     && a ∉ (grid[i,c][1] for c in 1:j-1)     
     && b ∉ (grid[i,c][2] for c in 1:j-1)     
     && (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 Walsh-Hadamard 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 Walsh-Hadamard 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{(element-wise 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 bit-wise 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 bit-wise 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 Walsh-Hadamard 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(N2-1, dmap(N2,i-N2)) * 2 + 1
    end
end

Now we can test whether we get the same sequence for 16 –

julia> dmap.(16,0:15)
16-element 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 “Walsh-Hadamard 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:N-1]
end

# N must be a power of 2
function walsh_hadamard_permutations(N)
    perms = []
    seq = vcat(1:N)
    for i in 0:(N-1)
        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 letter-digit 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 Walsh-Hadamard 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 “Walsh-Hadamard permutations”

I wanted to see whether the “Walsh-Hadamard permutations” can be extended to non-powers of 2. To get there, let’s look at the power-of-2 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 non-2 values, I used a generalization of the bit swap operation “1,2;2,1” that you get with a conventional Walsh-Hadamard, 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 Walsh-Hadamard 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 1-based index instead of a 0-based index.
    end 
    
    bn = N[end]
    N2 = N[1:end-1]
    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 Walsh-Hadamard 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:n-1
        for j in 0:n-1
            print(whp(i,j), " ")
        end
        println()
    end
end

One additional thing to note is that the indices are zero-based and not one-based … since we’re using modulo arithmetic here and zero-based is easier to work with when using modulo arithmetic. The entries of the permutation matrix are also therefore zero-based. Let’s see how to get the 8x8 Walsh-Hadamard 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 letter-digit matrix computed for non-power-of-2 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}^{k-1}{b_m}}\) where we take \(b_0 = 1\). If all our \(b_k = 2\), we can see that the “bit flip” operation is a digit-wise 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(N-1, 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:end-1]
    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 letter-digit 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:n-1 
        println(join([ld(i,j) for j in 0:n-1], " "))
    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:N-1
    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!

(22-02-2022) 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 unique-element 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 –

$$ (ad-bc)j \equiv (ad-bc)n\ \mod N $$

For that to imply \(j \equiv n\ \mod N\), we require \(\text{gcd}(ad-bc,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” :)



  1. 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! ↩︎

  2. 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? ↩︎