Undistorted rotated/offset lattices cropped to the tightest enclosing square, as multistart seeds for inflate-and-separate; fixed the structurally stuck n (120: 0.965 -> 0.99998).

+"""Equal circles in the unit square -- v6.
 
+Reduction: maximise the minimum pairwise distance m of n points in the unit square, then
+r = m / (2(1+m)) and centre_i = r + q_i (1 - 2r).
+
+v6 adds the initialiser that mattered.  v3 built hexagonal starts by forcing rows and
+columns to span [0,1] exactly, which distorts the lattice; the good starts are instead an
+undistorted triangular (or square) lattice, rotated and offset, cropped to the tightest
+axis-aligned square holding n of its points.  For n=120 that start alone is within 0.07%
+of the record, where the row-based family relaxed to 3.5% short.
+
+Search: multistart over the best few cropped lattices, then reinsertion / shake hill
+climbing on the incumbent, then one fine polish.
+"""
+
 import math
+import random
+import time
 
-def _radius(pts):
-    r = min(min(x, 1 - x, y, 1 - y) for x, y in pts)
-    n = len(pts)
+_S3 = 0.8660254037844386
+
+
+def _build(xs, ys, n, cell):
+    inv = 1.0 / cell
+    cells = {}
     for i in range(n):
-        xi, yi = pts[i]
-        for j in range(i + 1, n):
-            d = math.hypot(xi - pts[j][0], yi - pts[j][1]) / 2
-            if d < r:
-                r = d
-    return r
+        k = (int(xs[i] * inv), int(ys[i] * inv))
+        b = cells.get(k)
+        if b is None:
+            cells[k] = [i]
+        else:
+            b.append(i)
+    return cells, inv
 
-def _grid(n):
-    best = None
-    for rows in range(1, n + 1):
-        cols = math.ceil(n / rows)
-        m = max(rows, cols)
-        if best is None or m < best[0]:
-            best = (m, rows, cols)
-    m, rows, cols = best
-    step = 1.0 / m
-    return [((j + 0.5) * step, (i + 0.5) * step) for i in range(rows) for j in range(cols)][:n]
 
-def _hex(n):
-    best = None
-    for rows in range(1, n + 1):
-        for cols in range(1, n + 1):
-            count = rows * cols - rows // 2
-            if count < n:
-                continue
-            r = min(1.0 / (2 * cols), 1.0 / (2 + (rows - 1) * math.sqrt(3)))
-            if best is None or r > best[0]:
-                best = (r, rows, cols)
-    r, rows, cols = best
-    pts = []
-    for i in range(rows):
-        k = cols - (i % 2)
-        for j in range(k):
-            pts.append((r + (i % 2) * r + 2 * r * j, r + i * r * math.sqrt(3)))
-    return pts[:n]
+def _min_dist(xs, ys, n, cell):
+    cells, inv = _build(xs, ys, n, cell)
+    best = 9.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 > 8.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):
+    cells, inv = _build(xs, ys, n, cell)
+    best = 9.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 > 8.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
+                        xj = xs[j] + ux
+                        yj = ys[j] + uy
+                        xs[j] = 0.0 if xj < 0.0 else (1.0 if xj > 1.0 else xj)
+                        ys[j] = 0.0 if yj < 0.0 else (1.0 if yj > 1.0 else yj)
+                        xi -= ux
+                        yi -= uy
+                        if xi < 0.0:
+                            xi = 0.0
+                        elif xi > 1.0:
+                            xi = 1.0
+                        if yi < 0.0:
+                            yi = 0.0
+                        elif yi > 1.0:
+                            yi = 1.0
+            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 _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.5
+    bd = -1.0
+    rr = (-2, -1, 0, 1, 2)
+    for _ in range(tries):
+        px = rng.random()
+        py = rng.random()
+        cx = int(px * inv)
+        cy = int(py * inv)
+        d = 9.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):
+    """Undistorted lattices, rotated and offset, cropped to the tightest square.
+
+    For each basis, angle and offset the n points with the smallest max(|x|,|y|) are kept
+    and rescaled into the unit square, so the start has min distance exactly 1/L with L
+    the crop width.  Ranked by that, deduplicated, best `keep` returned.
+    """
+    out = []
+    seen = set()
+    for u1, u2, v1, v2, maxrad in (
+            (1.0, 0.0, 0.5, _S3, 0.5235987755982988),      # triangular, 30 deg domain
+            (1.0, 0.0, 0.0, 1.0, 0.7853981633974483)):     # square, 45 deg domain
+        area = abs(u1 * v2 - u2 * v1)
+        half = math.sqrt(n * area) * 0.5
+        lim = half * 1.55
+        rad = lim * 1.4143 + 1.0
+        ra = int(rad * 2.3) + 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
+        nang = 16
+        for ai in range(nang):
+            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(3):
+                for oj in range(3):
+                    dx = (oi / 3.0) * u1 + (oj / 3.0) * v1
+                    dy = (oi / 3.0) * u2 + (oj / 3.0) * 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
+                        k = px if px >= 0.0 else -px
+                        t = py if py >= 0.0 else -py
+                        if t > k:
+                            k = t
+                        if k <= lim:
+                            ap((k, px, py))
+                    if len(keys) < n:
+                        continue
+                    keys.sort()
+                    h = keys[n - 1][0]
+                    if h <= 1e-9:
+                        continue
+                    L = 2.0 * h
+                    m = 1.0 / L
+                    key = round(m, 10)
+                    if key in seen:
+                        continue
+                    seen.add(key)
+                    out.append((m, [((p[1] + h) / L, (p[2] + h) / L) for p in keys[:n]]))
+    out.sort(key=lambda t: -t[0])
+    return out[:keep]
+
+
+def _row_lattices(n):
+    """Hexagonal row splits forced to span the square; still best for some small n."""
+    out = []
+    for rows in range(2, int(2.2 * math.sqrt(n)) + 3):
+        dy = 1.0 / (rows - 1)
+        for longfirst in (True, False):
+            nlong = (rows + 1) // 2 if longfirst else rows // 2
+            nshort = rows - nlong
+            c = 2
+            while nlong * c + nshort * (c - 1) < n:
+                c += 1
+            for extra in (0, 1):
+                cc = c + extra
+                pts = []
+                for i in range(rows):
+                    islong = (i % 2 == 0) if longfirst else (i % 2 == 1)
+                    k = cc if islong else cc - 1
+                    y = i * dy
+                    for j in range(k):
+                        x = j / (cc - 1) if islong else (j + 0.5) / (cc - 1)
+                        pts.append((x, y))
+                if len(pts) >= n:
+                    out.append(pts)
+                    out.append([(b, a) for a, b in pts])
+    side = int(math.ceil(math.sqrt(n)))
+    if side > 1:
+        out.append([(j / (side - 1), i / (side - 1))
+                    for i in range(side) for j in range(side)])
+    return out
+
+
+def _starts(n, keep):
+    cell = max(1e-3, 0.7 / math.sqrt(n))
+    cands = []
+    seen = set()
+    for m, pts in _crop_lattices(n, keep):
+        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 _row_lattices(n):
+        p = pts[:n]
+        xs = [q[0] for q in p]
+        ys = [q[1] for q in p]
+        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 / (2.0 * (1.0 + m))
+    s = 1.0 - 2.0 * r
+    return [(r + xs[i] * s, r + ys[i] * s) for i in range(n)]
+
+
 def pack(n, time_budget, seed):
-    cands = [_grid(n), _hex(n)]
-    return max(cands, key=_radius)
+    t0 = time.perf_counter()
+    t_end = t0 + time_budget * 0.90
+    rng = random.Random((seed & 0xFFFFFFFF) * 1000003 + n)
 
+    if n == 1:
+        return [(0.5, 0.5)]
+    if n == 2:
+        return [(0.25, 0.25), (0.75, 0.75)]
+
+    sweeps = 10 if n > 90 else 16
+    tries = 300
+    cell0 = max(1e-3, 0.7 / math.sqrt(n))
+
+    lats = _starts(n, 8)
+
+    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.30 / math.sqrt(n)
+            for i in range(n):
+                x = xs[i] + rng.uniform(-j, j)
+                y = ys[i] + rng.uniform(-j, j)
+                xs[i] = 0.0 if x < 0.0 else (1.0 if x > 1.0 else x)
+                ys[i] = 0.0 if y < 0.0 else (1.0 if y > 1.0 else y)
+        else:
+            xs = [rng.random() for _ in range(n)]
+            ys = [rng.random() for _ in range(n)]
+        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)
+                x = xs[i] + rng.gauss(0.0, amp)
+                y = ys[i] + rng.gauss(0.0, amp)
+                xs[i] = 0.0 if x < 0.0 else (1.0 if x > 1.0 else x)
+                ys[i] = 0.0 if y < 0.0 else (1.0 if y > 1.0 else y)
+        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)
+