Spiral matrix

Today I found a very interesting set of 'simple' problems. It is a list of Berkeley University's "hardcore tech-interview style riddles and mathematical puzzles". I was specially interested in the Computer Science section (as you might have guessed) and, in particular, one problem caught my attention (probably because of the large matrix).

The problem

The problem reads as follows:

Write a program that will display a "spiral" of NxN numbers, using constant space (no arrays allowed). For example, here's what the spiral looks like for N=10:

 99    98    97    96    95    94    93    92    91    90
 64    63    62    61    60    59    58    57    56    89
 65    36    35    34    33    32    31    30    55    88
 66    37    16    15    14    13    12    29    54    87
 67    38    17     4     3     2    11    28    53    86
 68    39    18     5     0     1    10    27    52    85
 69    40    19     6     7     8     9    26    51    84
 70    41    20    21    22    23    24    25    50    83
 71    42    43    44    45    46    47    48    49    82
 72    73    74    75    76    77    78    79    80    81

Sounds simple enough.

The solution

As I've been reading a lot of "Real World Haskell" lately (and considering the language is great for the task at hand) I decided to code this program in Haskell. If you don't already know Haskell I recommend you to go and learn it on your free time. It is an excellent language with a very unorthodox view of programming (not only functional, but pure!) that will open your mind to a very new and interesting programming style. A good resource for learning Haskell (besides the book, which I highly recommend and will do so with further details in a future post about books) is the classical Gentle introduction to Haskell. But if you prefer a nice colourful website with cool drawings and good info to match the cool design (I know I do!) you can check Learn you a Haskell for great good.

Enough with the language. Into the problem.

To analyse the solution I enumerated some of the matrices, from which matrices of size 4 and 5 were the most useful

15, 14, 13, 12,
4,  3,  2, 11,
5,  0,  1, 10,
6,  7,  8,  9

24, 23, 22, 21, 20,
9,  8,  7,  6, 19,
10,  1,  0,  5, 18,
11,  2,  3,  4, 17,
12, 13, 14, 15, 16

A matrix of size n contains the numbers from n*n-1 to 0. So the next matrix will contain 2n-1 new elements and the rest will be the same elements in a different order. This gives us the idea that the algorithm can be solved recursively.
If we look with more detail into the two matrices outlined above, we can notice that the top row and rightmost column contains all the new elements and that, due to the spiral ordering described in the problem, they will always be located in these positions (there are 2n-1 slots there).

What's left is to figure out how the elements in the previous matrix map to those in the next one and we'll have ourselves a recursive algorithm. We can notice that the first element in the matrix of size 4 is the last element in the sub-matrix formed by removing the new elements from the matrix of size 5. Because of the spiral order the next element to the right in the old matrix becomes the next to the left in the new one and the element that was beneath a number is transposed above it. We can see it as a 90 degree rotation on the old matrix.

The code

Knowing this all that's left is to write the code.
We will need a function that given the position in the n-sized matrix will return the element
that belongs to that position.

spiral :: Int -> Int -> Int -> Int
spiral 0 j n = n*n - 1 - j
spiral i j n
   | j == (n-1) = n*n - n - i
   |otherwise   = spiral (n-1-i) (n-j-2) (n-1)

The first line of this function (after the signature) is one of our base cases: The top row of the matrix is formed by the numbers from (n*n)-1 to n*(n-1). The second definition our second base case: the rightmost column. In this case each element descends from the element in the top right corner (i.e.: n*(n-1) ) one number by row. The rest of the function (the otherwise clause) is our recursive case. On the matrix of size n-1, we are interested in the element that occupies the row n2-i where n2 is the size of the smaller matrix (i.e.: n-1) and the column n2-j-1. Note that the difference between the row and the column is that in the column case we need to substract an extra 1 which refers to the extra columns that appears in the larger matrix. If we don't do this, the transformation would be transposed to the right.

Extra

This code is enough to solve the problem. But to better calculate the results (without having to calculate each index at a time) I provided some extra functions.

-- claculates an index based on the previous matrix
spiral :: Int -> Int -> Int -> Int
spiral 0 j n = n*n - 1 - j
spiral i j n
   | j == (n-1) = n*n - n - i
   |otherwise   = spiral (n-1-i) (n-j-2) (n-1)

-- returns the numbers from 0 to n-1
range :: Int -> [Int]
range n = take n (iterate (1+) 0)

-- returns every possible index of an NxN matrix
matrix :: Int -> [(Int,Int)]
matrix n = [(i,j) | i<-range n, j<-range n]  -- a helper function for the map in printSpiral revSp :: Int -> (Int, Int) -> Int
revSp n (i,j) = spiral i j n

-- returns the result for each index
printSpiral :: Int -> [Int]
printSpiral n = map (revSp n) (matrix n)

-- calculates all spiral matrices from 0 to n
through :: Int -> [[Int]]
through n = map printSpiral (range (n+1))

Does anyone have a better (or different) solution?

Permutations

This last week I happened to have the necessity of generating all the possible permutations of an array. This seems to be a fairly easy task but many programmers won't come up with a simple solution to this problem very fast. There is at least one intuitive, recursive, solution to this problem, many more complex but *slightly* more efficient, iterative, solutions and no really time-efficient solution. After all, the problem is to generate n! orderings of an array so we, at least, have to visit each one once.

The well-known solution

The most widely spread solution to this problem is the Dijkstra algorithm for generating the next permutation. It is based on the assumption that there is a total order between the elements of the array and the elements of the array are unique (there are no repetitions). If these conditions are met, the algorithm always generates the next permutation in lexicographical order.

The algorithm is very straight forward and goes as follow:

def get_next(a):
    N = len(a)
    i = len(a)-1

    while a[i-1] >= a[i]:
        i -= 1

    j = N

    while a[j-1] <= a[i-1]:
        j -= 1

    a[i-1], a[j-1] = a[j-1], a[i-1]

    i += 1
    j = N

    while i < j:
        a[i-1], a[j-1] = a[j-1], a[i-1]
        i+=1
        j-=1

    return a

To use this algorithm to generate every possible permutation we just call it n! times:

def dijkstra(a):
    import operator
    fact = reduce(operator.mul, range(1,len(a)+1))
    the_list = [copy(a)]
    for i in range(fact-1):
        a = get_next(a)
        the_list.append(copy(a))

    return the_list

the only possibly uncommon thing here is line 3 which is a functional-style factorial function.

This works very well for generating the permutations of n (the permutations of the array [0..n-1]) but it won't work for arrays of arbitrary elements. So we explore other options.

The known recursive solution

There's a (not-so-well-)known recursive solution to solving the permutation problem. It is known as the HeapPermute algorithm or Heap's algorithm for permuting an array.
My python implementation of the algorithm goes as follow:

def heap_aux(a, n, the_list):
    if n == 1:
        the_list.append(copy(a))
    else:
        for i in range(n):
            heap_aux(a, n-1, the_list);

            if n % 2 == 1:
                a[0], a[n-1] = a[n-1], a[0]
            else:
                a[i], a[n-1] = a[n-1], a[i]


def heap(a):
    the_list = []
    heap_aux(a, len(a), the_list)
    return the_list

This algorithm works very well and it is very efficient too but it generates the permutations in an order that I find particularly disturbing.
Also it is not as easy to understand as I would like.

Enter my recursive solution
This is probably the most straight forward solution I could come up with. It is a clasical recursive design.
Want to permute an array of n elements? How would you do it if you already had a function to solve the problem for a smaller array?
Simple. Select each possible element of the array to be "the first element of the permutation" and then permute the rest of the array:
A permutation of 1,2,3,4 would be:
1 + permutation_of(2,3,4)
2 + permutation_of(1,3,4)
3 + permutation_of(2,1,4)
4 + permutation_of(2,3,1)

Simple, right?

Let's code it:

  for j in range(i, len(a)):
      a[i], a[j] = a[j], a[i]
      permutation(a, i+1, the_list)
      a[i], a[j] = a[j], a[i]

This loop captures the essence of the procedure described above. We have to remember to put the elements back where they were after the recursive call to allow next iteration to work. Not doing so would result in swapping the wrong element in the next iteration.

Puting it all together (and adding a small optimization):

def perm_aux(a, i, the_list):
    if i == len(a):
        the_list.append(copy(a))
    else:
        perm_aux(a, i+1, the_list)
        for j in range(i+1, len(a)):
            a[i], a[j] = a[j], a[i]
            perm_aux(a, i+1, the_list)
            a[i], a[j] = a[j], a[i]


def permutate(a):
    the_list = []
    perm_aux(a, 0, the_list)
    return the_list

This algorithm generates the list in a nicer order (though it's not lexicographical). And as its predecessor is quite efficient (taking into account the complexity of the problem).

Performance
Normally recursive algorithms are put down for performance reasons: They keep all local variables stored in the stack while the other recursive calls are made and have to incur in function calling overhead. In the case of these two algorithms the good thing is that the recurtion depth is not proportional to the size of the problem solution but to the size of the array (actually it is exactly as deep as the array's length). Also, if the list variables are passed by reference (this is programming language specific), only the reference to the variable needs to be stored in the stack and there is no overhead of coping the variables on every call.
It will probably never be as efficient as their iterative counterparts but for me readability (with small performance penalties) definitely pays off.

To measure the times I used the following python decorator:

def timing(func):
    def inside(*args):
        t1 = time.time()
        res = func(*args)
        t2 = time.time()
        print '%s took %0.3f ms' % (func.func_name, (t2-t1)*1000.0)
        return res
    return inside

Given the complexity of the problem, the graph grows rapidly out of scale. So a graph in logarithmic scale is in order



As you can see the algorithms behave very similarly. The largest difference (for an array of size 9) between the fastest (heap) and the slowest (permute) is 0.6 seconds. This of course will grow with the array size, but in a linear instead of exponential manner.


Have another solution to this problem? comment on! =)