I wondered if there was a way to automate the process of finding a way to tile a tile into a square.

The idea is to represent the tile with a matrix of $0$s for blank space and $1$s for filled spaces like in $\{\{0,1,1\}, \{1,0,0\}\}$.

The tile.

The square (or a general rectangular area) is represented by an array as well.

The function should then have as arguments two lists (representing the tile and the area) and should find the ways (if they exist) to do the tiling.

This is quite an hard problem for me, I have tried playing with Plus, so that if the result matrix has $2$ it means that the tile overlaps while if it is $0$ there is a gap but I can't find a way to tell Mathematica that they should fit the area.

  • 6
    $\begingroup$ Here you have it mathematica-journal.com/issue/v9i3/contents/polyominoes/… $\endgroup$ Jun 14, 2012 at 18:30
  • $\begingroup$ Man thanks so much!I'll read it now but at a first glance seems very useful. Thanks again! $\endgroup$ Jun 14, 2012 at 18:31
  • 1
    $\begingroup$ @belisarius I think this is nice and should be posted as an asnwer $\endgroup$ Jun 14, 2012 at 19:48
  • $\begingroup$ @VitaliyKaurov Ok, done. I usually don't post as an answer something that came up from a Google search ... but that also depends on how long it took me to find the answer. $\endgroup$ Jun 14, 2012 at 21:50

3 Answers 3


In this article the author solves the problem of tiling a rectangle by using pieces taken from a set of polyominoes, which are plane geometric figures formed by joining one or more equal squares edge to edge. For example, these are the pentaminoes, polyominoes formed by joining 5 squares:

enter image description here

Of course this problem is more difficult than the one you asked for, but it is also more interesting and ... there is Mathematica code in the article!

  • 5
    $\begingroup$ +1 Thank you for posting. I think that Mathematica Journal is great and people can find many gems there. I wish more people would read it. $\endgroup$ Jun 14, 2012 at 22:03
  • $\begingroup$ Some time ago i had found that great resource,the Mathematica Journal,but my skills in Mathematica were really poor(and still they are :D)and so i eventually forgot it.Anyway i'm getting better and better and i think i'll start giving a deeper look to it.Thank you all for remembering me this great source! $\endgroup$ Jun 15, 2012 at 8:20
  • $\begingroup$ +1 For letting me know about the existance of the Mathematica Journal. $\endgroup$
    – Red Banana
    Oct 20, 2012 at 3:58
  • 1
    $\begingroup$ @RHall I don't remember testing the code, and probably I didn't because I would posted a warning if it wasn't working. If you really want others looking into it I think a "what is missing in this Mma journal code?" question could attract a lot of attention in this site $\endgroup$ Dec 14, 2013 at 15:31
  • 1
    $\begingroup$ The Mathematica Journal article is no longer available through the link. I was able to get the notebook file by contacting the editor at [email protected]. $\endgroup$ Jul 21, 2020 at 1:22

With some diffidence (because there appears to be a Mathematica bug: see below), I would like to offer an answer in the spirit of the OP's original attempt to solve the problem algebraically.


This problem can be formulated as a binary integer linear program. The reformulation represents the square (or more generally, a rectangle as implemented below) as a vector. Its coefficients represent occupancy: $0$ for empty, $1$ for full. In this representation, every possible tile placement corresponds to a unique vector indicating the cells occupied by the tile.

Let's begin by generating the vectors corresponding to all possible placements of tiles in all possible orientations. It will be convenient to specify as little as possible, beginning with a set of unique tile shapes.

As a running example I use the L, C, T, and Z tiles:

tileSet = {{{1, 1, 1}, {1, 0, 0}}, {{1, 1, 1}, {1, 0, 1}}, {{1, 1, 1}, {0, 1, 0}}, {{0, 1, 1}, {1, 1, 0}}};
ArrayPlot[#, Frame -> False, ImageSize -> 10 Max[Length[First[#]]]] & /@ tileSet

Tile set

(There is no restriction that tiles have a common bounding box, but--to make sure they can be properly positioned--it is important that their matrix descriptions not have any zero padding along outside rows or columns.)

We need to rotate (and, if we wish, flip) all these tiles. An elegant way (perhaps not the simplest) is to generate the group action abstractly in terms of generators (alpha and beta).

apply[s_List, alpha] := Reverse /@ s;       (* Horizontal reflection *)
apply[s_List, beta] := Transpose[s];        (* Diagonal reflection *)
apply[s_List, g_List] := Fold[apply, s, g]; (* Group composition *)
group = FoldList[Append, {}, Riffle[ConstantArray[alpha, 4], beta]];

(Riffle and ConstantArray are just a trick, specific to Coxeter groups, to generate a unique set of group elements. If you don't want to allow reflections, then eliminate the four orientation-reversing group elements, which are the ones using odd numbers of generators.)

Apply the group action to the tiles and eliminate any duplicates:

tiles = Union[Flatten[Outer[apply[#1, #2] &, tileSet, group, 1], 1]];
ArrayPlot[#, Frame -> False, ImageSize -> 10 Max[Length[First[#]]]] & /@ tiles

All tiles

These tiles need to be converted from arrays (suitable for graphical display) into vectors (for input to the linear program). Because these vectors represent all possible positions (translations) of the oriented tiles into the square, here is where we specify the square's dimensions, $m$ rows by $n$ columns. I use $6$ by $9$ for this example (which is a significant size). makeAllTiles creates all possible positions of a single tile within a rectangle; mapping it over the set of all oriented tiles and flattening the result gives what we want. (No need any more to remove duplicates, because there is no possibility of an accidental match.)

makeAllTiles[tile_, m_Integer, n_Integer] := 
  With[{m0 = Length[tile], n0 = Length[First[tile]]},
   Flatten[Table[ArrayPad[tile, {{i, m - m0 - i}, {j, n - n0 - j}}], 
      {i, 0, m - m0}, {j, 0, n - n0}], 1]];
allTiles = Flatten[makeAllTiles[#, 6, 9] & /@ tiles, 1];
ArrayPlot[options = Transpose[Flatten /@ allTiles], ImageSize -> 800]

Tile positions

Each possible position of each possible orientation of each tile shape is represented by a column vector in this array. The rows correspond to cells in the rectangle.

It's now surprisingly straightforward to complete the setup: we seek a subset of columns summing to the vector representing a filled square: all the coefficients in this sum must be $1$'s. The subset can be picked out by an indicator vector $x$. We therefore require its entries to be integral, bounded below by $0$, and bounded above by $1$. A feasible vector, let's say, will be one where the corresponding sum fills the square, but sometimes fills its cells multiple times (due to overlap of tiles): that is, a lower bound $b$ for the result should be a $1$ in each cell of the square. Finally, we minimize the sum of the cell values: if a solution exists, this sum will equal the number of squares. The sum is a linear combination with unit coefficients, represented as a vector $c$.

b = ConstantArray[1, Length[options]]; (* Lower bounds on cell values *)
c = ConstantArray[1, Length[First[options]]]; (* Sum of cell values *)
lu = ConstantArray[{0, 1}, Length[First[options]]]; (* 0-1 constraints *)
x = Quiet[LinearProgramming[c, options, b, lu, Integers, Tolerance -> 0.0005], 
      {LinearProgramming::lpip, LinearProgramming::lpsnf}];
If[!ListQ[x] || Max[options.x] > 1, x = {}];
solution = allTiles[[Select[x Range[Length[x]], # > 0 &]]];
ArrayPlot[z = Sum[i solution[[i]], {i, 1, Length[solution]}], 
    PlotRange -> {0, Max[z]}, Mesh -> True, ColorFunction -> "Rainbow"]


Solution is a list of tiles (in their orientations and positions).


LinearProgramming (apparently) can return only one solution even when many are optimal. To explore all solutions, or a range of solutions, may require different methods. (There are still some tricks we can play: e.g., by permuting the order of columns in Options, we might wind up with a different solution.)

As with the systematic searching methods, when this problem gets big, the computing time grows rapidly. You can see this is so, just by considering the quadratic growth in the list of possible positions of the tiles as the square's dimensions increase. Computing time for the example is around one second (increasing as Tolerance decreases).

Finally, I had to fudge something to get this solution: LinearProgramming, with its default setting for Tolerance, returned a non-optimal solution! (options.x had some 2's in it, indicating some double overlaps.) By starting with Tolerance set to $0.5$ and gradually decreasing it, within a few steps a valid solution was obtained. (With the value set to 10^-1, though, I raised an Assert in the underlying C code :-(. That is a bug, at least in version

  • 1
    $\begingroup$ Wow! I love answers like this. It's humbling and enlightening all at once. $\endgroup$
    – Mr.Wizard
    Jun 16, 2012 at 15:29
  • $\begingroup$ Great answer! Perhaps your dihedral groups mastering may help me in this question on math.se related to this question in this site $\endgroup$ Jun 16, 2012 at 15:35
  • 1
    $\begingroup$ @belisarius Actually, I had been using these methods to study your peg solitaire game (and I have even simplified it using exactly the same LinearProgramming treatment described here). $\endgroup$
    – whuber
    Jun 16, 2012 at 18:43
  • $\begingroup$ I think the problem is that you are trying to minimize the number of pieces used to cover the square when you use c={1,1,1,1...}. If instead you tried to minimize the number of total squares by using c=Total@options it seems to work fine without the Tolerance specification. What do you think? $\endgroup$
    – Rojo
    Aug 8, 2012 at 1:51
  • $\begingroup$ @Rojo That's a good idea. What I was really seeking was any feasible solution. I could not find an option in LinearProgramming to do that and so, as a work-around, was forced to find a feasible solution by optimizing something. $\endgroup$
    – whuber
    Aug 8, 2012 at 13:01

So here's a clarification to the Mathematica Journal post that David Carraher has corrected and updated so it works in Mathematica 9 posted here. Now we all have working code for this cool Mathematica solution!


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.