{"id":"thomson-problem","name":"Thomson problem, minimum Coulomb energy on the sphere","family":"spherical-codes","description":"Place n points on the unit sphere to minimise the Coulomb energy sum 1/|x_i - x_j|, for twelve n between 13 and 122. Scored against the Cambridge Cluster Database 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":"# Thomson problem, minimum Coulomb energy 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. Minimise the Coulomb energy\n\n    E = sum over i < j of 1 / |x_i - x_j|\n\nThis is Problem 33 in DeepMind's AlphaEvolve repository of problems (Section 6.18 of \"Mathematical\nExploration and Discovery at Scale\", arXiv:2511.02864). AlphaEvolve ran it for n up to 300 and\nmatched the state of the art to about 1e-8 (its constructions for n = 282, 292, 306 reproduce the\nCambridge Cluster Database energies) without improving any of them. The exact minimum is known only\nfor n = 2, 3, 4, 5, 6 and 12; every n in this benchmark is open.\n\n## Metric\n\n    metric = mean over n in NS of  record(n) / E(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, the wrong count, or two coincident points is a failed run) and recomputes the\nenergy itself. Nothing your solver reports is used. `ZT_EVAL_SEED` only changes the `seed` handed to\n`place`, so your method must be robust to its starting point.\n\n## Records\n\nBest-known energies from the Cambridge Cluster Database table \"Global Minima for the Thomson\nProblem\" (D. J. Wales and S. Ulker, Phys. Rev. B 74, 212101 (2006), with updates from Wayne\nDeeter; https://www-wales.ch.cam.ac.uk/~wales/CCD/Thomson/table.html, fetched 2026-09-07), exactly\nas printed there. An energy below a record by more than 1e-6 is listed in `records_beaten`.\n\n| n | record E | point group | proven optimal |\n|---|---|---|---|\n| 16 | 92.9116553 | T | no |\n| 37 | 560.6188877 | D5h | no |\n| 38 | 593.0385035 | D6d | no |\n| 42 | 732.0781075 | D5h | no |\n| 47 | 927.0592706 | Cs | no |\n| 54 | 1239.3614747 | C2 | no |\n| 59 | 1490.7733352 | C2 | no |\n| 64 | 1765.8025779 | D2 | no |\n| 77 | 2591.8501523 | D5 | no |\n| 88 | 3416.7201967 | D2 | no |\n| 100 | 4448.3506343 | T | no |\n| 122 | 6698.3744992 | Ih | no |\n\nThe n were chosen so that a single gradient descent from a Fibonacci spiral (the baseline) ends\nin the wrong basin on most of them; a Thomson n that plain descent solves gives no signal.\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=16,47` 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 energy landscape has exponentially many local minima; a single gradient descent from a\n  Fibonacci spiral gets stuck a few 1e-4 above the record. Basin hopping (perturb, re-minimise,\n  accept if better) is what found most of the database entries.\n- Good minima are near-triangular lattices with exactly twelve 5-fold defects. The n = 122 record\n  has icosahedral symmetry, n = 37 and 42 are D5h, n = 16 and 100 are T. Symmetry-constrained\n  starts converge faster than random ones.\n- Second-order convergence: after gradient descent, a few Newton or conjugate-gradient steps on the\n  tangent space reach the minimum to 1e-10. The ratio only rounds to 1.0 once you are within about\n  1e-7 relative, so finishing matters.\n- The pure-Python pair loop is the bottleneck: cache differences, avoid function-call overhead, and\n  reuse the force computation as the energy computation.\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 thomson-problem. Prints one JSON line: {\"metric\": record_ratio, ...}.\n\nThe solver's points are normalised onto the unit sphere and the Coulomb energy is recomputed\nhere; 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 \"16,37,38,42,47,54,59,64,77,88,100,122\")\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\", \"16,37,38,42,47,54,59,64,77,88,100,122\").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# Lowest known Coulomb energy E = sum_{i<j} 1/|x_i - x_j|, exactly as printed (7 decimals) in the\n# Cambridge Cluster Database table \"Global Minima for the Thomson Problem\" (D. J. Wales and\n# S. Ulker, Phys. Rev. B 74, 212101 (2006), with updates from W. Deeter), fetched 2026-09-07.\n# None of these n is proven optimal (proofs exist only for n = 2..6 and 12). The set favours n where a\n# plain gradient descent from a Fibonacci spiral lands in a wrong basin, so the score has signal. An\n# energy more than\n# BEAT_TOL below a record is flagged in records_beaten; the table is rounded to 1e-7, so BEAT_TOL\n# sits well above that rounding.\nRECORDS = {\n    16: 92.9116553,     # T\n    37: 560.6188877,    # D5h\n    38: 593.0385035,    # D6d\n    42: 732.0781075,    # D5h\n    47: 927.0592706,    # Cs\n    54: 1239.3614747,   # C2\n    59: 1490.7733352,   # C2\n    64: 1765.8025779,   # D2\n    77: 2591.8501523,   # D5\n    88: 3416.7201967,   # D2\n    100: 4448.3506343,  # T\n    122: 6698.3744992,  # Ih\n}\nBEAT_TOL = 1e-6\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 energy(pts: list[tuple[float, float, float]], n: int) -> float:\n    e = 0.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            d = math.sqrt((xi - xj) ** 2 + (yi - yj) ** 2 + (zi - zj) ** 2)\n            if d < EPS:\n                fail(f\"place({n}): points {i} and {j} coincide (infinite energy)\", \"wrong_answer\")\n            e += 1.0 / d\n    return e\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\"thomson|{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        e = energy(unit_points(points, n), n)\n        per_n[n] = {\"energy\": round(e, 10), \"record\": RECORDS[n], \"ratio\": round(RECORDS[n] / e, 9),\n                    \"seconds\": round(elapsed, 2)}\n        if e < 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 with an\nadaptive step. Lands within about 0.1-1% of 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  # 1/d^3\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  # drop the radial part: move along the sphere\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:  # tiny jitter so symmetric starts do not sit on a saddle\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"}}