Problem 222: Sphere Packing

Detailed mathematical approach

Problem Summary

A pipe has an internal radius of \(50\) mm, and \(21\) balls, one of each radius \(30,31,\dots,50\) mm, must fit inside it. We must find the length of the shortest pipe that can fully contain all of them, in micrometres and rounded to the nearest integer; since \(1\) mm is \(1000\) micrometres, this is \(1000L\) rounded, where \(L\) is the shortest length in millimetres. ^[1] The statement gives no smaller example, so as a test case we use the one built into the C++ program: a pipe of radius \(8\) with balls of radii \(5\), \(6\), \(7\) and \(8\).

The whole solution in one sentence: the balls are too large to pass one another, so they stand in a row, the shortest pipe for a given row is the two end radii plus, for each pair of neighbours, a gap that depends only on the sum of their radii, and after proving that the two largest balls belong at the ends, a dynamic programme over the sets of balls still to be placed finds the best order of the other \(19\).

Mathematical Approach

We first work out the exact length that a given order of the balls needs, and then search for the best order: an exchange argument fixes the two end balls, and a dynamic programme does the rest. Lengths are in millimetres, and \(R=50\) is the radius of the pipe.

1. Two balls in the pipe

Take the balls of radii \(30\) and \(31\). Lined up on the axis of the pipe they need \(60+62=122\) mm, but less when they lean against opposite sides of the pipe. The centre of the smaller ball is then \(50-30=20\) mm from the axis and the centre of the other one \(50-31=19\) mm from it on the other side, so the centres are \(39\) mm apart across the pipe. Touching balls have centres \(30+31=61\) mm apart, so by Pythagoras the centres are \(\sqrt{61^2-39^2}=\sqrt{2200}\)\({}\approx46.904\) mm apart along the axis, and the pipe needs \(30+\sqrt{2200}+31\approx107.904\) mm: the second ball slides partly past the first one.

2. The gap between two neighbours

Now take any two balls of radii \(a,b\le R\) and look at the pipe along its axis. A ball of radius \(r\) fits only if its centre is at most \(R-r\) from the axis, so by the triangle inequality the two centres are at most \((R-a)+(R-b)=2R-a-b\) apart across the pipe. If they are \(h\) apart along the axis and \(t\) apart across it, the distance between the centres is \(\sqrt{h^2+t^2}\), ^[2] and the balls do not overlap exactly when this distance is at least \(a+b\). With \(t\le2R-a-b\) and the difference of two squares this gives

$$h^2\ge(a+b)^2-(2R-a-b)^2=4R\,(a+b-R).$$

So for \(a+b\ge R\) two balls need an axial distance of at least \(g(a+b)\), where

$$g(s)=2\sqrt{R\,(s-R)},$$

and this distance is reached exactly when the balls touch each other and touch the wall on opposite sides. The gap depends only on the sum \(s=a+b\). For \(R=50\) we get \(g(61)=\sqrt{2200}\approx46.904\), as in Section 1.

3. The balls stand in a row

Any two of our balls have \(a+b\ge30+31=61\gt R\), so \(g(a+b)\gt0\): no two balls can have their centres at the same position along the axis. Numbering the balls in the order of their centres along the axis gives a sequence of radii \(r_1,r_2,\dots,r_n\) with \(n=21\), and by Section 2 consecutive centres are at least \(g(r_i+r_{i+1})\) apart along the axis.

4. The length of one order

The first and the last ball reach \(r_1\) and \(r_n\) beyond their centres towards the two ends of the pipe, so a pipe that holds the balls in the order \(r_1,\dots,r_n\) has length at least

$$C(r_1,\dots,r_n)=r_1+\sum_{i=1}^{n-1}g(r_i+r_{i+1})+r_n.$$

This bound is reached. Put all centres in one plane through the axis, alternately on its two sides, each ball touching the wall, with consecutive centres exactly \(g(r_i+r_{i+1})\) apart; neighbours then touch, as in Section 2. Two balls that are not neighbours, \(r_i\) and \(r_j\) with \(j\ge i+2\), are at least \(g(r_i+r_{i+1})+g(r_{j-1}+r_j)\) apart along the axis. Every radius is at least \(30\), so \(g(r+c)\ge2\sqrt{50\,(r-20)}\) for any two radii \(r\) and \(c\), and this exceeds \(r\) for \(30\le r\le50\), because \(200\,(r-20)-r^2\)\({}=6000-(100-r)^2\gt0\). Hence the axial distance alone is larger than \(r_i+r_j\), and such balls never touch. No ball sticks out beyond an end ball either: every other centre is at least \(g(61)\approx46.9\) mm further inside, and two radii differ by at most \(20\) mm. So the shortest pipe has length

$$L=\min C(r_1,\dots,r_{21}),$$

taken over all \(21!\approx5.1\cdot10^{19}\) orders of the balls. Reversing an order does not change \(C\).

5. A toy pipe

In the test case of the C++ program, \(R=8\) and the radii are \(5\), \(6\), \(7\) and \(8\). The pair sums run from \(11\) to \(15\), and the gaps are \(g(11)=4\sqrt6\approx9.798\), \(g(12)=8\sqrt2\approx11.314\), \(g(13)=4\sqrt{10}\approx12.649\), \(g(14)=8\sqrt3\approx13.856\) and \(g(15)=4\sqrt{14}\approx14.967\). The conditions of Sections 3 and 4 hold here as well: every pair sum exceeds \(8\), and even the smallest gap is larger than every radius. Up to reversal, which keeps the length, the \(24\) orders reduce to \(12\), for example:

• \(8,6,5,7\): \(15+g(14)+g(11)+g(12)\)\({}\approx49.968\);
• \(6,5,7,8\): \(14+g(11)+g(12)+g(15)\)\({}\approx50.078\);
• \(8,5,6,7\): \(15+g(13)+g(11)+g(13)\)\({}\approx50.096\);
• \(5,8,7,6\): \(11+g(13)+g(15)+g(13)\)\({}\approx51.265\).

Trying all twelve shows that \(8,6,5,7\) is the shortest order, of length \(15+8\sqrt3+4\sqrt6+8\sqrt2\). The two largest balls stand at its ends. Putting them there is not enough by itself, though: \(6,5,7,8\) is shorter than \(8,5,6,7\).

6. A larger sum costs more than it saves

For \(R\le t\lt s\le2R\),

$$\begin{aligned}g(s)-g(t)&=2\sqrt R\,\bigl(\sqrt{s-R}-\sqrt{t-R}\bigr)\\&=\frac{2\sqrt R\,(s-t)}{\sqrt{s-R}+\sqrt{t-R}}\gt s-t,\end{aligned}$$

because both roots in the denominator are at most \(\sqrt R\) and the second one is smaller. In words: raising the sum of two neighbours by one millimetre lengthens their gap by more than one millimetre. All pair sums of our balls lie between \(61\) and \(99\), inside this range. In the toy pipe, \(g(15)-g(14)\approx1.110\gt1\).

7. The two largest balls stand at the ends

Take an order \(x_1,\dots,x_n\) and reverse its first \(k\) balls, for some \(2\le k\le n-1\): the new order is \(x_k,\dots,x_1,x_{k+1},\dots,x_n\). The gaps inside the reversed block stay the same, the end term \(x_1\) becomes \(x_k\), and the gap \(g(x_k+x_{k+1})\) becomes \(g(x_1+x_{k+1})\). So the length changes by

$$\Delta=(x_k-x_1)-\bigl(g(x_k+x_{k+1})-g(x_1+x_{k+1})\bigr),$$

and if \(x_k\gt x_1\), Section 6 with \(s=x_k+x_{k+1}\) and \(t=x_1+x_{k+1}\) gives \(\Delta\lt0\). For example, reversing the first three balls of \(6,5,7,8\) gives \(7,5,6,8\) with \(\Delta=1-(g(15)-g(14))\)\({}\approx-0.110\): the length drops from \(50.078\) to \(49.968\). Hence in a shortest order the first ball is larger than every ball strictly between the two ends, and reversing a final block instead of a first one shows the same for the last ball. So the ball of radius \(50\) cannot stand inside, since the first ball would have to be larger still, and neither can the ball of radius \(49\), since both end balls would have to be larger than \(49\).

8. A dynamic programme over sets of balls

Reversing if necessary, we put the ball of radius \(49\) first and the ball of radius \(50\) last. The \(19!\approx1.2\cdot10^{17}\) orders of the balls \(30,\dots,48\) between them are far too many to try, but the cost of the rest of a row depends only on the last ball placed and on the set of balls still to come. Let \(F(S,x)\) be the shortest length of the rest of the row when the row so far ends with the ball \(x\) and \(S\) is the set of middle balls still to be placed, counting the remaining gaps, the gap to the final ball \(50\) and its end term \(50\). Then

$$F(\varnothing,x)=g(x+50)+50,$$

$$F(S,x)=\min_{r\in S}\bigl(g(x+r)+F(S\setminus\{r\},r)\bigr),$$

$$L=49+F(\{30,31,\dots,48\},49).$$

The recurrence is proved by induction on the size of \(S\): every way to finish the row places some \(r\in S\) next, which costs the gap \(g(x+r)\), and after that the cheapest way to finish is \(F(S\setminus\{r\},r)\) by the induction hypothesis. This is the Held–Karp dynamic programme for the travelling salesman problem, with the two end terms added. ^[3] The ball \(49\) is the last placed ball only at the start; otherwise \(x\) is one of the \(19\) middle balls and \(S\) a nonempty set of the other \(18\), which gives \(19\,(2^{18}-1)+1=4\,980\,718\) values of \(F\) with \(S\ne\varnothing\).

9. The dynamic programme on the toy pipe

In the toy pipe the first ball is \(7\), the last ball is \(8\), and the middle balls are \(5\) and \(6\):

• \(F(\varnothing,5)=g(13)+8\approx20.649\) and \(F(\varnothing,6)=g(14)+8\approx21.856\);
• \(F(\{6\},5)=g(11)+F(\varnothing,6)\)\({}\approx31.654\) and \(F(\{5\},6)=g(11)+F(\varnothing,5)\)\({}\approx30.447\);
• \(F(\{5,6\},7)\) is the smaller of \(g(12)+31.654\approx42.968\) and \(g(13)+30.447\approx43.096\).

So \(L=7+F(\{5,6\},7)\)\({}\approx49.968\), reached by the order \(7,5,6,8\), the reverse of the shortest order of Section 5.

10. The shortest order is a pendulum

Recording which \(r\) attains each minimum and following these choices back gives the order \(49,47,\dots,31,30,32,\dots,48,50\). Read backwards, it is the pendulum order \(50,48,\dots,32,30,31,33,\dots,49\): every second radius from the largest down to the smallest, then the others back up, like the toy order \(8,6,5,7\). A theorem of Supnick explains this. Add an imaginary end point \(D\) with cost \(c(x,D)=x\) to each ball and cost \(c(x,y)=g(x+y)\) between two balls; a row of balls is then a round trip through \(D\) with cost \(C\). Number the balls by size as cities \(1,\dots,21\) and let \(D\) be city \(22\). If a symmetric cost matrix satisfies

$$c_{ir}+c_{js}\le c_{is}+c_{jr}$$

whenever \(i\lt j\), \(r\lt s\) and the four cities are distinct, then the round trip \(1,3,5,\dots,6,4,2\) is a shortest one. ^{[4] [5]} For four balls, with radii \(a_i\lt a_j\) and \(a_r\lt a_s\), the condition compares \(g\) at the smallest and the largest of the four sums with \(g\) at the two middle sums, which have the same total. For \(p\le m,m'\le q\) with \(m+m'=p+q\), squaring turns \(\sqrt p+\sqrt q\le\sqrt m+\sqrt{m'}\) into \(pq\le mm'\), which holds because \(mm'-pq=(m-p)(q-m)\ge0\); we apply this to the sums minus \(R\). If \(D\) is one of the four cities, the condition reads \(a_j-a_i\le g(a_j+a_r)-g(a_i+a_r)\) or the same with \(i,j\) and \(r,s\) exchanged, and this is Section 6. Cutting the round trip \(1,3,\dots,21,22,20,\dots,2\) at \(D=22\) leaves the path \(20,18,\dots,2,1,3,\dots,21\). City \(i\) has radius \(i+29\), so this is the order \(49,47,\dots,31,30,32,\dots,50\) found above, and its two edges to \(D\) are the end terms \(49+50\) of \(C\). Supnick's theorem therefore proves that the pendulum order is optimal, an independent check that the published programs do not use.

How the Code Works

The C++ program ^[8] reads its parameters into Options, by default pipe radius \(50\) and radii \(30\) to \(50\). parse_arguments accepts --skip-checkpoints, --pipe-radius=, --min-radius= and --max-radius=, whose digits parse_int_after_prefix reads without an overflow check; main returns \(1\) for an unknown argument or unless \(1\le\) minimum radius \(\le\) maximum radius \(\le\) pipe radius. The conditions of Sections 3 and 4 are not checked, so other parameters can give wrong lengths: with --pipe-radius=5 --min-radius=2 --max-radius=3 the program prints \(5000\), although the ball of radius \(3\) alone needs \(6000\) micrometres; a pair sum below the pipe radius makes distance_y take the square root of a negative number, and more than \(31\) balls would overflow the \(32\)-bit masks.

Unless --skip-checkpoints is given, run_checkpoints compares solve_optimized(8, 5, 8) with solve_bruteforce, which tries all \(24\) orders of the toy pipe with std::next_permutation and so also tests Section 7 there; if the two lengths differ by more than \(10^{-9}\), main returns \(2\).

distance_y returns \(\sqrt{(a+b)^2-(2R-a-b)^2}\), the gap of Section 2, and solve_optimized implements Section 8. A set \(S\) of middle balls is a bitmask mask whose bit index_of(r), that is \(r-30\), stands for the ball \(r\), and the table cache of \(2^{19}\cdot21\) doubles starts at \(-1\), an impossible length marking entries not yet computed. The recursive std::function dfs(mask, last_radius) returns \(F(S,x)\): for an empty mask the gap to the largest ball plus its radius, otherwise the entry at mask + (index_of(last_radius) << shift) with shift \(=19\). The index is unique because the bits of \(49\) and \(50\) are never set, so mask stays below \(2^{19}\). A new entry is the minimum, over the radii in the mask, of the gap plus dfs of the smaller mask. The start clears the bits of \(49\) and \(50\) from the full mask and returns \(49\) plus dfs(mask, 49), and main prints std::llround(1000.0 * best) as an unsigned 64-bit integer.

Python and Java run only the dynamic programme with the fixed parameters, without options or checkpoint: Python's solve with nested functions and a dictionary keyed by (mask, last_radius), computing shift without using it, and Java with the static methods distanceY, indexOf and dfs and a double[] table filled with \(-1\) by Arrays.fill. Masks are unsigned int in C++, int in Java and arbitrary-precision integers in Python; the recursion is only \(20\) calls deep, and only C++ and Java handle a single ball. At exact halves Python's round goes to the even integer, while std::llround and Java's Math.round go up for positive values; ^[7] this does not matter here, since \(1000L\) is far from a half-integer.

Complexity and Verification

For \(n\) balls the dynamic programme has about \((n-2)\,2^{n-3}\) states, each of which loops over the \(n\) radii, so the time is \(O(n^2\,2^n)\) and the memory \(O(n\,2^n)\), as for Held–Karp. ^[3] Here the \(4\,980\,718\) states try \(44\,826\,643\) next balls in total. The C++ and Java table has \(2^{19}\cdot21=11\,010\,048\) entries of \(8\) bytes, about \(88\) MB, of which \(45\%\) are used. In our tests on a desktop computer, the C++ and Java programs took under half a second and the Python program about \(18\) seconds, with about \(1\) GB of memory for its dictionary of \(4\,980\,718\) values.

The programs compute in IEEE 754 double precision, where square roots and additions are correctly rounded with relative error at most \(2^{-53}\approx1.1\cdot10^{-16}\), ^[6] and each gap is the square root of an exact integer. Since \(g(s)\le s\), because \((s-2R)^2\ge0\), no row is longer than the \(1680\) mm of all balls lined up on the axis, so each of the \(41\) rounded operations along a row, \(20\) square roots and \(21\) additions, errs by less than \(2\cdot10^{-13}\) mm. A minimum of candidates that are each within \(\varepsilon\) of their exact values is within \(\varepsilon\) of the exact minimum, even if rounding changes the winner, so by induction over the levels of the dynamic programme every computed value of \(F\) is within \(10^{-11}\) mm of the exact one; the final multiplication by \(1000\) adds less than \(2\cdot10^{-10}\) micrometres. A \(60\)-digit evaluation of the pendulum order differs from the computed minimum by \(2.5\cdot10^{-13}\) mm, and \(1000L\) lies \(0.38\) micrometres from the nearest half-integer, so the rounding is safe.

The C++ checkpoint passes, and the three programs print the same value. These further checks were made for this article and are not part of the published solutions:

• a dynamic programme with free end balls, over all \(2^{21}\) sets of balls, gives the same minimum, which confirms Section 7 for the real pipe;
• for all \(210\) runs of at least two consecutive radii between \(30\) and \(50\), both dynamic programmes and the pendulum order agree, and for the \(12\) runs of \(10\) radii so does a search through all \(10!=3\,628\,800\) orders;
• the Supnick condition of Section 10 holds for all \(43\,890\) index quadruples with \(i\lt j\), \(r\lt s\) and four distinct cities;
• in a placement of the \(21\) balls in space as in Section 4, neighbours touch, any two non-neighbours are farther apart along the axis than the sum of their radii by at least \(32.89\) mm, and the length equals the computed minimum.

Footnotes and References

The references below supply the background facts used above. The gap formula, the exchange argument for the end balls, the check of the Supnick condition and the rounding analysis are derived in this article.

1. Project Euler 222 — Sphere Packing. The official statement asks for the length of the shortest pipe of internal radius 50 mm that can fully contain 21 balls of radii 30 mm, 31 mm, …, 50 mm, in micrometres and rounded to the nearest integer.
2. Wikipedia — Pythagorean theorem. The squared distance between two points is the sum of the squares of their coordinate differences, in the plane and in space of any dimension; this gives the distance √(h² + t²) in Section 2.
3. Wikipedia — Held–Karp algorithm. The shortest path through all cities of a set S to an end city e is computed from the paths through smaller sets by minimising over the second-to-last city, in time Θ(2ⁿn²) and space Θ(n2ⁿ).
4. Wikipedia — Monge array. A Monge array satisfies A[i,j] + A[k,ℓ] ≤ A[i,ℓ] + A[k,j] for all i < k and j < ℓ; a symmetric square Monge matrix is a Supnick matrix, and a travelling salesman problem with a Supnick cost matrix has a predetermined optimal route.
5. V. G. Deĭneko, J. A. van der Veen, R. Rudolf and G. J. Woeginger — Three easy special cases of the euclidean travelling salesman problem, RAIRO Operations Research 31 (1997), 343–362. The paper calls a symmetric matrix a Supnick matrix if c_ir + c_js ≤ c_is + c_jr for i < j, r < s and {i, j} ∩ {r, s} = ∅, and states, citing F. Supnick (Annals of Mathematics, 1957), that for such matrices the tour (1, 3, 5, 7, …, 8, 6, 4, 2) is optimal.
6. Wikipedia — Double-precision floating-point format: a 53-bit significand and a maximum relative rounding error of 2⁻⁵³ ≈ 1.11 × 10⁻¹⁶. cppreference — std::sqrt and Java — Math.sqrt: the square root is correctly rounded, as IEEE 754 requires. Python — numeric types: floats are usually C doubles, and integers have unlimited precision.
7. cppreference — std::round, std::llround: halfway cases are rounded away from zero. Python — round: ties go to the even choice. Java — Math.round: the closest long, with ties rounding to positive infinity.
8. C++ — Euler222.cpp. The linked immutable revision is the C++ program described above.