{"id":"tammes-problem","name":"Tammes problem, maximum minimum distance on the sphere","family":"spherical-codes","description":"Place n points on the unit sphere to maximise the smallest pairwise distance (equivalently pack n equal spherical caps), for twelve n between 15 and 100. Scored against Sloane's spherical-code records; none of these n is proven optimal.","metric":"record_ratio","direction":"maximize","tolerance":0.05,"eval_timeout_seconds":300,"agent_timeout_seconds":1800,"mutable":["sphere.py"],"runtime":"python>=3.11, standard library only (math, random, itertools, functools, collections, heapq, time)","decomposable":true,"status":"active","captain":null,"parent_problem":null,"program_md":"# Tammes problem, maximum minimum distance on the sphere\n\n## Goal\n\n`sphere.py` exposes `place(n: int, time_budget: float, seed: int) -> list[tuple[float, float, float]]`:\n`n` points `(x, y, z)` on the unit sphere. Maximise the smallest pairwise Euclidean distance\n\n    D = min over i < j of |x_i - x_j|\n\nEquivalently: pack `n` equal spherical caps as large as possible, or find the best spherical code\nof size `n`. This is Problem 34 in DeepMind's AlphaEvolve repository of problems (Section 6.18 of\n\"Mathematical Exploration and Discovery at Scale\", arXiv:2511.02864). AlphaEvolve matched the\nrecords for n = 3, 7, 12 and 25 and fell slightly short at n = 32, 50, 100 and 200 (Table 5 of the\npaper), so the larger n here are where a good pure-Python search can still make a difference.\nOptimality is proven only for n <= 14 and n = 24; every n in this benchmark is open.\n\n## Metric\n\n    metric = mean over n in NS of  D(n) / record(n)\n\n1.0 means matching every record; above 1.0 means beating at least one. The eval projects every\nreturned point onto the unit sphere (any finite non-zero vector is accepted; a zero vector, a\nnon-finite coordinate or the wrong count is a failed run; coincident points simply score 0 for that\nn) and recomputes the minimum distance itself. Nothing your solver reports is used. `ZT_EVAL_SEED`\nonly changes the `seed` handed to `place`, so your method must be robust to its starting point.\n\n## Records\n\nBest-known configurations from N. J. A. Sloane, R. H. Hardin, W. D. Smith et al., \"Tables of\nSpherical Codes\" (NeilSloane.com/packings/, files `dim3/pack.3.<n>.txt`, fetched 2026-09-07). The\nrecord distance below was recomputed from those coordinates; each matches the table's angular\nseparation to its full 1e-7 degree precision. A distance above a record by more than 1e-8 is listed\nin `records_beaten`.\n\n| n | record D | angle (deg) | source | proven optimal |\n|---|---|---|---|---|\n| 15 | 0.902656188015 | 53.6578501 | Hardin, Sloane, Smith 1994 | no |\n| 17 | 0.862444879257 | 51.0903285 | Hardin, Sloane, Smith 1994 | no |\n| 19 | 0.808558114565 | 47.6919141 | Hardin, Sloane, Smith 1994 | no |\n| 21 | 0.775243921143 | 45.6132231 | Kottwitz, Acta Cryst. A47 (1991) | no |\n| 25 | 0.710776154955 | 41.6344612 | Hardin, Sloane, Smith 1994 | no |\n| 27 | 0.695141408884 | 40.6776007 | Kottwitz, Acta Cryst. A47 (1991) | no |\n| 32 | 0.642469275564 | 37.4752140 | Hardin, Sloane, Smith 1994 | no |\n| 33 | 0.622257802439 | 36.2545530 | Kottwitz, Acta Cryst. A47 (1991) | no |\n| 50 | 0.513472084621 | 29.7529564 | Hardin, Sloane, Smith 1994 | no |\n| 54 | 0.495975188171 | 28.7169205 | Kottwitz, Acta Cryst. A47 (1991) | no |\n| 64 | 0.453898297814 | 26.2350433 | Hardin, Sloane, Smith 1994 | no |\n| 100 | 0.365006496096 | 21.0312020 | Hardin, Sloane, Smith 1994 | no |\n\n## Constraints\n\n- Standard library only. No numpy, no scipy. The eval rejects other imports.\n- Respect `time_budget` (seconds, per call). The eval fails a call that runs past 1.25x + 3 s.\n- Deterministic given `seed`: use `random.Random(seed)`.\n\n## Iterating\n\n- `ZT_EVAL_NS=15,25` runs a subset of n; `ZT_EVAL_PER_N_SECONDS=2` shortens the per-n budget. The\n  default is all twelve n at 8 s each (about 100 s). Run `python eval.py` in your workspace.\n- The per-n detail in the eval output shows which n are furthest from their record.\n\n## Ideas that are known to matter (check the journal before repeating one)\n\n- The objective is a max-min, so plain gradient methods stall. Standard tricks: minimise a soft\n  energy sum d^-p with p ramped up (the baseline), or iterate \"find the active contact pairs, then\n  solve the LP that pushes them apart\" (an SLP), which converges to the true local optimum.\n- The record configurations have very specific contact graphs (Kottwitz's improvements at n = 21,\n  27, 33, 54 came from symmetric constructions). Search over starting symmetries, not just random\n  starts, and finish each candidate with a contact-graph polish.\n- Basin hopping: perturb the current best a little, re-optimise, accept if the minimum distance\n  did not drop. Many restarts beat one long run.\n- The pure-Python pair loop is the bottleneck: only the near-contact pairs matter once the\n  configuration is decent, so keep a neighbour list and refresh it occasionally.\n\nWrite one honest line in `NOTES.md`: the idea, and which `n` it helped.\n\nSimpler is better: all else equal prefer the shorter solver, and treat removing code for an\nequal score as a win. Log every experiment, including discards, in your results.tsv.\n","eval_py":"\"\"\"Eval for tammes-problem. Prints one JSON line: {\"metric\": record_ratio, ...}.\n\nThe solver's points are normalised onto the unit sphere and the minimum pairwise distance is\nrecomputed here; nothing the solver reports is trusted.\n\nEnv:\n  ZT_EVAL_SEED             seed handed to place()\n  ZT_EVAL_PER_N_SECONDS    time budget per n (default 8)\n  ZT_EVAL_NS               comma-separated n values (default \"15,17,19,21,25,27,32,33,50,54,64,100\")\n\"\"\"\n\nfrom __future__ import annotations\n\nimport ast\nimport json\nimport math\nimport os\nimport random\nimport sys\nimport time\nfrom pathlib import Path\n\nSEED = os.environ.get(\"ZT_EVAL_SEED\", \"dev-seed\")\nBUDGET = float(os.environ.get(\"ZT_EVAL_PER_N_SECONDS\", \"8\"))\nNS = [int(x) for x in os.environ.get(\"ZT_EVAL_NS\", \"15,17,19,21,25,27,32,33,50,54,64,100\").split(\",\")]\nSTDLIB_ALLOW = {\"math\", \"random\", \"itertools\", \"functools\", \"collections\", \"heapq\", \"time\", \"sys\", \"typing\", \"operator\"}\nEPS = 1e-9\nFORBIDDEN_NAMES = {\"__import__\", \"importlib\", \"builtins\", \"__builtins__\", \"open\", \"exec\", \"eval\", \"compile\",\n                   \"globals\", \"__loader__\", \"__spec__\", \"breakpoint\", \"input\", \"memoryview\", \"vars\"}\n\n# Best-known minimum Euclidean distance, computed (2026-09-07) from the coordinate files\n# pack.3.<n>.txt of N. J. A. Sloane, R. H. Hardin, W. D. Smith et al., \"Tables of Spherical Codes\",\n# NeilSloane.com/packings/. Each agrees with the table's angular separation to its full 1e-7 degree\n# precision (the angle is in the comment). n = 21, 27, 33, 54 are D. A. Kottwitz's improvements\n# (Acta Cryst. A47, 158-165 (1991)). None of these n is proven optimal (proofs cover n <= 14 and\n# n = 24). A distance more than BEAT_TOL above a record is flagged in records_beaten.\nRECORDS = {\n    15: 0.902656188015,   # 53.6578501 deg\n    17: 0.862444879257,   # 51.0903285 deg\n    19: 0.808558114565,   # 47.6919141 deg\n    21: 0.775243921143,   # 45.6132231 deg  (Kottwitz)\n    25: 0.710776154955,   # 41.6344612 deg\n    27: 0.695141408884,   # 40.6776007 deg  (Kottwitz)\n    32: 0.642469275564,   # 37.4752140 deg\n    33: 0.622257802439,   # 36.2545530 deg  (Kottwitz)\n    50: 0.513472084621,   # 29.7529564 deg\n    54: 0.495975188171,   # 28.7169205 deg  (Kottwitz)\n    64: 0.453898297814,   # 26.2350433 deg\n    100: 0.365006496096,  # 21.0312020 deg\n}\nBEAT_TOL = 1e-8\n\n\ndef fail(msg: str, kind: str = \"error\") -> None:\n    print(json.dumps({\"metric\": 0.0, \"error\": msg, \"kind\": kind}))\n    sys.exit(1)\n\n\ndef check_imports(path: Path) -> None:\n    try:\n        tree = ast.parse(path.read_text(encoding=\"utf-8\"))\n    except SyntaxError as e:\n        fail(f\"syntax error in pack.py: {e}\", \"compile_error\")\n    for node in ast.walk(tree):\n        names = []\n        if isinstance(node, ast.Import):\n            names = [a.name.split(\".\")[0] for a in node.names]\n        elif isinstance(node, ast.ImportFrom) and node.module:\n            names = [node.module.split(\".\")[0]]\n        for nm in names:\n            if nm not in STDLIB_ALLOW:\n                fail(f\"import of '{nm}' is not allowed (stdlib subset only: {sorted(STDLIB_ALLOW)})\", \"compile_error\")\n        # dynamic imports and raw file/process access are not part of the problem either\n        ident = node.id if isinstance(node, ast.Name) else node.attr if isinstance(node, ast.Attribute) else None\n        if ident in FORBIDDEN_NAMES:\n            fail(f\"use of '{ident}' is not allowed in a solver\", \"compile_error\")\n        if isinstance(node, ast.ImportFrom) and node.level:\n            fail(\"relative imports are not allowed in a solver\", \"compile_error\")\n\n\ndef unit_points(points: list, n: int) -> list[tuple[float, float, float]]:\n    \"\"\"n finite non-zero vectors, projected onto the unit sphere. Anything else is a wrong answer.\"\"\"\n    if not isinstance(points, (list, tuple)) or len(points) != n:\n        fail(f\"place({n}) must return {n} points\", \"wrong_answer\")\n    out = []\n    for p in points:\n        try:\n            x, y, z = float(p[0]), float(p[1]), float(p[2])\n        except Exception:\n            fail(f\"place({n}) returned a non-point {p!r}\", \"wrong_answer\")\n        if not all(map(math.isfinite, (x, y, z))):\n            fail(f\"place({n}) returned a non-finite coordinate\", \"wrong_answer\")\n        r = math.sqrt(x * x + y * y + z * z)\n        if r < EPS:\n            fail(f\"place({n}) returned a point at the origin, which has no direction\", \"wrong_answer\")\n        out.append((x / r, y / r, z / r))\n    return out\n\n\ndef min_distance(pts: list[tuple[float, float, float]], n: int) -> float:\n    best = 4.0\n    for i in range(n):\n        xi, yi, zi = pts[i]\n        for j in range(i + 1, n):\n            xj, yj, zj = pts[j]\n            d2 = (xi - xj) ** 2 + (yi - yj) ** 2 + (zi - zj) ** 2\n            if d2 < best:\n                best = d2\n    return math.sqrt(best)\n\n\ndef main() -> None:\n    here = Path(__file__).parent\n    check_imports(here / \"sphere.py\")\n    sys.path.insert(0, str(here))\n    try:\n        import sphere as cand  # noqa: E402\n    except SystemExit:\n        raise\n    except Exception as e:\n        fail(f\"import sphere.py failed: {e!r}\", \"compile_error\")\n    if not hasattr(cand, \"place\"):\n        fail(\"sphere.py must define place(n, time_budget, seed)\", \"compile_error\")\n\n    seed_int = random.Random(f\"tammes|{SEED}\").getrandbits(32)\n    per_n, beaten = {}, []\n    for n in NS:\n        if n not in RECORDS:\n            fail(f\"no record for n={n}\", \"error\")\n        t0 = time.perf_counter()\n        try:\n            points = cand.place(n, BUDGET, seed_int)\n        except SystemExit:\n            raise\n        except Exception as e:\n            fail(f\"place({n}) raised {e!r}\", \"runtime_error\")\n        elapsed = time.perf_counter() - t0\n        if elapsed > 1.25 * BUDGET + 3:\n            fail(f\"place({n}) took {elapsed:.1f}s against a {BUDGET:.0f}s budget\", \"timeout\")\n        d = min_distance(unit_points(points, n), n)\n        per_n[n] = {\"min_dist\": round(d, 12), \"record\": RECORDS[n], \"ratio\": round(d / RECORDS[n], 9),\n                    \"seconds\": round(elapsed, 2)}\n        if d > RECORDS[n] + BEAT_TOL:\n            beaten.append(n)\n    metric = sum(v[\"ratio\"] for v in per_n.values()) / len(per_n)\n    print(json.dumps({\"metric\": round(metric, 9), \"per_n\": per_n, \"records_beaten\": beaten}))\n\n\nif __name__ == \"__main__\":\n    main()\n","baseline":{"sphere.py":"\"\"\"Baseline: Fibonacci-sphere start, then projected gradient descent on the soft-min energy\nsum (s/d_ij)^p with p stepped up 8 -> 64, keeping the best minimum distance seen. Lands a few\npercent below the records. Beat it.\"\"\"\n\nimport math\nimport random\nimport time\n\n\ndef _fibonacci(n: int) -> list[list[float]]:\n    golden = math.pi * (3.0 - math.sqrt(5.0))\n    pts = []\n    for i in range(n):\n        z = 1.0 - (2.0 * i + 1.0) / n\n        r = math.sqrt(max(0.0, 1.0 - z * z))\n        pts.append([r * math.cos(golden * i), r * math.sin(golden * i), z])\n    return pts\n\n\ndef _min_dist(pts: list[list[float]]) -> float:\n    best = 4.0\n    for i in range(len(pts)):\n        xi, yi, zi = pts[i]\n        for j in range(i + 1, len(pts)):\n            xj, yj, zj = pts[j]\n            d2 = (xi - xj) ** 2 + (yi - yj) ** 2 + (zi - zj) ** 2\n            if d2 < best:\n                best = d2\n    return math.sqrt(best)\n\n\ndef _energy_and_forces(pts: list[list[float]], p: float, s: float) -> tuple[float, list[list[float]]]:\n    \"\"\"E = sum (s/d)^p; the force on i from j is p (s/d)^p / d^2 * (x_i - x_j).\"\"\"\n    n = len(pts)\n    e = 0.0\n    f = [[0.0, 0.0, 0.0] for _ in range(n)]\n    for i in range(n):\n        xi, yi, zi = pts[i]\n        fi = f[i]\n        for j in range(i + 1, n):\n            xj, yj, zj = pts[j]\n            dx, dy, dz = xi - xj, yi - yj, zi - zj\n            d2 = dx * dx + dy * dy + dz * dz\n            t = (s * s / d2) ** (p / 2)\n            e += t\n            g = p * t / d2\n            fi[0] += dx * g; fi[1] += dy * g; fi[2] += dz * g\n            fj = f[j]\n            fj[0] -= dx * g; fj[1] -= dy * g; fj[2] -= dz * g\n    return e, f\n\n\ndef _step(pts: list[list[float]], f: list[list[float]], lr: float) -> list[list[float]]:\n    out = []\n    for (x, y, z), (fx, fy, fz) in zip(pts, f):\n        rad = fx * x + fy * y + fz * z\n        nx, ny, nz = x + lr * (fx - rad * x), y + lr * (fy - rad * y), z + lr * (fz - rad * z)\n        r = math.sqrt(nx * nx + ny * ny + nz * nz)\n        out.append([nx / r, ny / r, nz / r])\n    return out\n\n\ndef place(n: int, time_budget: float, seed: int) -> list[tuple[float, float, float]]:\n    rng = random.Random(seed)\n    pts = _fibonacci(n)\n    for q in pts:\n        q[0] += rng.gauss(0, 1e-3); q[1] += rng.gauss(0, 1e-3); q[2] += rng.gauss(0, 1e-3)\n        r = math.sqrt(q[0] ** 2 + q[1] ** 2 + q[2] ** 2)\n        q[0] /= r; q[1] /= r; q[2] /= r\n    best, best_d = [list(q) for q in pts], _min_dist(pts)\n    start = time.perf_counter()\n    total = 0.85 * time_budget\n    powers = (8.0, 16.0, 32.0, 64.0)\n    for k, p in enumerate(powers):\n        deadline = start + total * (k + 1) / len(powers)\n        s = _min_dist(pts)  # scale so the largest term is about 1 and nothing overflows\n        e, f = _energy_and_forces(pts, p, s)\n        lr = 0.05 * s / max(1.0, max(math.sqrt(fx * fx + fy * fy + fz * fz) for fx, fy, fz in f))\n        while time.perf_counter() < deadline:\n            trial = _step(pts, f, lr)\n            e2, f2 = _energy_and_forces(trial, p, s)\n            if e2 < e:\n                pts, e, f = trial, e2, f2\n                lr *= 1.2\n                d = _min_dist(pts)\n                if d > best_d:\n                    best, best_d = [list(q) for q in pts], d\n            else:\n                lr *= 0.5\n                if lr < 1e-14:\n                    break\n    return [tuple(q) for q in best]\n"}}