Container-general reduction r = m*RHO/(2*RHO+m) with gauge-cropped undistorted lattice starts; radius prefilter (gauge >= |p|) pays for a fine angle/offset sweep.

-"""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 regular heptagon.
 
-import itertools
+Reduction.  The container is a regular polygon of circumradius 1 (inradius RHO), or a
+disc.  Eroding it by r gives a similar copy scaled by (RHO - r) / RHO, so packing n equal
+circles is equivalent to spreading n points in the container 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.  (For the unit square this specialises to the familiar r = m / (2(1+m)).)
+
+Method.  Inflate-and-separate: hold a target separation d, run Gauss-Seidel overlap-removal
+sweeps over a uniform grid with cell d (expected O(n) per sweep), and ratchet d up whenever
+the achieved minimum distance improves.  Starts come from undistorted triangular and square
+lattices, rotated and offset, cropped by the container's own gauge to the tightest scaled
+copy holding n points -- an undistorted lattice beats one forced to span the container.
+Then reinsertion / shake hill climbing, then one fine polish.
+"""
+
 import math
+import random
+import time
 
-DIM = 2
+SIDES = 7                       # 0 = disc
+RHO = 0.9009688679024191                           # inradius when the circumradius is 1
+_NRM = tuple((math.cos((2 * j + 1) * math.pi / SIDES),
+              math.sin((2 * j + 1) * math.pi / SIDES)) for j in range(SIDES)) if SIDES else ()
 
 
-def _boundary(c):
-    return min(1.0 * math.cos(math.pi / 7) - c[0] * math.cos((2 * j + 1) * math.pi / 7) - c[1] * math.sin((2 * j + 1) * math.pi / 7) for j in range(7))
+def _gauge(x, y):
+    """1.0 exactly on the boundary, <1 inside; scaling a point by 1/gauge lands it on it."""
+    if not SIDES:
+        return math.sqrt(x * x + y * y) / RHO
+    g = -9.0
+    for nx, ny in _NRM:
+        v = x * nx + y * ny
+        if v > g:
+            g = v
+    return g / RHO
 
 
-def _weight(i):
-    return 1.0
+def _proj(x, y):
+    """Nearest-ish point of the container; always returns something feasible."""
+    if not SIDES:
+        d = math.sqrt(x * x + y * y)
+        if d <= RHO:
+            return x, y
+        f = RHO / d
+        return x * f, y * f
+    for _ in range(3):
+        worst = 0.0
+        wx = wy = 0.0
+        for nx, ny in _NRM:
+            e = x * nx + y * ny - RHO
+            if e > worst:
+                worst = e
+                wx = nx
+                wy = ny
+        if worst <= 0.0:
+            return x, y
+        x -= worst * wx
+        y -= worst * wy
+    g = _gauge(x, y)
+    if g > 1.0:
+        f = 1.0 / g
+        x *= f
+        y *= f
+    return x, y
 
 
-def _lattice(n, r):
-    """Cubic lattice points at spacing 2r whose distance to the boundary is at least r."""
-    lo, hi = ((-1.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_dist(xs, ys, n, cell):
+    inv = 1.0 / cell
+    cells = {}
+    for i in range(n):
+        k = (int(xs[i] * inv), int(ys[i] * inv))
+        b = cells.get(k)
+        if b is None:
+            cells[k] = [i]
+        else:
+            b.append(i)
+    best = 99.0
+    for i in range(n):
+        xi = xs[i]
+        yi = ys[i]
+        cx = int(xi * inv)
+        cy = int(yi * inv)
+        for ox in (-1, 0, 1):
+            for oy in (-1, 0, 1):
+                lst = cells.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
+                    dd = dx * dx + dy * dy
+                    if dd < best:
+                        best = dd
+    if best > 98.0:
+        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
+                dd = dx * dx + dy * dy
+                if dd < best:
+                    best = dd
+    return math.sqrt(best)
 
 
+def _min_pair(xs, ys, n, cell):
+    inv = 1.0 / cell
+    cells = {}
+    for i in range(n):
+        k = (int(xs[i] * inv), int(ys[i] * inv))
+        b = cells.get(k)
+        if b is None:
+            cells[k] = [i]
+        else:
+            b.append(i)
+    best = 99.0
+    bi, bj = 0, 1
+    for i in range(n):
+        xi = xs[i]
+        yi = ys[i]
+        cx = int(xi * inv)
+        cy = int(yi * inv)
+        for ox in (-1, 0, 1):
+            for oy in (-1, 0, 1):
+                lst = cells.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
+                    dd = dx * dx + dy * dy
+                    if dd < best:
+                        best = dd
+                        bi, bj = i, j
+    if best > 98.0:
+        for i in range(n - 1):
+            for j in range(i + 1, n):
+                dx = xs[j] - xs[i]
+                dy = ys[j] - ys[i]
+                dd = dx * dx + dy * dy
+                if dd < best:
+                    best = dd
+                    bi, bj = i, j
+    return bi, bj
+
+
+def _separate(xs, ys, n, d, sweeps, rng):
+    d2 = d * d
+    inv = 1.0 / d
+    for _ in range(sweeps):
+        cells = {}
+        for i in range(n):
+            k = (int(xs[i] * inv), int(ys[i] * inv))
+            b = cells.get(k)
+            if b is None:
+                cells[k] = [i]
+            else:
+                b.append(i)
+        worst = 0.0
+        for i in range(n):
+            xi = xs[i]
+            yi = ys[i]
+            cx = int(xi * inv)
+            cy = int(yi * inv)
+            for ox in (-1, 0, 1):
+                for oy in (-1, 0, 1):
+                    lst = cells.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
+                        dd = dx * dx + dy * dy
+                        if dd >= d2:
+                            continue
+                        if dd < 1e-20:
+                            a = rng.random() * 6.283185307179586
+                            dx = math.cos(a) * 1e-7
+                            dy = math.sin(a) * 1e-7
+                            dd = 1e-14
+                        dist = math.sqrt(dd)
+                        gap = d - 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)
+                        xs[j] = px
+                        ys[j] = py
+                        xi, yi = _proj(xi - ux, yi - uy)
+            xs[i] = xi
+            ys[i] = yi
+        if worst <= 1e-15:
+            return True
+    return False
+
+
+def _ratchet(xs, ys, n, m, t_end, rng, step, sweeps, floor):
+    while step > floor:
+        if time.perf_counter() >= t_end:
+            break
+        d = m * (1.0 + step)
+        sx = xs[:]
+        sy = ys[:]
+        _separate(xs, ys, n, d, sweeps, rng)
+        m2 = _min_dist(xs, ys, n, d)
+        if m2 > m:
+            m = m2
+            step *= 1.25
+        else:
+            xs[:] = sx
+            ys[:] = sy
+            step *= 0.55
+    return m
+
+
+def _rand_pt(rng):
+    while True:
+        x = rng.uniform(-1.0, 1.0)
+        y = rng.uniform(-1.0, 1.0)
+        if _gauge(x, y) <= 1.0:
+            return x, y
+
+
+def _hole(xs, ys, n, skip, m, rng, tries):
+    inv = 1.0 / m
+    cells = {}
+    for i in range(n):
+        if i == skip:
+            continue
+        k = (int(xs[i] * inv), int(ys[i] * inv))
+        b = cells.get(k)
+        if b is None:
+            cells[k] = [i]
+        else:
+            b.append(i)
+    bx = by = 0.0
+    bd = -1.0
+    rr = (-2, -1, 0, 1, 2)
+    for _ in range(tries):
+        px, py = _rand_pt(rng)
+        cx = int(px * inv)
+        cy = int(py * inv)
+        d = 99.0
+        for ox in rr:
+            for oy in rr:
+                lst = cells.get((cx + ox, cy + oy))
+                if lst is None:
+                    continue
+                for j in lst:
+                    dx = xs[j] - px
+                    dy = ys[j] - py
+                    dd = dx * dx + dy * dy
+                    if dd < d:
+                        d = dd
+        if d > bd:
+            bd = d
+            bx = px
+            by = py
+    return bx, by
+
+
+def _crop_lattices(n, keep, t_end, nang, noff):
+    """Undistorted lattices cropped by the container gauge to the tightest scaled copy.
+
+    gauge(p) >= |p| always (the gauge equals the radius only where p points at a vertex),
+    so radius-filtering the raw lattice to the estimated crop scale keeps every point that
+    can make the cut.  That is what makes a fine angle/offset sweep affordable.
+    """
+    out = []
+    seen = set()
+    area_pp = (0.8660254037844386, 1.0)
+    a1 = math.pi if not SIDES else 0.5 * SIDES * math.sin(2.0 * math.pi / SIDES)
+    for bi, (u1, u2, v1, v2, maxrad) in enumerate((
+            (1.0, 0.0, 0.5, 0.8660254037844386, 1.0471975511965976),
+            (1.0, 0.0, 0.0, 1.0, 1.5707963267948966))):
+        g_est = math.sqrt(n * area_pp[bi] / a1)
+        rad = g_est * 1.18 + 1.7
+        ra = int(rad * 2.4) + 3
+        raw = []
+        for a in range(-ra, ra + 1):
+            for b in range(-ra, ra + 1):
+                x = a * u1 + b * v1
+                y = a * u2 + b * v2
+                if x * x + y * y <= rad * rad:
+                    raw.append((x, y))
+        if len(raw) < n:
+            continue
+        for ai in range(nang):
+            if time.perf_counter() >= t_end:
+                break
+            th = maxrad * ai / nang
+            c = math.cos(th)
+            s = math.sin(th)
+            rot = [(x * c - y * s, x * s + y * c) for x, y in raw]
+            for oi in range(noff):
+                for oj in range(noff):
+                    dx = (oi / noff) * u1 + (oj / noff) * v1
+                    dy = (oi / noff) * u2 + (oj / noff) * v2
+                    sx = dx * c - dy * s
+                    sy = dx * s + dy * c
+                    keys = []
+                    ap = keys.append
+                    for x, y in rot:
+                        px = x + sx
+                        py = y + sy
+                        ap((_gauge(px, py), px, py))
+                    if len(keys) < n:
+                        continue
+                    keys.sort()
+                    g = keys[n - 1][0]
+                    if g <= 1e-9:
+                        continue
+                    m = 1.0 / g
+                    key = round(m, 10)
+                    if key in seen:
+                        continue
+                    seen.add(key)
+                    out.append((m, [(p[1] * m, p[2] * m) for p in keys[:n]]))
+    out.sort(key=lambda t: -t[0])
+    return out[:keep]
+
+
+def _rings(n, rng):
+    """Concentric shells: the natural competitor to a lattice in a round container."""
+    out = []
+    for outer in range(max(3, int(math.sqrt(n))), int(3.4 * math.sqrt(n)) + 2):
+        if outer > n:
+            break
+        pts = []
+        left = n - outer
+        for k in range(outer):
+            a = 2.0 * math.pi * k / outer
+            pts.append((math.cos(a) * RHO, math.sin(a) * RHO))
+        ring = 1
+        while left > 0:
+            cnt = max(1, int(outer * max(0.0, 1.0 - ring * 1.0 / (0.5 * outer))))
+            cnt = min(cnt, left)
+            rad = RHO * max(0.0, 1.0 - ring * 1.0 / (0.5 * outer))
+            if cnt == 1 or rad <= 1e-6:
+                pts.append((0.0, 0.0))
+                left -= 1
+                break
+            off = rng.random() * 2.0 * math.pi
+            for k in range(cnt):
+                a = off + 2.0 * math.pi * k / cnt
+                pts.append((math.cos(a) * rad, math.sin(a) * rad))
+            left -= cnt
+            ring += 1
+        pts = [_proj(x, y) for x, y in pts[:n]]
+        while len(pts) < n:
+            pts.append(_rand_pt(rng))
+        out.append(pts)
+    return out
+
+
+def _starts(n, keep, rng, t_end, nang, noff):
+    cell = max(1e-3, 1.4 / math.sqrt(n))
+    cands = []
+    seen = set()
+    for m, pts in _crop_lattices(n, keep, t_end, nang, noff):
+        key = round(m, 9)
+        if key not in seen:
+            seen.add(key)
+            cands.append((m, [p[0] for p in pts], [p[1] for p in pts]))
+    for pts in _rings(n, rng):
+        xs = [p[0] for p in pts]
+        ys = [p[1] for p in pts]
+        m = _min_dist(xs, ys, n, cell)
+        key = round(m, 9)
+        if key not in seen:
+            seen.add(key)
+            cands.append((m, xs, ys))
+    cands.sort(key=lambda t: -t[0])
+    return [(c[1], c[2]) for c in cands[:keep]]
+
+
+def _to_centres(xs, ys, n, m):
+    r = m * RHO / (2.0 * RHO + m)
+    s = (RHO - r) / RHO
+    return [(xs[i] * s, ys[i] * s) 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)
+
+    if n == 1:
+        return [(0.0, 0.0)]
+
+    sweeps = 10 if n > 90 else 16
+    tries = 300
+    cell0 = max(1e-3, 1.4 / math.sqrt(n))
+
+    nang = 30 if n >= 60 else 18
+    noff = 5 if n >= 60 else 4
+    lats = _starts(n, 8, rng, t0 + time_budget * 0.16, nang, noff)
+
+    best_xs = best_ys = None
+    best_m = -1.0
+    starts = 0
+    t_multi = t0 + time_budget * 0.45
+    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.08))
+        if starts < len(lats):
+            xs = lats[starts][0][:]
+            ys = lats[starts][1][:]
+        elif lats and starts % 3:
+            base = lats[starts % len(lats)]
+            xs = base[0][:]
+            ys = base[1][:]
+            j = 0.60 / math.sqrt(n)
+            for i in range(n):
+                xs[i], ys[i] = _proj(xs[i] + rng.uniform(-j, j), ys[i] + rng.uniform(-j, j))
         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]]
+            pts = [_rand_pt(rng) for _ in range(n)]
+            xs = [p[0] for p in pts]
+            ys = [p[1] for p in pts]
+        m = _min_dist(xs, ys, n, cell0)
+        if m < 1e-9:
+            m = 1e-6
+        m = _ratchet(xs, ys, n, m, slice_end, rng, 0.02, sweeps, 1e-7)
+        if m > best_m:
+            best_m, best_xs, best_ys = m, xs[:], ys[:]
+        starts += 1
 
+    t_moves = t0 + time_budget * 0.80
+    move_cap = time_budget * 0.05
+    while time.perf_counter() < t_moves:
+        xs = best_xs[:]
+        ys = best_ys[:]
+        if rng.random() < 0.65:
+            i, j = _min_pair(xs, ys, n, best_m)
+            k = i if rng.random() < 0.5 else j
+            hx, hy = _hole(xs, ys, n, k, best_m, rng, tries)
+            xs[k] = hx
+            ys[k] = hy
+        else:
+            amp = best_m * rng.uniform(0.05, 0.35)
+            for _ in range(max(1, int(n * rng.uniform(0.03, 0.20)))):
+                i = rng.randrange(n)
+                xs[i], ys[i] = _proj(xs[i] + rng.gauss(0.0, amp), ys[i] + rng.gauss(0.0, amp))
+        m = _min_dist(xs, ys, n, best_m)
+        if m < 1e-9:
+            continue
+        m = _ratchet(xs, ys, n, m,
+                     min(t_moves, time.perf_counter() + move_cap),
+                     rng, 0.004, sweeps, 1e-7)
+        if m > best_m:
+            best_m, best_xs, best_ys = m, xs[:], ys[:]
+
+    best_m = _ratchet(best_xs, best_ys, n, best_m, t_end, rng, 1e-4, sweeps, 1e-13)
+    return _to_centres(best_xs, best_ys, n, best_m)
+