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 using planex2d.

Open this example in the playground, the interactive code editor in the browser, to change the code and run it.

import numpy as np
import planex as px

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:

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))

21 squares, side 112 | sizes [2, 4, 6, 7, 8, 9, 11, 15, 16, 17, 18, 19, 24, 25, 27, 29, 33, 35, 37, 42, 50]

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):

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))

all different True | sum of areas 12544 | side² 12544
union: faces 1 | area 12544.0 | box ((0.0, 0.0), (112.0, 112.0))

Figure 1: 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:

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))

11 nodes, 21 resistors
largest difference between a current and its square's side 2.842170943040401e-14
current from the top 112.00000000000006 A at 112 V

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).

(a) The squares

(b) The network

Figure 2: The tiling and its electrical network, one resistor per square

See Many faces for px.Unify, and Square packing for squares packed into the smallest square.