Project Euler

Sphere Packing

What is the minimum length of a tube of internal radius R = 50 mm that can contain 21 spheres of radii r = 30, 31,..., 50 mm? Give the answer in micrometers (mu m), rounded to the nearest integer.

Source sync Apr 19, 2026

Problem #0222

Level Level 10

Solved By 2,380

Languages C++, Python

Answer 1590933

Length 485 words

dynamic_programmingoptimizationgeometry

What is the length of the shortest pipe, of internal radius $\SI {50}{\milli \meter }$, that can fully contain $21$ balls of radii $\SI {30}{\milli \meter }, \SI {31}{\milli \meter }, \dots , \SI {50}{\milli \meter }$?

Give your answer in micrometres ($10^{-6} m$) rounded to the nearest integer.

Problem 222: Sphere Packing

Mathematical Foundation

Theorem 1 (Sphere Center Offset). A sphere of radius $r$ placed inside a cylinder of internal radius $R$ (with $r < R$ ) and touching the cylinder wall has its center at distance $R - r$ from the cylinder axis.

Proof. The sphere touches the cylinder wall at a point where the cylinder’s inner surface has distance $R$ from the axis. The center of the sphere is at distance $R - r$ from the axis (the sphere radius measured inward from the contact point). $\square$

Theorem 2 (Vertical Gap Formula). Two tangent spheres of radii $r_i$ and $r_j$ , both touching the cylinder wall of radius $R$ , placed on opposite sides of the cylinder axis, have a vertical separation of

\Delta z(r_i, r_j) = 2\sqrt{R(r_i + r_j - R)}

provided $r_i + r_j > R$ .

Proof. The centers are at distances $R - r_i$ and $R - r_j$ from the axis, on opposite sides ( $\theta = \pi$ ). The horizontal distance between centers is

d_h = (R - r_i) + (R - r_j) = 2R - r_i - r_j.

Since the spheres are tangent, the center-to-center distance is $r_i + r_j$ . By the Pythagorean theorem:

\Delta z^2 = (r_i + r_j)^2 - d_h^2 = (r_i + r_j)^2 - (2R - r_i - r_j)^2.

Let $s = r_i + r_j$ :

\Delta z^2 = s^2 - (2R - s)^2 = s^2 - 4R^2 + 4Rs - s^2 = 4R(s - R).

Hence $\Delta z = 2\sqrt{R(r_i + r_j - R)}$ , valid when $s > R$ . $\square$

Lemma 1 (Validity). For our problem, $r_i + r_j \geq 30 + 31 = 61 > 50 = R$ , so the formula is always valid.

Proof. Immediate from the minimum sphere radii. $\square$

Theorem 3 (Total Tube Length). For a permutation $\sigma$ of the 21 spheres, the minimum tube length is

L(\sigma) = r_{\sigma(1)} + r_{\sigma(21)} + \sum_{k=1}^{20} \Delta z(r_{\sigma(k)}, r_{\sigma(k+1)}).

Proof. The first sphere requires $r_{\sigma(1)}$ clearance from one end, the last requires $r_{\sigma(21)}$ from the other end. Between consecutive spheres, the vertical space consumed is $\Delta z$ . $\square$

Theorem 4 (Optimality via Bitmask DP). The minimum of $L(\sigma)$ over all $21!$ permutations can be computed exactly by dynamic programming on subsets.

Proof. Define $\text{dp}[S][j]$ as the minimum partial tube length when the set $S$ of spheres has been placed and $j \in S$ is the last sphere. Initialize $\text{dp}[\{j\}][j] = r_j$ . The recurrence

\text{dp}[S \cup \{k\}][k] = \min_{j \in S}\bigl(\text{dp}[S][j] + \Delta z(r_j, r_k)\bigr)

correctly computes the optimum by the principle of optimality (Bellman). The final answer is $L^* = \min_j(\text{dp}[\text{all}][j] + r_j)$ . $\square$

Editorial

Key insight: spheres alternate sides in the tube. The vertical gap between consecutive spheres i,j (on opposite sides) is 2sqrt(R(ri+rj-R)). We use bitmask DP to find the optimal ordering. Note: 2^21 * 21 ~ 44M states, feasible in C++ but slow in Python. For Python, we use a greedy/heuristic approach validated against the known answer, or we optimize the DP. Actually, with 21 spheres the DP has 2^21 = 2M masks * 21 = ~44M states. In Python this is too slow with dicts, so we use a known mathematical insight: the optimal arrangement alternates large and small spheres, and we can verify with a reduced search. The optimal ordering can be found by noting that we should interleave small and large radii. Specifically, arrange as: 30, 50, 31, 49, 32, 48, …, 39, 41, 40 But one sphere (the smallest) might be better dropped to the end. We’ll try all arrangements of the “interleaved” type and pick the best. We iterate over j in S. Finally, iterate over k not in S.

Pseudocode

for j in S
for k not in S

Complexity Analysis

Time: $O(2^n \cdot n^2)$ where $n = 21$ . This gives $2^{21} \times 21^2 \approx 9.2 \times 10^8$ operations.
Space: $O(2^n \cdot n)$ for the DP table, i.e., $2^{21} \times 21 \approx 4.4 \times 10^7$ entries.

Answer

\boxed{1590933}

C++ project_euler/problem_222/solution.cpp

#include <bits/stdc++.h>
using namespace std;

int main(){
    // Sphere packing in a tube of radius R=50mm.
    // 21 spheres with radii 30..50 mm.
    // Bitmask DP to find optimal ordering.

    const int N = 21;
    const double R = 50.0;
    double radii[N];
    for(int i = 0; i < N; i++) radii[i] = 30.0 + i;

    // Precompute vertical gaps
    double dz[N][N];
    for(int i = 0; i < N; i++)
        for(int j = 0; j < N; j++){
            double s = radii[i] + radii[j];
            double val = 4.0 * R * (s - R);
            dz[i][j] = (val > 0) ? sqrt(val) : 0.0;
        }

    // dp[mask][last] = min partial length
    int FULL = (1 << N);
    vector<vector<double>> dp(FULL, vector<double>(N, 1e18));

    // Initialize: place one sphere
    for(int i = 0; i < N; i++)
        dp[1 << i][i] = radii[i];

    for(int mask = 1; mask < FULL; mask++){
        for(int j = 0; j < N; j++){
            if(!(mask & (1 << j))) continue;
            if(dp[mask][j] >= 1e17) continue;
            // Try adding sphere k
            for(int k = 0; k < N; k++){
                if(mask & (1 << k)) continue;
                int nmask = mask | (1 << k);
                double cost = dp[mask][j] + dz[j][k];
                if(cost < dp[nmask][k])
                    dp[nmask][k] = cost;
            }
        }
    }

    double best = 1e18;
    for(int j = 0; j < N; j++)
        best = min(best, dp[FULL-1][j] + radii[j]);

    // Convert mm to micrometers
    long long ans = llround(best * 1000.0);
    cout << ans << endl;

    return 0;
}

Python project_euler/problem_222/solution.py

"""
Problem 222: Sphere Packing

Find the minimum length tube (radius 50mm) to contain 21 spheres of radii 30-50mm.

Key insight: spheres alternate sides in the tube. The vertical gap between
consecutive spheres i,j (on opposite sides) is 2*sqrt(R*(ri+rj-R)).

We use bitmask DP to find the optimal ordering.
Note: 2^21 * 21 ~ 44M states, feasible in C++ but slow in Python.
For Python, we use a greedy/heuristic approach validated against the known answer,
or we optimize the DP.

Actually, with 21 spheres the DP has 2^21 = 2M masks * 21 = ~44M states.
In Python this is too slow with dicts, so we use a known mathematical insight:
the optimal arrangement alternates large and small spheres, and we can verify
with a reduced search.

The optimal ordering can be found by noting that we should interleave
small and large radii. Specifically, arrange as:
30, 50, 31, 49, 32, 48, ..., 39, 41, 40
But one sphere (the smallest) might be better dropped to the end.

We'll try all arrangements of the "interleaved" type and pick the best.
"""

import math
from itertools import permutations

def tube_length(order, R=50.0):
    """Compute tube length for a given ordering of sphere radii."""
    total = order[0] + order[-1]  # first and last sphere touch the ends
    for i in range(len(order) - 1):
        s = order[i] + order[i+1]
        total += 2 * math.sqrt(R * (s - R))
    return total

def solve():
    R = 50.0
    radii = list(range(30, 51))  # 30 to 50, 21 spheres
    N = len(radii)

    # Optimal strategy: interleave small and large.
    # Try different interleaving patterns.
    # Pattern: pick from both ends alternately.

    best = float('inf')
    best_order = None

    # Try all 2^21 isn't feasible in Python, but we can try systematic interleavings.
    # The key insight: split into two groups and interleave.
    # Group A (small): 30-39 (10 spheres), Group B (large): 40-50 (11 spheres)
    # Interleave: B, A, B, A, ..., B

    # Actually, let's just try: for each possible "excluded from alternation" sphere,
    # try placing it at each position.

    # Better: the known optimal for this problem uses a specific interleaving.
    # Let's try the approach: split into odds and evens by index.

    # Approach: try all ways to split 21 spheres into two groups and interleave.
    # The DP approach with bitmask is the correct one. Let's implement it
    # efficiently using numpy-style or just accept it'll be slow.

    # For Python, let's use a simplified bitmask DP with pruning.
    # Actually, with careful implementation using arrays, 44M states might work.

    # Alternative fast approach for Python: since we know the answer,
    # let's use the greedy nearest-neighbor heuristic with verification.

    # Actually the cleanest Python approach: use the known structure.
    # The optimal arrangement places spheres in a zigzag pattern.
    # All spheres with r_i + r_j > R (always true here) go on alternating sides.
    # We need to minimize sum of 2*sqrt(R*(r_i + r_{i+1} - R)).
    # This is a TSP-like problem. For 21 cities, we can solve it with DP.

    # Let's implement bitmask DP in Python, but with bit operations.
    # 2^21 = 2097152, times 21 = ~44M. We'll use a flat array.

    FULL = 1 << N
    INF = float('inf')

    # dp[mask * N + last] = min cost
    dp = [INF] * (FULL * N)

    # Precompute dz
    dz = [[0.0]*N for _ in range(N)]
    for i in range(N):
        for j in range(N):
            s = radii[i] + radii[j]
            dz[i][j] = 2.0 * math.sqrt(R * (s - R))

    # Initialize
    for i in range(N):
        dp[(1 << i) * N + i] = radii[i]

    for mask in range(1, FULL):
        for j in range(N):
            if not (mask & (1 << j)):
                continue
            val = dp[mask * N + j]
            if val >= INF:
                continue
            for k in range(N):
                if mask & (1 << k):
                    continue
                nmask = mask | (1 << k)
                nval = val + dz[j][k]
                idx = nmask * N + k
                if nval < dp[idx]:
                    dp[idx] = nval

    best = INF
    full_base = (FULL - 1) * N
    for j in range(N):
        val = dp[full_base + j] + radii[j]
        if val < best:
            best = val

    ans = round(best * 1000)
    print(ans)

if __name__ == "__main__":
    solve()