{"id":"sphere-max-volume","name":"Points on a sphere, maximum convex-hull volume","family":"spherical-codes","description":"Place n points on the unit sphere so that the volume of their convex hull is as large as possible (the Fejes Toth problem), for twelve n between 9 and 50. Scored against Sloane's maximal-volume 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":"# Points on a sphere, maximum convex-hull volume\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 volume of their convex hull, i.e. find the\nlargest polyhedron with `n` vertices inscribed in the unit sphere (Fejes Toth's 1964 problem).\n\nThis is Problem 41 in DeepMind's AlphaEvolve repository of problems (Section 6.24 of \"Mathematical\nExploration and Discovery at Scale\", arXiv:2511.02864). AlphaEvolve matched the first ~60 entries\nof Sloane's table to all 13 printed digits and improved none of them. Optimality is proven only for\nn <= 8 (Berman and Hanes 1970); Mutoh (2003) found numerical candidates for n <= 30 and Sloane's\ntable extends to n = 130. Every n in this benchmark is open.\n\n## Metric\n\n    metric = mean over n in NS of  V(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 or coplanar points just lose\nvolume) and computes the hull volume itself with an incremental convex hull, summing the origin\ntetrahedra over the hull's faces. Nothing your solver reports is used. `ZT_EVAL_SEED` only changes\nthe `seed` handed to `place`, so your method must be robust to its starting point.\n\n## Records\n\nBest-known volumes from R. H. Hardin, N. J. A. Sloane and W. D. Smith, \"Maximal Volume Spherical\nCodes\" (NeilSloane.com/maxvolumes/, 1994; fetched 2026-09-07), exactly as printed there. The eval's\nhull code reproduces each from the site's coordinate files (`dim3/maxvol.3.<n>.txt`) to 5e-13. A\nvolume above a record by more than 1e-9 is listed in `records_beaten`.\n\n| n | record volume | source | proven optimal |\n|---|---|---|---|\n| 9 | 2.043750115900 | Hardin, Sloane, Smith 1994 | no |\n| 10 | 2.218711131545 | Hardin, Sloane, Smith 1994 | no |\n| 11 | 2.354634495069 | Hardin, Sloane, Smith 1994 | no |\n| 13 | 2.612834152060 | Hardin, Sloane, Smith 1994 | no |\n| 14 | 2.720977899349 | Hardin, Sloane, Smith 1994 | no |\n| 16 | 2.886455392275 | Hardin, Sloane, Smith 1994 | no |\n| 18 | 3.009613252523 | Hardin, Sloane, Smith 1994 | no |\n| 20 | 3.118538793195 | Hardin, Sloane, Smith 1994 | no |\n| 24 | 3.283995205283 | Hardin, Sloane, Smith 1994 | no |\n| 30 | 3.455125752062 | Hardin, Sloane, Smith 1994 | no |\n| 40 | 3.634130342837 | Hardin, Sloane, Smith 1994 | no |\n| 50 | 3.740940879707 | 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=9,13` 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 volume has a clean gradient: for a triangulated hull, the derivative of the volume with\n  respect to vertex `v` is one sixth of the sum of the cross products of the two other vertices over\n  the faces containing `v` (the \"area vector\" of its link). Build the hull once per step (the eval's\n  incremental hull is a template you may copy), project the gradient onto the tangent plane, and\n  ascend. The baseline only does Coulomb repulsion, which is not the same objective.\n- Berman and Hanes' necessary condition: at an optimum every vertex is the normalised sum of the\n  area vectors of its incident faces, and all faces are triangles. Use it as a stopping test and a\n  polish step.\n- The optimal polyhedra are \"medial\" (vertex degrees 5 and 6 only, like Thomson minima but not the\n  same configurations). Start from a good Thomson or Tammes configuration and ascend the volume.\n- Basin hopping over the combinatorial type: small perturbations that flip an edge, then re-ascend.\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 sphere-max-volume. Prints one JSON line: {\"metric\": record_ratio, ...}.\n\nThe solver's points are normalised onto the unit sphere and the volume of their convex hull is\nrecomputed here with an incremental hull; 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 \"9,10,11,13,14,16,18,20,24,30,40,50\")\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\", \"9,10,11,13,14,16,18,20,24,30,40,50\").split(\",\")]\nSTDLIB_ALLOW = {\"math\", \"random\", \"itertools\", \"functools\", \"collections\", \"heapq\", \"time\", \"sys\", \"typing\", \"operator\"}\nEPS = 1e-9\nHULL_EPS = 1e-10\nFORBIDDEN_NAMES = {\"__import__\", \"importlib\", \"builtins\", \"__builtins__\", \"open\", \"exec\", \"eval\", \"compile\",\n                   \"globals\", \"__loader__\", \"__spec__\", \"breakpoint\", \"input\", \"memoryview\", \"vars\"}\n\n# Largest known convex-hull volume, exactly as printed (12 decimals) in R. H. Hardin, N. J. A. Sloane\n# and W. D. Smith, \"Maximal Volume Spherical Codes\", NeilSloane.com/maxvolumes/ (1994; fetched\n# 2026-09-07). The hull code below reproduces every value from the site's coordinate files\n# (dim3/maxvol.3.<n>.txt) to 5e-13. Optimality is proven only for n <= 8 (Berman and Hanes 1970);\n# none of these n is proven. A volume more than BEAT_TOL above a record is flagged in records_beaten.\nRECORDS = {\n    9: 2.043750115900,\n    10: 2.218711131545,\n    11: 2.354634495069,\n    13: 2.612834152060,\n    14: 2.720977899349,\n    16: 2.886455392275,\n    18: 3.009613252523,\n    20: 3.118538793195,\n    24: 3.283995205283,\n    30: 3.455125752062,\n    40: 3.634130342837,\n    50: 3.740940879707,\n}\nBEAT_TOL = 1e-9\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 _sub(a, b):\n    return (a[0] - b[0], a[1] - b[1], a[2] - b[2])\n\n\ndef _cross(a, b):\n    return (a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0])\n\n\ndef _dot(a, b):\n    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]\n\n\ndef _hull_faces(pts: list[tuple[float, float, float]]) -> set[tuple[int, int, int]] | None:\n    \"\"\"Incremental 3-D convex hull. Returns outward-oriented triangles, or None if the mesh it built\n    is not a closed surface containing every point (a degeneracy the caller retries after rotating).\"\"\"\n    n = len(pts)\n    i0 = 0\n    i1 = max(range(n), key=lambda i: _dot(_sub(pts[i], pts[i0]), _sub(pts[i], pts[i0])))\n    d = _sub(pts[i1], pts[i0])\n    i2 = max(range(n), key=lambda i: _dot(_cross(d, _sub(pts[i], pts[i0])), _cross(d, _sub(pts[i], pts[i0]))))\n    nrm = _cross(d, _sub(pts[i2], pts[i0]))\n    i3 = max(range(n), key=lambda i: abs(_dot(nrm, _sub(pts[i], pts[i0]))))\n    if _dot(nrm, nrm) < HULL_EPS or abs(_dot(nrm, _sub(pts[i3], pts[i0]))) < HULL_EPS:\n        return set()  # flat: zero volume\n    cen = tuple(sum(pts[i][k] for i in (i0, i1, i2, i3)) / 4 for k in range(3))\n\n    def plane(f):\n        a, b, c = f\n        nn = _cross(_sub(pts[b], pts[a]), _sub(pts[c], pts[a]))\n        length = math.sqrt(_dot(nn, nn))\n        nn = (nn[0] / length, nn[1] / length, nn[2] / length)\n        return nn, _dot(nn, pts[a])\n\n    def orient(a, b, c):\n        nn = _cross(_sub(pts[b], pts[a]), _sub(pts[c], pts[a]))\n        return (a, b, c) if _dot(nn, _sub(pts[a], cen)) > 0 else (a, c, b)\n\n    faces = {orient(i0, i1, i2), orient(i0, i1, i3), orient(i0, i2, i3), orient(i1, i2, i3)}\n    planes = {f: plane(f) for f in faces}\n    for p in range(n):\n        if p in (i0, i1, i2, i3):\n            continue\n        P = pts[p]\n        visible = [f for f in faces if _dot(planes[f][0], P) - planes[f][1] > HULL_EPS]\n        if not visible:\n            continue\n        edges: dict[tuple[int, int], int] = {}\n        for (a, b, c) in visible:\n            for e in ((a, b), (b, c), (c, a)):\n                edges[e] = edges.get(e, 0) + 1\n        horizon = [e for e in edges if (e[1], e[0]) not in edges]\n        for f in visible:\n            faces.discard(f)\n            planes.pop(f)\n        for (a, b) in horizon:\n            f = (a, b, p)\n            faces.add(f)\n            planes[f] = plane(f)\n    # closed surface: every directed edge has its reverse exactly once\n    edges = {}\n    for (a, b, c) in faces:\n        for e in ((a, b), (b, c), (c, a)):\n            edges[e] = edges.get(e, 0) + 1\n    if any(v != 1 or (e[1], e[0]) not in edges for e, v in edges.items()):\n        return None\n    # every point is inside or on every face plane\n    for f, (nn, off) in planes.items():\n        for P in pts:\n            if _dot(nn, P) - off > 1e-7:\n                return None\n    return faces\n\n\ndef hull_volume(pts: list[tuple[float, float, float]], n: int) -> float:\n    for attempt in range(4):\n        faces = _hull_faces(pts)\n        if faces is not None:\n            return sum(_dot(pts[a], _cross(pts[b], pts[c])) for (a, b, c) in faces) / 6.0\n        # degenerate mesh: rotate everything by a fixed small angle and try again (volume is invariant)\n        th = 0.37 * (attempt + 1)\n        c, s = math.cos(th), math.sin(th)\n        pts = [(x * c - y * s, x * s + y * c, z) for (x, y, z) in pts]\n        pts = [(x, y * c - z * s, y * s + z * c) for (x, y, z) in pts]\n    fail(f\"place({n}): convex hull construction failed on a degenerate configuration\", \"error\")\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\"maxvol|{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        v = hull_volume(unit_points(points, n), n)\n        per_n[n] = {\"volume\": round(v, 12), \"record\": RECORDS[n], \"ratio\": round(v / RECORDS[n], 9),\n                    \"seconds\": round(elapsed, 2)}\n        if v > 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 Coulomb energy. A\nThomson minimiser is not a volume maximiser, so this sits about 0.5-1% 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 _energy_and_forces(pts: list[list[float]]) -> tuple[float, list[list[float]]]:\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            inv = 1.0 / math.sqrt(d2)\n            e += inv\n            g = inv / 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 p in pts:\n        p[0] += rng.gauss(0, 1e-3); p[1] += rng.gauss(0, 1e-3); p[2] += rng.gauss(0, 1e-3)\n        r = math.sqrt(p[0] ** 2 + p[1] ** 2 + p[2] ** 2)\n        p[0] /= r; p[1] /= r; p[2] /= r\n    deadline = time.perf_counter() + 0.85 * time_budget\n    lr = 0.5 / (n * math.sqrt(n))\n    e, f = _energy_and_forces(pts)\n    while time.perf_counter() < deadline:\n        trial = _step(pts, f, lr)\n        e2, f2 = _energy_and_forces(trial)\n        if e2 < e:\n            pts, e, f = trial, e2, f2\n            lr *= 1.2\n        else:\n            lr *= 0.5\n            if lr < 1e-12:\n                break\n    return [tuple(p) for p in pts]\n"}}