Friday, February 27, 2015

Fun and Games with CRV: The N-Queens Problem

It's been quite a while since we've solved the zebra puzzle using SystemVerilog. In this post we'll look at another oldie, but a goldie called the n-queens problem. This problem first appeared in a more specific form as the eight queens puzzle, first published in 1848. In this puzzle, the player is asked to place eight queens on a chessboard in such a way that they don't threaten each other. It has fascinated mathematicians (including the great Carl Friedrich Gauss) ever since to find solutions to it. Edsger Dijkstra used it to illustrate the power of backtracking. More recently, Team Specman also solved this problem using constraint programming and published it in this post. The e code they employed is short and sweet, but it's not so easy to digest, though.

In pure "re-inventing the wheel" fashion that we engineers love, I'm going to do my own solution, but in SystemVerilog.

We'll model the chess board as an n by n array of bit values. A 1 will mean that a queen is present on that square, while a 0 will indicate that the square is empty.

class n_queens_solver #(int unsigned n = 8);
  rand bit board[n][n];
  
  // ...
endclass

Modeling the board like this will (hopefully) make it easier to solve the puzzle in a way closer to how a human might solve it on a real board. We'll want to set constraints on the positions of the queens based on the their ranges of motion. Concretely, we'll want to say that only one queen can occupy a row, column or diagonal.

Before we start, however, we'll want to make sure the language supports a couple of things. First, we want to make sure that we can constrain a vector (a one-dimensional array) to contain only a single 1. We can do this using the classical double-for approach:

class singular_on_line #(int unsigned len = 8);
  rand bit line[len];
  
  constraint double_for {
    // the '1' can only be in one place
    foreach (line[i])
      foreach (line[j])
        (i != j) -> ((line[i] == 1) -> (line[j] == 0));
      
    // the '1' must exist in the array
    1 inside { line };
  }
endclass

The two foreach loops say that if one location of the vector contains the 1, then no other location can hold it. This doesn't ensure, however, that the 1 exists inside the vector, hence the need for the inside constraint. This approach works, but it's pretty long and verbose.

Luckily, starting with the 2012 version of the standard we can use array reduction methods in constraints. A good candidate for this is the sum() method. We want the sum of all elements in the array to be 1:

class singular_on_line #(int unsigned len = 8);
  constraint sum {
    line.sum() == 1;
  }
endclass

When testing this constraint, something strange happens. In most cases we won't get a single 1 inside the array. We will, however, always get an odd number of 1s. What's happening here? Because the elements of line are bits, the result is also interpreted as a one bit value, leading to truncation when adding up all the elements. The proper way to do it is to use an explicit cast to force summation to an integer value:

class singular_on_line #(int unsigned len = 8);
  constraint sum {
    line.sum() with ( int'(item) ) == 1;
  }
endclass

It's a shame we can't use the array locator methods inside constraints, as that would have made everything even more expressive. Maybe in a future release of the standard...

The second thing we should look at before diving into the full problem is how to construct the diagonals. This is going to be a bit more funky, since not all diagonals have the same length. For an array of n x n elements we'll have 2*n - 1 diagonals. Since the diagonals are of different lengths, it makes the most sense to store them as dynamic arrays:

class array_of_diags #(int unsigned n = 8);
  rand bit[2:0] array[n][n];
  
  rand bit[2:0] diags[n*2 - 1][];
endclass

We need to establish a convention on how we'll number the diagonals. Let's look at a 3 x 3 array:

+-------+-------+-------+
|       |       |       |
| (0,0) | (0,1) | (0,2) |
|       |       |       |
+-------+-------+-------+
|       |       |       |
| (1,0) | (1,1) | (1,2) |
|       |       |       |
+-------+-------+-------+
|       |       |       |
| (2,0) | (2,1) | (2,2) |
|       |       |       |
+-------+-------+-------+

We'll number the diagonals going horizontally from left to right and vertically from top to bottom. We'll insert elements going also from left to right, but from bottom to top. This is a bit difficult to explain in words (and I'm not good enough yet with HTML to draw it for you), but it should be easy to understand by example. In our case, the diagonals will contain the following elements:

  • diags[0] = { (0,0) }
  • diags[1] = { (1,0), (0,1) }
  • diags[2] = { (2,0), (1,1), (0,2) }
  • diags[3] = { (2,1), (1,2) }
  • diags[4] = { (2,2) }

With a little observation and mathematical induction, we can figure out that the constraint to create the diagonals looks like this:

class array_of_diags #(int unsigned n = 8);
  constraint create_diags {
    foreach (diags[i,j])
      if (i < n)
        diags[i][j] == array[i - j][j];
      else
        diags[i][j] == array[(n - 1) - j][i + j - (n - 1)];
  }
endclass

What's very important, though, when working with constraints on dynamic arrays is to make sure that their sizes are set up correctly. We can easily do this inside the pre_randomize() function:

class array_of_diags #(int unsigned n = 8);
  function void pre_randomize();
    foreach (diags[i])
      diags[i] = new[get_len_of_diag(i)];
  endfunction
  
  
  function int unsigned get_len_of_diag(int unsigned idx);
    if (idx < n)
      return idx + 1;
    
    return 2*n - (idx + 1);
  endfunction
endclass

The formula to get the length of a diagonal can be easily worked out by observation. In fact, that's how I figured out most of these array constraints. I just took an example array and tried to relate different things to an element's row/column index and the size of the array.

Just to be on the safe side, I've also made sure that any constraints we apply on the diagonals' elements will also propagate back to the initial array we've constructed them from. We don't want to get any issues with unidirectional constraints. For brevity, I won't show that code here.

With these topics sorted out we can get started. We'll use some helper variables to store our rows, columns and diagonals:

class n_queens_solver #(int unsigned n = 8);
  rand bit rows[n][n];
  rand bit cols[n][n];
  rand bit main_diags[2*n - 1][];
  rand bit anti_diags[2*n - 1][];
endclass

Unfortunately, SystemVerilog doesn't have the concept of pointers, so these extra arrays will use some extra memory, but this shouldn't decrease generation performance. Don't hold me to this last statement, but this is a hunch of mine, since the constraints we use to relate these variables to the board don't increase our randomization state space (they are 1:1 mappings).

It's not really necessary to use a separate variable for the rows (as rows are easily indexable from an array), but it makes the code more expressive (plus, I'm a bit of a sucker for symmetry and it would feel unbalanced to treat the rows differently). Here's the constraint to create the rows:

class n_queens_solver #(int unsigned n = 8);
  constraint create_rows {
    foreach (rows[i, j])
      rows[i][j] == board[i][j];
  }
endclass

And here's how to create the columns:

class n_queens_solver #(int unsigned n = 8);
  constraint create_cols {
    foreach (cols[i, j])
      cols[i][j] == board[j][i];
  }
endclass

The main diagonals are the diagonals that sweep the array from left to right and from top to bottom:

main_diags

Here's the constraint to create them:

class n_queens_solver #(int unsigned n = 8);
  constraint create_main_diags {
    foreach (main_diags[i,j])
      if (i < n)
        main_diags[i][j] == board[j][(n - 1) - i + j];
      else
        main_diags[i][j] == board[i - (n - 1) + j][j];
  }
endclass

The anti diagonals sweep the array from right to left and from top to bottom:

anti_diags

And here's the constraint to create them too:

class n_queens_solver #(int unsigned n = 8);
  constraint create_anti_diags {
    foreach (anti_diags[i,j])
      if (i < n)
        anti_diags[i][j] == board[i - j][j];
      else
        anti_diags[i][j] == board[(n - 1) - j][i + j - (n - 1)];
  }
endclass

I would love to see a concept where these fields can get created without occupying extra memory, since they're basically just copies of other fields. Until then, we'll have to make the most of what we have.

After creating the lines of our chess board, it's time to finally start writing the constraints to solve the puzzle. Since a queen has an unlimited horizontal range of motion, it must occupy an entire row by itself:

class n_queens_solver #(int unsigned n = 8);
  constraint singular_on_row {
    foreach (rows[i])
      rows[i].sum() with ( int'(item) ) == 1;
  }
endclass

Similarly, a queen must occupy an entire column by itself:

class n_queens_solver #(int unsigned n = 8);
  constraint singular_on_col {
    foreach (cols[i])
      cols[i].sum() with ( int'(item) ) == 1;
  }
endclass

Queens also have unlimited range on any diagonals (both main and anti) they occupy. This means that no two queens can occupy the same diagonal. My first idea was to put a constraint similar to the ones above on each diagonal:

class n_queens_solver #(int unsigned n = 8);
  constraint singular_on_main_diag {
    foreach (main_diags[i])
      main_diags[i].sum() with ( int'(item) ) == 1;
  }
endclass

After adding this constraint the constraint solver started failing. I was scratching my head in wonder and was just about ready to cry foul on the solver, when I realized that I had been wrong all along. It's true that there can only be one queen on a certain diagonal, but the key thing to note is that there doesn't have to be a queen on each diagonal. Of course, this makes sense, since as we saw above there are 2*n - 1 diagonals and we can only place n queens.

What we can do instead is only constrain those diagonals that cross a position where a queen is located:

class n_queens_solver #(int unsigned n = 8);  
  constraint singular_on_main_diag {
    foreach (cols[j,i])
      if (cols[j][i] == 1)
        main_diags[(n - 1) - j + i].sum() with ( int'(item) ) == 1;
  }
endclass

Adding a similar constraint to the anti diagonals will solve the puzzle. The constraint doesn't look very nice, though. It's looks pretty complicated and it won't be easy to understand (even I'm having problems it while writing this post). It seems to me that we're doing too much work for the solver!

As we saw above, we can either have one queen or no queens on a diagonal. Why not write that as a constraint? This leads to much less code:

class n_queens_solver #(int unsigned n = 8);
  constraint singular_on_main_diag {
    foreach (main_diags[i])
      main_diags[i].sum() with ( int'(item) ) inside {0, 1};
  }
  
  constraint singular_on_anti_diag {
    foreach (anti_diags[i])
      anti_diags[i].sum() with ( int'(item) ) inside {0, 1};
  }
endclass

And that's all there is to it! We've written a lot of code just to set up our board's lines, but the constraints we wrote on them are pretty straightforward and easy to understand.

I've tried running the solution on my modest machine and it works pretty fast for an n equal to 8. It's a bit slower for 9 and even slower for 10. It kind of runs out of steam for greater values, as Vitaly predicted in his post. Intelligen might be faster for this particular problem (I haven't tried it out myself), but whether it's faster in all practical situations for usual constraints I can't say. Cadence might be stretching the truth here, or rather put its best foot forward, for marketing purposes. I mean, how often do you write this style of constraints in production code?

You can find the full code on SourceForge. Feel free to download it and see how large an n you can handle.

Also, while I was writing the post I had the idea that instead of using vectors for the rows, columns and diagonals, maybe we could use packed arrays. We might be able to set $onehot constraints on them, though I'm not sure if this is supported. If any of you try it out, please share your experience in the comments section.

4 comments:

  1. Very prettily done - thanks.

    I've used $onehot constraints but I'm not 100% sure that all the main tools support that.

    Did you, at any point in your experiments, try using an auxiliary variable (one for each column, and possibly another for each row) that contains the location of the queen? That might make some of the constraints easier to write, but perhaps would give the solver a harder time.

    When doing "interesting" things with SV constraints - whether for intellectual pastime like this, or for the day job - I've often found myself hitting incomprehensible performance problems, simply because I have no idea how the tools' solvers are doing their job. Each tool seems to have its own "no go areas", constraint idioms that cause solver performance to drop horribly - but it's really hard to get any kind of grip on why that happens, and what you can do to avoid it. That's a slightly simpler problem for Specman users because they have only one implementation to worry about (two, I suppose, if you count both Intelligen and pgen?) and the Specman documentation says quite a lot about how things work behind the curtain. Constraint solver performance in SV, and how to avoid hobbling it, would be a good topic for a future post or series of posts, if anyone has any expertise in that area.

    ReplyDelete
    Replies
    1. I didn't want to use the method you describe (with the auxiliary variable per row/column) because that's how Vitaly did it in his Specman post and I didn't want my post to just be a simple translation from e to SV. I also wanted to experiment with the array reduction method in constraints (as that's something new in the 2012 standard).

      With respect to solver performance, I've tried it out in 2 simulators and I have noticed considerably worse performance on one vs. the other. It's as you say, pretty much impossible to know why. I usually try to stay away from such topics as simulator performance as I don't want to get slapped with a lawsuit. Even ignoring that, I would have no idea where to start, since what applies to one simulator won't apply to another.

      Delete
  2. That's a pity. Only 1 of 3 simulators I tried supports $onehot in constraints. Using an auxiliary "which bit should be set" variable works fine, and does not seem to cause any trouble when combined with other constraints on the vector:

    rand bit [N-1:0] v;
    rand int position;
    constraint v_onequeen {
    position inside {[0:N-1]};
    v == (1 << position);
    }

    ReplyDelete
  3. class n_queen_puzzle#(int n=8);
    rand bit board[n][n];
    rand bit board_tx[n][n];

    constraint unique_row_c {
    foreach(board[i,])
    board[i].sum() with (int'(item)) == 1;
    }

    constraint board_tx_c {
    foreach (board[i,j])
    board_tx[j][i] == board[i][j];

    solve board before board_tx;
    }

    constraint unique_col_c {
    foreach(board_tx[i,])
    board_tx[i].sum() with (int'(item)) == 1;


    }

    constraint board_c2 {
    foreach(board[i, j]) {
    if(board[i][j]) {
    foreach(board[k, l]) {
    if(i != k && k+l == i+j) board[k][l] != 1;
    if(i != k && l-k == j-i) board[k][l] != 1;
    }
    }
    }
    }
    endclass

    ReplyDelete