Spiral Matrices with Julia

Julia gives you [almost] C-like speed with Python-like learning curve and easy syntax. In an attempt to solidify my Julia-lang skills I was working through the exercises on Exercism¹ which, by the way, is an amazing resource to practice programming languages. That's where I came across the Spiral Matrix problem.

As you can see, the problem requires you to write a function that takes in an integer n, and outputs a spiral matrix of dimensions n x n (which implies that it's also a square matrix) with n^2 elements sorted from smallest to largest in a [clockwise] spiral order².

I couldn't think of a way to accomplish this task right off the top of my head. But I refused to look up solutions and decided to figure it out entirely on my own. Three days of working on paper, in code, and hunting down repeating patterns which could be expressed in code through mathematical formulae - finally let me obtain the ever-desired spiral matrix. My hard-work paid off.

The purpose of this article is to walk the reader through my entire thought process from the moment I began brainstorming about ways of solving this problem all the way to the moment when I finally acquired the desired result. We'll come across some pattern-hunting, a bit of math, a few code snippets in Julia trying to achieve the individual parts that are finally put together to acquire the solution to the exercise. A basic understanding of linear algebra (vectors and matrices) will be helpful albeit not necessary.

If anything in the Julia syntax looks confusing or odd, make sure to check out the things you should know about Julia at the bottom of this article.

---

How I did it

The basic approach I took was to take an array containing the integers from 1 through n^2 and map them to a two-dimensional array (a n x n matrix) in a spiral manner. Every n in the one-dimensional array can be paired with some tuple (i,j) describing the location it is mapped to in the two dimensional array.³

Thus we obtain a new array of tuples.

Once we've obtained that array of tuples containing the coordinates of another two-dimensional array where each consecutive n in the range of 1 to n^2 is to be mapped to in order to obtain a spiral matrix, the remaining task is pretty straight forward. Here's how it looks (the linear_map variable is set for n=4)

# initiate a spiral matrix of n x n dimensions filled with zeros
spiral_mat = zeros(n, n)
display(spiral_mat)
#=
0 0 0 0
0 0 0 0
0 0 0 0
0 0 0 0
=#

# the array containing the coordinates
linear_map = [(1,1), (1,2), (1,3), (1,4), (2,4), (3,4), (4,4), (4,3), (4,2), (4,1), (3,1), (2,1), (2,2), (2,3), (3,3), (3,2)]

for n in 1:length(linear_map)
  # set n as the value of spiral_mat[i,j]
  spiral_mat[linear_map[n][1], linear_map[n][2]] = n
end

display(spiral_mat)
#=
 1  2  3 4
12 13 14 5
11 16 15 6
10  9  8 7
=#

Note: We're using display() instead of print() or println() because the latter two print the output as inline and thus don't have the pretty matrix-like formatting like display() does.

Okay, that was easy enough, just a couple of lines of code. But how do we obtain the linear_map? That's the fun part.. and also the relatively harder one 🙄.

Let's go find some patterns! 🙌

---

Hunting for Patterns

The first thing that came to my mind was to write down the is and js side by side and see if there's a specific incremental and/or decremental pattern that they followed.

Here are the hand-written linear_maps for n=2, n=3, and n=4 (2x2, 3x3, and 4x4 matrices).

2x2    3x3    4x4
1 1    1 1    1 1
1 2    1 2    1 2
2 2    1 3    1 3
2 1    2 3    1 4
       3 3    2 4
       3 2    3 4
       3 1    4 4
       2 1    4 3
       2 2    4 2
              4 1
              3 1
              2 1
              2 2
              2 3
              3 3
              3 2

Observe that in all three cases,

1. the right column increments while the left column remains at the same value
2. the left column increments while the right column remains at the same value
3. the right column decrements while the left column remains at the same value
4. the left column decrements while the right column remains at the same value

These four events take place in the order I've listed them and keep repeating till the center of the spiral matrix is reached.⁴ I'm going to call this list the actions and each of the elements, an action, through out the rest of the article for the ease of addressing.

We'll discuss shortly about how the algorithm (albeit a pretty abstract one so far) identifies the center.

To make it easier for you to remember, here's the same list in a more compact form. Keeping this in memory is going to prove helpful in a bit.

1. right increments, left stays same
2. left increments, right stays same
3. right decrements, left stays same
4. left decrements, right stays same

Alright. Now, how does the algorithm know when it's reached the center? The short answer is that it doesn't. I've left out a constraint from the algorithm which I'm going to add in now. Before each of the events take place, a check must be made whether the same configuration of i and j has been went over before. Namely, you can imagine adding (i,j) pairs to an array as you go through the steps mentioned above; but before adding the tuple to the array you check whether it's already in the array or not - you carry out the step if it isn't in the array, pass on to the next step if it is already in the array, and stop the algorithm if you find yourself in the situation where you skip to the next step because you found a certain pair, say, (a,b) in the array, but the pair generated from the next step, say, (c,d), is also in the array. That is to say that if you find two consecutively generated pairs to be in the array, that's when you stop.

Now we make another observation - exactly how many times does each of the columns increment or decrement before moving to the next step in the actions? Here's a Julia repl you can play with and generate the i and j columns for any given n to verify the claim I'm about to make.

Here's the claim, or the observation, rather. If you observe enough of those columns, you would notice a general pattern. Here's what it looks like. It might be a little long, so bear with me. We take i=1 and j=1.

i   j
i   j+1
i   j+2
.   .
.   .
.   .
i   n
i+1 n
i+2 n
i+3 n
.   .
.   .
.   .
n   n
n   n-1
n   n-2
n   n-3
.   .
.   .
.   .
n   j
n-1 j
n-2 j
n-3 j
.   .
.   .
.   .
i+1 j
i+1 j+1
i+1 j+2
i+1 j+3
.   .
.   .
.   .
i+1 n-1
i+2 n-1
i+3 n-1
i+4 n-1
.   .
.   .
.   .
n-1 n-1
n-1 n-2
n-1 n-3
n-1 n-4
.   .
.   .
.   .
n-1 j+1
n-2 j+1
n-3 j+1
n-4 j+1
.   .
.   .
.   .
i+2 j+1
i+2 j+2
i+2 j+3
.   .
.   .
.   .
i+2 n-2
i+3 n-2
i+4 n-2
.   .
.   .
.   .
n-2 n-2
n-2 n-3
n-2 n-4
.   .
.   .
.   .
n-2 j+2
n-3 j+2
n-4 j+2
.   .
.   .
.   .
i+3 j+2
i+3 j+3
i+3 j+4
.   .
.   .
.   .
i+3 n-3
i+4 n-3
i+5 n-3
.   .
.   .
.   .
n-3 n-3

Oookay, I suppose that was sufficient to show the pretty vague pattern that I spotted by observing a number of (i,j) pair columns for hours.

We can encode this pattern in a piece of pseudo-code. For a given n,

while j != n, j++;
while i != n, i++;
while j != 1, j--;
while i != 2, i--;
while j != n-1, j++;
while i != n-1, i++;
while j != 2, j--;
while i != 3, i--;
while j != n-2, j++;
while i != n-2, i++;
while j != 3, j--;
while i != 4, i--;
while j != n-3, j++;
while i != n-3, i++;

Note: This piece of pseudo-code, although still not dynamic, does the same thing as the general pattern showed earlier.

At this point, we've got ourselves a new array to work with -

[n, n, 1, 2, n-1, n-1, 2, 3, n-2, n-2, 3, 4, n-3, n-3, ...]

Let's recall what this array actually contains. The right column (the j column) increments till it reaches n, then the i increments till it reaches n, then the j decrements till it reaches 1, then the i decrements till it reaches 2, and so on. In other words, our new array contains the upper-bounds and the lower-bounds up to and down to which i and j increment and decrement based on the action. Naturally enough, I'm going to address this array as the bounds for the rest of the article.

Now let's implement the algorithm in the pseudo-code in julia syntax and find the linear_map for a value of n, say, n=2.

# initialize the linear map array
linear_map = []

# initialize i and j and set them equal to 1
i=j=1

# set a value for n
n = 2

while j != n
    push!(linear_map, (i, j))
    j += 1
end

while i != n
    push!(linear_map, (i, j))
    i += 1
end

while j != 1
    push!(linear_map, (i, j))
    j -= 1
end

push!(linear_map, (i,j))
#= we're pushing (i,j) first and then incrementing/decrementing;
that is why we need to push in the final values
of i and j after the loop is over
=#


display(linear_map)
#=
[(1,1), (1,2), (2,2), (2,1)]
=#

The capability of stopping the algorithm when the center is reached can not be deduced from this array. However, we'll find out soon that we won't need to stop when the center is reached, we can make use of another observation to select a part of this array, the bounds, to work with and neglect the rest.⁵

The next task is to establish a one-to-one correspondence between the set of natural numbers (1, 2, 3, ...) and the bounds. Although it looks like a tricky job at the face of it, we're going to figure out a clever work around to simplify the problem in hand.

---

Order within the Chaos

Let's begin with establishing a one-to-one correspondence between the set of natural numbers and the bounds.

Apologies: Please pardon my messy and lousy illustrations, working with graphics isn't quite my thing.

Now as you can see, there's no easy function that can be used to map each of the elements of the set of natural numbers to the bounds. Hence we're not going to work with the set of natural numbers directly. Instead, we're going to break it up into two subsets - one containing the natural numbers that are mapped to the elements in the bounds that have an n in them, and the other containing the natural numbers being mapped to the constants in the bounds. Then again, we'll make two more subsets from the two subsets of bounds - one containing the odd numbers of its parent set and the other containing the even numbers of the same. Like so...

We have thus acquired four simple arrays/sets - A, B, C, and D - whose general formulae can now be easily deduced unlike earlier. The general formulae are as follows

• A(k) = 4k - 3
• B(k) = 4k - 2
• C(k) = 4k - 1
• D(k) = 4k

Now, by noticing the mapping from the set of natural numbers to the bounds once again, we observe that the general formulae of the sets to which each of A, B, C, and D are being mapped to are

• A'(k) = n - (k-1)
• B'(k) = n - (k-1)
• C'(k) = k
• D'(k) = k + 1

where A', B', C', and D' are the respective ranges of the functions A, B, C, and D.

Now is a good time to recall what we're working towards here. We want to be able to map each of the elements of the set of natural numbers to the bounds. Alright, let's get back to doing that. We took bounds apart, forming A, B, C, and D. This made our task of finding the general formulae easier. Now it's time to put them back together to obtain the original bounds array. But that's not all, we also need to encode the natural numbers that are mapped to each element in bounds. To do that, we'll simply replace each element of bounds with a tuple whose first element will be the natural number mapped to the index of the element and second element will be the original element of bounds. We'll not be doing that in the bounds itself, though. We'll be doing it in the broken down arrays - A, B, C, and D. In fact, that is exactly why we broke down bounds, to make the task just mentioned easier to accomplish. Here's how that's going to look like.

[(1,n)  (2,n)  (3,1)  (4,2)  (5,n-1)  (6,n-1)  (7,2)  (8,3)  (9,n-2)  (10,n-2)  (11,3)  (12,4)  (13,n-3)  (14,n-3)  ... ]

Label: We'll call this one mapped_bounds.

Here's how to do it in code

# broken down arrays with tuples
a = [(4*k - 3, n - (k-1)) for k in 1:n]
b = [(4*k - 2, n - (k-1)) for k in 1:n]
c = [(4*k - 1, k) for k in 1:n]
d = [(4*k, k+1) for k in 1:n]

# put the array containing 'n' together from its odd and even number subsets
ens = cat(a, b, dims=(1, 1))

# do the same for the array of the constants
nums = cat(c, d, dims=(1, 1))

#= Note: We're keeping 'ens' and 'nums' in variables because
later we'll also need them individually =#

# finally, put bounds together by concatenating ens and nums
bounds = sort(cat(ens, nums, dims=(1, 1)))

#= Note: We're sorting bounds by the first elements of the
tuples (by default) which puts the natural numbers in
progressive order =#

Now that we have a dynamic piece of code (with respect to n, that is) to generate the bounds array, we're finally ready to dynamically generate the linear_map and execute the few lines of code we saw toward the beginning of the article, finally acquiring our ever-desired spiral matrix.

---

Putting it all together

Alright, let's get the linear_map. Here's what we gotta do - we gotta loop through bounds, check whether the iterant (which is a tuple) belongs to the ens or the nums, check whether the first element of the iterant is even or odd, then do something with i or j on the basis of the previous checks, and finally add (push) the tuple (i,j) to the linear_map.

Here's how the code looks like

#initialize linear_map variable
linear_map = []

#= notice that the first elements of the tuples are A(k), B(k),
C(k), and D(k), while the second elements are A'(k), B'(k), C'(k),
and D'(k) =#

a = [(4*k - 3, n - (k-1)) for k in 1:n]
b = [(4*k - 2, n - (k-1)) for k in 1:n]
c = [(4*k - 1, k) for k in 1:n]
d = [(4*k, k+1) for k in 1:n]

ens = cat(a, b, dims=(1, 1))
nums = cat(c, d, dims=(1, 1))
bounds = sort(cat(ens, nums, dims=(1, 1)))

i = j = 1

for e in bounds[1:2*n - 1]
    # why the slicing? check note 5 in footer
    if e in ens
        if e[1] % 2 == 1
            while j != e[2]
                push!(linear_map, (i, j))
                j += 1
            end
        elseif e[1] % 2 == 0
            while i != e[2]
                push!(linear_map, (i, j))
                i += 1
            end
        end

    elseif e in nums
        if e[1] % 2 == 1
            while j != e[2]
                push!(linear_map, (i, j))
                j -= 1
            end
        elseif e[1] % 2 == 0
            while i != e[2]
                push!(linear_map, (i, j))
                i -= 1
            end
        end        
    end
end

# don't forget to push the final values of i and j after the loops
push!(linear_map, (i,j))

Now that we have linear_map, the rest is as easy as cake. Remember what we did earlier? That's exactly what we're going to do now.

# create a n x n matrix of zeros
spiral_mat = zeroes(n, n)

for k in 1:length(linear_map)
    spiral_mat[linear_map[k][1], linear_map[k][2]] = k
end

display(spiral_mat)

And there you have it, a spiral matrix of n x n dimensions and n^2 elements. Let's put this whole thing in the form of a function now.

function spiral(n::Int)
    n == 0 && return Matrix{Int}(undef,0,0)

    lin_map = [] 
    ens = sort(cat([(4*k - 3, n - (k-1)) for k ∈ 1:n], [(4*k - 2, n - (k-1)) for k ∈ 1:n], dims=(1, 1)))
    nums = sort(cat([(4*k - 1, k) for k ∈ 1:n], [(4*k, k+1) for k ∈ 1:n], dims=(1, 1)))

    i = j = 1

    for e ∈ sort(cat(ens, nums, dims=(1, 1)))[1:2*n - 1]
        if e ∈ ens
            if e[1] % 2 == 1
                while j != e[2]
                    push!(lin_map, (i, j))
                    j += 1
                end
            elseif e[1] % 2 == 0
                while i != e[2]
                    push!(lin_map, (i, j))
                    i += 1
                end
            end

        elseif e ∈ nums
            if e[1] % 2 == 1
                while j != e[2]
                    push!(lin_map, (i, j))
                    j -= 1
                end
            elseif e[1] % 2 == 0
                while i != e[2]
                    push!(lin_map, (i, j))
                    i -= 1
                end
            end        
        end
    end

    push!(lin_map, (i, j))

    spiral_mat = zeros(Int64, n, n)

    for k ∈ 1:length(lin_map)
        spiral_mat[lin_map[k][1], lin_map[k][2]] = k
    end

    return spiral_mat 
end

End

---

Things you should know about Julia

• Julia, unlike most other languages you might be familiar with, is not 0-indexed.

arr = ["a", "b", "c"]
println(arr[0]) # ERROR: BoundsError: attempt to access 3-element Array{Int64,1} at index [0]
println(arr[1]) # a

• Julia supports plenty of ASCII characters as part of its syntax. The only one used in this article is ∈, which is the same as in.

back to article

---

Closing Remarks

Any mistakes pointed out, corrections recommended, and suggestions made, will be highly appreciated. Thank you for taking the time to read this article 🙏

---

1. Exercism is a brilliant resource if you're looking to learn new languages or just want to solidify your skills in an existing one. It has plenty of exercises ranged from easy to difficult for you to take on, solve, and improve your efficiency in the language that interests you the most. Although it doesn't teach you the highly advanced stuff, it's still a very good resource to widen your grasp on the basics. ↩

2. The counter-clockwise implementation of it can also be figured out using a similar approach to this one. The linear map of n=2, for example, would then be [(1,2),(1,1),(2,1),(2,2)] ↩

3. Note that we don't ever explicitly do the pairing. Instead, since Julia is 1-indexed instead of being 0-indexed, the very indexe of each of the elements in the linear_map is the natural number that the element is paired with. ↩

4. There's no explicit center when n is odd. To generalize, it can be referred to as the innermost element of the spiral. ↩

5. The number of iterations over the elements of bounds needed to reach the center of the matrix is equal to 2n-1. That is, when n=2, you'll reach the center after completing the 3rd step; when n=3, the 5th step; when n=7, the 13th step; and so on. ↩

Comments

[Jaheim Archibald]

An alternative solution using recursion:

def main():
    n = 4
    Matrix = generateEmptyMatrix(n)

    generateMatrix(Matrix, n, 0, 0, n, n)

def generateMatrix(matrix, dimension, startRow, startColumn, endRow, endColumn, counter=10, delta=0):
    if dimension == 1:
        matrix[startRow][startColumn] = counter
        return displayMatrix(matrix)
    if dimension == 2:
        # GENERATE FIRST ROW
        for i in range(startRow, endRow):
            matrix[startRow][i] = counter
            counter += 1
        # GENERATE LAST ROW
        for i in range(2):
            matrix[startRow+1][startRow+1 -i] = counter
            counter += 1
        return displayMatrix(matrix)
    else:
        # GENERATE FIRST ROW
        for i in range(startRow, endRow):
            matrix[startRow][i] = counter
            counter += 1
        # GENERATE LAST COLUMN
        for i in range(startColumn+1, endColumn):
            matrix[i][endColumn - 1] = counter
            counter += 1
        # GENERATE LAST ROW
        for i in range(startRow + 1, endRow):
            matrix[endRow - 1][endColumn - 1 - i + delta] = counter
            counter += 1
        # GENERATE FIRST COLUMN
        for i in range(startRow + 1, endColumn - 1):
            matrix[endRow - 1 - i + delta][startColumn] = counter
            counter += 1
        return generateMatrix(matrix, dimension - 2, startRow + 1, startColumn + 1, endRow - 1, endColumn - 1, counter, delta + 1)

Output when n = 4 (I started my counter at 10 instead of 1):

[10, 11, 12, 13]
[21, 22, 23, 14]
[20, 25, 24, 15]
[19, 18, 17, 16]

[Warisul Imam]

If it's okay with you, may I add it to the article (with credits to you for the section, of course) so that people can have two different ways of solving the problem in the same place?

[Jaheim Archibald]

Yes, it's ok with me

[Warisul Imam]

Could you check your DM, please?