{"id":"spherical-codes-4d","name":"Spherical codes in 4 dimensions, maximum minimum distance on S^3","family":"spherical-codes","description":"Place n points on the unit sphere in R^4 to maximise the smallest pairwise distance (the best spherical code of size n), for thirteen n between 15 and 100. Scored against Sloane's tables of putatively optimal packings; 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":"# Spherical codes in 4 dimensions, maximum minimum distance on S^3\n\n## Goal\n\n`sphere.py` exposes `place(n: int, time_budget: float, seed: int) -> list[tuple[float, ...]]`:\n`n` points, each a tuple of 4 coordinates, on the unit sphere S^3 in R^4. Maximise the smallest\npairwise Euclidean distance\n\n    D = min over i < j of |x_i - x_j|\n\nEquivalently: find the best spherical code of size `n` in 4 dimensions, or pack `n` equal caps on\nS^3. This is the 4-dimensional case of the Tammes problem (`tammes-problem` is the same\nbenchmark on S^2). Sloane's tables of putatively optimal packings cover n = 5, ..., 130 in\n4 dimensions; the thirteen `n` here were chosen away from every proven case and from the regular\npolytopes and their subsets, so every instance is open and a better configuration is a new record.\nProven optimality in R^4 covers n <= 8 (the simplex, then Rankin's bound: for d + 1 < n <= 2d the\ncross-polytope's 90 degrees is optimal), n = 10 (Bachoc and Vallentin 2009, the Petersen code, via\nsemidefinite programming) and n = 120 (the 600-cell, Boroczky 1978; also universally optimal by Cohn\nand Kumar). The 24-cell at n = 24 is only conjectured optimal (Cohn's 24-cell conjecture is open) and\nn = 23 and n = 113..120 are subsets of the 24-cell and 600-cell, so they are excluded as well.\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 with exactly 4 coordinates is\naccepted; a zero vector, a non-finite coordinate, the wrong dimension or the wrong count is a\nfailed run; coincident points simply score 0 for that `n`) and recomputes the minimum distance\nitself. 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 configurations from N. J. A. Sloane, R. H. Hardin, W. D. Smith et al., \"Tables of\nSpherical Codes\" (NeilSloane.com/packings/, files `dim4/pack.4.<n>.txt`, fetched 2026-09-07; every\nentry credited on the site to Hardin, Sloane and Smith, 1994). The record distance below was\nrecomputed from those coordinates (the files are double precision, so the recomputed distance is\nexact far beyond the 12 decimals kept); the angle is `2 asin(D / 2)` in degrees, quoted to the\ntables' 1e-7 degree precision. Running the files themselves through `eval.py` scores exactly 1.0 on\nevery `n`. A distance above a record by more than 1e-8 is listed in `records_beaten`.\n\n| n | record D | angle (deg) | source | proven optimal |\n|---|---|---|---|---|\n| 15 | 1.139304896784 | 69.4519870 | Hardin, Sloane, Smith 1994 | no |\n| 18 | 1.074412011970 | 64.9872827 | Hardin, Sloane, Smith 1994 | no |\n| 22 | 1.002113635232 | 60.1398863 | Hardin, Sloane, Smith 1994 | no |\n| 26 | 0.958342717702 | 57.2625923 | Hardin, Sloane, Smith 1994 | no |\n| 30 | 0.911885065192 | 54.2511897 | Hardin, Sloane, Smith 1994 | no |\n| 35 | 0.870029899428 | 51.5724903 | Hardin, Sloane, Smith 1994 | no |\n| 40 | 0.836064816081 | 49.4208547 | Hardin, Sloane, Smith 1994 | no |\n| 48 | 0.787615150888 | 46.3832527 | Hardin, Sloane, Smith 1994 | no |\n| 56 | 0.749397826979 | 44.0114102 | Hardin, Sloane, Smith 1994 | no |\n| 64 | 0.721722892927 | 42.3062196 | Hardin, Sloane, Smith 1994 | no |\n| 72 | 0.695669254241 | 40.7098566 | Hardin, Sloane, Smith 1994 | no |\n| 84 | 0.659737516843 | 38.5216197 | Hardin, Sloane, Smith 1994 | no |\n| 100 | 0.627382237223 | 36.5636318 | 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,26` runs a subset of n; `ZT_EVAL_PER_N_SECONDS=2` shortens the per-n budget. The\n  default is all thirteen n at 8 s each (about 110 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- In 4 dimensions there is no Fibonacci-style spiral start: use many random starts, or seed from\n  structured sets (subsets and unions of orbits of the cross-polytope, the 24-cell, or a\n  product of lower-dimensional codes) and let the optimiser break the symmetry.\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; the record configurations often have rigid\n  contact graphs that a single descent rarely finds.\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 spherical-codes-4d. Prints one JSON line: {\"metric\": record_ratio, ...}.\n\nThe solver's points are normalised onto the unit sphere S^3 in R^4 and the minimum pairwise\ndistance is recomputed 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,18,22,26,30,35,40,48,56,64,72,84,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,18,22,26,30,35,40,48,56,64,72,84,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\nDIM = 4\n\n# Best-known minimum Euclidean distance between n points on the unit sphere in R^4, computed\n# (2026-09-07) from the coordinate files pack.4.<n>.txt of N. J. A. Sloane, R. H. Hardin, W. D. Smith\n# et al., \"Tables of Spherical Codes\", NeilSloane.com/packings/ (all found in 1994 by Hardin, Sloane\n# and Smith; the minimal angle in degrees is in the comment). None of these n is proven optimal\n# (in R^4 the proven cases are n <= 8 (simplex and Rankin's 90-degree bound), n = 10 (Bachoc and Vallentin 2009) and n = 120 (the 600-cell); the 24-cell at n = 24 is only conjectured optimal). A distance more than BEAT_TOL above a record is flagged in records_beaten.\nRECORDS = {\n    15: 1.139304896784,   # 69.4519870 deg\n    18: 1.074412011970,   # 64.9872827 deg\n    22: 1.002113635232,   # 60.1398863 deg\n    26: 0.958342717702,   # 57.2625923 deg\n    30: 0.911885065192,   # 54.2511897 deg\n    35: 0.870029899428,   # 51.5724903 deg\n    40: 0.836064816081,   # 49.4208547 deg\n    48: 0.787615150888,   # 46.3832527 deg\n    56: 0.749397826979,   # 44.0114102 deg\n    64: 0.721722892927,   # 42.3062196 deg\n    72: 0.695669254241,   # 40.7098566 deg\n    84: 0.659737516843,   # 38.5216197 deg\n    100: 0.627382237223,  # 36.5636318 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, ...]]:\n    \"\"\"n finite non-zero vectors in R^DIM, 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            if len(p) != DIM:\n                fail(f\"place({n}) returned a point with {len(p)} coordinates, expected {DIM}\", \"wrong_answer\")\n            v = tuple(float(c) for c in p)\n        except SystemExit:\n            raise\n        except Exception:\n            fail(f\"place({n}) returned a non-point {p!r}\", \"wrong_answer\")\n        if not all(map(math.isfinite, v)):\n            fail(f\"place({n}) returned a non-finite coordinate\", \"wrong_answer\")\n        r = math.sqrt(sum(c * c for c in v))\n        if r < EPS:\n            fail(f\"place({n}) returned a point at the origin, which has no direction\", \"wrong_answer\")\n        out.append(tuple(c / r for c in v))\n    return out\n\n\ndef min_distance(pts: list[tuple[float, ...]], n: int) -> float:\n    best = 4.0\n    for i in range(n):\n        pi = pts[i]\n        for j in range(i + 1, n):\n            pj = pts[j]\n            d2 = 0.0\n            for a, b in zip(pi, pj):\n                d2 += (a - b) * (a - b)\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\"sphcode4|{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: seeded random start on S^3, 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\nDIM = 4\n\n\ndef _normalise(v: list[float]) -> list[float]:\n    r = math.sqrt(sum(c * c for c in v))\n    return [c / r for c in v]\n\n\ndef _random_start(n: int, rng: random.Random) -> list[list[float]]:\n    return [_normalise([rng.gauss(0.0, 1.0) for _ in range(DIM)]) for _ in range(n)]\n\n\ndef _min_dist(pts: list[list[float]]) -> float:\n    best = 4.0\n    for i in range(len(pts)):\n        pi = pts[i]\n        for j in range(i + 1, len(pts)):\n            pj = pts[j]\n            d2 = 0.0\n            for a, b in zip(pi, pj):\n                d2 += (a - b) * (a - b)\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] * DIM for _ in range(n)]\n    for i in range(n):\n        pi = pts[i]\n        fi = f[i]\n        for j in range(i + 1, n):\n            pj = pts[j]\n            diff = [a - b for a, b in zip(pi, pj)]\n            d2 = sum(c * c for c in diff)\n            t = (s * s / d2) ** (p / 2)\n            e += t\n            g = p * t / d2\n            fj = f[j]\n            for k in range(DIM):\n                fi[k] += diff[k] * g\n                fj[k] -= diff[k] * 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, fx in zip(pts, f):\n        rad = sum(a * b for a, b in zip(fx, x))\n        out.append(_normalise([a + lr * (b - rad * a) for a, b in zip(x, fx)]))\n    return out\n\n\ndef place(n: int, time_budget: float, seed: int) -> list[tuple[float, ...]]:\n    rng = random.Random(seed)\n    pts = _random_start(n, rng)\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(sum(c * c for c in fx)) for fx 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"}}