Same tangential-container reduction in d dimensions; giving the axis-aligned lattice its own offset sweep removed the variance that capped the first port (0.9732 -> 0.9869, four n exact).

-"""Baseline: a lattice of candidate points inside the container, spacing found by bisection so
-that at least n fit. Deliberately naive; scores well below the records. Beat it."""
+"""Equal spheres in the unit cube [0,1]^3.
 
-import itertools
+Reduction.  The container is tangential (a cube, or a ball), so eroding it by r gives a
+similar copy scaled by (RHO - r)/RHO with RHO its inradius.  Packing n equal objects is
+therefore equivalent to spreading n points in the container so as to maximise the minimum
+pairwise distance m, with
+
+    r = m RHO / (2 RHO + m)          centre_i = q_i (RHO - r) / RHO
+
+That removes the boundary term: the container becomes a projection, and everything below
+maximises m.  This is the same reduction I used for the 2-D containers, which holds in any
+dimension because it only needs the erosion to be similar.
+
+Method.  Inflate-and-separate: hold a target separation d, run Gauss-Seidel overlap-removal
+sweeps, and ratchet d up whenever the achieved minimum distance improves.  Starts are
+undistorted lattices (Z^d and the checkerboard D_d, which is the densest lattice packing in
+3, 4 and 5 dimensions), randomly rotated and offset, cropped by the container gauge to the
+tightest scaled copy holding n points.  Then reinsertion / shake hill climbing, then polish.
+
+Neighbour search switches on dimension: a uniform grid below 4-D, where 3^d
+neighbouring cells is cheaper than n^2/2 pairs, and plain all-pairs above it, where it is not.
+"""
+
 import math
+import random
+import time
 
 DIM = 3
+BALL = False                     # True: unit ball. False: unit cube [0,1]^DIM
+RHO = 0.5                       # inradius; 1.0 for the ball, 0.5 for the centred cube
+SHIFT = 0.5                   # added to every coordinate on output
+_R2 = RHO * RHO
+_USE_GRID = DIM <= 4
+_OFFS = None
 
 
-def _boundary(c):
-    return min(c[0], 1 - c[0], c[1], 1 - c[1], c[2], 1 - c[2])
+def _offsets():
+    global _OFFS
+    if _OFFS is None:
+        acc = [()]
+        for _ in range(DIM):
+            acc = [t + (o,) for t in acc for o in (-1, 0, 1)]
+        _OFFS = acc
+    return _OFFS
 
 
-def _weight(i):
-    return 1.0
+def _gauge(P, b):
+    if BALL:
+        s = 0.0
+        for k in range(DIM):
+            v = P[b + k]
+            s += v * v
+        return math.sqrt(s) / RHO
+    g = 0.0
+    for k in range(DIM):
+        v = P[b + k]
+        if v < 0.0:
+            v = -v
+        if v > g:
+            g = v
+    return g / RHO
 
 
-def _lattice(n, r):
-    """Cubic lattice points at spacing 2r whose distance to the boundary is at least r."""
-    lo, hi = ((0.0,) * DIM, (1.0,) * DIM)
-    step = 2.0 * r
-    axes = []
-    for d in range(DIM):
-        k = int((hi[d] - lo[d]) / step) + 1
-        axes.append([lo[d] + r + i * step for i in range(k)])
-    pts = [p for p in itertools.product(*axes) if _boundary(p) >= r]
-    return pts
+def _proj(P, b):
+    """Project point at offset b back into the container."""
+    if BALL:
+        s = 0.0
+        for k in range(DIM):
+            v = P[b + k]
+            s += v * v
+        if s > _R2:
+            f = RHO / math.sqrt(s)
+            for k in range(DIM):
+                P[b + k] *= f
+    else:
+        for k in range(DIM):
+            v = P[b + k]
+            if v < -RHO:
+                P[b + k] = -RHO
+            elif v > RHO:
+                P[b + k] = RHO
 
 
+def _rand_pt(rng, out):
+    if BALL:
+        while True:
+            s = 0.0
+            v = []
+            for _ in range(DIM):
+                x = rng.uniform(-RHO, RHO)
+                v.append(x)
+                s += x * x
+            if s <= _R2:
+                out.extend(v)
+                return
+    else:
+        for _ in range(DIM):
+            out.append(rng.uniform(-RHO, RHO))
+
+
+def _cells(P, n, cell):
+    inv = 1.0 / cell
+    g = {}
+    for i in range(n):
+        b = i * DIM
+        key = tuple(int(math.floor(P[b + k] * inv)) for k in range(DIM))
+        e = g.get(key)
+        if e is None:
+            g[key] = [i]
+        else:
+            e.append(i)
+    return g, inv
+
+
+def _min_dist(P, n, cell):
+    best = 1e18
+    if _USE_GRID and cell > 1e-12:
+        g, inv = _cells(P, n, cell)
+        offs = _offsets()
+        for i in range(n):
+            bi = i * DIM
+            key = tuple(int(math.floor(P[bi + k] * inv)) for k in range(DIM))
+            for off in offs:
+                lst = g.get(tuple(key[k] + off[k] for k in range(DIM)))
+                if lst is None:
+                    continue
+                for j in lst:
+                    if j <= i:
+                        continue
+                    bj = j * DIM
+                    s = 0.0
+                    for k in range(DIM):
+                        t = P[bj + k] - P[bi + k]
+                        s += t * t
+                    if s < best:
+                        best = s
+        if best < 1e17:
+            return math.sqrt(best)
+    for i in range(n - 1):
+        bi = i * DIM
+        for j in range(i + 1, n):
+            bj = j * DIM
+            s = 0.0
+            for k in range(DIM):
+                t = P[bj + k] - P[bi + k]
+                s += t * t
+            if s < best:
+                best = s
+    return math.sqrt(best)
+
+
+def _min_pair(P, n, cell):
+    best = 1e18
+    bi_, bj_ = 0, 1
+    for i in range(n - 1):
+        bi = i * DIM
+        for j in range(i + 1, n):
+            bj = j * DIM
+            s = 0.0
+            for k in range(DIM):
+                t = P[bj + k] - P[bi + k]
+                s += t * t
+            if s < best:
+                best = s
+                bi_, bj_ = i, j
+    return bi_, bj_
+
+
+def _push(P, i, j, dd, d2, rng):
+    bi = i * DIM
+    bj = j * DIM
+    s = 0.0
+    for k in range(DIM):
+        t = P[bj + k] - P[bi + k]
+        s += t * t
+    if s >= d2:
+        return 0.0
+    if s < 1e-24:
+        for k in range(DIM):
+            P[bj + k] += rng.uniform(-1e-7, 1e-7)
+        return dd
+    dist = math.sqrt(s)
+    gap = dd - dist
+    f = 0.5 * gap / dist
+    for k in range(DIM):
+        t = (P[bj + k] - P[bi + k]) * f
+        P[bj + k] += t
+        P[bi + k] -= t
+    _proj(P, bi)
+    _proj(P, bj)
+    return gap
+
+
+def _separate(P, n, dd, sweeps, rng):
+    d2 = dd * dd
+    offs = _offsets() if _USE_GRID else None
+    for _ in range(sweeps):
+        worst = 0.0
+        if _USE_GRID:
+            g, inv = _cells(P, n, dd)
+            for i in range(n):
+                bi = i * DIM
+                key = tuple(int(math.floor(P[bi + k] * inv)) for k in range(DIM))
+                for off in offs:
+                    lst = g.get(tuple(key[k] + off[k] for k in range(DIM)))
+                    if lst is None:
+                        continue
+                    for j in lst:
+                        if j > i:
+                            w = _push(P, i, j, dd, d2, rng)
+                            if w > worst:
+                                worst = w
+        else:
+            for i in range(n - 1):
+                for j in range(i + 1, n):
+                    w = _push(P, i, j, dd, d2, rng)
+                    if w > worst:
+                        worst = w
+        if worst <= 1e-15:
+            return True
+    return False
+
+
+def _hole(P, n, skip, rng, tries):
+    """Emptiest sampled spot: the reinsertion target. Dropping the capped point somewhere
+    random wastes the move; the 2-D solver showed the search for a real hole is what pays."""
+    bestv = -1.0
+    best = None
+    for _ in range(tries):
+        c = []
+        _rand_pt(rng, c)
+        d = 1e18
+        for j in range(n):
+            if j == skip:
+                continue
+            b = j * DIM
+            s = 0.0
+            for k in range(DIM):
+                t = P[b + k] - c[k]
+                s += t * t
+                if s >= d:
+                    break
+            if s < d:
+                d = s
+        if d > bestv:
+            bestv = d
+            best = c
+    return best
+
+
+def _ratchet(P, n, m, t_end, rng, step, sweeps, floor):
+    while step > floor:
+        if time.perf_counter() >= t_end:
+            break
+        dd = m * (1.0 + step)
+        save = P[:]
+        _separate(P, n, dd, sweeps, rng)
+        m2 = _min_dist(P, n, dd)
+        if m2 > m:
+            m = m2
+            step *= 1.25
+        else:
+            P[:] = save
+            step *= 0.55
+    return m
+
+
+def _lattice_pts(kind, want):
+    """Integer lattice points inside a ball, enough of them, minimal distance normalised."""
+    R = max(2.0, (want ** (1.0 / DIM)))
+    while True:
+        pts = []
+        cur = [0] * DIM
+
+        def rec(k, rem):
+            if k == DIM:
+                if kind == 1:                       # D_d: even coordinate sum
+                    if sum(cur) & 1:
+                        return
+                pts.append(tuple(cur))
+                return
+            lim = int(math.sqrt(rem))
+            for v in range(-lim, lim + 1):
+                cur[k] = v
+                rec(k + 1, rem - v * v)
+            cur[k] = 0
+
+        rec(0, R * R)
+        if len(pts) >= want or R > 40:
+            scale = 1.0 / (math.sqrt(2.0) if kind == 1 else 1.0)
+            return [tuple(c * scale for c in p) for p in pts]
+        R *= 1.3
+
+
+def _rot(rng):
+    """Random orthonormal matrix by Gram-Schmidt on Gaussian columns."""
+    rows = []
+    for _ in range(DIM):
+        while True:
+            v = [rng.gauss(0.0, 1.0) for _ in range(DIM)]
+            for u in rows:
+                d = sum(a * b for a, b in zip(u, v))
+                for k in range(DIM):
+                    v[k] -= d * u[k]
+            nrm = math.sqrt(sum(a * a for a in v))
+            if nrm > 1e-9:
+                rows.append([a / nrm for a in v])
+                break
+    return rows
+
+
+def _crop_starts(n, keep, t_end, rng, tries):
+    """Rotated, offset lattices cropped by the container gauge to n points."""
+    out = []
+    seen = set()
+    for kind in (1, 0):                              # D_d first: denser for DIM 3,4,5
+        base = _lattice_pts(kind, max(3 * n, 60))
+        if len(base) < n:
+            continue
+        for t in range(tries):
+            if time.perf_counter() >= t_end:
+                break
+            # The ball gauge is rotation invariant, so only offsets matter there.  For a
+            # cube, the axis-aligned lattice is a strong candidate and needs its own offset
+            # sweep: rotating every try but the first meant the aligned family was only ever
+            # sampled at zero offset, which is the one offset least likely to crop well.
+            Q = None if (BALL or t < tries // 2) else _rot(rng)
+            sh = [0.0] * DIM
+            if t:
+                for k in range(DIM):
+                    sh[k] = rng.uniform(-0.5, 0.5)
+            flat = []
+            for p in base:
+                if Q is None:
+                    for k in range(DIM):
+                        flat.append(p[k] + sh[k])
+                else:
+                    for k in range(DIM):
+                        flat.append(sum(Q[k][c] * p[c] for c in range(DIM)) + sh[k])
+            cnt = len(base)
+            keyed = sorted(range(cnt), key=lambda i: _gauge(flat, i * DIM))
+            g = _gauge(flat, keyed[n - 1] * DIM)
+            if g <= 1e-9:
+                continue
+            m = 1.0 / g
+            kk = round(m, 10)
+            if kk in seen:
+                continue
+            seen.add(kk)
+            P = []
+            for i in keyed[:n]:
+                b = i * DIM
+                for k in range(DIM):
+                    P.append(flat[b + k] * m)
+            out.append((m, P))
+    out.sort(key=lambda t: -t[0])
+    return out[:keep]
+
+
+def _to_centres(P, n, m):
+    r = m * RHO / (2.0 * RHO + m)
+    s = (RHO - r) / RHO
+    return [tuple(P[i * DIM + k] * s + SHIFT for k in range(DIM)) for i in range(n)]
+
+
 def pack(n, time_budget, seed):
-    # treat every object as the largest one when choosing the lattice spacing
-    wmax = max(_weight(i + 1) for i in range(n))
-    # find a feasible spacing by halving, then bisect between it and the last infeasible one
-    a = 1.0
-    while len(_lattice(n, a * wmax)) < n and a > 1e-9:
-        a /= 2
-    b = 2 * a
-    for _ in range(40):
-        m = (a + b) / 2
-        if len(_lattice(n, m * wmax)) >= n:
-            a = m
+    t0 = time.perf_counter()
+    t_end = t0 + time_budget * 0.90
+    rng = random.Random((seed & 0xFFFFFFFF) * 1000003 + n * 7919 + DIM)
+
+    if n == 1:
+        return [tuple([SHIFT] * DIM)]
+
+    sweeps = 8 if n > 120 else 14
+    tries = 40 if n > 120 else 80
+    cell0 = max(1e-6, 2.0 * RHO / (n ** (1.0 / DIM)))
+
+    lats = _crop_starts(n, 6, t0 + time_budget * 0.20, rng, 26)
+
+    best = None
+    best_m = -1.0
+    starts = 0
+    t_multi = t0 + time_budget * 0.50
+    while True:
+        now = time.perf_counter()
+        if starts and now >= t_multi:
+            break
+        slice_end = min(t_end, now + max(0.25, time_budget * 0.10))
+        if starts < len(lats):
+            P = lats[starts][1][:]
+        elif lats and starts % 3:
+            P = lats[starts % len(lats)][1][:]
+            j = 0.35 * lats[starts % len(lats)][0]
+            for i in range(n):
+                b = i * DIM
+                for k in range(DIM):
+                    P[b + k] += rng.uniform(-j, j)
+                _proj(P, b)
         else:
-            b = m
-    pts = _lattice(n, a * wmax)
-    pts.sort(key=lambda p: -_boundary(p))   # keep the most interior points
-    return [tuple(p) for p in pts[:n]]
+            P = []
+            for _ in range(n):
+                _rand_pt(rng, P)
+        m = _min_dist(P, n, cell0)
+        if m < 1e-9:
+            m = 1e-6
+        m = _ratchet(P, n, m, slice_end, rng, 0.02, sweeps, 1e-7)
+        if m > best_m:
+            best_m, best = m, P[:]
+        starts += 1
 
+    t_moves = t0 + time_budget * 0.82
+    cap = time_budget * 0.06
+    while time.perf_counter() < t_moves:
+        P = best[:]
+        i, j = _min_pair(P, n, best_m)
+        k = i if rng.random() < 0.5 else j
+        b = k * DIM
+        if rng.random() < 0.6:
+            hole = _hole(P, n, k, rng, tries)
+            for c in range(DIM):
+                P[b + c] = hole[c]
+        else:
+            amp = best_m * rng.uniform(0.1, 0.5)
+            for c in range(DIM):
+                P[b + c] += rng.gauss(0.0, amp)
+            _proj(P, b)
+        m = _min_dist(P, n, best_m)
+        if m < 1e-9:
+            continue
+        m = _ratchet(P, n, m, min(t_moves, time.perf_counter() + cap), rng, 0.004, sweeps, 1e-7)
+        if m > best_m:
+            best_m, best = m, P[:]
+
+    best_m = _ratchet(best, n, best_m, t_end, rng, 1e-4, sweeps, 1e-13)
+    return _to_centres(best, n, best_m)
+