Rectangle erosion is not similar, so ratchet on r directly instead of the max-min-distance reduction; fitting each start into the eroded box first took n=101 from 0.359 to 0.990.

-"""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 circles in the rectangle [0,1] x [0,0.4].
 
-import itertools
+Why this one is different.  For a square, a disc or a regular polygon, eroding the container
+by r gives a similar copy, so packing reduces to spreading points to maximise the minimum
+pairwise distance and the radius follows in closed form.  A rectangle is not similar to its
+erosion: [0,1]x[0,H] erodes to a box of size (1-2r) x (H-2r), whose aspect ratio moves with r.
+So there is no fixed target to maximise, and the solver optimises r itself.
+
+Method.  Ratchet on r directly.  The achieved radius of a configuration is
+
+    r = min( wall clearance , half the smallest pairwise distance )
+
+Hold a target r', run Gauss-Seidel sweeps that push every pair out to 2r' and clamp every
+centre into the eroded box [r', 1-r'] x [r', H-r'], then recompute the achieved r and keep the
+move if it improved.  Same inflate-and-separate skeleton as my other entries, but the
+container shrinks as the target grows, which is what the non-similar erosion costs.
+
+Starts: row structures (the natural family for a thin box) and triangular / square lattices
+cropped to the container aspect, both ranked by achieved r before any relaxation.
+"""
+
 import math
+import random
+import time
 
-DIM = 2
+H = 0.4
 
 
-def _boundary(c):
-    return min(c[0], 1 - c[0], c[1], 0.4 - c[1])
+def _achieved_r(xs, ys, n, cell):
+    """min(wall clearance, half the closest pair), grid-accelerated."""
+    r = 1.0
+    for i in range(n):
+        x = xs[i]
+        y = ys[i]
+        v = x if x < 1.0 - x else 1.0 - x
+        w = y if y < H - y else H - y
+        if w < v:
+            v = w
+        if v < r:
+            r = v
+    if cell <= 1e-12:
+        cell = 1e-12
+    inv = 1.0 / cell
+    g = {}
+    for i in range(n):
+        k = (int(math.floor(xs[i] * inv)), int(math.floor(ys[i] * inv)))
+        b = g.get(k)
+        if b is None:
+            g[k] = [i]
+        else:
+            b.append(i)
+    best = 1e18
+    for i in range(n):
+        xi = xs[i]
+        yi = ys[i]
+        cx = int(math.floor(xi * inv))
+        cy = int(math.floor(yi * inv))
+        for ox in (-1, 0, 1):
+            for oy in (-1, 0, 1):
+                lst = g.get((cx + ox, cy + oy))
+                if lst is None:
+                    continue
+                for j in lst:
+                    if j <= i:
+                        continue
+                    dx = xs[j] - xi
+                    dy = ys[j] - yi
+                    d = dx * dx + dy * dy
+                    if d < best:
+                        best = d
+    if best > 1e17:
+        for i in range(n - 1):
+            for j in range(i + 1, n):
+                dx = xs[j] - xs[i]
+                dy = ys[j] - ys[i]
+                d = dx * dx + dy * dy
+                if d < best:
+                    best = d
+    h = math.sqrt(best) * 0.5
+    return h if h < r else r
 
 
-def _weight(i):
-    return 1.0
+def _fit(xs, ys, n, cell):
+    """Shrink a start into the box eroded by r, choosing the largest r it supports.
 
+    The row and lattice starts span [0,1]x[0,H] exactly, so their wall clearance is 0 and
+    their achieved radius is 0 however good the structure is.  Left uncorrected the ratchet
+    has to climb from ~0, which small n survives and large n does not.  Bisect instead: map
+    the start into [r,1-r]x[r,H-r] and keep the largest r for which the mapped configuration
+    still has every pair at least 2r apart.
+    """
+    lo = 0.0
+    hi = 0.5 * (H if H < 1.0 else 1.0)
+    bx = by = None
+    br = 0.0
+    for _ in range(28):
+        r = 0.5 * (lo + hi)
+        sx = 1.0 - 2.0 * r
+        sy = (H - 2.0 * r) / H
+        if sx <= 0.0 or sy <= 0.0:
+            hi = r
+            continue
+        mx = [r + v * sx for v in xs]
+        my = [r + v * sy for v in ys]
+        got = _achieved_r(mx, my, n, max(1e-9, 2.0 * r) if r > 0 else cell)
+        if got >= r:
+            br = got
+            bx, by = mx, my
+            lo = r
+        else:
+            hi = r
+    if bx is None:
+        return xs[:], ys[:], _achieved_r(xs, ys, n, cell)
+    return bx, by, br
 
-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 _min_pair(xs, ys, n):
+    best = 1e18
+    bi, bj = 0, 1
+    for i in range(n - 1):
+        xi = xs[i]
+        yi = ys[i]
+        for j in range(i + 1, n):
+            dx = xs[j] - xi
+            dy = ys[j] - yi
+            d = dx * dx + dy * dy
+            if d < best:
+                best = d
+                bi, bj = i, j
+    return bi, bj
 
+
+def _separate(xs, ys, n, r, sweeps, rng):
+    """Push every pair out to 2r and keep every centre inside the box eroded by r."""
+    dd = 2.0 * r
+    d2 = dd * dd
+    lo = r
+    hix = 1.0 - r
+    hiy = H - r
+    if hix < lo or hiy < lo:
+        return False
+    inv = 1.0 / dd
+    for _ in range(sweeps):
+        g = {}
+        for i in range(n):
+            k = (int(math.floor(xs[i] * inv)), int(math.floor(ys[i] * inv)))
+            b = g.get(k)
+            if b is None:
+                g[k] = [i]
+            else:
+                b.append(i)
+        worst = 0.0
+        for i in range(n):
+            xi = xs[i]
+            yi = ys[i]
+            cx = int(math.floor(xi * inv))
+            cy = int(math.floor(yi * inv))
+            for ox in (-1, 0, 1):
+                for oy in (-1, 0, 1):
+                    lst = g.get((cx + ox, cy + oy))
+                    if lst is None:
+                        continue
+                    for j in lst:
+                        if j <= i:
+                            continue
+                        dx = xs[j] - xi
+                        dy = ys[j] - yi
+                        s = dx * dx + dy * dy
+                        if s >= d2:
+                            continue
+                        if s < 1e-24:
+                            a = rng.random() * 6.283185307179586
+                            dx = math.cos(a) * 1e-7
+                            dy = math.sin(a) * 1e-7
+                            s = 1e-14
+                        dist = math.sqrt(s)
+                        gap = dd - dist
+                        if gap > worst:
+                            worst = gap
+                        f = 0.5 * gap / dist
+                        ux = dx * f
+                        uy = dy * f
+                        xj = xs[j] + ux
+                        yj = ys[j] + uy
+                        xs[j] = lo if xj < lo else (hix if xj > hix else xj)
+                        ys[j] = lo if yj < lo else (hiy if yj > hiy else yj)
+                        xi -= ux
+                        yi -= uy
+                        if xi < lo:
+                            xi = lo
+                        elif xi > hix:
+                            xi = hix
+                        if yi < lo:
+                            yi = lo
+                        elif yi > hiy:
+                            yi = hiy
+            xs[i] = xi
+            ys[i] = yi
+        if worst <= 1e-15:
+            return True
+    return False
+
+
+def _ratchet(xs, ys, n, r, t_end, rng, step, sweeps, floor):
+    while step > floor:
+        if time.perf_counter() >= t_end:
+            break
+        tgt = r * (1.0 + step)
+        sx = xs[:]
+        sy = ys[:]
+        _separate(xs, ys, n, tgt, sweeps, rng)
+        r2 = _achieved_r(xs, ys, n, 2.0 * tgt)
+        if r2 > r:
+            r = r2
+            step *= 1.25
+        else:
+            xs[:] = sx
+            ys[:] = sy
+            step *= 0.55
+    return r
+
+
+def _rows_starts(n):
+    """Rows of points spanning the box: the natural family for a thin rectangle."""
+    out = []
+    top = int(math.sqrt(2.0 * n / max(H, 1e-9))) + 3
+    for R in range(1, min(n, top) + 1):
+        base = n // R
+        extra = n % R
+        for stagger in (0, 1):
+            xs = []
+            ys = []
+            for i in range(R):
+                cnt = base + (1 if i < extra else 0)
+                if cnt <= 0:
+                    continue
+                y = H * (i + 0.5) / R
+                for j in range(cnt):
+                    if cnt == 1:
+                        x = 0.5
+                    else:
+                        x = (j + 0.5 * (i % 2) * stagger) / (cnt - 1 + 0.5 * (i % 2) * stagger)
+                        x = min(1.0, max(0.0, x))
+                    xs.append(x)
+                    ys.append(y)
+            if len(xs) == n:
+                out.append((xs, ys))
+    return out
+
+
+def _crop_starts(n, tries, rng):
+    """Triangular and square lattices cropped to a box of the container aspect."""
+    out = []
+    ax = 0.5
+    ay = 0.5 * H
+    for u1, u2, v1, v2, maxrad, app in ((1.0, 0.0, 0.5, 0.8660254037844386,
+                                         1.0471975511965976, 0.8660254037844386),
+                                        (1.0, 0.0, 0.0, 1.0, 1.5707963267948966, 1.0)):
+        est = math.sqrt(n * app / (4.0 * ax * ay))
+        rad = est * max(ax, ay) * 2.6 + 3.0
+        ra = int(rad * 2.4) + 3
+        raw = [(a * u1 + b * v1, a * u2 + b * v2)
+               for a in range(-ra, ra + 1) for b in range(-ra, ra + 1)
+               if (a * u1 + b * v1) ** 2 + (a * u2 + b * v2) ** 2 <= rad * rad]
+        if len(raw) < n:
+            continue
+        for t in range(tries):
+            th = maxrad * (t % 8) / 8.0
+            c = math.cos(th)
+            s = math.sin(th)
+            sx = rng.uniform(-0.5, 0.5) if t else 0.0
+            sy = rng.uniform(-0.5, 0.5) if t else 0.0
+            keyed = []
+            for x, y in raw:
+                px = x * c - y * s + sx
+                py = x * s + y * c + sy
+                gx = px / ax
+                gy = py / ay
+                if gx < 0.0:
+                    gx = -gx
+                if gy < 0.0:
+                    gy = -gy
+                keyed.append((gx if gx > gy else gy, px, py))
+            if len(keyed) < n:
+                continue
+            keyed.sort()
+            g = keyed[n - 1][0]
+            if g <= 1e-9:
+                continue
+            f = 1.0 / g
+            xs = [p[1] * f + ax for p in keyed[:n]]
+            ys = [p[2] * f + ay for p in keyed[:n]]
+            out.append((xs, ys))
+    return out
+
+
 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 * 131 + int(H * 1000))
+
+    if n == 1:
+        return [(0.5, H * 0.5)]
+
+    sweeps = 12 if n > 90 else 18
+    cell0 = max(1e-4, math.sqrt(H / n))
+
+    cands = _rows_starts(n) + _crop_starts(n, 16, rng)
+    scored = []
+    seen = set()
+    for xs0, ys0 in cands:
+        xs, ys, r = _fit(xs0, ys0, n, cell0)
+        k = round(r, 10)
+        if k in seen:
+            continue
+        seen.add(k)
+        scored.append((r, xs, ys))
+    scored.sort(key=lambda t: -t[0])
+    scored = scored[:8]
+
+    best_xs = best_ys = None
+    best_r = -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.09))
+        if starts < len(scored):
+            xs = scored[starts][1][:]
+            ys = scored[starts][2][:]
+        elif scored and starts % 3:
+            src = scored[starts % len(scored)]
+            j = 0.35 * max(src[0], 1e-6)
+            xs = [min(1.0, max(0.0, v + rng.uniform(-j, j))) for v in src[1]]
+            ys = [min(H, max(0.0, v + rng.uniform(-j, j))) for v in src[2]]
         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]]
+            xs = [rng.random() for _ in range(n)]
+            ys = [rng.uniform(0.0, H) for _ in range(n)]
+        r = _achieved_r(xs, ys, n, cell0)
+        if r < 1e-9:
+            r = 1e-7
+        r = _ratchet(xs, ys, n, r, slice_end, rng, 0.02, sweeps, 1e-7)
+        if r > best_r:
+            best_r, best_xs, best_ys = r, xs[:], ys[:]
+        starts += 1
 
+    t_moves = t0 + time_budget * 0.82
+    cap = time_budget * 0.06
+    while time.perf_counter() < t_moves:
+        xs = best_xs[:]
+        ys = best_ys[:]
+        i, j = _min_pair(xs, ys, n)
+        k = i if rng.random() < 0.5 else j
+        if rng.random() < 0.6:
+            bx = by = 0.0
+            bd = -1.0
+            for _ in range(200):
+                px = rng.random()
+                py = rng.uniform(0.0, H)
+                d = min(px, 1.0 - px, py, H - py)
+                dd = d * d * 4.0
+                for q in range(n):
+                    if q == k:
+                        continue
+                    ex = xs[q] - px
+                    ey = ys[q] - py
+                    e = ex * ex + ey * ey
+                    if e < dd:
+                        dd = e
+                        if dd < bd:
+                            break
+                if dd > bd:
+                    bd = dd
+                    bx = px
+                    by = py
+            xs[k] = bx
+            ys[k] = by
+        else:
+            amp = best_r * rng.uniform(0.2, 1.0)
+            xs[k] = min(1.0, max(0.0, xs[k] + rng.gauss(0.0, amp)))
+            ys[k] = min(H, max(0.0, ys[k] + rng.gauss(0.0, amp)))
+        r = _achieved_r(xs, ys, n, max(1e-6, 2.0 * best_r))
+        if r < 1e-9:
+            continue
+        r = _ratchet(xs, ys, n, r, min(t_moves, time.perf_counter() + cap),
+                     rng, 0.004, sweeps, 1e-7)
+        if r > best_r:
+            best_r, best_xs, best_ys = r, xs[:], ys[:]
+
+    _ratchet(best_xs, best_ys, n, best_r, t_end, rng, 1e-4, sweeps, 1e-13)
+    return [(best_xs[i], best_ys[i]) for i in range(n)]
+