Erosion is not similar here, so ratchet on r directly with the container as half planes plus a disc; starts must be fitted about the incentre or their achieved r is 0.

-"""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 a right triangle.
 
-import itertools
+Like the rectangle family and unlike the square, the disc and the regular polygons, this
+container is not similar to its own erosion, so there is no fixed max-min-distance target to
+optimise and the solver works on r directly.
+
+The container is described once, as unit-normal half planes plus an optional bounding disc:
+
+    boundary(p) = min over half planes of (c - n.p), and (R - |p|) if there is a disc
+
+Eroding by r lowers every c by r and R by r, which is all the ratchet needs.
+
+Method.  Hold a target r', push every pair out to 2r', project every centre into the eroded
+container, recompute r = min(boundary, half the closest pair), keep the move if it improved.
+Starts are lattices filled into the container at the largest spacing that still fits n points,
+then scaled about the incentre to the factor that maximises achieved r before any relaxation.
+"""
+
 import math
+import random
+import time
 
-DIM = 2
+HP = [(-1.0, 0.0, 0.0), (0.0, -1.0, 0.0), (0.7071067811865475, 0.7071067811865475, 0.7071067811865475)]                      # (nx, ny, c) with unit normal: n.p <= c
+DISK = 0.0                  # bounding radius, or 0.0 for none
 
 
-def _boundary(c):
-    return min(c[0], c[1], (1 - c[0] - c[1]) / math.sqrt(2))
+def _bd(x, y):
+    b = 1e18
+    for nx, ny, c in HP:
+        v = c - nx * x - ny * y
+        if v < b:
+            b = v
+    if DISK:
+        v = DISK - math.sqrt(x * x + y * y)
+        if v < b:
+            b = v
+    return b
 
 
-def _weight(i):
-    return 1.0
+def _bbox():
+    lo = -DISK if DISK else -2.0
+    hi = DISK if DISK else 2.0
+    xs = [lo, hi]
+    ys = [lo, hi]
+    return min(xs), max(xs), min(ys), max(ys)
 
 
-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 _incentre():
+    """Deepest interior point, found numerically so the container spec stays the only input."""
+    x0, x1, y0, y1 = _bbox()
+    bx = by = 0.0
+    bb = -1e18
+    step = (x1 - x0) / 40.0
+    for i in range(41):
+        for j in range(41):
+            x = x0 + (x1 - x0) * i / 40.0
+            y = y0 + (y1 - y0) * j / 40.0
+            v = _bd(x, y)
+            if v > bb:
+                bb, bx, by = v, x, y
+    for _ in range(400):
+        improved = False
+        for dx, dy in ((step, 0.0), (-step, 0.0), (0.0, step), (0.0, -step),
+                       (step, step), (step, -step), (-step, step), (-step, -step)):
+            v = _bd(bx + dx, by + dy)
+            if v > bb:
+                bb, bx, by = v, bx + dx, by + dy
+                improved = True
+                break
+        if not improved:
+            step *= 0.5
+            if step < 1e-15:
+                break
+    return bx, by, bb
 
 
+CX, CY, RHO = _incentre()
+
+
+def _proj(x, y, r):
+    """Nearest-ish point of the container eroded by r; always returns something feasible."""
+    for _ in range(6):
+        worst = 0.0
+        wx = wy = 0.0
+        disk = False
+        for nx, ny, c in HP:
+            e = nx * x + ny * y - (c - r)
+            if e > worst:
+                worst, wx, wy, disk = e, nx, ny, False
+        if DISK:
+            d = math.sqrt(x * x + y * y)
+            e = d - (DISK - r)
+            if e > worst:
+                worst, disk = e, True
+        if worst <= 0.0:
+            return x, y
+        if disk:
+            d = math.sqrt(x * x + y * y)
+            if d > 1e-15:
+                f = (DISK - r) / d
+                x *= f
+                y *= f
+        else:
+            x -= worst * wx
+            y -= worst * wy
+    # guaranteed fallback: walk toward the incentre until feasible
+    for _ in range(80):
+        if _bd(x, y) >= r - 1e-15:
+            break
+        x = 0.5 * (x + CX)
+        y = 0.5 * (y + CY)
+    return x, y
+
+
+def _rand_pt(rng):
+    x0, x1, y0, y1 = _bbox()
+    for _ in range(2000):
+        x = rng.uniform(x0, x1)
+        y = rng.uniform(y0, y1)
+        if _bd(x, y) > 0.0:
+            return x, y
+    return CX, CY
+
+
+def _achieved_r(xs, ys, n, cell):
+    r = 1e18
+    for i in range(n):
+        v = _bd(xs[i], ys[i])
+        if v < r:
+            r = v
+    if r < 0.0:
+        return -1.0
+    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 _fit(xs, ys, n, cell):
+    """Scale a start about the incentre to the factor that maximises achieved r."""
+    bx = by = None
+    br = -1.0
+    f = 1.0
+    for _ in range(46):
+        mx = [CX + (v - CX) * f for v in xs]
+        my = [CY + (v - CY) * f for v in ys]
+        r = _achieved_r(mx, my, n, cell)
+        if r > br:
+            br, bx, by = r, mx, my
+        f *= 0.97
+    return (bx, by, br) if bx is not None else (xs[:], ys[:], -1.0)
+
+
+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):
+    dd = 2.0 * r
+    d2 = dd * dd
+    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
+                        px, py = _proj(xs[j] + ux, ys[j] + uy, r)
+                        xs[j] = px
+                        ys[j] = py
+                        xi, yi = _proj(xi - ux, yi - uy, r)
+            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 _fill(n, rng, tries):
+    """Lattices laid over the container at the largest spacing that still fits n points."""
+    out = []
+    x0, x1, y0, y1 = _bbox()
+    for u1, u2, v1, v2, dom in ((1.0, 0.0, 0.5, 0.8660254037844386, 1.0471975511965976),
+                                (1.0, 0.0, 0.0, 1.0, 1.5707963267948966)):
+        for t in range(tries):
+            th = dom * (t % 6) / 6.0
+            c = math.cos(th)
+            s = math.sin(th)
+            ox = rng.uniform(-0.5, 0.5) if t else 0.0
+            oy = rng.uniform(-0.5, 0.5) if t else 0.0
+            lo, hi = 1e-3, 2.0 * (x1 - x0)
+            best = None
+            for _ in range(34):
+                sp = 0.5 * (lo + hi)
+                pts = []
+                rng_a = int((x1 - x0) / sp) + 3
+                rng_b = int((y1 - y0) / sp) + 3
+                for a in range(-rng_a, rng_a + 1):
+                    for b in range(-rng_b, rng_b + 1):
+                        px = (a * u1 + b * v1 + ox) * sp
+                        py = (a * u2 + b * v2 + oy) * sp
+                        qx = px * c - py * s
+                        qy = px * s + py * c
+                        if _bd(qx, qy) >= 0.0:
+                            pts.append((qx, qy))
+                            if len(pts) > n + 60:
+                                break
+                    if len(pts) > n + 60:
+                        break
+                if len(pts) >= n:
+                    best = pts
+                    lo = sp
+                else:
+                    hi = sp
+            if best:
+                best.sort(key=lambda p: -_bd(p[0], p[1]))
+                out.append(([p[0] for p in best[:n]], [p[1] for p in best[:n]]))
+    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 * 271)
+
+    if n == 1:
+        return [(CX, CY)]
+
+    sweeps = 12 if n > 90 else 18
+    cell0 = max(1e-4, 2.0 * RHO / math.sqrt(n))
+
+    cands = _fill(n, rng, 8)
+    scored = []
+    seen = set()
+    for xs0, ys0 in cands:
+        xs, ys, r = _fit(xs0, ys0, n, cell0)
+        k = round(r, 10)
+        if k in seen or r <= 0.0:
+            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 = []
+            ys = []
+            for i in range(n):
+                px, py = _proj(src[1][i] + rng.uniform(-j, j),
+                               src[2][i] + rng.uniform(-j, j), 0.0)
+                xs.append(px)
+                ys.append(py)
         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 = []
+            ys = []
+            for _ in range(n):
+                px, py = _rand_pt(rng)
+                xs.append(px)
+                ys.append(py)
+        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 = CX, CY
+            bx, by = CX, CY
+            bd = -1.0
+            for _ in range(160):
+                px, py = _rand_pt(rng)
+                dd = _bd(px, py) ** 2 * 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, bx, by = dd, px, py
+            xs[k] = bx
+            ys[k] = by
+        else:
+            amp = best_r * rng.uniform(0.2, 1.0)
+            px, py = _proj(xs[k] + rng.gauss(0.0, amp), ys[k] + rng.gauss(0.0, amp), 0.0)
+            xs[k] = px
+            ys[k] = py
+        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)]
+