---
jupyter:
  description: 'Squaring the square with planex2d: Duijvestijn''s 21 different squares built from their Bouwkamp code,
    checked to fill the square of side 112 without gap or overlap, and solved as the electrical network of Brooks,
    Smith, Stone and Tutte, whose currents are the sides of the squares.'
---

# Squaring the square

Can a square be cut into smaller squares that are all different? It was once thought impossible, until Brooks, Smith,
Stone and Tutte turned the question into one about electrical networks and found the first such squares. The fewest
squares that can do it is 21. A. J. W. Duijvestijn found that tiling in a computer search, and it is the only one of its
order. Its squares have whole sides and fill a square of side 112.

A tiling like this is written down as its Bouwkamp code, which lists the sides of the squares row by row, each row in
brackets. The example builds Duijvestijn's square from its code, checks that the 21 squares fill it without a gap or an
overlap, and then solves it as the network of Brooks, Smith, Stone and Tutte, in which the currents are the sides of the
squares.

This page is part of a collection of [geometry and design examples](geometry-and-design.md) using planex2d.

<!-- #raw raw_mimetype="text/restructuredtext" -->
.. rst-class:: px-example-links

`Open this example in the playground <https://planex2d.org/playground/>`__, the interactive code editor in the browser, to change the code and run it.
<!-- #endraw -->

![Duijvestijn's square of side 112, cut into 21 different squares](../_static/gallery/squaring-the-square.png){.result}

```python
import numpy as np
import planex as px
```

```python nbsphinx="hidden"
from docs_figures import Row
```

## From the code to the squares

The code reads from the top. Each bracket is a row of squares, placed side by side from the left into the highest free
gap, or into the leftmost one if several are equally high. The squares of a bracket fill that gap's width exactly. The
free gaps are kept in a skyline, which lists the levels reached so far across the width:

```python
code = [(50, 35, 27), (8, 19), (15, 17, 11), (6, 24), (29, 25, 9, 2), (7, 18), (16,), (42,), (4, 37), (33,)]
side = sum(code[0])


def squares_from(code):
    """the squares as (left, top, size), from the top down, with the depth measured downwards"""
    skyline = [[0, side, 0]]  # gaps as [left, width, depth reached]
    placed = []
    for row in code:
        k = min(range(len(skyline)), key=lambda i: (skyline[i][2], skyline[i][0]))
        left, width, depth = skyline[k]
        assert sum(row) == width, f"the row {row} does not fit the gap of width {width}"
        new = []
        for size in row:
            placed.append((left, depth, size))
            new.append([left, size, depth + size])
            left += size
        skyline[k:k + 1] = new
        merged = [skyline[0]]
        for gap in skyline[1:]:
            if gap[2] == merged[-1][2]:
                merged[-1][1] += gap[1]
            else:
                merged.append(gap)
        skyline = merged
    return placed


placed = squares_from(code)
squares = [px.QuickSketch((x, side - y - s)).Rectangle(s, s).Face() for x, y, s in placed]
print(len(squares), "squares, side", side, "| sizes", sorted(s for _, _, s in placed))
```

## The checks

The squares are all different, their areas add up to the square's, and their union is one face of side 112 with that
same area. Had two squares overlapped, the union would be smaller than the sum, and a gap would leave a hole or a
smaller outline (Figure 1):

```python
sizes = [s for _, _, s in placed]
union = px.Unify(squares)
print("all different", len(set(sizes)) == len(sizes), "| sum of areas", sum(s * s for s in sizes), "| side²", side**2)
print("union: faces", len(union.faces), "| area", union.mass, "| box", union.bounding_box)
# px.Plot(px.Compound(squares))
```

```python nbsphinx="hide-input"
Row((px.Compound(squares), "21 squares of side 2 to 50",
     {"labels": [(str(s), (x + s / 2, side - y - s / 2)) for x, y, s in placed if s >= 8], "lw": 1.2}),
    caption="Duijvestijn's square of side 112, with the sides of the squares from 8 up")
```

## The electrical network

Brooks, Smith, Stone and Tutte read a tiling of squares as an electrical network. Each horizontal line in it, made of
the tops and bottoms of squares that touch, is a node at the voltage of its height. Each square is a resistor of one ohm
between the node of its top and the node of its bottom. The voltage across a square is its side, so the current through
it is its side too, and Kirchhoff's laws hold because the squares fit. The square as a whole is a network of one ohm
between its top and its bottom, so 112 volts across it drive 112 amperes.

The check runs the other way. The nodes and resistors are read from the geometry, the network is solved with 112
volts between top and bottom, and the currents come out as the sides of the squares:

```python
def nodes(placed):
    """the horizontal lines of the tiling: their height, their ends and the squares below and above them"""
    pieces = sorted([(y, x, x + s, "below", k) for k, (x, y, s) in enumerate(placed)]
                    + [(y + s, x, x + s, "above", k) for k, (x, y, s) in enumerate(placed)])
    lines = []
    for height, left, right, role, k in pieces:
        if not lines or lines[-1]["height"] != height or left > lines[-1]["right"]:
            lines.append({"height": height, "left": left, "right": right, "below": [], "above": []})
        lines[-1]["right"] = max(lines[-1]["right"], right)
        lines[-1][role].append(k)
    return lines


lines = nodes(placed)
top_of = {k: n for n, line in enumerate(lines) for k in line["below"]}
bottom_of = {k: n for n, line in enumerate(lines) for k in line["above"]}
L = np.zeros((len(lines), len(lines)))
for k in range(len(placed)):
    a, b = top_of[k], bottom_of[k]
    L[a, a] += 1
    L[b, b] += 1
    L[a, b] -= 1
    L[b, a] -= 1
top = next(n for n, line in enumerate(lines) if line["height"] == 0)
bottom = next(n for n, line in enumerate(lines) if line["height"] == side)
inner = [n for n in range(len(lines)) if n not in (top, bottom)]
voltage = np.zeros(len(lines))
voltage[top] = side
voltage[inner] = np.linalg.solve(L[np.ix_(inner, inner)], -L[np.ix_(inner, [top])][:, 0] * side)
currents = [voltage[top_of[k]] - voltage[bottom_of[k]] for k in range(len(placed))]
print(len(lines), "nodes,", len(placed), "resistors")
print("largest difference between a current and its square's side", max(abs(i - s) for i, s in zip(currents, sizes)))
print("current from the top", sum(currents[k] for k in lines[top]["below"]), "A at", side, "V")
middle = [((line["left"] + line["right"]) / 2, side - line["height"]) for line in lines]
wires = [px.FromPolyline([middle[top_of[k]], (x + s / 2, side - y - s / 2), middle[bottom_of[k]]])
         for k, (x, y, s) in enumerate(placed)]
# px.Plot(px.Compound(squares + wires))
```

Drawn over the tiling, the nodes sit in the middle of their horizontal lines, and each resistor runs from the node
above its square through the square's middle to the node below (Figure 2).

```python nbsphinx="hide-input"
Row((px.Compound(squares), "The squares", {"fill": False, "lw": 1.0}),
    (px.Compound(wires), "The network", {"dots": middle, "dashed": squares, "lw": 1.2}),
    caption="The tiling and its electrical network, one resistor per square")
```

See [Many faces](../core/booleans/many-faces.md) for `px.Unify`, and [Square packing](square-packing.md) for squares
packed into the smallest square.

```python nbsphinx="hidden"
import sys
import matplotlib.pyplot as plt
sys.path.insert(0, "..")
from docs_figures import FACE, INK, MAP
from gallery import save_figure

fig, ax = plt.subplots(figsize=(6, 4))
for face, s in zip(squares, sizes):
    shade = MAP(0.15 + 0.7 * (s - min(sizes)) / (max(sizes) - min(sizes)))
    face.Plot(ax=ax, face_col=shade, line_col=INK, legend=False)
for artist in ax.lines + ax.collections + ax.patches:
    artist.set_linewidth(1.4)
for x, y, s in placed:
    ax.text(x + s / 2, side - y - s / 2, str(s), ha="center", va="center", color=INK, fontsize=6 + s / 5)
ax.set_xlim(-28, side + 28)
ax.set_ylim(-1.5, side + 1.5)
ax.set_aspect("equal", adjustable="box")
ax.axis("off")
save_figure(fig, "squaring-the-square")
plt.close(fig)
```