Count the ways to tile an MxN rectangle with 1x2 dominos.

Additional material

[ Editor: from a November 1993 thread on picking a random tiling of
  the MxN rectangle with dominos, all tilings equally probable.  The C
  program in the posting is omitted. ]

In article <2d53mt$nf7@news.mantis.co.uk> Tony Lezard
(tony@mantis.co.uk) writes:
>I think the idea of working out the appropriate biases for each
>placing of a domino when constructing a random layout has something
>going for it.

Here ya go: a program that generates a specified tiling of the R by C
rectangle by calculating the permanent of the adjacency matrix at each
step.  But the problem is, it's much slower than the other program I
wrote which uses brute force.

If you type "domrand" without an argument, it counts the number of
solutions of tiling the R by C rectangle.  For the 6x6 case, this
takes 17 seconds on my SPARC 10 (compared with a negligible amount of
time for "domino.c").  If you type "domrand 12", it finds solution
number 12 out of 6728 [Editor: "6278" in the original] (for the 6x6
case).  So, to generate a random tiling of the 6x6 plane where all
tilings are equally probable, you simply choose a random number x from
the uniform distribution in [1,6728] and type "domrand x".  As I say,
it's much faster to use "domino x" instead of "domrand x".

The bottleneck is in calculate_permanent().  Is there a faster way to
evaluate permanents than doing Cramer's rule?  Notice that I store the
adjacency matrix in sparse form.  For the 6x6 rectangle, the matrix
that we have to calculate the permanent of is 18x18, and even though
it's sparse, the C code to evaluate the permanent is still slow.

-- James A. Cherry <jac@doe.carleton.ca>

[ Editor: 6728 = 2 * 58^2, twice a square as exercise 2 above requires
  for M = N = 6, while 6278 = 2 * 3139 is not; the posting itself also
  uses [1,6728].  The answer to Cherry's question is the main entry
  above: the permanent can be turned into a determinant. ]

Here is the best domino-tiler so far.  It converts the permanent
adjacency matrix to the determinant one and finds the determinant.  It
finds the number of different ways to tile a 24x24 rectangle with 2x1
dominos as

7435775679122922965141930394167684120105535297322308613569173527105536

in just over two minutes on my SPARC 10.

I had to represent all integers as strings [...].  The determinant-calculating routine ignores the
sign of the determinant, swapping rows freely without keeping track of
signs.  The matrix is stored as a sparse matrix, and when necessary,
each row is divided by the gcd of all the elements to keep the numbers
small.  Still, in the 24x24 case, the largest number the program has to
deal with is 463 digits long.  Perhaps by improving the algorithm to
find the determinant, this could be reduced.

Again, it's not too hard to modify the program to generate a tiling so
that one tiling is no more likely than another.  [...] So, you can cover
the square in the upper left with a domino horizontally and find the
number of boards that can be generated starting from this position
[ Editor: s1 ] [...].  Then, cover the first square with a domino
vertically and count the number of boards [ Editor: s2 ] [...].  You
will know there are s1 solutions with the
domino placed horizontally and s2 solutions with it placed vertically.
Thus, place the domino horizontally with probability s1/(s1+s2).
Continue this with each domino until the board is filled.

-- James A. Cherry <jac@doe.carleton.ca>, December 1993

[ Editor: the posting included the program, domdet.c (about 700 lines
of C), which is not reproduced here; [...] marks passages about its
internals that are omitted. ]

[ Editor: more from the November 1993 thread "Picking a tiling of the
  M*N rectangle with dominos", in which Tony Lezard asked for an
  algorithm to select one of the tilings at random, all equally
  likely.  The C programs in the postings (domino.c and domrand.c) are
  omitted. ]

In article <2d01n6$srj@salmon.maths.tcd.ie> Nicholas Gray
(asdfghjk@maths.tcd.ie) writes:
> I am interested in this problem because I need to generate a random
> layout for the domino puzzle, where the 28 double-6 dominoes are
> arranged in a 7 by 8 rectangle, the score on each unit square given,
> and the object is to find a possible layout which corresponds to the
> scores.

Tony suggested one strategy in his original article: for the 7x8
rectangle, we can calculate that there are 1,292,697 possible tilings.
I have written a program that will generate all the tilings in under
20 seconds on a SPARC 10, compiling the code using

   gcc -O2 -o domino domino.c

To generate solution number 123,456, for example, type

   domino 123456

at the command line.  The program's output will be:

>Searching for solution 123456...
> 1  1  2  2  3  3  4  4
> 5  6  7  7  8  8  9  9
> 5  6 10 11 12 13 14 14
>15 15 10 11 12 13 16 16
>17 18 19 19 20 20 21 21
>17 18 22 22 23 23 24 24
>25 25 26 26 27 27 28 28
>Total solutions=1292697

-- James A. Cherry <jac@doe.carleton.ca>

I managed to improve the algorithm in "domrand.c".  I was doing it a
silly way before; the new way has indices into the sparse matrix so I
know where each new row begins, and thus I greatly reduce the time to
calculate the permanent.

The new program is still inferior by a wide margin: "domino.c" only
took 0:17 to find all 7x8 solutions by brute force, while this new
program takes 10:36.

-- James A. Cherry <jac@doe.carleton.ca>

tony@mantis.co.uk (Tony Lezard) wrote:
>Perhaps if you have a quick determinant finder (which can easily be
>modified to find permanents, as far as I can tell)
and followed up:
> On second thoughts, perhaps it isn't so easy....

Indeed, Valiant showed that finding the permanent is "#P-complete"
[Garey and Johnson, problem AN13], which implies it's NP-hard.  So
it's very, very unlikely to be that easy.

> It doesn't matter though.  Set elements in the adjacency matrix that
> correspond to _vertical_ adjacency to i instead of 1 and the
> determinant of this matrix is the permanent of the original one,
> times some power of i (which depends on the dimensions of the grid
> and the way you number the squares of the grid).

I don't think this approach works, but I can't be sure.  Have you, for
instance, verified this for fairly large adjacency matrices?  Have you
a proof?  I would you would need to relate permutation parity to the
subset of vertical dominoes, and I doubt there is a relation there.

[ Editor: the main entry above describes exactly this transformation
  (the second "key idea") as one that works for rectangles. ]

I do, however, have an approach that will greatly speed up the
exhaustive-search calculation of the number of tilings of large
regions:

1. Cut the region to be tiled approximately in half with a straight
cut along the edges of squares.  (I advise you to cut the short way if
there is a difference).  The two parts, while not necessarily equal,
will be called the `halves' below.

2. In every domino tiling, each edge of this straight cut will either
bisect a domino, or will lie along the edge of the domino.  Choose
each subset of the edges that will bisect dominoes (you can just about
cut the work in half by choosing only the subsets of the correct
cardinality modulo 2.)

3. For each subset, separate the two halves of the region minus those
dominoes chosen in step 2.  We are going to compute the number of
tilings of each half.

4. If either half has forced dominoes (squares with only one neighbor
left) take them out at this time.  If isolated squares show up in
either half, that half will be untileable, and the number of tilings
of the opposite half can be ignored.

5. Count the number of tilings of each half (using this divide-and-
conquer technique recursively if the halves are large enough, or brute
force if not.  Exactly where the cutoff is depends on your program.)
The number of tilings of the entire region is the sum of the products
of the number of half-tilings over all the subsets of step 2.

Adapting the algorithm to finding the N'th tiling (using an ordering
appropriate to the algorithm) is straightforward.  I should warn James
A. Cherry (and anyone else who has a hack attack over this one) that I
believe you will have to use gmp to make effective use of this
technique.  It gives you the power to count tilings in the range where
longs (and probably double floats) are just not big enough.

-- Dan Hoey <Hoey@AIC.NRL.Navy.Mil>

Counting tilings by computer, and generating random tilings:

I implemented Dan Hoey's divide-and-conquer suggestion in C.  He was
right: the resulting program is much faster than my previous two
efforts.  This program can count the 53,060,477,521,960,000 tilings of
the 12x12 rectangle in just under two minutes on my SPARC 10.  (I
checked that this is the correct number of tilings of the 12x12 using
the formula Tony Lezard gave in his original article.)

Because the program is so quick, it can calculate tilings that are
beyond the range of 32-bit long integers.  As a result, it stores the
number of solutions in a string and performs additions and
multiplications with special routines add() and multiply().

Dan Hoey was right: it makes a _big_ difference if you divide the board
along the longer dimension.  This is because each recursive step takes
time O(2^n) where n is the maximum number of dominos that can cross a
given division.  Clearly, then, it is essential to keep "n" low.

For Tony Lezard, who asked how to generate a particular tiling with
equal probability, the problem of modifying the program to do this is
left as an exercise. :-) There's an easy, inelegant way to do it: for
each domino, place it on the board and call recursive_count() twice,
once placing the domino horizontally and once placing it vertically.
If the first way has X tilings and the second way has Y tilings, then
place the domino horizontally with probability X/(X+Y) or vertically
otherwise.  Continue this with each domino until the board is full.
Note that a board square (i,j) is occupied if the b[i][j] contains a 0,
and it is empty if b[i][j] is non-zero.  As well, some dominos might
only fit one direction, in which case there's no need to count the
tilings.

The more elegant way would be to go into the main loop which permutes
the dominos crossing row "divrow" and find how many tilings there are
with each possible set of crossing dominos.  Suppose that over all
crossings there are Z tilings, and that over one particular crossing,
there are W tilings: then, use that crossing with probability W/Z.
Then, divide and conquer and continue the method.  This is a little
tricky, though, given the recursive nature of the tiling method.

[ Editor: the C program domdiv.c (about 400 lines; at most 12 rows and
  columns and 19-digit counts) that followed is omitted. ]

-- James A. Cherry <jac@doe.carleton.ca>, November 1993

[ Editor: in a follow-up he wrote: ]

>Dan Hoey was right: it makes a _big_ difference if you divide the board
>along the longer dimension.

I meant "shorter dimension" in the second line, not longer.  Sorry.

-- James A. Cherry <jac@doe.carleton.ca>

This thread seems to have been sparked by Noam Elkies' entry in
the rec.puzzles FAQ.  Counting domino tilings of planar regions has
gotten some attention lately from me and some other mathematicians.
So I'm happy to relate the state of the art in this topic, including
a pretty good answer to the original question in this thread:
How do you make a random tiling of an M x N rectangle?

James A Cherry has already described a way to make a random tiling that
would work well, if you had a fast way to count the number of tilings
of an arbitrary region in the plane.  Happily, the
*permanent-determinant* method fits the bill.

As Tony Lezard has already related to the net, you begin by coloring
squares in the plane black and white in alternation, and then you make
a matrix whose rows are the black squares whose columns are the white
squares.  An entry of the matrix is 1 if the row and column containing
it are adjacent squares and 0 if they are not.  The permanent is then
the number of tilings.

By permanent-determinant method, in a planar situation such as this
one, you can change the signs of the entries of the matrix to change
the permanent to a determinant.  This is an amazing trick due to the
mathematical physicist P. W. Kasteleyn.  And to answer the doubts
raised by Dan Hoey, yes it has all been proved, and the proofs aren't
even all that hard once you know what to prove.  The proper way to
alter the signs of the matrix is a little bit tricky in the general
case, but there is a particularly simple rule when the region is a
rectangle or any other region obtained from a rectangle by deleting
dominos.  Namely, a given non-zero entry in the matrix corresponds to
a white square adjacent to a black square.  If the white square is
above the black square, make the entry a -1.  If the white square is
below, to the left, or to the right of the black square, put in a -1.
[ Editor: with -1 in all four cases every entry would change sign and
  the determinant would still not count tilings; presumably the second
  case should read +1, which gives each 2x2 block exactly one -1. ]
As before, if the squares are not adjacent, put in a zero.  The terms
in the determinant of this matrix are all either 1 or -1, therefore the
determinant is, up to sign, the number of tilings.

To answer another question, finding the determinant at the end is much
faster than finding a permanent or any other method suggested here to
count the number of tilings.  The general way to find the determinant
of a matrix is to row reduce it.  The matrix whose determinant we want
here is sparse, which means that there are relatively few non-zero
entries.  One obvious trick is to store only the non-zero entries of
the matrix as you perform the row reduction.  To minimize the number of
non-zero entries during the row-reduction process, recall that the rows
of our matrix are the black squares of the region under consideration.
Order the squares left-to-right, top-to-bottom, and order the rows of
the matrix the same way.  Now row reduce by reducing the first row of
the matrix, then the first two rows, then the first three rows, and so
on.

Finally, row reduction involves division in various places, which is
most conveniently done with rational arithmetic if you are starting
with an integral matrix.  If you don't have rational arithmetic, you
can fake it with integral arithmetic and clearing denominators.

Computer challenge:  If you followed everything I said, you should
be able to count the number of domino tilings of the region

[][]
[]

where each "[]" is a 40 x 40 array of unit squares.

-- Greg Kuperberg <gk00@midway.uchicago.edu>

[ Editor: further postings from the November 1993 thread "Picking a
tiling of the M*N rectangle with dominos", not already in the entry. ]

In article <2d29gi$hks@news.mantis.co.uk> Tony Lezard (tony@mantis.co.uk)
writes:
>Nice program.  Seems to work too (its calculation of the number of
>tilings of the 5x6 rectangle agrees with my calculation (by hand) of
>1183).  Takes more like a minute on the Sparc here to generate all
>solutions.  Nevertheless, it seems a pity that this "brute force"
>solution is the only one anyone has been able to come up with.

I'm sure there are parts of the code that could be improved a bit to
increase the speed.  However, I was thinking of another way to solve
the problem:

Going back to the other person's solution (I've forgotten his name,
sorry), he said he picks a random orientation for the first piece: half
the time it's horizontal, and half the time it's vertical.  So in the
2x3 case, he would start one of the following ways:

aa*  (Solution: aab )  a**  (Solutions: abb abc)
***  (          ccb )  a**  (           acc abc)

But as Tony points out, there is only _one_ solution to the puzzle if
you pick the first way, and there are _two_ solutions to the puzzle if
you pick the second way, as shown above.  Thus, to have the program
pick one of the three possible tilings with equal probability, you
should place the "aa" piece horizontally with probability 1/3 and place
it vertically with probability 2/3.

So, my proposal is the following: given an m by n board, is there a
_quick_ way to work out the number of tilings of a given
non-rectangular field?  So, consider a 3x4 board: I want to put down
the first piece either horizontally or vertically:

aa**  a***
****  a***
****  ****

For each configuration of asterisks, what is the formula for the number
of tilings of those regions?  If I knew it, then I would know how to
weight the probability of the placement of my first piece.

(There are 7 possible tilings of the region of asterisks on the left
and four of the region on the right.  If I knew this, I would then
place the "aa" piece horizontally with probability 7/11.)

Is there some matrix method for calculating the number of tilings of an
irregularly-shaped region with 2x1 dominos?

-- James A. Cherry <jac@kira.doe.carleton.ca>

I think this can be done with an adjacency matrix idea (I described
these in a previous post).  Placing a domino somewhere is like crossing
out a row and a column on one of these things and you end up with a
smaller matrix representing the reduced field.  You can produce an
adjacency matrix for any field shape and the number of solutions will
once again be the permanent of this matrix.  Perhaps if you have a
quick determinant finder (which can easily be modified to find
permanents, as far as I can tell) this may be an efficient decision
procedure.  The problem is you have to run it quite a few times to get
what you want.

How fast can determinants be evaluated anyway?  Does it help to know
that the matrix is sparse (no more than 4 1's (or i's - see the
archive) in any row or column)?

-- Tony Lezard <tony@mantis.co.uk>

On second thoughts, perhaps it isn't so easy to find permanents as they
are not conserved through elementary row operations, which would be
used to quickly find determinants.

It doesn't matter though.  Set elements in the adjacency matrix that
correspond to _vertical_ adjacency to i instead of 1 and the
determinant of this matrix is the permanent of the original one, times
some power of i (which depends on the dimensions of the grid and the
way you number the squares of the grid).

I think the idea of working out the appropriate biases for each placing
of a domino when constructing a random layout has something going for
it.

-- Tony Lezard <tony@mantis.co.uk>

In article <CHF33B.86K@cunews.carleton.ca>
jac@doe.carleton.ca (James A. Cherry) writes:
>OK, OK, here is the best domino-tiler so far.  It converts the permanent
>adjacency matrix to the determinant one and finds the determinant.  It
>finds the number of different ways to tile a 24x24 rectangle with 2x1
>dominos as
>
>7435775679122922965141930394167684120105535297322308613569173527105536
>
>in just over two minutes on my SPARC 10.
>
>I had to represent all integers as strings of length NUM_LEN, and there
>are routines to handle all of them: negate(), abs(), compare(), add(),
>multiply(), divide(), quotient(), remainder(), and gcd(). [...]

>the gcd of all the elements to keep the numbers small.
>Still, in the 24x24 case, the largest number the program has to deal with
>is 463 digits long.  Perhaps by improving the algorithm to find the
>determinant, this could be reduced. [......]

One standard trick around this is to compute the determinant modulo a
bunch of small primes (which lets you divide at will, though you do
have to keep track of the sign of the determinant mod p), then
Chinese-remainder the results to recover the determinant itself.  This
avoids most of the multiprecision arithmetic, and never requires you to
deal with numbers 7 times as long as the final answer (provided you
have a good guess on how large the answer will be to begin with).  The
first 64 primes, from 2 through 311, suffice to distinguish all
positive 125-digit integers.  Or you could use the largest 26 primes
less than 2^16 for the same effect.

You could also implement one of the known algorithms (Wiedemann,
Lanczos,...) for computing determinants of sparse matrices; for a
matrix of order 288 with no more than 4 nonzero entries in each row and
column you're already bound to get a speedup of more than an order of
magnitude.

-- Noam D. Elkies <elkies@zariski.harvard.edu>, December 1993