Circle packing in a square

How large can \(n\) equal circles be when they must fit into a unit square without overlapping? The question is old and hard. The best packings are proven optimal for every \(n\) up to 33, and beyond that the best known ones are collected, as on Eckard Specht’s Packomania. For two circles, the answer lies on the diagonal, \(r = 1/(2 + \sqrt{2})\), four sit in a grid of \(r = 1/4\), and five put the fifth in the middle, \(r = (\sqrt{2} - 1)/2\).

Here, SciPy searches for the packings, and planex2d checks them with exact discs: the discs must lie inside the square, and their union must have the area of \(n\) separate discs, so none of them overlaps another. Their radii are compared with the proven ones.

The best packings of 5 to 10 equal circles in a square

[1]:

import math

import numpy as np
import planex as px
from scipy.optimize import minimize

The search

Packing circles of radius \(r\) in the unit square is the same as spreading \(n\) points in a square of side \(1 - 2r\) as far apart as possible: points at least \(d\) apart in the unit square give circles of radius \(r = d / (2(1 + d))\) once the square is scaled down by \(1 - 2r\). spread maximizes the smallest distance \(d\) with SciPy’s SLSQP from 100 random starts, since the problem has many local optima, and returns the circles’ radius and centers:

[3]:

def spread(n, starts=100, seed=1):
    rng = np.random.default_rng(seed)
    apart = [{"type": "ineq", "fun": lambda x, i=i, j=j: np.sum((x[2 * i:2 * i + 2] - x[2 * j:2 * j + 2]) ** 2) - x[-1] ** 2}
             for i in range(n) for j in range(i + 1, n)]
    best = None
    for _ in range(starts):
        found = minimize(lambda x: -x[-1], np.r_[rng.random(2 * n), 0.1], method="SLSQP", constraints=apart,
                         bounds=[(0, 1)] * (2 * n) + [(0, 2)], options={"ftol": 1e-15, "maxiter": 500})
        if found.success and (best is None or found.x[-1] > best.x[-1]):
            best = found
    d = best.x[-1]
    r = d / (2 * (1 + d))
    return r, r + best.x[:-1].reshape(n, 2) * (1 - 2 * r)

The checks

For each \(n\), the discs are drawn exactly, a hair smaller than the radius found so that circles that touch stay apart. Nothing of them may lie outside the square, and their union must have the area \(n\pi r^2\). The radii agree with the proven optimal ones, as Packomania lists them, to the twelve digits printed:

[4]:

proven = {2: 0.292893218813, 3: 0.254333095030, 4: 0.25, 5: 0.207106781187, 6: 0.187680601147,
          7: 0.174457630187, 8: 0.170540688701, 9: 0.166666666667, 10: 0.148204322565}
square = px.QuickSketch().Rectangle(1, 1).Face()
packings = {}
for n, best in proven.items():
    r, centers = spread(n)
    discs = px.Unify([px.Circle(tuple(c), r * (1 - 1e-12)).Face() for c in centers])
    packings[n] = discs
    print("n", n, "| radius", round(r, 12), "| proven", best, "| outside", (discs - square).mass,
          "| union against n π r²", discs.mass / (n * math.pi * (r * (1 - 1e-12)) ** 2) - 1)

n 2 | radius 0.292893218813 | proven 0.292893218813 | outside 0.0 | union against n π r² 2.220446049250313e-16
n 3 | radius 0.25433309503 | proven 0.25433309503 | outside 0.0 | union against n π r² -2.220446049250313e-16
n 4 | radius 0.25 | proven 0.25 | outside 0.0 | union against n π r² 0.0
n 5 | radius 0.207106781187 | proven 0.207106781187 | outside 0.0 | union against n π r² 0.0
n 6 | radius 0.187680601147 | proven 0.187680601147 | outside 0.0 | union against n π r² -3.3306690738754696e-16
n 7 | radius 0.174457630187 | proven 0.174457630187 | outside 0.0 | union against n π r² 0.0
n 8 | radius 0.170540688701 | proven 0.170540688701 | outside 0.0 | union against n π r² 0.0
n 9 | radius 0.166666666667 | proven 0.166666666667 | outside 0.0 | union against n π r² 4.440892098500626e-16
n 10 | radius 0.148204322565 | proven 0.148204322565 | outside 0.0 | union against n π r² -4.440892098500626e-16

5 circles — Figure 1: The best packings of 5, 7 and 10 circles

7 circles — Figure 1: The best packings of 5, 7 and 10 circles