I wrote a tangram solver.
It's harder than you think.

Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Morbi non arcu risus quis varius quam quisque id diam.

Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Morbi non arcu risus quis varius quam quisque id diam.

Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Morbi non arcu risus quis varius quam quisque id diam. Fermentum dui faucibus in ornare. Suscipit adipiscing bibendum est ultricies integer quis. Amet est placerat in egestas erat imperdiet sed euismod nisi. Id faucibus nisl tincidunt eget nullam non nisi est sit. Nisl vel pretium lectus quam. Nullam eget felis eget nunc lobortis mattis. Sed viverra tellus in hac habitasse platea dictumst vestibulum. Egestas purus viverra accumsan in nisl nisi scelerisque. Pharetra magna ac placerat vestibulum lectus mauris ultrices eros in.

Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Morbi non arcu risus quis varius quam quisque id diam. Fermentum dui faucibus in ornare. Suscipit adipiscing bibendum est ultricies integer quis. Amet est placerat in egestas erat imperdiet sed euismod nisi. Id faucibus nisl tincidunt eget nullam non nisi est sit. Nisl vel pretium lectus quam. Nullam eget felis eget nunc lobortis mattis. Sed viverra tellus in hac habitasse platea dictumst vestibulum. Egestas purus viverra accumsan in nisl nisi scelerisque. Pharetra magna ac placerat vestibulum lectus mauris ultrices eros in.

1.  The Search for Solutions

1. Show a two piece puzzle. 2. Draw the search tree. 3. Explain recursive DFS idea.

Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Morbi non arcu risus quis varius quam quisque id diam.

Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Morbi non arcu risus quis varius quam quisque id diam.

Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Morbi non arcu risus quis varius quam quisque id diam.

Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Morbi non arcu risus quis varius quam quisque id diam.

Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Morbi non arcu risus quis varius quam quisque id diam.

Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Morbi non arcu risus quis varius quam quisque id diam.

Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Morbi non arcu risus quis varius quam quisque id diam.

2.  Floating-Point Precision Problems

Floating point precision is the #1 reason tangram solving is hard. Under the hood, computers represent non-integer numbers like 1.2 or 5.123 using a representation standardized under the name IEEE-754. IEEE-754 floats are great at representing small numbers like 6.626e-34 or huge numbers like 6.022e+23, but they can't represent every number exactly.

For example, if we try printing the IEEE-754 double precision representation of 1/3, we get 0.3333333333333333148296162... This is a serious problem if you want to answer exact geometric questions like "is this point exactly on that line?".

From our macro perspective, it's easy to say "of course it's on the line!", but your computer will zoom in to the 20th decimal place and say "nope" .

The standard way to resolve these issues is to change the question. Instead of asking the computer "Is this point exactly on that line?", you ask: "Is this point extremely close to that line?" Sometimes, you might have trouble defining what you mean by extremely close, but usually this hack works just fine.

With tangrams though, we really do need exact answers to questions like "Is this point inside that polygon? Or just on its edge?"

We can try asking approximate questions instead like "Is this point within 1e-16 units of that line?". Things will even appear to work. But once you add few pieces to your puzzle, those errors in the nth decimal place start to add up.

Eventually, something will exceed your 1e-16 threshold when it shouldn't. Or by coincidence, something else will come in below your threshold when it shouldn't.

My first attempt at tangram solving ran into these problems immediately, and no amount of fiddling with tolerances could fix these nasty bugs.

3.  Squintegers

Squintegers ( square-root-integers ) are a new number system I invented to solve tangrams using exact math. I noticed that all edge lengths and coordinates in a tangram puzzle can be represented as (a + bs), where a and b are integer values and s is the square root of two.

// Squints are numbers of the form (a+bs), where
// a and b are integers and s is the square root of two.
struct Squint {
  static const int s2 = 2;
  long int a;
  long int b;
};

Equality of Squintegers is much simpler than equality between floats. To compare two squintegers, it's sufficient to compare their a and b components. Two Squints are identical if and only if their a and b components are equal.

$$(a+bs) \stackrel{?}{=} (c+ds) \iff (a \stackrel{?}{=} c) \land (b \stackrel{?}{=} d)$$

Addition and Subtraction of Squintegers is also straightforward. The components of the result are equal to the sum or difference of the components of the original values.

$$(a+bs) + (c+ds) = (a+c) + (b+d)s \\ (a+bs) - (c+ds) = (a-c) + (b-d)s$$

Multiplication is interesting! Multiplying two Squintegers is very similar to multiplying complex numbers. Importantly, the result of multiplying two Squintegers is also a Squinteger.

$$(a+bs)(c+ds) = ac + ads + bcs + bds^2 = (ac+2bd) + (ad+bc)s$$

Division is a problem. What is the solution to (1+0s)/(2+0s)? The result, 1/2 cannot be represented as a Squinteger!

$$\frac{(a+bs)}{(c+ds)} = \; ??$$

Squinteger coordinates solve the exact comparison problems we ran into with floating point coordinates, but they don't allow us to do all of the math we need to solve interesting geometry problems using them.

4.  Sqratios

Sqratios ( squinteger-ratios, pronounced: skray-shee-owz ) are my solution to the division of Squintegers. Instead of representing coordinates as single Squintegers, we'll represent coordinates as fractions composed of Squintegers.

struct Sqratio {
  Squint num;
  Squint den;
};

Similar rational datatypes arise frequently in computational geometry when solving similar (but more important) problems to tangrams. The rules for manipulating Sqratios are the same rules you learned in grade school for manipulating fractions.

Equality

$$\frac{a_{num}}{a_{den}} \stackrel{?}{=} \frac{b_{num}}{b_{den}} \iff a_{num}b_{den} \stackrel{?}{=} a_{den}b_{num}$$

Comparisons

$$\frac{a_{num}}{a_{den}} \stackrel{?}{<} \frac{b_{num}}{b_{den}} \iff a_{num}b_{den} \stackrel{?}{<} a_{den}b_{num}$$

Addition

$$\frac{a_{num}}{a_{den}} + \frac{b_{num}}{b_{den}} = \frac{ a_{num} b_{den} + a_{den} b_{num} }{ a_{den} b_{den} }$$

Subtraction

$$\frac{a_{num}}{a_{den}} - \frac{b_{num}}{b_{den}} = \frac{ a_{num} b_{den} - a_{den} b_{num} }{ a_{den} b_{den} }$$

Multiplication

$$\frac{a_{num}}{a_{den}} \cdot \frac{b_{num}}{b_{den}} = \frac{ a_{num} b_{num} }{ a_{den} b_{den} }$$

Division

$$\frac{a}{b} = a \cdot \frac{1}{b} = \frac{a_{num}}{a_{den}} \cdot \frac{b_{den}}{b_{num}}$$

Simplification

$$\frac{a}{b} = \frac{ \frac{a}{\mathrm{gcd}(a,b)} }{ \frac{b}{\mathrm{gcd}(a,b)} } = \frac{a'}{b'}$$

Simplifying Sqratios is similar to simplifying fractions. Calculate the greatest-common-divisor of the numerator and denominator, and divide it out of the numerator and the denominator. Sqratios are simplified after every addition, subtraction, multiplication, or division to prevent overflowing the underlying Squintegers.

Sign Normalization

$$\frac{-a}{-b} \rightarrow \frac{ a }{ b } \qquad \frac{a}{-b} \rightarrow \frac{ -a }{ b }$$

Sign normalization is an important implementation detail. For cross-multiplication to work correctly, it is necessary to hold the sign of a Sqratio in the numerator. Sign normalization flips the sign of the numerator and denominator if the denominator is negative. Sqratios are normalized after every addition, subtraction, multiplication, or division to prevent incorrect comparisons.

5.  Points and Polygons

Compared to Squintegers and Sqratios, Points and Polygons are mercifully simple.

Points are pairs of (x, y) coordinates, where x and y are Sqratios. Often, I find it helpful to think of Points as 2-D vectors, so I can do things like compute dot products or add two points together.

struct Point {
  Sqratio x, y;
};

Polygons are simply lists of Points with implicit edges between adjacent points. Algorithms operating on these polygons must take care to remember the implicit edge connecting the last point in the polygon to the first.

struct Poly {
  std::vector<Point> verts;
};

In order to distinguish between filled polygons (islands) and unfilled polygons (holes), I use counterclockwise winding order to signify a filled polygon. Almost all polygons that arise during tangram solving are represented by counterclockwise polygons.

6.  Transforming Polygons

Placing a tangram piece in a puzzle requires rotating, translating, and (sometimes) flipping it. Translation is straightforward. I simply add the translation vector to each vertex in the piece and return the result.

Poly translate( const Poly& poly, const Point& t )
{
    Poly newPoly = poly;
    for( auto& v : newPoly ) {
      v += t;
    }
    return newPoly;
}

Rotating pieces is slightly harder because sqratios can't represent arbitrary real numbers. The only valid rotations are increments of 45°. If we tried to rotate a sqratio-valued coordinate by say, 30°, the result would no longer be representable using sqratios.

To make things simple, this routine accepts an integer n, representing the number of times you would like to rotate the piece 45° counterclockwise. I look up the sine and cosine of the given angle and use a 2D rotation matrix to rotate each vertex of the piece.

Poly rotate_ccw_45deg( const Poly& poly, int n )
{
    Sqratio cos[8], sin[8];
    cos[0] = Sqratio(  1, 1 );    sin[0] = Sqratio(  0, 1 );
    cos[1] = Sqratio(  s, 2 );    sin[1] = Sqratio(  s, 2 );
    cos[2] = Sqratio(  0, 1 );    sin[2] = Sqratio(  1, 1 );
    cos[3] = Sqratio( -s, 2 );    sin[3] = Sqratio(  s, 2 );
    cos[4] = Sqratio( -1, 1 );    sin[4] = Sqratio(  0, 1 );
    cos[5] = Sqratio( -s, 2 );    sin[5] = Sqratio( -s, 2 );
    cos[6] = Sqratio(  0, 1 );    sin[6] = Sqratio( -1, 1 );
    cos[7] = Sqratio(  s, 2 );    sin[7] = Sqratio( -s, 2 );

    const Sqratio& c = cos[n%8]; 
    const Sqratio& s = sin[n%8]; 

    Poly newPoly = poly;
    for( auto& v : newPoly ) {
      Sqratio x = c*v.x - s*v.y;
      Sqratio y = s*v.x + c*v.y;
      v.x = x; v.y = y;
    }
    return newPoly;
}

The last transformation we need is the flip. Flipping a piece is the easiest transformation of all—I simply invert the x coordinate of every vertex in the piece.

Poly flip( const Poly& poly )
    {
        Poly newPoly = poly;
        for( auto& v : newPoly ) {
          v.x *= -1;
        }
        return newPoly;
    }

Of the seven pieces in the standard tangram set, only the parallelogram requires flipping. The triangles and square are symmetric, so flipping them would only cause the solver to waste time re-exploring symmetric possibilities. In the solver, every piece includes a boolean to indicate whether or not flipping is required.

7.  Polygon Area

One way to decide if a piece fits in a puzzle is to compare the puzzle's area before you remove a piece to the puzzle's area after you remove that piece. If the area of the puzzle decreases by exactly the area of the piece, then you know that the piece fit entirely inside the puzzle.

A simple way to calculate the area of a 2-D polygon (or an N-D polytope) is to divide the polygon into triangles and sum the signed area of each triangle.

First, pick an arbitrary point. I chose the origin to reduce the amount of math required. Then, create a triangle using each edge of the polygon and the origin. The result will be a set of triangles, some with counterclockwise winding order (green) and some with clockwise winding order (red).

By convention, my polygons are stored in counterclockwise order. Therefore the green triangles contribute positive area, and the red triangles contribute negative area. The sum of the signed triangle areas is equal to the area of the counterclockwise polygon.

This implementation uses the 2D determinant formula to calculate (twice) the signed area of each triangle. These values are summed together, divided by two, and returned.

Sqratio polygon_area( const Poly& poly )
{
    Sqratio a(0);
    const auto n = poly.size();

    for(int i=0; i<n; ++i) {
      const Point& vc = poly[i];           // current vertex
      const Point& vn = poly[(i+1)%n];     // next vertex
      a += ((vc.x * vn.y) - (vc.y * vn.x));
    }
    return a / Squint(2);
}

Factoring out the divide-by-two saves us a few operations. Normally I wouldn't worry about it, but since we are computing polygon areas in the inner loop of the solver, and since division of Sqratios implies expensive simplification and sign fixup, I figured we might as well.

8.  Subtracting Polygons

1. Polygon subtraction is surprisingly hard! Literal edge cases. 2. Here's the basic strategy. 3. Unit testing needed! 4. Creating self-intersecting polygons. Yuck!

9.  Resolving Self-Intersections

1. The fast pointer, slow pointer algorithm. 2. Find first intersection, start walking. 3. Finally! All of the pieces of the solver are working

10.  Solving Tangrams

Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Morbi non arcu risus quis varius quam quisque id diam.

PieceVector solve(puzzle, in_pieces, out_pieces)
{
  // Recursive base case: the puzzle is solved!
  if( out_pieces.empty() and puzzle.area() == 0 ) {
    return in_pieces;
  }

  // Sort pieces by area from largest to smallest.
  out_pieces.sort_by_area();

  // Pick the largest piece to insert next.
  next_piece = out_pieces[0];

  // Keep track of the remaining pieces.
  next_out_pieces = out_pieces[1:];

  // For each polygon in the puzzle...
  for( Poly& puzzle_part : puzzle ) {
  {
    // For each vertex in the polygon...
    for( puzzle_vertex : puzzle_part )
    {
      // For each flip of the piece...
      for( int f=0; f<2; ++f )
      {
        // For each vertex in the piece we are inserting...
        for( piece_vert : next_piece.verts )
        {
          // For each rotation of the piece...
          for( int r=0; r<8; ++r )
          {
            // Flip, rotate, and translate the next_piece to place it
            // in the puzzle in the next orientation with its piece_vert
            // on top of the chosen puzzle_vert.

            Point t = puzzle_vert - piece_vert;
            Piece placed_piece = transform_piece(next_piece, f, r, t);
               
            // Remove placed_piece from all puzzle parts.
            next_puzzle = remove_piece_from_puzzle(puzzle, placed_piece);

            // If removing the piece did not reduce the puzzle area by
            // exactly the area of the piece, it didn't fit in the puzzle.
            if( puzzle.area() != next_puzzle.area() + placed_piece.area() ) {
                continue;
            }

            // Otherwise, we found a place where next_piece fits.
            // Let's remove it from the puzzle and recurse.
            subpuzzle_solution = solve(next_puzzle,
              next_in_pieces, next_out_pieces);

            // If the subpuzzle was unsolvable, forget this
            // piece placement and go on to the next one.
            if( !subpuzzle_solution.valid ) { continue; }

            // Otherwise, add placed_piece to
            // the subpuzzle solution and return it!
            subpuzzle_solution.push_back( placed_piece );
            return subpuzzle_solution;
          }
        }
      }
    }
  } 
}

11.  Results

1. The solver can handle standard tangrams nooo problem. 2. How long does it take to solve? 3. What is the effect of placing small pieces first? 4. How big can we go? Where to from now? Can you solve a giant floor mosaic?

Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Morbi non arcu risus quis varius quam quisque id diam.

This is the result! My algorithm searches for valid piece placements until it solves the puzzle.

Lorem ipsum dolor sit amet, consectetur adipiscing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua. Morbi non arcu risus quis varius quam quisque id diam.

12.  Optimization

The solver works, but can we make it faster? Of course we can! I'd love to hook up a profiler and optimize every instruction, but it's Sunday night, and I'm supposed to be writing a dissertation. Instead, let's pick some low-hanging fruit.

Right now we're using polygon subtraction to decide if a piece fits into the puzzle. Polygon subtraction is fairly involved, and we're doing it in the inner loop of the solver. Let's see if we can cut down the number of polygon subtractions we do.

Here's my idea: instead of performing a full subtraction to decide if a piece fits in the puzzle, let's check if all of the piece's vertices fit in the puzzle first. Since the point-in-polygon test is a lot simpler than the polygon subtraction, this could be a cheap way to disqualify piece placements that don't fit inside the puzzle.

// Check if every vert of the placed piece is inside of the same part of the
// puzzle. If not, go on to the next piece orientation. This skips a ton
// of polygon subtractions and leads to something like a 10x speedup!
{
  bool abort = true;
  for( const auto& part : puzzle.polys ) { 
    bool all_verts_in_this_part = true;
    for( const auto& vtx : placed_piece.poly.verts ) {
      if( point_in_polygon(part, vtx) > 0 ) {
        all_verts_in_this_part = false;
        break;
      }
    }
    if( all_verts_in_this_part ) {
      abort = false;
      break;
    }
  }
  if( abort ) { continue; }
}

This chunk of code tests each vertex of the puzzle against each polygon component of the current puzzle. If all piece vertices fit inside the same puzzle component, we go ahead and perform the subtraction, just like before. If any piece vertex falls outside of the puzzle, we skip this piece placement and move on to the next one.

After including this change, the solver is nearly eleven times faster! Before, the solver required 12.5 seconds to solve the classic square puzzle. After this change, the solver takes only 1.16 seconds, an easy 10.7x win! I suspect we could extract another 100x performance gain if we really needed to, but I'm happy with that for now.

13.  Conclusions

Someday, computers will take over the world, but they won't do it without a lot of hand-holding from engineers like us. With this work, I've advanced the state-of-the-art in computerized tangram solving. I may have hastened our inevitable subjugation by machines, but at least I learned a lot along the way!

In order to circumvent floating-point precision issues, I invented Squintegers, a stupid way for computers to represent tangram coordinates with exact math. No, I can't think of any other use for Squintegers, but if you can, send me an email at jsford94@gmail.com. I'd love to talk to you.

In the end, I'm happy to have rid the world of another useless, time wasting puzzle. Next to Deep Mind and AlphaGo, I'm sure this blog post will be added to the archive of computer game-solving. Yet another point scored for the machines!