Verification notebook

Finite-field Kakeya

Proofs and exact checks for a new infinite Kakeya family, a 53-point set in F₃⁵, higher-dimensional families, and 25 benchmarks.

This notebook is a self-contained mathematical and computational companion to the Finite-field Kakeya section in the paper.

The task template, including the prompt given to the agents, is available here.

The outline is as follows:

The required packages are:

  • NumPy

1. Problem Definitions

A finite-field Kakeya set contains a complete affine line in every direction. The problem is to construct such a set with as few points as possible; exact lower-order improvements remain meaningful even though the leading asymptotic density is known.

Let pp be a prime. An affine line in Fpd\mathbb F_p^d is

L(w,v)={w+tv:tFp},v0. L(w,v)=\{w+t v:t\in\mathbb F_p\},\qquad v\ne0.

Two nonzero vectors determine the same projective direction if one is a nonzero scalar multiple of the other. We use the unique representative whose first nonzero coordinate equals 11. There are

1+p++pd1=pd1p1 1+p+\cdots+p^{d-1}=\frac{p^d-1}{p-1}

such representatives. A set is Kakeya if it contains an affine line in every projective direction.

The next cell defines the shared routines used below for enumerating projective directions, constructing affine lines, and checking finite Kakeya certificates.

Show code

Code cell 3 · In [1]

from __future__ import annotations

from collections import Counter
from fractions import Fraction
from itertools import combinations, product
from math import comb, log
from pathlib import Path
import hashlib
import json

import numpy as np


ARTIFACT_DIR = Path.cwd()


def canonical_directions(p: int, d: int):
    """Unique projective representatives: first nonzero coordinate is 1."""
    for pivot in range(d):
        prefix = (0,) * pivot + (1,)
        for tail in product(range(p), repeat=d - pivot - 1):
            yield prefix + tail


def line_points(anchor, direction, p: int):
    return {
        tuple((int(a) + t * int(v)) % p for a, v in zip(anchor, direction))
        for t in range(p)
    }


def find_kakeya_witnesses(points, p: int, d: int):
    point_set = {tuple(map(int, point)) for point in points}
    witnesses = {}
    for direction in canonical_directions(p, d):
        # Every point on a contained line can serve as its anchor.
        for anchor in point_set:
            if line_points(anchor, direction, p) <= point_set:
                witnesses[direction] = anchor
                break
        else:
            raise AssertionError(f'No contained line in direction {direction}')
    return witnesses


def assert_point_certificate(array: np.ndarray, p: int, d: int, size: int):
    assert array.shape == (size, d)
    assert np.issubdtype(array.dtype, np.integer)
    assert np.all((0 <= array) & (array < p))
    assert len({tuple(map(int, row)) for row in array}) == size
    witnesses = find_kakeya_witnesses(array, p, d)
    expected = (p**d - 1) // (p - 1)
    assert len(witnesses) == expected
    return witnesses

2. S1. A new infinite family in d=3d = 3 for p3(mod4)p \equiv 3 \pmod 4.

2.1 Construction and exact size

Let Q={u2:uFp}Q=\{u^2:u\in\mathbb F_p\}, including zero, and put S=Q=(p+1)/2S=|Q|=(p+1)/2. Define

Bp={(x,y,z):x2+4yQ, x2+4zQ} B_p=\{(x,y,z):x^2+4y\in Q,\ x^2+4z\in Q\}

and, in the plane x=0x=0,

Up=c1{(0,t,ct+c/(c1)):tFp}{(0,t,t):tFp}{(0,0,t):tFp}. U_p=\bigcup_{c\ne1}\{(0,t,ct+c/(c-1)):t\in\mathbb F_p\} \cup\{(0,t,t):t\in\mathbb F_p\} \cup\{(0,0,t):t\in\mathbb F_p\}.

Set Kp=BpUpK_p=B_p\cup U_p.

Theorem 2.1 (coverage and exact size). For every odd prime pp, KpK_p is Kakeya and

Kp={(2p3+7p21)/8,p1(mod4),(2p3+7p2+3)/8,p3(mod4). |K_p|= \begin{cases} (2p^3+7p^2-1)/8,&p\equiv1\pmod4,\\ (2p^3+7p^2+3)/8,&p\equiv3\pmod4. \end{cases}

Proof of coverage. For each graph direction (1,a,b)(1,a,b), the line

(t,at+a2,bt+b2),tFp, (t,at+a^2,bt+b^2),\qquad t\in\mathbb F_p,

lies in BpB_p, because t2+4(at+a2)=(t+2a)2t^2+4(at+a^2)=(t+2a)^2, and similarly in the last coordinate. The displayed boundary lines cover every direction (0,1,c)(0,1,c) with c1c\ne1, the diagonal covers (0,1,1)(0,1,1), and the vertical line covers (0,0,1)(0,0,1). These are all canonical directions. \square

Proof of the size formula. It remains to count Up(Q×Q)U_p\setminus(Q\times Q) in the plane x=0x=0. Write χ\chi for the quadratic character, extended by χ(0)=0\chi(0)=0, and ε=χ(1)\varepsilon=\chi(-1).

Lemma 2.2 (finite-line selector). Put

Δ(y,z)=(1yz)24yz. \Delta(y,z)=(1-y-z)^2-4yz.

The finite lines with c1c\ne1, together with the vertical line, cover exactly the pairs for which Δ(y,z)Q\Delta(y,z)\in Q.

Proof. Eliminating cc from z=cy+c/(c1)z=cy+c/(c-1) gives

yc2+(1yz)c+z=0. yc^2+(1-y-z)c+z=0.

The value c=1c=1 is never a root, since substitution gives 11. Thus for y0y\ne0 there is an allowed root exactly when the discriminant Δ\Delta is a square. At y=0y=0, the finite intercept map cc/(c1)c\mapsto c/(c-1) covers every value except 11; the vertical line restores (0,1)(0,1) and covers the whole row, while Δ(0,z)=(1z)2\Delta(0,z)=(1-z)^2. \square

Lemma 2.3 (square-grid selector count).

T:=#{(y,z)Q2:Δ(y,z)Q}=p2+8p18. T:=\#\{(y,z)\in Q^2:\Delta(y,z)\in Q\}=\frac{p^2+8p-1}{8}.

Proof. The exact indicator of QQ is (1+χ(u)+[u=0])/2(1+\chi(u)+[u=0])/2. Hence

2T=S2+y,zQχ(Δ(y,z))+#{(y,z)Q2:Δ=0}. 2T=S^2+\sum_{y,z\in Q}\chi(\Delta(y,z))+\#\{(y,z)\in Q^2:\Delta=0\}.

For y=a2,z=b2y=a^2,z=b^2,

Δ(a2,b2)=[1(a+b)2][1(ab)2]. \Delta(a^2,b^2)=[1-(a+b)^2][1-(a-b)^2].

The change (a,b)(u,v)=(a+b,ab)(a,b)\mapsto(u,v)=(a+b,a-b) is bijective. The standard quadratic-character identity

uFpχ(1u2)=ε \sum_{u\in\mathbb F_p}\chi(1-u^2)=-\varepsilon

therefore makes the double root-lift sum equal to 11. Correcting for the two axes gives

y,zQχ(Δ(y,z))=(p1)/2. \sum_{y,z\in Q}\chi(\Delta(y,z))=(p-1)/2.

If G(y,z)=χ(Δ(y,z))G(y,z)=\chi(\Delta(y,z)), then

a,bG(a2,b2)=4y,zQG(y,z)2zQG(0,z)2yQG(y,0)+G(0,0). \sum_{a,b}G(a^2,b^2) =4\sum_{y,z\in Q}G(y,z)-2\sum_{z\in Q}G(0,z)-2\sum_{y\in Q}G(y,0)+G(0,0).

Each axis sum is S1=(p1)/2S-1=(p-1)/2, and G(0,0)=1G(0,0)=1.

For the zero count, the same root lift counts the union of u=±1u=\pm1 and v=±1v=\pm1, which has 4p44p-4 ordered pairs. The two axis zero counts are each 11, and the origin is not a zero. The identical correction formula therefore gives exactly pp zeros in Q2Q^2. Substitution yields the stated value of TT. The one-variable character identity itself follows by counting the p1p-1 solutions of (ru)(r+u)=1(r-u)(r+u)=-1 and comparing them with u(1+χ(1u2))\sum_u(1+\chi(1-u^2)). \square

Lemma 2.4 (diagonal correction). The diagonal adds

A=#{yQ:χ(14y)=1}=p2+ε4 A=\#\{y\in Q:\chi(1-4y)=-1\}=\frac{p-2+\varepsilon}{4}

new square-grid points.

Proof. On z=yz=y the selector discriminant is 14y1-4y. Moreover

yQχ(14y)=1ε2, \sum_{y\in Q}\chi(1-4y)=\frac{1-\varepsilon}{2},

obtained by lifting y=a2y=a^2 and using aχ(14a2)=ε\sum_a\chi(1-4a^2)=-\varepsilon; y=1/4y=1/4 is the unique zero. Solving for the number of negative values gives AA. \square

By Lemmas 2.2–2.4,

UpQ2=T+A=p2+10p5+2ε8. |U_p\cap Q^2|=T+A=\frac{p^2+10p-5+2\varepsilon}{8}.

Lemma 2.5 (total boundary-union count). Let FF be the union of the finite lines with c1c\ne1, let D={(y,y):yFp}D=\{(y,y):y\in\mathbb F_p\}, and let V={(0,z):zFp}V=\{(0,z):z\in\mathbb F_p\}. Then

F=p(p1)2+(p1),DF=p12,V(FD)=1. |F|=\frac{p(p-1)}2+(p-1),\qquad |D\setminus F|=\frac{p-1}2,\qquad |V\setminus(F\cup D)|=1.

Proof. On the row with first coordinate yy, the finite lines cover the image of

ϕy(c)=cy+cc1,cFp{1}. \phi_y(c)=cy+\frac{c}{c-1},\qquad c\in\mathbb F_p\setminus\{1\}.

For y=0y=0, this fractional-linear map is a bijection from Fp{1}\mathbb F_p\setminus\{1\} onto itself, so the row has p1p-1 points. Assume y0y\ne0. For distinct c,dc,d,

ϕy(c)ϕy(d)=(cd)(y1(c1)(d1)). \phi_y(c)-\phi_y(d)=(c-d)\left(y-\frac1{(c-1)(d-1)}\right).

Thus a collision occurs exactly when (c1)(d1)=1/y(c-1)(d-1)=1/y. Every fiber has size at most two, because ϕy(c)=z\phi_y(c)=z is the quadratic from Lemma 2.2. Writing u=c1u=c-1 and v=d1v=d-1, there are p1p-1 ordered solutions of uv=1/yuv=1/y. The diagonal solutions satisfy u2=1/yu^2=1/y, so there are two when yy is a square and none when it is a nonsquare. Hence the number of unordered distinct collisions is (p3)/2(p-3)/2 or (p1)/2(p-1)/2, respectively, and

imϕy=p+χ(y)2. |\operatorname{im}\phi_y|=\frac{p+\chi(y)}2.

Summing over y0y\ne0 and using y0χ(y)=0\sum_{y\ne0}\chi(y)=0 gives p(p1)/2p(p-1)/2 off-axis points. This proves the formula for F|F|.

For y0y\ne0, Lemma 2.2 says that (y,z)F(y,z)\in F exactly when Δ(y,z)\Delta(y,z) is a square. On the diagonal z=yz=y, this becomes 14y1-4y. As y14yy\mapsto1-4y is a bijection, precisely (p1)/2(p-1)/2 diagonal points have nonsquare discriminant and are new; the point at y=0y=0 already belongs to FF. Finally the row y=0y=0 of FF contains every point except (0,1)(0,1), so the vertical line adds exactly that one point. \square

By Lemma 2.5,

Up=p2+2p12. |U_p|=\frac{p^2+2p-1}{2}.

Consequently

UpQ2=3p22p+12ε8. |U_p\setminus Q^2|=\frac{3p^2-2p+1-2\varepsilon}{8}.

Finally, every xx-slice of BpB_p contains S2S^2 points, so Bp=pS2|B_p|=pS^2, and Bp{x=0}=Q2B_p\cap\{x=0\}=Q^2. Adding the penalty gives Theorem 2.1. \square

The next cell constructs KpK_p directly, checks the displayed size formula, and verifies an explicit line in every projective direction for every odd prime through 4747. It also defines the finite-field helpers used in later sections.

Show code

Code cell 5 · In [2]

def is_prime(n: int) -> bool:
    if n < 2:
        return False
    if n % 2 == 0:
        return n == 2
    return all(n % q for q in range(3, int(n**0.5) + 1, 2))


def quadratic_residues(p: int):
    return {x * x % p for x in range(p)}


def legendre(a: int, p: int) -> int:
    a %= p
    if a == 0:
        return 0
    return 1 if pow(a, (p - 1) // 2, p) == 1 else -1


def one_pole_construction(p: int):
    Q = quadratic_residues(p)
    body = {
        (x, y, z)
        for x, y, z in product(range(p), repeat=3)
        if (x * x + 4 * y) % p in Q and (x * x + 4 * z) % p in Q
    }
    boundary = set()
    for c in range(p):
        if c == 1:
            continue
        intercept = c * pow(c - 1, -1, p) % p
        boundary.update((0, t, (c * t + intercept) % p) for t in range(p))
    boundary.update((0, t, t) for t in range(p))
    boundary.update((0, 0, t) for t in range(p))
    return body | boundary


def one_pole_formula(p: int) -> int:
    constant = -1 if p % 4 == 1 else 3
    return (2 * p**3 + 7 * p**2 + constant) // 8


def explicit_one_pole_witnesses(p: int):
    witnesses = {}
    for a, b in product(range(p), repeat=2):
        witnesses[(1, a, b)] = (0, a * a % p, b * b % p)
    for c in range(p):
        if c != 1:
            witnesses[(0, 1, c)] = (0, 0, c * pow(c - 1, -1, p) % p)
    witnesses[(0, 1, 1)] = (0, 0, 0)
    witnesses[(0, 0, 1)] = (0, 0, 0)
    return witnesses


one_pole_rows = []
for p in [q for q in range(3, 48, 2) if is_prime(q)]:
    K = one_pole_construction(p)
    W = explicit_one_pole_witnesses(p)
    assert len(K) == one_pole_formula(p)
    assert set(W) == set(canonical_directions(p, 3))
    assert all(line_points(w, v, p) <= K for v, w in W.items())
    Q = quadratic_residues(p)
    delta_count = sum(
        ((1 - y - z) ** 2 - 4 * y * z) % p in Q for y in Q for z in Q
    )
    assert delta_count == (p * p + 8 * p - 1) // 8
    one_pole_rows.append((p, len(K), '1 mod 4' if p % 4 == 1 else '3 mod 4'))

print('Verified exact construction rows (p, size, class):')
print(one_pole_rows)

Saved output 1

Verified exact construction rows (p, size, class):
[(3, 15, '3 mod 4'), (5, 53, '1 mod 4'), (7, 129, '3 mod 4'), (11, 439, '3 mod 4'), (13, 697, '1 mod 4'), (17, 1481, '1 mod 4'), (19, 2031, '3 mod 4'), (23, 3505, '3 mod 4'), (29, 6833, '1 mod 4'), (31, 8289, '3 mod 4'), (37, 13861, '1 mod 4'), (41, 18701, '1 mod 4'), (43, 21495, '3 mod 4'), (47, 27889, '3 mod 4')]

2.2 Equivalence with the p1(mod4)p \equiv 1 \pmod 4 AlphaEvolve construction

We use the set defined in AlphaEvolve’s commit-pinned kakeya2.tex. Let g=(p1)/4=1/4g=(p-1)/4=-1/4 in Fp\mathbb F_p and define KAE=KAK0BK0CK^{\rm AE}=K_A\cup K_{0B}\cup K_{0C} by

KA={(x,q1+q22x2g,q1q22):xFp, q1,q2Q}, K_A=\left\{\left(x,\frac{q_1+q_2}{2}-x^2-g,\frac{q_1-q_2}{2}\right):x\in\mathbb F_p,\ q_1,q_2\in Q\right\},
K0B={(0,y,z):y+z2Q},K0C={(0,y,0):yFp}. K_{0B}=\{(0,y,z):y+z^2\in Q\},\qquad K_{0C}=\{(0,y,0):y\in\mathbb F_p\}.

Theorem 2.6 (affine equivalence). If p1(mod4)p\equiv1\pmod4, then KAEK^{\rm AE} is affinely equivalent to the set KpK_p of Theorem 2.1.

Proof. Choose uFpu\in\mathbb F_p^* with u2=4u^2=-4 and define

Φ(x,y,z)=(ux,u24(y+z+g),u24(yz+g)). \Phi(x,y,z)=\left(ux,\frac{u^2}{4}(y+z+g),\frac{u^2}{4}(y-z+g)\right).

Its linear part has determinant 2u0-2u\ne0, so it is an invertible affine map. If a point lies in KAK_A, then its image (X,Y,Z)(X,Y,Z) satisfies

X2+4Y=u2q1,X2+4Z=u2q2. X^2+4Y=u^2q_1,\qquad X^2+4Z=u^2q_2.

Since u2u^2 is a nonzero square, this maps KAK_A bijectively onto BpB_p. On the plane x=0x=0, direct substitution, using u2=4u^2=-4 and 1+4g=01+4g=0, gives

Δ(Φ(0,y,z))=4(y+z2). \Delta(\Phi(0,y,z))=4(y+z^2).

Therefore Φ(K0B)\Phi(K_{0B}) is exactly the selector set {(0,Y,Z):Δ(Y,Z)Q}\{(0,Y,Z):\Delta(Y,Z)\in Q\}. Lemma 2.2 identifies this set with the finite one-pole lines together with the vertical line. Moreover Φ(K0C)={(0,t,t):tFp}\Phi(K_{0C})=\{(0,t,t):t\in\mathbb F_p\}, the diagonal line. Thus Φ(KAE)=BpUp=Kp\Phi(K^{\rm AE})=B_p\cup U_p=K_p. \square

The next cell reconstructs both sets at p=5,13,17,29p=5,13,17,29, applies the affine map above, and checks exact equality of the resulting point sets.

Show code

Code cell 7 · In [3]

def alphaevolve_kakeya2(p):
    assert p % 4 == 1
    Q = quadratic_residues(p)
    inv2 = pow(2, -1, p)
    g = (p - 1) // 4
    body = {
        (x, ((q1 + q2) * inv2 - x * x - g) % p,
         ((q1 - q2) * inv2) % p)
        for x, q1, q2 in product(range(p), Q, Q)
    }
    conic_boundary = {
        (0, y, z) for y, z in product(range(p), repeat=2)
        if (y + z * z) % p in Q
    }
    diagonal_source = {(0, y, 0) for y in range(p)}
    return body | conic_boundary | diagonal_source


def alphaevolve_to_one_pole(point, p):
    x, y, z = map(int, point)
    g = (p - 1) // 4
    u = next(a for a in range(1, p) if a * a % p == (-4) % p)
    inv4 = pow(4, -1, p)
    return (
        u * x % p,
        u * u * (y + z + g) * inv4 % p,
        u * u * (y - z + g) * inv4 % p,
    )


affine_equivalence_rows = []
for p in (5, 13, 17, 29):
    alpha_set = alphaevolve_kakeya2(p)
    mapped_set = {alphaevolve_to_one_pole(point, p) for point in alpha_set}
    station_set = one_pole_construction(p)
    assert mapped_set == station_set
    u = next(a for a in range(1, p) if a * a % p == (-4) % p)
    affine_equivalence_rows.append((p, u, len(mapped_set)))

print({
    'affine_equivalence_checks_(p,u,size)': affine_equivalence_rows,
})

Saved output 1

{'affine_equivalence_checks_(p,u,size)': [(5, 1, 53), (13, 3, 697), (17, 8, 1481), (29, 5, 6833)]}

Station Reference. The exact one-pole formula, corrected y=0y=0 selector, and affine-equivalence calculation are based on Archive #160, message archive_160-3, and Eval #2414, with supporting checks in Evals #2399, #2406, and #2410.

Related Work. Dvir (2009) proved the finite-field Kakeya conjecture, while Bukh and Chao (2021) determined the sharp leading density and developed a general framework for completing almost-Kakeya sets using lower-dimensional Kakeya sets. The infinite family above is related to that framework; its explicitly analyzed one-pole completion yields the displayed exact size formulas. Georgiev et al. (2025) reported the AlphaEvolve families used for comparison. The p1(mod4)p\equiv1\pmod4 branch is an independent rediscovery up to affine equivalence, whereas the p3(mod4)p\equiv3\pmod4 branch is new and saves (p3)/4(p-3)/4 points over AlphaEvolve’s applicable infinite family.

3. S2. Finite improvements and a 53-point Kakeya set in F35\mathbb{F}_3^5.

3.1 The 53-point construction

Theorem 3.1. If k5k_5 denotes the minimum cardinality of a Kakeya set in F35\mathbb F_3^5, then

k553. k_5\le53.

Proof. The point array kakeya_F3_d5_53.npy contains 5353 distinct points. The exact verification below shows that it contains an affine line in each of the 121121 projective directions of F35\mathbb F_3^5, so it is the required Kakeya set. \square

The next cell authenticates the stored 5353-point array and independently searches for a contained line in all 121121 directions. It also counts every contained affine line, providing a stronger check than replaying stored witnesses.

Show code

Code cell 10 · In [4]

DATA_HASHES = {'kakeya_F3_d3_13.npy': '8b1cac563e00cc53ecb5c0d38e23078fbbde36e6b6796f6c375769fbfa393737', 'kakeya_F3_d4_27.npy': 'dae8bf4b72750146aadad1de525be39fbfa49e78473202581363bacc82c14f2d', 'kakeya_F3_d5_53.npy': '4ddbba43b2cff772e4eb1ffe42064fb05e95a03972a2c448472acdcf803da812'}

path_53 = ARTIFACT_DIR / 'kakeya_F3_d5_53.npy'
if not path_53.is_file():
    raise FileNotFoundError(f'Required certificate is missing: {path_53}')
assert hashlib.sha256(path_53.read_bytes()).hexdigest() == DATA_HASHES[path_53.name]

K_53 = np.load(path_53, allow_pickle=False)
witnesses_53 = assert_point_certificate(K_53, p=3, d=5, size=53)

contained_line_counts = {}
point_set_53 = {tuple(map(int, row)) for row in K_53}
for direction in canonical_directions(3, 5):
    lines = {
        frozenset(line_points(anchor, direction, 3))
        for anchor in point_set_53
        if line_points(anchor, direction, 3) <= point_set_53
    }
    contained_line_counts[direction] = len(lines)

assert sum(contained_line_counts.values()) == 131
print({
    'points': len(point_set_53),
    'canonical_directions': len(witnesses_53),
    'contained_affine_lines': sum(contained_line_counts.values()),
    'minimum_lines_per_direction': min(contained_line_counts.values()),
    'maximum_lines_per_direction': max(contained_line_counts.values()),
    'sha256': DATA_HASHES[path_53.name],
})

Saved output 1

{'points': 53, 'canonical_directions': 121, 'contained_affine_lines': 131, 'minimum_lines_per_direction': 1, 'maximum_lines_per_direction': 2, 'sha256': '4ddbba43b2cff772e4eb1ffe42064fb05e95a03972a2c448472acdcf803da812'}

Station Reference. The 53 points originated in Eval #2520. The verification reconstructs witnesses independently from the 5353 points.

3.2 The complete 25-pair Station result

Table 3.1 (verified Station constructions).

Dimension ddPrime ppNumber of points
3313
3553
37128
311437
313697
3192,030
3233,504
3296,833
3318,288
33713,861
34118,701
34321,495
34727,889
35339,677
4327
45161
47527
4112,684
4134,962
41713,509
41920,579
5353
55490
572,135
51116,288

Theorem 3.2. For every triple (d,p,N)(d,p,N) listed in Table 3.1, there is a Kakeya set of NN points in Fpd\mathbb F_p^d.

Verification. For each row, the bundled certificate supplies a point set of cardinality NN and a witness line in every projective direction. The next cell checks the point count, distinctness, coordinate range, complete direction list, and containment of every witness line. A successful run verifies all 2525 existence claims. \square

The comparison with AlphaEvolve and the earlier literature is deferred to Section 7.

Show code

Code cell 12 · In [5]

STATION_MANIFEST_SHA256 = 'f19d7b03cfdd5956386f2ef4ef2aa41009783ba0b5137d8b54d0ec360b22afad'
STATION_BUNDLE_SHA256 = '105a642a701a4d26f1e360ac041504731c2d7d7bc13d1c185dc5b5013c21dfd5'

station_manifest_path = ARTIFACT_DIR / 'station_point_certificates.json'
assert hashlib.sha256(station_manifest_path.read_bytes()).hexdigest() == STATION_MANIFEST_SHA256
station_manifest = json.loads(station_manifest_path.read_text())

station_bundle_path = ARTIFACT_DIR / station_manifest['bundle']
assert hashlib.sha256(station_bundle_path.read_bytes()).hexdigest() == STATION_BUNDLE_SHA256
assert station_manifest['bundle_sha256'] == STATION_BUNDLE_SHA256

station_certificate_checks = []
with np.load(station_bundle_path, allow_pickle=False) as station_bundle:
    for record in sorted(station_manifest['records'], key=lambda row: (row['d'], row['p'])):
        d, p, size = record['d'], record['p'], record['size']
        point_array = station_bundle[record['points_key']]
        witness_array = station_bundle[record['witnesses_key']]
        directions = list(canonical_directions(p, d))

        assert point_array.shape == (size, d)
        assert witness_array.shape == (len(directions), d)
        assert np.all((0 <= point_array) & (point_array < p))
        points = {tuple(map(int, point)) for point in point_array}
        assert len(points) == size

        for direction, anchor in zip(directions, witness_array):
            assert line_points(anchor, direction, p) <= points
        station_certificate_checks.append((d, p, size, len(directions)))

assert [(d, p, size) for d, p, size, _ in station_certificate_checks] == [
    (3, 3, 13), (3, 5, 53), (3, 7, 128), (3, 11, 437),
    (3, 13, 697), (3, 19, 2030), (3, 23, 3504), (3, 29, 6833),
    (3, 31, 8288), (3, 37, 13861), (3, 41, 18701), (3, 43, 21495),
    (3, 47, 27889), (3, 53, 39677), (4, 3, 27), (4, 5, 161),
    (4, 7, 527), (4, 11, 2684), (4, 13, 4962), (4, 17, 13509),
    (4, 19, 20579), (5, 3, 53), (5, 5, 490), (5, 7, 2135),
    (5, 11, 16288),
]

print('Verified Station constructions:', len(station_certificate_checks))
print('Verified projective directions:', sum(row[3] for row in station_certificate_checks))

Saved output 1

Verified Station constructions: 25
Verified projective directions: 50066

Station Reference. The 25 Station point certificates are reconstructed from Evals #3, #12, #42, #52, #690, #730, #975, #2053, #2520, and #2529. The bundled point sets and witness lines are checked above.

Related Work. Lev (2009) recorded k3=13k_3=13 and k427k_4\le27 and proposed the recurrence that predicts k5=53k_5=53. The previous bound k563k_5\le63 is attained by the missing-digit construction of Kopparty et al. (2011) and was also AlphaEvolve’s value. The 53-point construction improves that bound and matches Lev’s prediction.

4. S3. Structural analysis of the new infinite family.

Fix λFp\lambda\in\mathbb F_p^* and put E0=(λQ)2E_0=(\lambda Q)^2. Choose A,B,r,u,vFpA,B,r,u,v\in\mathbb F_p, set D=Ar+BD=Ar+B, and assume D0D\ne0. For crc\ne r, define

βc=Ac+Bcr=A+Dcr,Lc={(y,cy+βc):yFp}. \beta_c=\frac{Ac+B}{c-r}=A+\frac{D}{c-r},\qquad L_c=\{(y,cy+\beta_c):y\in\mathbb F_p\}.

Add Lr(u)={(y,ry+u)}L_r(u)=\{(y,ry+u)\} and the vertical line L(v)={(v,z)}L_\infty(v)=\{(v,z)\}. Let UU be the union of these p+1p+1 lines and define the boundary penalty by pen(p)=UE0\operatorname{pen}(p)=|U\setminus E_0|.

Theorem 4.1 (uniform full-family penalty). Let pp be an odd prime. For every λFp\lambda\in\mathbb F_p^* and every A,B,r,u,vFpA,B,r,u,v\in\mathbb F_p satisfying Ar+B0Ar+B\ne0,

pen(p)3p28<5p. \left|\operatorname{pen}(p)-\frac{3p^2}{8}\right|<5p.

Proof. Let Uf=crLcU_f=\bigcup_{c\ne r}L_c. Eliminating cc gives

Δ(y,z)=(zAry)24Dy. \Delta(y,z)=(z-A-ry)^2-4Dy.

For y0y\ne0, membership in UfU_f is equivalent to Δ(y,z)Q\Delta(y,z)\in Q. As a polynomial in zz, Δ\Delta is monic with discriminant 16Dy16Dy. Therefore zχ(Δ)=1\sum_z\chi(\Delta)=-1 and the number of its zeros is 1+χ(Dy)1+\chi(Dy), so the number of covered zz is (p+χ(Dy))/2(p+\chi(Dy))/2. Summation over y0y\ne0 gives p(p1)/2p(p-1)/2 points. At y=0y=0, the values βc=A+D/(cr)\beta_c=A+D/(c-r) run bijectively through Fp{A}\mathbb F_p\setminus\{A\}. Hence

Uf=p(p1)2+(p1).(1) |U_f|=\frac{p(p-1)}2+(p-1).\tag{1}

It remains to bound O=UfE0O=|U_f\cap E_0|. Put

Y=χ(y/λ),Z=χ(z/λ),X=χ(Δ), Y=\chi(y/\lambda),\quad Z=\chi(z/\lambda),\quad X=\chi(\Delta),

and write ey,ez,eΔe_y,e_z,e_\Delta for the indicators of y=0,z=0,Δ=0y=0,z=0,\Delta=0. The exact square-set indicator is (1+Y+ey)/2(1+Y+e_y)/2. The selector (1+X+eΔ)/2(1+X+e_\Delta)/2 agrees with UfU_f off the row y=0y=0; on that row it additionally selects the omitted point (0,A)(0,A). Thus, if OO^* denotes the selector overlap,

8O=y,z(1+Y+ey)(1+Z+ez)(1+X+eΔ),O=O[AλQ].(2) 8O^*=\sum_{y,z}(1+Y+e_y)(1+Z+e_z)(1+X+e_\Delta),\qquad O=O^*-[A\in\lambda Q].\tag{2}

The part without XX or eΔe_\Delta is (p+1)2(p+1)^2. We next bound the XX-block. The monic-quadratic identity used above gives

X=YX=0. \sum X=\sum YX=0.

For S01=χ(z)χ(Δ)S_{01}=\sum\chi(z)\chi(\Delta), substitute x=zAryx=z-A-ry. If r0r\ne0, the inner sum over yy is a correlation of two linear characters. Their roots coincide precisely when

rx2+4Dx+4DA=0. rx^2+4Dx+4DA=0.

If this quadratic has N{0,1,2}N\in\{0,1,2\} roots, the correlation identity yχ((ay+b)(cy+d))=χ(ac)\sum_y\chi((a y+b)(c y+d))=-\chi(ac) for distinct roots and (p1)χ(ac)(p-1)\chi(ac) for equal roots gives

S01=p(N1)χ(4Dr). S_{01}=p(N-1)\chi(-4Dr).

Thus S01p|S_{01}|\le p; when r=0r=0 the inner sum is zero. Multiplication by the harmless factor χ(λ)\chi(\lambda) proves ZXp|\sum ZX|\le p.

For S11=χ(yz)χ(Δ)S_{11}=\sum\chi(yz)\chi(\Delta), restrict to y0y\ne0 and set t=z/yt=z/y. The inner polynomial in yy is

qt(y)=(tr)2y22(A(tr)+2D)y+A2, q_t(y)=(t-r)^2y^2-2(A(t-r)+2D)y+A^2,

whose discriminant is 16D(At+B)16D(At+B). If A0A\ne0, the complete character sum is 1-1 except that it is p1p-1 at t=B/At=-B/A and 00 at t=rt=r; subtracting the y=0y=0 term and summing with weight χ(t)\chi(t) gives

S11=pχ(B/A)+χ(r). S_{11}=p\chi(-B/A)+\chi(r).

If A=0A=0, the same calculation gives S11=χ(r)S_{11}=\chi(r). Hence YZXp+1|\sum YZX|\le p+1. The remaining zero-coordinate XX-terms group as

z(1+Z+ez)X(0,z),y(1+Y)X(y,0), \sum_z(1+Z+e_z)X(0,z),\qquad \sum_y(1+Y)X(y,0),

with absolute values at most p+1p+1 and pp. The entire XX-block therefore has absolute value at most 4p+24p+2.

Finally Δ=0\Delta=0 has y(1+χ(Dy))=p\sum_y(1+\chi(Dy))=p affine solutions. The eΔe_\Delta-block is exactly four times the number of these solutions in E0E_0, hence is at most 4p4p. From (2),

8Op2(2p+1)+(4p+2)+4p=10p+3. |8O^*-p^2|\le(2p+1)+(4p+2)+4p=10p+3.

Removing the possible extra point changes the numerator by at most 88, so, for every odd pp,

8Op210p+11<16p,Op2/8<2p.(3) |8O-p^2|\le10p+11<16p,\qquad |O-p^2/8|<2p.\tag{3}

Let EE be the new contribution outside E0UfE_0\cup U_f from the pole and vertical lines. Since each has pp points, 0E2p0\le E\le2p. By (1) and (3),

pen(p)3p28=(p21)(Op28)+E, \operatorname{pen}(p)-\frac{3p^2}{8}=\left(\frac p2-1\right)-\left(O-\frac{p^2}{8}\right)+E,

and therefore pen(p)3p2/8<(p/21)+2p+2p<5p|\operatorname{pen}(p)-3p^2/8|<(p/2-1)+2p+2p<5p. \square

The next cell exhausts all admissible parameters at p=3p=3 and p=5p=5. For each choice it checks the exact finite-line union count and the bound pen(p)3p2/8<5p|\operatorname{pen}(p)-3p^2/8|<5p proved above.

Show code

Code cell 15 · In [6]

def full_one_pole_boundary(p, lam, A, B, r, u, v):
    Q = quadratic_residues(p)
    D = (A * r + B) % p
    assert lam % p != 0 and D != 0
    finite = set()
    for c in range(p):
        if c == r:
            continue
        beta = (A * c + B) * pow(c - r, -1, p) % p
        finite.update((y, (c * y + beta) % p) for y in range(p))
    boundary = finite | {(y, (r * y + u) % p) for y in range(p)}
    boundary |= {(v, z) for z in range(p)}
    footprint = {(lam * q % p, lam * s % p) for q in Q for s in Q}
    return finite, boundary, len(boundary - footprint)


full_family_checks = 0
for p in (3, 5):
    for lam in range(1, p):
        for A, B, r, u, v in product(range(p), repeat=5):
            if (A * r + B) % p == 0:
                continue
            finite, boundary, penalty = full_one_pole_boundary(
                p, lam, A, B, r, u, v
            )
            assert len(finite) == p * (p - 1) // 2 + (p - 1)
            assert abs(penalty - Fraction(3 * p * p, 8)) < 5 * p
            full_family_checks += 1

print({'full_one_pole_parameter_checks': full_family_checks})

Saved output 1

{'full_one_pole_parameter_checks': 10324}

Station Reference. The uniform full-family penalty theorem is based on Archive #160 and the supporting checks in Evals #2399, #2406, and #2410.

Related Work. Bukh and Chao (2021) established the sharp leading density and a construction framework that determines the first two asymptotic terms. AlphaEvolve, as reported by Georgiev et al. (2025), gives particular completion families rather than a uniform bound for the full one-pole class. The theorem above proves the class-wide 3p2/8+O(p)3p^2/8+O(p) boundary penalty for every nondegenerate one-pole completion.

5. Higher-dimensional uniform families

The infinite families for d=4d=4 and d=5d=5 described in the paper’s limitations are presented here.

5.1 The four-dimensional prefix-sum family

Let

Q={u2:uFp},S=Q=p+12. \mathcal Q=\{u^2:u\in\mathbb F_p\},\qquad S=|\mathcal Q|=\frac{p+1}{2}.

For n1n\ge1, let MnM_n be the lower-triangular matrix

(Mn)ij={1,ji,0,j>i. (M_n)_{ij}=\begin{cases}1,&j\le i,\\0,&j>i.\end{cases}

Thus Mn(a1,,an)T=(a1,a1+a2,,a1++an)TM_n(a_1,\ldots,a_n)^T=(a_1,a_1+a_2,\ldots,a_1+\cdots+a_n)^T. Define

Fn(a)=Mn1((Mna)12(Mna)n2). F_n(a)=-M_n^{-1} \begin{pmatrix} (M_na)_1^2\\ \vdots\\ (M_na)_n^2 \end{pmatrix}.

Set K1=FpK_1=\mathbb F_p. Recursively, for n1n\ge1, define

Tn+1={(t,Fn(a)+ta):tFp, aFpn},Kn+1=Tn+1({0}×Kn).(4) T_{n+1}=\{(t,F_n(a)+ta):t\in\mathbb F_p,\ a\in\mathbb F_p^n\}, \qquad K_{n+1}=T_{n+1}\cup(\{0\}\times K_n). \tag{4}

Theorem 5.1. For every odd prime p5p\ge5, the set K4K_4 defined in (4) is Kakeya in Fp4\mathbb F_p^4 and has size

K4={18p4+1932p3+2532p2332p1332,p1(mod4),18p4+1932p3+2532p2732p4932,p3(mod4). |K_4|= \begin{cases} \dfrac18p^4+\dfrac{19}{32}p^3+\dfrac{25}{32}p^2-\dfrac3{32}p-\dfrac{13}{32},&p\equiv1\pmod4,\\[4pt] \dfrac18p^4+\dfrac{19}{32}p^3+\dfrac{25}{32}p^2-\dfrac7{32}p-\dfrac{49}{32},&p\equiv3\pmod4. \end{cases}

Proof. For each aFpna\in\mathbb F_p^n, the line

t(t,Fn(a)+ta) t\longmapsto(t,F_n(a)+ta)

lies in Tn+1T_{n+1} and has direction (1,a)(1,a). By induction, {0}×Kn\{0\}\times K_n contains a line in every direction whose first coordinate is zero. Hence Kn+1K_{n+1} is Kakeya.

Fix t0t\ne0. Applying MnM_n to the transverse slice gives

Mn(Fn(a)+ta)=((Mna)i2+t(Mna)i)i=1n. M_n(F_n(a)+ta)= \bigl(-(M_na)_i^2+t(M_na)_i\bigr)_{i=1}^n.

Each coordinate ranges over a set of SS values, so the slice has SnS^n points. The nonzero slices of K4K_4 therefore contain (p1)S3(p-1)S^3 points.

At t=0t=0, the top slice is

W3=M31((Q)3), W_3=M_3^{-1}((-\mathcal Q)^3),

and the zero slice of K4K_4 is W3K3W_3\cup K_3. For δ0\delta\ne0, put

N(δ)=Q(Q+δ)=p+1+χ(δ)+χ(δ)4. N(\delta)=|\mathcal Q\cap(\mathcal Q+\delta)| =\frac{p+1+\chi(\delta)+\chi(-\delta)}4.

Partitioning W3K3W_3\cap K_3 by its first coordinate gives

W3K3=2S2+xQ{0,4}N(xx24)2. |W_3\cap K_3| =2S^2+\sum_{x\in-\mathcal Q\setminus\{0,-4\}} N\left(-x-\frac{x^2}{4}\right)^2.

Substituting the displayed formula for NN and summing the quadratic characters yields

W3K3={11p3+25p2+5p932,p1(mod4),11p3+25p2+p4532,p3(mod4). |W_3\cup K_3|= \begin{cases} \dfrac{11p^3+25p^2+5p-9}{32},&p\equiv1\pmod4,\\[4pt] \dfrac{11p^3+25p^2+p-45}{32},&p\equiv3\pmod4. \end{cases}

Adding the nonzero slices gives the stated formula. \square

5.2 The five-dimensional nested tangent family

Write χ\chi for the quadratic character and retain Q={u2:uFp}\mathcal Q=\{u^2:u\in\mathbb F_p\} and S=(p+1)/2S=(p+1)/2. Choose αFp\alpha\in\mathbb F_p^* with χ(α)=χ(1)\chi(\alpha)=\chi(-1), and choose eQe\in\mathcal Q so that

χ(αe2+e)=1if p1(mod4),αe2+e0if p3(mod4). \chi(\alpha e^2+e)=-1\quad\text{if }p\equiv1\pmod4, \qquad \alpha e^2+e\ne0\quad\text{if }p\equiv3\pmod4.

Define the planar tangent completion

P2=mFp{(y,z):y=mz(m1)24α}    {(y,e):yFp}. P_2= \bigcup_{m\in\mathbb F_p} \left\{(y,z):y=mz-\frac{(m-1)^2}{4\alpha}\right\} \;\cup\; \{(y,e):y\in\mathbb F_p\}.

For d3d\ge3, define

Ad={(x0,,xd1):xi+x02Q for 1i<d}. A_d=\{(x_0,\ldots,x_{d-1}):x_i+x_0^2\in\mathcal Q \text{ for }1\le i<d\}.

Set

K3=A3({0}×P2),Kd=Ad({0}×Kd1)(d4).(5) K_3=A_3\cup(\{0\}\times P_2), \qquad K_d=A_d\cup(\{0\}\times K_{d-1})\quad(d\ge4). \tag{5}

Lemma 5.2 (planar base calculation). A permitted value of ee exists. The set P2P_2 is Kakeya and

P2=pS+S1, |P_2|=pS+S-1,

while

P2Q2=p2+8p+78. |P_2\cap\mathcal Q^2|=\frac{p^2+8p+7}{8}.

Consequently,

K3=pS2+P2Q2=2p3+7p2+2p118. |K_3|=pS^2+|P_2\setminus\mathcal Q^2| =\frac{2p^3+7p^2+2p-11}{8}.

Proof. Let TT be the union of the pp sloped lines in P2P_2 and put c(z)=αz2+zc(z)=\alpha z^2+z. Eliminating mm gives

m2(2+4αz)m+(1+4αy)=0, m^2-(2+4\alpha z)m+(1+4\alpha y)=0,

whose discriminant is 16α(c(z)y)16\alpha(c(z)-y). Thus (y,z)T(y,z)\in T exactly when α(c(z)y)Q\alpha(c(z)-y)\in\mathcal Q. Every zz-row of TT has SS points, so T=pS|T|=pS; the horizontal line z=ez=e adds S1S-1 points. The sloped lines cover the directions (m,1)(m,1) and the horizontal line covers (1,0)(1,0), so P2P_2 is Kakeya.

Write ε=χ(1)=χ(α)\varepsilon=\chi(-1)=\chi(\alpha). If ε=1\varepsilon=-1, choose a nonzero square e1/αe\ne-1/\alpha; such a choice exists for p5p\ge5, and c(e)0c(e)\ne0. If ε=1\varepsilon=1, multiplication by α\alpha permutes the squares, and

xQ{0}χ(1+x)=1. \sum_{x\in\mathcal Q\setminus\{0\}}\chi(1+x)=-1.

This identity follows by expanding the square indicator and using the standard quadratic sum xχ(x(x+1))=1\sum_x\chi(x(x+1))=-1. Hence some nonzero square xx has χ(1+x)=1\chi(1+x)=-1; taking e=x/αe=x/\alpha gives χ(c(e))=1\chi(c(e))=-1.

For c0c\ne0, expanding the same indicator gives

#{yQ:cyQ}=p+1+2χ(c)4(ε=1) \#\{y\in\mathcal Q:c-y\in\mathcal Q\} =\frac{p+1+2\chi(c)}4\qquad(\varepsilon=1)

and

#{yQ:cy(FpQ){0}}=p+14(ε=1); \#\{y\in\mathcal Q:c-y\in(\mathbb F_p\setminus\mathcal Q)\cup\{0\}\} =\frac{p+1}4\qquad(\varepsilon=-1);

at c=0c=0 the corresponding count is SS. The polynomial c(z)=z(αz+1)c(z)=z(\alpha z+1) has exactly two zeros in Q\mathcal Q. Summing over zQz\in\mathcal Q gives

TQ2={(p2+6p+1)/8,ε=1,(p2+6p+5)/8,ε=1. |T\cap\mathcal Q^2|= \begin{cases} (p^2+6p+1)/8,&\varepsilon=1,\\ (p^2+6p+5)/8,&\varepsilon=-1. \end{cases}

The horizontal line adds respectively (p+3)/4(p+3)/4 or (p+1)/4(p+1)/4 new points of Q2\mathcal Q^2. Therefore P2Q2=(p2+8p+7)/8|P_2\cap\mathcal Q^2|=(p^2+8p+7)/8 in both cases, and the remaining formulas follow. \square

Put

r=p+2+χ(1)4. r=\frac{p+2+\chi(-1)}4.

For n3n\ge3, inclusion–exclusion in (5) gives

Kn+1=pSn+KnLn,Ln=Sn1+(S1)rn1.(6) |K_{n+1}|=pS^n+|K_n|-L_n, \qquad L_n=S^{n-1}+(S-1)r^{n-1}. \tag{6}

Theorem 5.3. For every odd prime p5p\ge5, the set K5K_5 defined by (5) is Kakeya in Fp5\mathbb F_p^5 and has size

K5={116p5+47128p4+2532p3+3764p21732p161128,p1(mod4),116p5+47128p4+5364p3+2732p22564p219128,p3(mod4). |K_5|= \begin{cases} \dfrac1{16}p^5+\dfrac{47}{128}p^4+\dfrac{25}{32}p^3+\dfrac{37}{64}p^2-\dfrac{17}{32}p-\dfrac{161}{128},&p\equiv1\pmod4,\\[4pt] \dfrac1{16}p^5+\dfrac{47}{128}p^4+\dfrac{53}{64}p^3+\dfrac{27}{32}p^2-\dfrac{25}{64}p-\dfrac{219}{128},&p\equiv3\pmod4. \end{cases}

Proof. At each lift, the new body has pSnpS^n points and supplies all graph directions; the embedded copy of KnK_n supplies the directions in the coordinate hyperplane. Their intersection has size LnL_n. Indeed, the overlap splits into Sn1S^{n-1} points from the zero parameter and (S1)rn1(S-1)r^{n-1} points from the nonzero square parameters, where the standard quadratic-character count gives

r=Sqp(Sqp+1)=p+2+χ(1)4. r=|\operatorname{Sq}_p\cap(\operatorname{Sq}_p+1)| =\frac{p+2+\chi(-1)}4.

Inclusion–exclusion proves (6). Substituting the dimension-three size and expanding the cases χ(1)=1\chi(-1)=1 and χ(1)=1\chi(-1)=-1 through dimensions four and five gives the displayed formulas. \square

The next cell reconstructs both higher-dimensional families at three primes, checks their exact sizes, expands the size laws, and checks the third-coefficient comparison with AlphaEvolve.

Show code

Code cell 19 · In [7]

def poly_add(*polynomials):
    degree = max(map(len, polynomials))
    return tuple(
        sum((poly[i] if i < len(poly) else Fraction(0)) for poly in polynomials)
        for i in range(degree)
    )


def poly_scale(poly, scalar):
    scalar = Fraction(scalar)
    return tuple(scalar * coefficient for coefficient in poly)


def poly_mul(left, right):
    result = [Fraction(0)] * (len(left) + len(right) - 1)
    for i, a in enumerate(left):
        for j, b in enumerate(right):
            result[i + j] += a * b
    return tuple(result)


def poly_pow(poly, exponent):
    result = (Fraction(1),)
    for _ in range(exponent):
        result = poly_mul(result, poly)
    return result


p_poly = (Fraction(0), Fraction(1))
one_poly = (Fraction(1),)
S_poly = poly_scale(poly_add(p_poly, one_poly), Fraction(1, 2))

prefix_expected = {
    1: (Fraction(-13, 32), Fraction(-3, 32), Fraction(25, 32),
        Fraction(19, 32), Fraction(1, 8)),
    -1: (Fraction(-49, 32), Fraction(-7, 32), Fraction(25, 32),
         Fraction(19, 32), Fraction(1, 8)),
}
prefix_zero_slice = {
    1: (Fraction(-9, 32), Fraction(5, 32), Fraction(25, 32),
        Fraction(11, 32)),
    -1: (Fraction(-45, 32), Fraction(1, 32), Fraction(25, 32),
         Fraction(11, 32)),
}
prefix_body = poly_mul((-1, 1), poly_pow(S_poly, 3))
for character in (1, -1):
    assert poly_add(prefix_body, prefix_zero_slice[character]) == prefix_expected[character]


def prefix_sum_kakeya(p, dimension):
    if dimension == 1:
        return {(t,) for t in range(p)}
    n = dimension - 1
    top = set()
    for slope in product(range(p), repeat=n):
        cumulative = []
        total = 0
        for value in slope:
            total = (total + value) % p
            cumulative.append(total)
        squared = [(-value * value) % p for value in cumulative]
        anchor = tuple(
            (squared[i] - (squared[i - 1] if i else 0)) % p
            for i in range(n)
        )
        for t in range(p):
            top.add((t,) + tuple(
                (anchor[i] + t * slope[i]) % p for i in range(n)
            ))
    boundary = {(0,) + point for point in prefix_sum_kakeya(p, n)}
    return top | boundary


def poly_evaluate(poly, value):
    return sum(coefficient * value**degree
               for degree, coefficient in enumerate(poly))


prefix_construction_sizes = {}
for prime in (5, 7, 11):
    character = 1 if prime % 4 == 1 else -1
    points = prefix_sum_kakeya(prime, 4)
    expected_size = poly_evaluate(prefix_expected[character], prime)
    assert expected_size.denominator == 1
    assert len(points) == expected_size
    prefix_construction_sizes[prime] = len(points)

tangent_three = (
    Fraction(-11, 8), Fraction(1, 4), Fraction(7, 8), Fraction(1, 4)
)


def tangent_lift(size_poly, dimension, character):
    r_poly = poly_scale(poly_add(p_poly, (2 + character,)), Fraction(1, 4))
    overlap = poly_add(
        poly_pow(S_poly, dimension - 1),
        poly_mul(poly_add(S_poly, (-1,)), poly_pow(r_poly, dimension - 1)),
    )
    body = poly_mul(p_poly, poly_pow(S_poly, dimension))
    return poly_add(body, size_poly, poly_scale(overlap, -1))


tangent_five_expected = {
    1: (Fraction(-161, 128), Fraction(-17, 32), Fraction(37, 64),
        Fraction(25, 32), Fraction(47, 128), Fraction(1, 16)),
    -1: (Fraction(-219, 128), Fraction(-25, 64), Fraction(27, 32),
         Fraction(53, 64), Fraction(47, 128), Fraction(1, 16)),
}
tangent_five = {}
for character in (1, -1):
    tangent_four = tangent_lift(tangent_three, 3, character)
    tangent_five[character] = tangent_lift(tangent_four, 4, character)
    assert tangent_five[character] == tangent_five_expected[character]

def quadratic_character(value, p):
    value %= p
    if value == 0:
        return 0
    return 1 if pow(value, (p - 1) // 2, p) == 1 else -1


def tangent_kakeya(p, dimension):
    squares = quadratic_residues(p)
    character = quadratic_character(-1, p)
    alpha = next(a for a in range(1, p)
                 if quadratic_character(a, p) == character)
    if character == 1:
        e = next(x for x in sorted(squares)
                 if quadratic_character(alpha * x * x + x, p) == -1)
    else:
        e = next(x for x in sorted(squares)
                 if (alpha * x * x + x) % p != 0)

    inverse = pow(4 * alpha, -1, p)
    completion = {
        ((m * z - (m - 1) ** 2 * inverse) % p, z)
        for m in range(p)
        for z in range(p)
    }
    completion |= {(y, e) for y in range(p)}

    def leading_block(dimension):
        return {
            (x0,) + tuple((q - x0 * x0) % p for q in tail)
            for x0 in range(p)
            for tail in product(squares, repeat=dimension - 1)
        }

    points = leading_block(3) | {(0,) + point for point in completion}
    for current_dimension in range(4, dimension + 1):
        points = leading_block(current_dimension) | {
            (0,) + point for point in points
        }
    return points


tangent_construction_sizes = {}
for prime in (5, 7, 11):
    character = 1 if prime % 4 == 1 else -1
    points = tangent_kakeya(prime, 5)
    expected_size = poly_evaluate(tangent_five_expected[character], prime)
    assert expected_size.denominator == 1
    assert len(points) == expected_size
    tangent_construction_sizes[prime] = len(points)


d4_third_gap = prefix_expected[1][2] - Fraction(11, 16)
d5_third_gap = tangent_five[1][3] - Fraction(177, 256)
assert d4_third_gap == Fraction(3, 32) > 0
assert d5_third_gap == Fraction(23, 256) > 0

print('Exact coefficient gaps over AlphaEvolve:', {
    'd=4, coefficient of p^2': d4_third_gap,
    'd=5, coefficient of p^3': d5_third_gap,
})

Saved output 1

Exact coefficient gaps over AlphaEvolve: {'d=4, coefficient of p^2': Fraction(3, 32), 'd=5, coefficient of p^3': Fraction(23, 256)}

Station Reference. The nested tangent construction and its overlap identity are from Archive #50. The four-dimensional prefix-sum construction and its zero-slice decomposition are from Archive #71.

Related Work. On p1(mod4)p\equiv1\pmod4, the uniform constructions reported by Georgiev et al. (2025) have sizes

18p4+1932p3+1116p2+O(p3/2) \frac18p^4+\frac{19}{32}p^3+\frac{11}{16}p^2+O(p^{3/2})

in dimension four and

116p5+47128p4+177256p3+O(p5/2) \frac1{16}p^5+\frac{47}{128}p^4+\frac{177}{256}p^3+O(p^{5/2})

in dimension five. The Station formulas have the same first two coefficients but exceed these third coefficients by 3/323/32 in dimension four and 23/25623/256 in dimension five. Since smaller Kakeya sets are better, both Station families are asymptotically worse than the corresponding AlphaEvolve families.

6. Additional findings

These are additional results that may be interesting but are not included in the spotlight.

6.1 Exact polynomial evaluation-code distance

Dvir’s polynomial method detects a small Kakeya set through a nonzero polynomial of degree at most p1p-1 that vanishes on the entire set. We determine exactly how close the classical quadratic-residue construction CdC_d comes to admitting such a certificate: every such polynomial that is not identically zero on CdC_d remains nonzero on at least p+1p+1 points when d=3d=3, and at least 2((p+1)/2)d22((p+1)/2)^{d-2} points when d4d\ge4, with equality in both cases. Thus the theorem gives the exact obstruction to applying Dvir’s vanishing argument to this construction, while also determining the minimum distance of the associated punctured Reed–Muller code. We prove the upper bounds with an explicit product of p1p-1 hyperplanes and the lower bounds by a weighted footprint argument, supplemented in dimension three by a shadow code and body-kernel decomposition.

The set and its evaluation code

Let pp be an odd prime, let

S=p+12,Q={u2:uFp}, S=\frac{p+1}{2},\qquad \mathcal Q=\{u^2:u\in\mathbb F_p\},

where 0Q0\in\mathcal Q, and define the quadratic-residue cone

Cd=PB, C_d=P\sqcup B,

where

P={(0,x1,,xd1):xiFp} P=\{(0,x_1,\ldots,x_{d-1}):x_i\in\mathbb F_p\}

and

B={(x0,,xd1):x00, xi+x02Q for every i1}. B=\{(x_0,\ldots,x_{d-1}):x_0\ne0,\ x_i+x_0^2\in\mathcal Q\text{ for every }i\ge1\}.

For a polynomial ff of total degree at most p1p-1, its evaluation word is the function fCdf|_{C_d}. Its Hamming weight is

wtCd(f)=#{xCd:f(x)0}. \operatorname{wt}_{C_d}(f)=\#\{x\in C_d:f(x)\ne0\}.

The minimum distance μ(Cd)\mu(C_d) is the least positive weight of such an evaluation word.

Theorem 6.1 (exact distance). For every odd prime pp,

μ(Cd)=2Sd2(d4). \mu(C_d)=2S^{d-2}\quad(d\ge4).

For d=3d=3,

μ(C3)=p+1(p5). \mu(C_3)=p+1\quad(p\ge5).

The restriction p5p\ge5 in dimension three is necessary: at p=3p=3 the distance is 33, not 44.

An explicit polynomial certificate

Lemma 6.2 (explicit upper certificate). For every d3d\ge3, a polynomial of degree p1p-1 has weight exactly 2Sd22S^{d-2} on CdC_d.

Proof. Take

F(x)=x0xd1cFp{1/2,1/2}(x0+cxd114c). F^*(x)=x_0x_{d-1} \prod_{c\in\mathbb F_p^*\setminus\{1/2,-1/2\}} \left(x_0+c x_{d-1}-\frac1{4c}\right).

It has 2+(p3)=p12+(p-3)=p-1 linear factors. The first factor makes it zero on PP, and the second covers the part of BB with xd1=0x_{d-1}=0. At a remaining body point, one of the other factors vanishes precisely when

4xd1c2+4x0c1=0. 4x_{d-1}c^2+4x_0c-1=0.

Its discriminant is 16(x02+xd1)16(x_0^2+x_{d-1}), a square by the definition of BB. If the two roots are distinct, they cannot both be the excluded pair {1/2,1/2}\{1/2,-1/2\}: their sum is x0/xd10-x_0/x_{d-1}\ne0, whereas the excluded pair has sum zero. The only possible failure is therefore a double root. This requires

xd1=x02, x_{d-1}=-x_0^2,

and the root is c=1/(2x0)c=1/(2x_0). It is excluded exactly when x0=±1x_0=\pm1. Hence the nonzero set of FF^* on CdC_d is

{(ε,q11,,qd21,1):ε{1,1}, qiQ}, \{(\varepsilon,q_1-1,\ldots,q_{d-2}-1,-1): \varepsilon\in\{1,-1\},\ q_i\in\mathcal Q\},

which has 2Sd22S^{d-2} points. \square

Two elementary polynomial bounds

We use the following finite-grid form of the footprint bound.

Lemma 6.3 (Cartesian-grid footprint). Let AiFpA_i\subseteq\mathbb F_p have NiN_i elements. Suppose a nonzero polynomial RR, reduced so that degXiR<Ni\deg_{X_i}R<N_i, has lexicographic leading monomial X1a1XmamX_1^{a_1}\cdots X_m^{a_m}. Then RR is nonzero at at least

i=1m(Niai) \prod_{i=1}^m(N_i-a_i)

points of A1××AmA_1\times\cdots\times A_m.

Proof. Induct on mm. View RR as a polynomial in X1X_1. Its leading coefficient, as a polynomial in the remaining variables, has leading exponents a2,,ama_2,\ldots,a_m, so by induction it is nonzero at at least i=2m(Niai)\prod_{i=2}^m(N_i-a_i) choices of those variables. At each such choice the resulting nonzero polynomial in X1X_1 has degree a1a_1, hence at most a1a_1 roots in A1A_1. \square

Lemma 6.4 (full-grid Reed–Muller bound). A nonzero polynomial in rr variables, of total degree at most p1p-1, is nonzero at at least pr1p^{r-1} points of Fpr\mathbb F_p^r.

Proof. Apply Lemma 6.3 to the full grid. If the leading exponents sum to at most p1p-1, then

i=1r(pai)pr1(piai)pr1. \prod_{i=1}^r(p-a_i)\ge p^{r-1}(p-\textstyle\sum_i a_i)\ge p^{r-1}.

The product inequality follows successively from (pa)(pb)p(pab)(p-a)(p-b)\ge p(p-a-b). \square

The lower bound when x0x_0 divides the polynomial

Lemma 6.5 (body footprint). If f=x0hf=x_0h, degfp1\deg f\le p-1, and fCd0f|_{C_d}\ne0, then

wtCd(f)2Sd2. \operatorname{wt}_{C_d}(f)\ge2S^{d-2}.

Proof. We have deghp2\deg h\le p-2, and ff vanishes on PP. On BB, introduce

t=x0,yi=xi+x02. t=x_0,\qquad y_i=x_i+x_0^2.

This is a bijection from BB to

Γ=Fp×Qd1. \Gamma=\mathbb F_p^*\times\mathcal Q^{d-1}.

Give tt weight 11 and every yiy_i weight 22. After substituting xi=yit2x_i=y_i-t^2, the transformed polynomial HH has weighted degree at most 2(p2)=4S62(p-2)=4S-6. Reduce it by

tp11,qQ(yiq). t^{p-1}-1,\qquad \prod_{q\in\mathcal Q}(y_i-q).

These reductions preserve its values on Γ\Gamma and do not increase weighted degree. The normal form is nonzero: f=x0hf=x_0h vanishes on PP, so fCd0f|_{C_d}\ne0 forces hh, and hence HH, to be nonzero somewhere on Γ\Gamma. Its leading monomial in the same lexicographic order as Lemma 6.3 is

ta0y1a1yd1ad1 t^{a_0}y_1^{a_1}\cdots y_{d-1}^{a_{d-1}}

with

0a02S3,0aiS1,a0+2i=1d1ai4S6. 0\le a_0\le2S-3,\quad 0\le a_i\le S-1, \quad a_0+2\sum_{i=1}^{d-1}a_i\le4S-6.

Lemma 6.3 gives at least

(2S2a0)i=1d1(Sai) (2S-2-a_0)\prod_{i=1}^{d-1}(S-a_i)

nonzero body values. We now minimize this product under the displayed constraints. At fixed i1ai\sum_{i\ge1}a_i, moving two transverse exponents farther apart cannot increase (Sai)(Saj)(S-a_i)(S-a_j). Thus at a minimum all but at most one transverse exponent lie in {0,S1}\{0,S-1\}. The degree budget permits at most one exponent equal to S1S-1.

If none equals S1S-1 and there is no interior exponent, the product is at least Sd12Sd2S^{d-1}\ge2S^{d-2}. If there is one interior exponent, maximizing it together with a0a_0 gives at least

12Sd2=2Sd2. 1\cdot2\cdot S^{d-2}=2S^{d-2}.

If one exponent equals S1S-1, write the possible other interior exponent as g[0,S2]g\in[0,S-2]. The budget gives a0+2g2S4a_0+2g\le2S-4, and the product is at least

(2+2g)(Sg)Sd3=(2S+2g(S1g))Sd32Sd2. (2+2g)(S-g)S^{d-3} =\bigl(2S+2g(S-1-g)\bigr)S^{d-3} \ge2S^{d-2}.

Therefore ff has at least 2Sd22S^{d-2} nonzero values. \square

Proof of Theorem 6.1 for d4d\ge4. For d4d\ge4, the other case is immediate. If x0fx_0\nmid f, then f(0,x1,,xd1)f(0,x_1,\ldots,x_{d-1}) is a nonzero polynomial of degree at most p1p-1. Lemma 6.4 gives at least pd2p^{d-2} nonzero values on PP, and

pd22Sd2(d4). p^{d-2}\ge2S^{d-2}\qquad(d\ge4).

Indeed 2p/(p+1)3/22p/(p+1)\ge3/2, and its (d2)(d-2)-nd power is at least 9/4>29/4>2. Lemmas 6.2 and 6.5 now prove Theorem 6.1 for every d4d\ge4. \square

Dimension three is the sole delicate case: Lemma 6.4 gives only pp plane points, one short of p+1p+1. The rest of this section closes that gap.

The one-coordinate shadow code

Assume henceforth that p=2n+15p=2n+1\ge5, so S=n+1S=n+1. Let H=Q{0}H=\mathcal Q\setminus\{0\}, the multiplicative group of nonzero squares.

The elementary finite-field facts used below are as follows. The squaring map on Fp\mathbb F_p^* has kernel {1,1}\{1,-1\}, so H=n|H|=n and Q=n+1|\mathcal Q|=n+1. If gg generates the cyclic group Fp\mathbb F_p^*, then (gk)n=(1)k(g^k)^n=(-1)^k. Thus sn=1s^n=1 for sHs\in H, while rn=1r^n=-1 for a nonsquare rr. It follows in particular that yn+1=yy^{n+1}=y for every yQy\in\mathcal Q.

Give Fp[τ,x]\mathbb F_p[\tau,x] the weights

wdeg(τ)=2,wdeg(x)=1, \operatorname{wdeg}(\tau)=2,\qquad\operatorname{wdeg}(x)=1,

and set

D={(τ,x):τH, x+τQ}. D=\{(\tau,x):\tau\in H,\ x+\tau\in\mathcal Q\}.

Define the shadow code UFpFpU\subseteq\mathbb F_p^{\mathbb F_p} by

U={xF(0,x):FD=0, wdegF2n}. U=\{x\mapsto F(0,x):F|_D=0,\ \operatorname{wdeg}F\le2n\}.

For any linear code LL, write d(L)d(L) for the least Hamming weight of a nonzero word of LL.

The next lemma determines its dimension without assuming any unproved interpolation statement.

Lemma 6.6 (triangle unisolvence). The monomials

T={τjxa:2j+a2n1} \mathcal T=\{\tau^jx^a:2j+a\le2n-1\}

form an interpolation basis on DD. Consequently, dimU=n+1=S\dim U=n+1=S.

Proof. Both T\mathcal T and DD have n(n+1)n(n+1) elements. Put y=x+τy=x+\tau, so D=H×QD=H\times\mathcal Q. Values on the n+1n+1 distinct points of Q\mathcal Q and coefficients in the basis 1,y,,yn1,y,\ldots,y^n are related by an invertible Vandermonde matrix. We may therefore work with those coefficients.

We will also use the multiplicative Fourier transform on HH. For a function g:HFpg:H\to\mathbb F_p, set

g^(r)=τHg(τ)τr,0r<n. \widehat g(r)=\sum_{\tau\in H}g(\tau)\tau^{-r}, \qquad0\le r<n.

It is invertible, with

g(τ)=1nr=0n1g^(r)τr. g(\tau)=\frac1n\sum_{r=0}^{n-1}\widehat g(r)\tau^r.

Indeed, after writing HH as a cyclic group, the finite geometric sum τHτm\sum_{\tau\in H}\tau^m equals nn when nmn\mid m and 00 otherwise. Substitution in the proposed inverse leaves only the term indexed by the same residue.

Here is the algebra producing the Fourier blocks. For a column τjxa\tau^jx^a,

τj(yτ)a=k=0a(1)ak(ak)ykτj+ak. \tau^j(y-\tau)^a =\sum_{k=0}^a(-1)^{a-k}\binom ak y^k\tau^{j+a-k}.

The frequency-rr Fourier component retains exactly the terms satisfying

kj+ar(modn). k\equiv j+a-r\pmod n.

Let cc be the residue of j+aj+a; equivalently, a transformed entry can be nonzero only when

cj+ar+d(modn). c\equiv j+a\equiv r+d\pmod n.

Write j+a=c+tnj+a=c+tn. The degree constraint forces t{0,1}t\in\{0,1\}. Since j=c+tnaj=c+tn-a, the two inequalities j0j\ge0 and 2j+a2n12j+a\le2n-1 give

t=0:0ac,t=1:2c+1an+c. t=0:\quad0\le a\le c, \qquad t=1:\quad2c+1\le a\le n+c.

This derives the row set below. Since a2n1a\le2n-1, each nonzero residue 1dn11\le d\le n-1 has at most the two exponents k=d,d+nk=d,d+n, and reduction by yn+1=yy^{n+1}=y identifies yd+ny^{d+n} with ydy^d on Q\mathcal Q. At residue zero, the exponents k=0,nk=0,n remain the distinct columns 1,yn1,y^n, because yn1y^n\ne1 at y=0y=0. Multiply the monomial side by (1)a(-1)^a and the coefficient side by (1)d(-1)^d, then apply the determinant-one shear

C0C0+(1)nCn. C_0\longleftarrow C_0+(-1)^nC_n.

The resulting uniform block indexed by c{0,,n1}c\in\{0,\ldots,n-1\} has columns d=0,,nd=0,\ldots,n, rows

Ac={0,,c}{2c+1,,n+c}, A_c=\{0,\ldots,c\}\cup\{2c+1,\ldots,n+c\},

and entries

Bc(a,d)=(ad)+(1)n(ad+n), B_c(a,d)=\binom ad+(-1)^n\binom a{d+n},

where (a2n)=0\binom a{2n}=0 in the column d=nd=n.

For every affected row a=n+ra=n+r, the finite binomial identity

(n+rd+n)=k=0r(n+rk1n1)(kd) \binom{n+r}{d+n} =\sum_{k=0}^{r}\binom{n+r-k-1}{n-1}\binom{k}{d}

follows by taking the coefficient of zrz^r in

(1z)nzd(1z)d+1=zd(1z)n+d+1. (1-z)^{-n}\frac{z^d}{(1-z)^{d+1}} =\frac{z^d}{(1-z)^{n+d+1}}.

The convolution on the left gives the displayed sum, while the coefficient on the right is (n+rrd)=(n+rn+d)\binom{n+r}{r-d}=\binom{n+r}{n+d}. The identity therefore allows the determinant-one row operation

Rn+rRn+r(1)nk=0r(n+rk1n1)Rk R_{n+r}\longleftarrow R_{n+r} -(-1)^n\sum_{k=0}^{r} \binom{n+r-k-1}{n-1}R_k

to remove the second term. The required lower rows 0,,r0,\ldots,r all belong to AcA_c. Thus

detBc=det[(ad)]aAc,0dn=a<a(aa)d=0nd!. \det B_c =\det\left[\binom ad\right]_{a\in A_c,\,0\le d\le n} =\frac{\prod_{a<a'}(a'-a)}{\prod_{d=0}^{n}d!}.

The last equality is the Vandermonde determinant, since (xd)\binom xd is a degree-dd polynomial with leading coefficient 1/d!1/d!. Every nonzero difference between nodes of Ac{0,,2n1}A_c\subseteq\{0,\ldots,2n-1\}, and every denominator factor, has magnitude less than p=2n+1p=2n+1. Hence every block determinant is nonzero in Fp\mathbb F_p, proving unisolvence.

The space of monomials of weighted degree at most 2n2n has dimension (n+1)2(n+1)^2, while D=n(n+1)|D|=n(n+1). Since T\mathcal T already evaluates bijectively, the kernel at weighted degree 2n2n has dimension n+1n+1. The slice map FF(0,x)F\mapsto F(0,x) is injective on that kernel: if the slice is zero as a function on Fp\mathbb F_p, then its degree is at most 2n=p12n=p-1, so one-variable root counting makes it the zero polynomial. Hence F=τGF=\tau G. Because τ0\tau\ne0 on DD, GG also vanishes on DD, but its weighted degree is at most 2n22n-2, contradicting unisolvence unless G=0G=0. Its image is UU, so dimU=n+1\dim U=n+1. \square

An explicit description of UU

Define functions on Fp\mathbb F_p

h(x)=x221xn+,1n, h_\ell(x)=x^\ell-\frac{2\ell}{2\ell-1}x^{n+\ell}, \qquad1\le\ell\le n,

and let

V=span{1,h1,,hn}. V=\operatorname{span}\{1,h_1,\ldots,h_n\}.

The denominators are nonzero because 121p21\le2\ell-1\le p-2. On HH, h(s)=s/(21)h_\ell(s)=-s^\ell/(2\ell-1). If a linear combination of the displayed generators vanishes, evaluation at 00 first kills its constant coefficient; evaluation on the cyclic group HH then kills the coefficients of the nn distinct characters ss^\ell, 1n1\le\ell\le n. Thus dimV=n+1\dim V=n+1.

Lemma 6.7 (explicit lifts). Every generator of VV belongs to UU. Hence U=VU=V.

Proof. We give the lift, including the coefficient identity that makes it work. For 1n1\le\ell\le n, put

H0(r)=(1+r)n+1(1+r) H_0(r)=(1+r)^{n+1}-(1+r)

and

G(r)=j=01(1)j(1)j(n+1)j(12)jj!rj, G_\ell(r)=\sum_{j=0}^{\ell-1} \frac{(-1)^j(1-\ell)_j(n-\ell+1)_j} {(1-2\ell)_j\,j!}r^j,

where (a)j=a(a+1)(a+j1)(a)_j=a(a+1)\cdots(a+j-1) and (a)0=1(a)_0=1. All denominators are products of nonzero integers of magnitude below pp.

Write H0G=afaraH_0G_\ell=\sum_a f_ar^a. We now perform the coefficient extraction explicitly. We first need the following terminating Pfaff identity: for an integer L0L\ge0,

j=0L(L)j(b)j(c)jj!zj=(1z)Lj=0L(L)j(cb)j(c)jj!(zz1)j. \sum_{j=0}^{L}\frac{(-L)_j(b)_j}{(c)_j\,j!}z^j =(1-z)^L\sum_{j=0}^{L} \frac{(-L)_j(c-b)_j}{(c)_j\,j!} \left(\frac{z}{z-1}\right)^j.

This is a polynomial identity, so the apparent denominator on the right cancels. Indeed, the coefficient of zmz^m on its right side is

(1)mj=0m(L)j(cb)j(c)jj!(Ljmj)=(L)m(b)m(c)mm!; (-1)^m\sum_{j=0}^{m} \frac{(-L)_j(c-b)_j}{(c)_j\,j!} \binom{L-j}{m-j} =\frac{(-L)_m(b)_m}{(c)_m\,m!};

after cancelling (L)m/m!(-L)_m/m!, the equality is the finite binomial identity proved in the next paragraph, with the parameters relabelled. Thus the coefficients agree term by term.

Apply this identity with L=1L=\ell-1, b=n+1b=n-\ell+1, c=12c=1-2\ell, and z=rz=-r. It gives

G(r)=(1+r)1j=01(1)j(n)j(12)jj!(r1+r)j. G_\ell(r) =(1+r)^{\ell-1} \sum_{j=0}^{\ell-1} \frac{(1-\ell)_j(-n-\ell)_j} {(1-2\ell)_j\,j!} \left(\frac r{1+r}\right)^j.

For a+1a\ge\ell+1, the polynomial (1+r)G(1+r)G_\ell has degree at most \ell, so it contributes nothing to faf_a. Consequently,

fa=[ra](1+r)n+1G(r)=j=01(1)j(n)j(12)jj!(n+jaj)=(n+a)j=01(1)j(a)j(12)jj!, \begin{aligned} f_a=[r^a](1+r)^{n+1}G_\ell(r) &=\sum_{j=0}^{\ell-1} \frac{(1-\ell)_j(-n-\ell)_j} {(1-2\ell)_j\,j!} \binom{n+\ell-j}{a-j}\\ &=\binom{n+\ell}{a} \sum_{j=0}^{\ell-1} \frac{(1-\ell)_j(-a)_j} {(1-2\ell)_j\,j!}, \end{aligned}

where the second equality uses

(n+jaj)=(n+a)(a)j(n)j. \binom{n+\ell-j}{a-j} =\binom{n+\ell}{a}\frac{(-a)_j}{(-n-\ell)_j}.

The temporary denominator is legitimate: for 0j10\le j\le\ell-1, none of the factors in (n)j(-n-\ell)_j is zero modulo p=2n+1p=2n+1.

It remains to evaluate the finite sum. The required binomial identity is

j=01(1)j(a)j(12)jj!=(12+a)1(12)1. \sum_{j=0}^{\ell-1} \frac{(1-\ell)_j(-a)_j}{(1-2\ell)_j\,j!} =\frac{(1-2\ell+a)_{\ell-1}}{(1-2\ell)_{\ell-1}}\,.

It gives, for a+1a\ge\ell+1,

fa=(n+a)(12+a)1(12)1. f_a=\binom{n+\ell}{a} \frac{(1-2\ell+a)_{\ell-1}}{(1-2\ell)_{\ell-1}}.

For an integer L0L\ge0,

j=0L(L)j(b)j(c)jj!=(cb)L(c)L. \sum_{j=0}^{L}\frac{(-L)_j(b)_j}{(c)_j\,j!} =\frac{(c-b)_L}{(c)_L}.

It follows by induction on LL: split (L+1j)=(Lj)+(Lj1)\binom{L+1}{j}=\binom Lj+\binom L{j-1}, apply the induction hypothesis to the two resulting sums, and simplify

(cb)L(c)Lbc(cb)L(c+1)L=(cb)L+1(c)L+1. \frac{(c-b)_L}{(c)_L} -\frac bc\,\frac{(c-b)_L}{(c+1)_L} =\frac{(c-b)_{L+1}}{(c)_{L+1}}.

Taking L=1L=\ell-1, b=ab=-a, and c=12c=1-2\ell proves the displayed coefficient formula. For <a<2\ell<a<2\ell, its numerator contains a zero factor, so fa=0f_a=0. At the endpoints, direct extraction gives

fn+f=221. \frac{f_{n+\ell}}{f_\ell}=-\frac{2\ell}{2\ell-1}.

Let P=(n+1)1P=(n-\ell+1)_{\ell-1} and E=(1)1(1)!(n+)E=(-1)^{\ell-1}(\ell-1)!\binom{n+\ell}{\ell}. Before normalization, fn+=P/(12)1f_{n+\ell}=P/(1-2\ell)_{\ell-1} and f=(EP)/(12)1f_\ell=(E-P)/(1-2\ell)_{\ell-1}. Since n=1/2n=-1/2 in Fp\mathbb F_p,

P=(1)135(21)21,E=(1)113(21)2, P=\frac{(-1)^{\ell-1}3\cdot5\cdots(2\ell-1)}{2^{\ell-1}}, \qquad E=\frac{(-1)^{\ell-1}1\cdot3\cdots(2\ell-1)}{2^\ell\ell},

so P=2EP=2\ell E. Here and below an empty product is 11. All displayed factors in EE are nonzero modulo pp, and

f=EP(12)1=E(12)(12)10. f_\ell=\frac{E-P}{(1-2\ell)_{\ell-1}} =\frac{E(1-2\ell)}{(1-2\ell)_{\ell-1}}\ne0.

The endpoint ratio is therefore the one claimed. Rescale GG_\ell, and hence every faf_a, so that f=1f_\ell=1.

Now define

F(τ,x)=a=1faτaxa+a=2n+faτn+axa. F_\ell(\tau,x)= \sum_{a=1}^{\ell}f_a\tau^{\ell-a}x^a +\sum_{a=2\ell}^{n+\ell}f_a\tau^{n+\ell-a}x^a.

We have f0=0f_0=0, because H0(0)=0H_0(0)=0, and deg(H0G)n+\deg(H_0G_\ell)\le n+\ell. Thus the two displayed intervals, together with the vanished gap <a<2\ell<a<2\ell, account for every coefficient of H0GH_0G_\ell. The polynomial FF_\ell has weighted degree at most 2n2n. On DD, put w=x/τw=x/\tau. Since τn=1\tau^n=1, its two displayed diagonals combine to give

F(τ,x)=τ(H0G)(w). F_\ell(\tau,x)=\tau^\ell(H_0G_\ell)(w).

Moreover 1+w=(x+τ)/τQ1+w=(x+\tau)/\tau\in\mathcal Q, so (1+w)n+1=1+w(1+w)^{n+1}=1+w, and the right-hand side is zero. At τ=0\tau=0, only its two endpoint terms remain, giving F(0,x)=h(x)F_\ell(0,x)=h_\ell(x). Finally, 1τn1-\tau^n lifts the constant function. Thus VUV\subseteq U; equal dimensions give U=VU=V. \square

The one-coordinate distance

Every word of VV is determined by a polynomial AA of degree at most nn: its values are AA on {0}H\{0\}\cup H, and A4xAA-4xA' on the nonsquares N=FpHN=\mathbb F_p^*\setminus H. Indeed, on HH the two terms of hh_\ell combine with coefficient 1/(21)-1/(2\ell-1), whereas on NN they combine with coefficient (41)/(21)(4\ell-1)/(2\ell-1); multiplication of the square-side coefficient by 141-4\ell gives the nonsquare-side coefficient.

For clarity, we also derive the matrix below. The interpolation nodes {0}H\{0\}\cup H are the roots of x(xn1)x(x^n-1), so

A(x)=E(1xn)+sHusx(xn1)n(xs). A(x)=E(1-x^n)+\sum_{s\in H}u_s \frac{x(x^n-1)}{n(x-s)}.

Apply 14xd/dx1-4x\,d/dx and evaluate at rNr\in N. Using n=1/2n=-1/2 and rn=1r^n=-1, the first term becomes 4E4E, while the ss-th interpolation term becomes 8r(r+s)/(rs)28r(r+s)/(r-s)^2. Therefore

fN=Tu+4E1,Tr,s=8r(r+s)(rs)2, f|_N=Tu+4E\mathbf1,\qquad T_{r,s}=\frac{8r(r+s)}{(r-s)^2},

where E=A(0)E=A(0), u=(A(s))sHu=(A(s))_{s\in H}, rNr\in N, and sHs\in H. Here 1\mathbf1 denotes the all-ones vector indexed by NN.

Lemma 6.8 (distance of the shadow code). For every prime p5p\ge5, d(U)=d(V)3d(U)=d(V)\ge3.

Proof. The weight is

1E0+suppu+wt(Tu+4E1). \mathbf1_{E\ne0}+|\operatorname{supp}u| +\operatorname{wt}(Tu+4E\mathbf1).

Only inputs of weight one or two require attention.

If u=0u=0 and E0E\ne0, the last block is a nonzero constant, giving weight n+1n+1. If E=0E=0 and uu is supported at one square ss, its nonsquare values are multiples of

K(t)=8t(t+1)(t1)2,t=r/s, K(t)=\frac{8t(t+1)}{(t-1)^2},\qquad t=r/s,

which vanishes only at t=1t=-1. Thus the total weight is at least 1+(n1)=n1+(n-1)=n, sufficient for p7p\ge7; for p=5p=5, 1-1 is a square and no nonsquare ratio is lost, giving weight 33. If E0E\ne0 and uu is one-sparse, the equation K(t)=cK(t)=c is a nonzero polynomial of degree at most two and has at most two solutions. The two input coordinates plus at least n2n-2 output coordinates give total weight at least nn, sufficient for p7p\ge7.

Finally, if E=0E=0 and uu is supported at distinct s1,s2s_1,s_2, clearing denominators shows that its zero nonsquare coordinates are roots of the nonzero cubic

R(X)=u1(X+s1)(Xs2)2+u2(X+s2)(Xs1)2. R(X)=u_1(X+s_1)(X-s_2)^2 +u_2(X+s_2)(X-s_1)^2.

It is nonzero because R(s1)=2u1s1(s1s2)20R(s_1)=2u_1s_1(s_1-s_2)^2\ne0. It has at most three nonsquare roots, so the two input coordinates plus at least n3n-3 output coordinates give weight n13n-1\ge3 for p11p\ge11. The two smaller primes are closed directly, without a search: with rows indexed by nonsquares and columns by nonzero squares,

T5=(3443)(H=(1,4), N=(2,3)), T_{5}=\begin{pmatrix}3&4\\4&3\end{pmatrix} \quad\bigl(H=(1,4),\ N=(2,3)\bigr),

and

T7=(310103031)(H=(1,2,4), N=(3,5,6)). T_{7}=\begin{pmatrix}3&1&0\\1&0&3\\0&3&1\end{pmatrix} \quad\bigl(H=(1,2,4),\ N=(3,5,6)\bigr).

For p=5p=5, each nonconstant column settles the E0E\ne0, one-sparse-uu case, and detT50\det T_5\ne0 settles the two-sparse case. For p=7p=7, each pair of columns is independent, settling the remaining two-sparse case. These observations close precisely the cases not covered by the general root counts. Thus every nonzero word has weight at least three. \square

From one coordinate to the whole plane

Let WW be the vector space of total-degree-p1=2n\le p-1=2n polynomials in x0,x1,x2x_0,x_1,x_2 that vanish on the body BB, and let

ρ(f)(x1,x2)=f(0,x1,x2). \rho(f)(x_1,x_2)=f(0,x_1,x_2).

Lemma 6.9 (body-kernel decomposition).

ρ(W)=U(x1)+U(x2). \rho(W)=U(x_1)+U(x_2).

Proof. The lifts in Lemma 6.7, embedded separately in x1x_1 and x2x_2, show U(x1)+U(x2)ρ(W)U(x_1)+U(x_2)\subseteq\rho(W). The two copies meet exactly in the constants, so this subspace has dimension 2(n+1)1=p2(n+1)-1=p. It remains to prove dimWp\dim W\le p.

Write τ=x02\tau=x_0^2 and decompose uniquely

f=E(τ,x1,x2)+x0O(τ,x1,x2), f=E(\tau,x_1,x_2)+x_0O(\tau,x_1,x_2),

where EE has weighted degree at most 2n2n and OO at most 2n12n-1, with weights (2,1,1)(2,1,1). Every τH\tau\in H has the two square roots ±x0\pm x_0; because pp is odd, vanishing at both signs separates the even and odd parts. We use the following rank claim.

Weighted-grid rank claim. Evaluation on

{(τ,x1,x2):τH, xi+τQ} \{(\tau,x_1,x_2):\tau\in H,\ x_i+\tau\in\mathcal Q\}

is injective through weighted degree 2n12n-1, and its kernel through weighted degree 2n2n has dimension at most 2n+1=p2n+1=p.

To prove the claim, put zi=(xi+τ)/τQz_i=(x_i+\tau)/\tau\in\mathcal Q and wi=zi1w_i=z_i-1. Let MaM_a be the coefficient vector of the remainder of waw^a modulo

Qn(w)=(w+1)n+1(w+1). Q_n(w)=(w+1)^{n+1}-(w+1).

Thus Ma=eaM_a=e_a for 0an0\le a\le n, where eae_a is the aa-th standard coordinate vector. Evaluation on the n+1n+1 distinct roots of QnQ_n and the coefficient vector of the remainder modulo QnQ_n are related by an invertible Vandermonde matrix. A monomial τjx1ax2b\tau^jx_1^ax_2^b evaluates as τj+a+bw1aw2b\tau^{j+a+b}w_1^aw_2^b. The invertible Fourier transform on HH splits the rows by the residue ej+a+b(modn)e\equiv j+a+b\pmod n, and the row in that block is MaMbM_a\otimes M_b.

For fixed e{0,,n1}e\in\{0,\ldots,n-1\}, put q=a+bq=a+b and write j+q=e+tnj+q=e+tn. Below weighted degree 2n2n, the inequalities j0j\ge0 and 2j+q2n12j+q\le2n-1 give

t=0:0qe,t=1:2e+1qn+e. t=0:\quad0\le q\le e, \qquad t=1:\quad2e+1\le q\le n+e.

No other tt is possible. Thus the totals a+ba+b arising below weight 2n2n are exactly

Ae={0,,e}{2e+1,,n+e}. A_e=\{0,\ldots,e\}\cup\{2e+1,\ldots,n+e\}.

At weighted degree exactly 2n2n, the additional rows with a,b>0a,b>0 have total q=2eq=2e for e>0e>0, and q=2nq=2n for e=0e=0. We add these mixed top-weight rows. Rows with a,bna,b\le n are coordinate vectors eaebe_a\otimes e_b and pivot on their matching columns; call them core rows. Every remaining row has exactly one exponent larger than nn. Call it horizontal when it is (a,n+δ)(a,n+\delta), grouping rows with the same low exponent aa into a horizontal fiber; call its transpose vertical.

We now enumerate every non-core row. When e=0e=0, only the added total 2n2n can be non-core, so a+n+δ=2na+n+\delta=2n gives Δa={na}\Delta_a=\{n-a\} for 1an11\le a\le n-1. When e>0e>0, the upper bound qn+eq\le n+e gives δea\delta\le e-a. The ordinary high totals begin at 2e+12e+1, while the added mixed total is 2e2e. The latter is unavailable when a=0a=0. Therefore the possible values of δ\delta are exactly the sets displayed next.

For e=0e=0, a horizontal fiber has Δa={na}\Delta_a=\{n-a\}, 1an11\le a\le n-1, and we select the column (a,n)(a,n). For e>0e>0, its possible δ\delta’s are

Δ0={max(1,2e+1n),,e}, \Delta_0=\{\max(1,2e+1-n),\ldots,e\},

and, for 1ae11\le a\le e-1,

Δa={max(1,2ena),,ea}. \Delta_a=\{\max(1,2e-n-a),\ldots,e-a\}.

After using the core pivots, a column (a,r)(a,r) is still unused precisely when its total a+ra+r is not one of the core-row totals. For a=0a=0, this leaves totals e+1,,2ee+1,\ldots,2e; for a>0a>0, the mixed top total 2e2e is already used and this leaves e+1,,2e1e+1,\ldots,2e-1. To keep horizontal selections disjoint from their vertical transposes, also require r>ar>a. Hence the eligible second-coordinate interval GaG_a is

{e+1,,min(n,2e)}(a=0) \{e+1,\ldots,\min(n,2e)\}\quad(a=0)

or

{max(a+1,e+1a),,min(n,2e1a)}(a>0). \{\max(a+1,e+1-a),\ldots,\min(n,2e-1-a)\}\quad(a>0).

A direct comparison of endpoints gives GaΔa|G_a|\ge|\Delta_a|. Explicitly, Δ0=G0=min(e,ne)|\Delta_0|=|G_0|=\min(e,n-e). For a>0a>0, write L=max(1,2ena)L=\max(1,2e-n-a) and k=eaL+1=Δak=e-a-L+1=|\Delta_a|. If 2e1an2e-1-a\le n, then L=1L=1 and maxGak+1=e\max G_a-k+1=e; if 2e1a>n2e-1-a>n, then L=2enaL=2e-n-a and the same equality follows with maxGa=n\max G_a=n. Thus GaG_a contains at least kk elements. Choose RaR_a to be the Δa|\Delta_a| largest elements of GaG_a, and select the columns (a,r)(a,r) with rRar\in R_a. It follows that minRa=e\min R_a=e for a>0a>0, while minR0=e+1\min R_0=e+1. Hence every active low coordinate satisfies ae1<ra\le e-1<r, and also δ<r\delta<r. Vertical rows use the transposed columns (r,a)(r,a).

For each non-core row we selected exactly one column, and different same-orientation fibers have different fixed low coordinates. The selected horizontal and vertical sets are disjoint. When e=0e=0, every low coordinate and every δ=na\delta=n-a is less than the selected terminal coordinate nn. When e>0e>0, every low coordinate of every active fiber is at most e1e-1, whereas every selected terminal coordinate is at least ee. A horizontal row has first factor Ma=eaM_a=e_a, so its entry on every selected vertical column is zero. The same global inequality, with the coordinates transposed, kills the other cross block after eliminating the core pivots. Together with the core pivot columns, these choices therefore form a square block-triangular minor; it remains only to prove that each diagonal arm block is nonsingular.

For each remaining arm block,

Mn+δ[r]=(1)δ(n+δr)(n+δr1δ1)=(1)δδ(n+δδ)(nr)n+δr. M_{n+\delta}[r] =(-1)^\delta\binom{n+\delta}{r} \binom{n+\delta-r-1}{\delta-1} =\frac{(-1)^\delta\delta\binom{n+\delta}{\delta}\binom nr} {n+\delta-r}.

To obtain this identity, first divide wn+δw^{n+\delta} by (w+1)n+1(w+1)^{n+1}. Replacing that divisor by Qn(w)Q_n(w) changes the remainder only in degrees at most δ\delta, so coefficients with r>δr>\delta are unchanged. Taylor expansion at w=1w=-1, followed by

k=0m(1)k(Nk)=(1)m(N1m), \sum_{k=0}^{m}(-1)^k\binom Nk=(-1)^m\binom{N-1}{m},

gives the first expression; cancellation of factorials gives the second. The row scale (1)δδ(n+δδ)(-1)^\delta\delta\binom{n+\delta}{\delta} and column scale (nr)\binom nr are nonzero: 1δn11\le\delta\le n-1, n+δp2n+\delta\le p-2, and 0rn0\le r\le n. After these nonzero scalings the block is the Cauchy matrix 1/(δ+y)1/(\delta+y), y=nry=n-r. Its determinant is

i<j(δjδi)i<j(yjyi)i,j(δi+yj). \frac{\prod_{i<j}(\delta_j-\delta_i) \prod_{i<j}(y_j-y_i)}{\prod_{i,j}(\delta_i+y_j)}\,.

Here is a proof of the determinant formula. After multiplication by i,j(δi+yj)\prod_{i,j}(\delta_i+y_j), the determinant is a polynomial alternating separately in the δi\delta_i’s and in the yjy_j’s. It is therefore divisible by both Vandermonde products in the numerator. The degrees are equal, so the quotient is constant; comparing successive highest-degree coefficients in δm,δm1,\delta_m,\delta_{m-1},\ldots gives constant 11.

The determinant is nonzero: all labels are distinct and every denominator lies in {1,,n1}\{1,\ldots,n-1\}, hence is nonzero modulo pp. The selected minor is therefore nonsingular. This proves independence of all rows through weight 2n12n-1, even after adding the mixed top-weight rows.

The only top-weight monomials omitted from this independent augmented family have a=0a=0 or b=0b=0. There are (n+1)+(n+1)1=2n+1=p(n+1)+(n+1)-1=2n+1=p of them. Thus the odd part is injective and the even kernel has dimension at most pp, proving the weighted-grid rank claim and dimWp\dim W\le p.

We already exhibited a pp-dimensional subspace of ρ(W)\rho(W). (The intersection of the two copies of UU consists only of constants: a function of x1x_1 equal to a function of x2x_2 for all pairs must be constant.) Hence equality holds. \square

Lemma 6.10 (two-coordinate sum code). If X=N|X|=N and a linear code LFXL\subseteq\mathbb F^X over a finite field F\mathbb F contains the constant functions, then

d({a(x)+b(y):a,bL})=Nd(L). d(\{a(x)+b(y):a,b\in L\})=N\,d(L).

Proof. Let a(x)+b(y)a(x)+b(y) be a nonzero sum-code word. If aa is nonconstant, every fiber a1(t)a^{-1}(t) has at most Nd(L)N-d(L) elements because ata-t is a nonzero word of LL. Hence the number of zeros of a(x)+b(y)a(x)+b(y) is

ta1(t)b1(t)N(Nd(L)). \sum_t|a^{-1}(t)|\,|b^{-1}(-t)| \le N(N-d(L)).

If aa is constant, then a+ba+b is a nonzero word of LL, repeated NN times, and obeys the same lower bound. A minimum word depending only on xx attains Nd(L)Nd(L). \square

Lemmas 6.8–6.10 give

d(ρ(W))=pd(U)3p. d(\rho(W))=p\,d(U)\ge3p.

Completion of the dimension-three proof

Proof of Theorem 6.1 for d=3d=3. Let ff have degree at most p1p-1, with x0fx_0\nmid f. Its plane restriction g=f(0,x1,x2)g=f(0,x_1,x_2) is nonzero, so Lemma 6.4 gives wtP(g)p\operatorname{wt}_P(g)\ge p. If wtC3(f)p\operatorname{wt}_{C_3}(f)\le p, equality must hold and ff must vanish on all of BB. Thus fWf\in W, so gρ(W)g\in\rho(W); but every nonzero word in ρ(W)\rho(W) has weight at least 3p3p, a contradiction. Hence every such ff has weight at least p+1p+1. Together with Lemmas 6.2 and 6.5, this proves

μ(C3)=p+1(p5). \mu(C_3)=p+1\qquad(p\ge5).

At p=3p=3, Lemma 6.5 gives weight at least 44 when x0fx_0\mid f, while Lemma 6.4 gives at least 33 nonzero plane values when x0fx_0\nmid f. The polynomial x1(x1+1)x_1(x_1+1) vanishes on the body because there x1{0,1}x_1\in\{0,-1\}, and on the plane it is nonzero exactly when x1=1x_1=1, at the three choices of x2x_2. Its weight is therefore 33, proving μ(C3)=3\mu(C_3)=3. \square

The next cell checks the finite-field linear algebra used in this proof. It verifies the shadow code, the body-kernel decomposition, the explicit support, and the exceptional p=3p=3 case at the stated sample primes.

Show code

Code cell 23 · In [8]

def gf_rref(matrix, p):
    A = np.array(matrix, dtype=np.int64) % p
    rows, cols = A.shape
    pivots = []
    r = 0
    for c in range(cols):
        pivot = next((i for i in range(r, rows) if A[i, c] % p), None)
        if pivot is None:
            continue
        if pivot != r:
            A[[r, pivot]] = A[[pivot, r]]
        A[r] = A[r] * pow(int(A[r, c]), -1, p) % p
        for i in range(rows):
            if i != r and A[i, c] % p:
                A[i] = (A[i] - int(A[i, c]) * A[r]) % p
        pivots.append(c)
        r += 1
        if r == rows:
            break
    return A, pivots


def gf_rank(matrix, p):
    A = np.array(matrix, dtype=np.int64)
    if A.size == 0:
        return 0
    return len(gf_rref(A, p)[1])


def gf_nullspace(matrix, p):
    R, pivots = gf_rref(matrix, p)
    cols = R.shape[1]
    free = [c for c in range(cols) if c not in pivots]
    basis = []
    for c in free:
        v = np.zeros(cols, dtype=np.int64)
        v[c] = 1
        for row, pivot in enumerate(pivots):
            v[pivot] = -R[row, c] % p
        basis.append(v)
    return np.array(basis, dtype=np.int64).reshape(len(free), cols)


def row_spaces_equal(A, B, p):
    A = np.array(A, dtype=np.int64) % p
    B = np.array(B, dtype=np.int64) % p
    return gf_rank(A, p) == gf_rank(B, p) == gf_rank(np.vstack([A, B]), p)


def total_degree_exponents(variables, degree):
    return [
        exponent
        for exponent in product(range(degree + 1), repeat=variables)
        if sum(exponent) <= degree
    ]


def evaluate_monomials(points, exponents, p):
    return np.array([
        [
            np.prod([pow(int(x), int(a), p) for x, a in zip(point, exponent)],
                    dtype=np.int64) % p
            for exponent in exponents
        ]
        for point in points
    ], dtype=np.int64)


def pochhammer_mod(a, length, p):
    out = 1
    for j in range(length):
        out = out * (a + j) % p
    return out


def shadow_lift_coefficients(p, ell):
    n = (p - 1) // 2
    G = []
    for j in range(ell):
        numerator = ((-1) ** j
                     * pochhammer_mod(1 - ell, j, p)
                     * pochhammer_mod(n - ell + 1, j, p))
        denominator = pochhammer_mod(1 - 2 * ell, j, p)
        factorial = 1
        for factor in range(1, j + 1):
            factorial = factorial * factor % p
        denominator = denominator * factorial % p
        G.append(numerator * pow(int(denominator), -1, p) % p)

    H = [comb(n + 1, a) % p for a in range(n + 2)]
    H[0] = (H[0] - 1) % p
    H[1] = (H[1] - 1) % p
    raw = [0] * (len(H) + len(G) - 1)
    for i, left in enumerate(H):
        for j, right in enumerate(G):
            raw[i + j] = (raw[i + j] + left * right) % p
    scale = pow(int(raw[ell]), -1, p)
    return [(value * scale) % p for value in raw]


def shadow_code_audit(p, exhaustive=True):
    n = (p - 1) // 2
    Q0 = sorted(quadratic_residues(p))
    H = [q for q in Q0 if q]
    N = [x for x in range(1, p) if x not in H]

    D = [(tau, (q - tau) % p) for tau in H for q in Q0]
    weighted = [
        (j, a)
        for j in range(n + 1)
        for a in range(2 * n - 2 * j + 1)
    ]
    D_eval = evaluate_monomials(D, weighted, p)
    kernel = gf_nullspace(D_eval, p)
    slice_eval = evaluate_monomials([(0, x) for x in range(p)], weighted, p)
    U = kernel @ slice_eval.T % p
    triangle_columns = [
        index for index, (j, a) in enumerate(weighted) if 2 * j + a <= 2 * n - 1
    ]
    assert len(triangle_columns) == n * (n + 1)
    assert gf_rank(D_eval[:, triangle_columns], p) == n * (n + 1)
    assert gf_rank(U, p) == n + 1

    V = [np.ones(p, dtype=np.int64)]
    for ell in range(1, n + 1):
        coefficient = 2 * ell * pow(2 * ell - 1, -1, p) % p
        V.append(np.array([
            (pow(x, ell, p) - coefficient * pow(x, n + ell, p)) % p
            for x in range(p)
        ], dtype=np.int64))
    V = np.array(V, dtype=np.int64)
    assert gf_rank(V, p) == n + 1
    assert row_spaces_equal(U, V, p)

    for ell in range(1, n + 1):
        coeff = shadow_lift_coefficients(p, ell)
        assert all(coeff[a] == 0 for a in range(ell + 1, 2 * ell))
        assert coeff[n + ell] == (-2 * ell * pow(2 * ell - 1, -1, p)) % p
        for tau, x in D:
            value = 0
            for a in range(1, ell + 1):
                value += coeff[a] * pow(tau, ell - a, p) * pow(x, a, p)
            for a in range(2 * ell, n + ell + 1):
                value += coeff[a] * pow(tau, n + ell - a, p) * pow(x, a, p)
            assert value % p == 0

    minimum = None
    if exhaustive:
        minimum = p
        for coefficients in product(range(p), repeat=n + 1):
            if not any(coefficients):
                continue
            word = np.array(coefficients, dtype=np.int64) @ V % p
            minimum = min(minimum, int(np.count_nonzero(word)))
        assert minimum >= 3

    T = np.array([
        [8 * r * (r + s) * pow((r - s) % p, -2, p) % p for s in H]
        for r in N
    ], dtype=np.int64)
    if p == 5:
        assert np.array_equal(T, np.array([[3, 4], [4, 3]]))
    if p == 7:
        assert np.array_equal(T, np.array([[3, 1, 0], [1, 0, 3], [0, 3, 1]]))

    return U, minimum


def body_decomposition_audit(p, U):
    n = (p - 1) // 2
    Q0 = sorted(quadratic_residues(p))
    body = [
        (x0, (q1 - x0 * x0) % p, (q2 - x0 * x0) % p)
        for x0 in range(1, p)
        for q1, q2 in product(Q0, repeat=2)
    ]
    exponents = total_degree_exponents(3, 2 * n)
    body_eval = evaluate_monomials(body, exponents, p)
    W = gf_nullspace(body_eval, p)
    assert W.shape[0] == p

    plane = list(product(range(p), repeat=2))
    plane_eval = evaluate_monomials([(0, x1, x2) for x1, x2 in plane],
                                    exponents, p)
    rho = W @ plane_eval.T % p
    sum_code = np.array(
        [[word[x1] for x1, x2 in plane] for word in U]
        + [[word[x2] for x1, x2 in plane] for word in U],
        dtype=np.int64,
    )
    assert gf_rank(rho, p) == p
    assert gf_rank(sum_code, p) == p
    assert row_spaces_equal(rho, sum_code, p)
    return len(body), len(exponents), gf_rank(rho, p)


def certificate_support(p, d):
    Q0 = sorted(quadratic_residues(p))
    plane = [(0,) + tail for tail in product(range(p), repeat=d - 1)]
    body = [
        (x0,) + tuple((q - x0 * x0) % p for q in qs)
        for x0 in range(1, p)
        for qs in product(Q0, repeat=d - 1)
    ]
    support = set()
    excluded = {pow(2, -1, p), (-pow(2, -1, p)) % p}
    for point in plane + body:
        x0, last = point[0], point[-1]
        value = x0 * last % p
        for c in range(1, p):
            if c in excluded:
                continue
            value = value * (
                x0 + c * last - pow(4 * c % p, -1, p)
            ) % p
        if value:
            support.add(point)
    return support


audit_rows = []
for p in [5, 7]:
    U, dU = shadow_code_audit(p)
    body_size, polynomial_dimension, rho_dimension = body_decomposition_audit(p, U)
    support = certificate_support(p, 3)
    assert len(support) == p + 1
    audit_rows.append((p, dU, body_size, polynomial_dimension, rho_dimension,
                       len(support)))

for p in [11]:
    U, dU = shadow_code_audit(p, exhaustive=False)
    support = certificate_support(p, 3)
    assert len(support) == p + 1
    audit_rows.append((p, dU, 'not expanded', 'not expanded', 'not expanded',
                       len(support)))

for p in [3, 5, 7]:
    for d in [3, 4, 5]:
        S = (p + 1) // 2
        candidates = []
        for a0 in range(p - 1):
            for transverse in product(range(S), repeat=d - 1):
                if a0 + 2 * sum(transverse) <= 2 * (p - 2):
                    candidates.append(
                        (p - 1 - a0)
                        * np.prod([S - a for a in transverse], dtype=np.int64)
                    )
        assert min(candidates) == 2 * S ** (d - 2)

p3_cone = (
    [(0, x1, x2) for x1, x2 in product(range(3), repeat=2)]
    + [
        (x0, (q1 - x0 * x0) % 3, (q2 - x0 * x0) % 3)
        for x0 in [1, 2]
        for q1, q2 in product(sorted(quadratic_residues(3)), repeat=2)
    ]
)
assert sum((x1 * (x1 + 1)) % 3 != 0 for _, x1, _ in p3_cone) == 3

print('Polynomial-distance audits (p, d(U), |B|, polynomial dimension, '
      'dim rho(W), certificate weight):')
print(audit_rows)

Saved output 1

Polynomial-distance audits (p, d(U), |B|, polynomial dimension, dim rho(W), certificate weight):
[(5, 3, 36, 35, 5, 6), (7, 3, 96, 84, 7, 8), (11, None, 'not expanded', 'not expanded', 'not expanded', 12)]

Station Reference. The explicit certificate and body-footprint proof originate in Archives #124 and #126. The dimension-three proof is reconstructed from Archives #127, #128, #130, #132, #142, #145, and #147. These identifiers record the provenance of the proof above.

Related Work. The Cartesian-grid footprint bound and the full-grid Schwartz–Zippel bound are classical; see Geil and Høholdt (2000) and Dvir (2009). Here these bounds are combined with an explicit upper certificate and a separate dimension-three analysis to obtain equality. The exact distance formula for this quadratic-residue support appears to be new.

6.2 Exact inversion-map shear profile

We determine the image size of every shear xx1+txx\mapsto x^{-1}+tx and use it to compute the normalized second moment exactly. For p1(mod4)p\equiv1\pmod4, all nonzero shears have the same image size; for p3(mod4)p\equiv3\pmod4, their image sizes differ by one according to the quadratic character of tt. The proof reduces the image count to a quadratic discriminant calculation, revealing the precise role of the secant geometry of the inversion conic in this one-pole completion.

Let f(0)=0f(0)=0 and f(x)=x1f(x)=x^{-1} for x0x\ne0. Put Pt=im(f+tid)P_t=|\operatorname{im}(f+t\,\mathrm{id})|, S=(p+1)/2S=(p+1)/2, δt=PtS\delta_t=P_t-S, E=tδtE=\sum_t\delta_t, and, when E>0E>0,

ρ(f)=tδt22(S1)E. \rho(f)=\frac{\sum_t\delta_t^2}{2(S-1)E}.

Theorem 6.11.

ρ(x1)={(p+3)/(2(p1)),p3(mod4),1/2,p1(mod4). \rho(x^{-1})= \begin{cases} (p+3)/(2(p-1)),&p\equiv3\pmod4,\\ 1/2,&p\equiv1\pmod4. \end{cases}

Proof. At t=0t=0, inversion is a bijection, so P0=pP_0=p and δ0=S1\delta_0=S-1. For t0t\ne0 and x0x\ne0, a value vv is attained precisely when

tx2vx+1=0 tx^2-vx+1=0

has a solution. Its discriminant is v24tv^2-4t. The identity

vχ(v24t)=1 \sum_v\chi(v^2-4t)=-1

follows by counting the p1p-1 solutions of (vu)(v+u)=4t(v-u)(v+u)=4t. Therefore the image from nonzero xx has size (p+χ(t))/2(p+\chi(t))/2. The separately supplied value gt(0)=0g_t(0)=0 is already present exactly when x2=1/tx^2=-1/t is soluble.

If χ(1)=1\chi(-1)=-1, square tt gives δt=1\delta_t=1 and nonsquare tt gives δt=1\delta_t=-1. Hence E=S1E=S-1 and

tδt2=(S1)2+(p1)=(p1)(p+3)4. \sum_t\delta_t^2=(S-1)^2+(p-1)=\frac{(p-1)(p+3)}4.

If χ(1)=1\chi(-1)=1, every nonzero tt has δt=0\delta_t=0, so E=S1E=S-1 and tδt2=(S1)2\sum_t\delta_t^2=(S-1)^2. Substitution gives the two formulas. \square

The next cell evaluates the shear profile directly at every odd prime from 55 through 5353 and checks the two exact formulas according to pmod4p\bmod4.

Show code

Code cell 26 · In [9]

def inversion_rho(p: int):
    inv = [0] + [pow(x, -1, p) for x in range(1, p)]
    S = (p + 1) // 2
    deltas = []
    for t in range(p):
        P_t = len({(inv[x] + t * x) % p for x in range(p)})
        deltas.append(P_t - S)
    E = sum(deltas)
    return Fraction(sum(delta * delta for delta in deltas), 2 * (S - 1) * E), deltas


inversion_rows = []
for p in [q for q in range(5, 54, 2) if is_prime(q)]:
    rho, deltas = inversion_rho(p)
    expected = Fraction(p + 3, 2 * (p - 1)) if p % 4 == 3 else Fraction(1, 2)
    assert rho == expected
    inversion_rows.append((p, str(rho), deltas[0], sorted(set(deltas[1:]))))
print('Inversion checks (p, rho, delta_0, nonzero-slope delta values):')
print(inversion_rows)

Saved output 1

Inversion checks (p, rho, delta_0, nonzero-slope delta values):
[(5, '1/2', 2, [0]), (7, '5/6', 3, [-1, 1]), (11, '7/10', 5, [-1, 1]), (13, '1/2', 6, [0]), (17, '1/2', 8, [0]), (19, '11/18', 9, [-1, 1]), (23, '13/22', 11, [-1, 1]), (29, '1/2', 14, [0]), (31, '17/30', 15, [-1, 1]), (37, '1/2', 18, [0]), (41, '1/2', 20, [0]), (43, '23/42', 21, [-1, 1]), (47, '25/46', 23, [-1, 1]), (53, '1/2', 26, [0])]

Station Reference. The inversion-map theorem and its statistic are based on Archive #72 and Eval #967.

Related Work. The maps xtx+x1x\mapsto tx+x^{-1} are scalar multiples of the classical x+a/xx+a/x maps associated with Dickson polynomials; see Chou, Gómez-Calderón, and Mullen (1988). We combine their classical value-set structure with the shear statistic ρ\rho to derive the displayed closed formula.

7. Appendix verification

This section reconstructs the 25-entry literature comparison from explicit constructions and witness lines.

Write Sqp={u2:uFp}\operatorname{Sq}_p=\{u^2:u\in\mathbb F_p\} and S=(p+1)/2S=(p+1)/2.

7.1 Saraf–Sudan quadratic-residue construction

This is the quadratic-residue construction of Mockenhaupt and Tao (2004), in the form recorded by Saraf and Sudan (2008).

Theorem 7.1. For odd prime pp, define

Dd={(t,y1,,yd1):yi+t2Sqp for every i},H={0}×Fpd1. D_d=\{(t,y_1,\ldots,y_{d-1}):y_i+t^2\in\operatorname{Sq}_p\text{ for every }i\},\qquad H=\{0\}\times\mathbb F_p^{d-1}.

Then Cd=DdHC_d=D_d\cup H is Kakeya and

Cd=(p1)Sd1+pd1. |C_d|=(p-1)S^{d-1}+p^{d-1}.

Proof. The hyperplane HH contains a line in every direction whose first coordinate is zero. For direction (1,a1,,ad1)(1,a_1,\ldots,a_{d-1}), put bi=ai/2b_i=a_i/2; the line

t(t,a1t+b12,,ad1t+bd12) t\longmapsto(t,a_1t+b_1^2,\ldots,a_{d-1}t+b_{d-1}^2)

lies in DdD_d because ait+bi2+t2=(t+bi)2a_it+b_i^2+t^2=(t+b_i)^2. There are pSd1pS^{d-1} points in DdD_d, pd1p^{d-1} in HH, and Sd1S^{d-1} in their intersection. Inclusion–exclusion gives the formula. Thus the often-quoted sum pSd1+pd1pS^{d-1}+p^{d-1} is a valid upper bound but is not the exact size of the displayed union. \square

7.2 Bukh–Chao recursion with full translations

The recursion is due to Bukh and Chao (2021).

Let

Gd={(t,y1,,yd1):yi+t2/4Sqp for every i}. G_d=\{(t,y_1,\ldots,y_{d-1}):y_i+t^2/4\in\operatorname{Sq}_p\text{ for every }i\}.

Starting with R1=FpR_1=\mathbb F_p, recursively set

Rd=Gd((cd,hd)+({0}×Rd1)),(cd,hd)Fp×Fpd1.(7) R_d=G_d\cup\left((c_d,h_d)+(\{0\}\times R_{d-1})\right),\qquad(c_d,h_d)\in\mathbb F_p\times\mathbb F_p^{d-1}. \tag{7}

Every translation is allowed. The certificate uses two explicitly recorded special cases: complete transverse histories with every cd=0c_d=0, and slice-placement histories with every hd=0h_d=0. Together they improve twelve cells of the earlier column.

Theorem 7.2. Let d2d\ge2. For every sequence ((cj,hj))j=2d((c_j,h_j))_{j=2}^d with cjFpc_j\in\mathbb F_p and hjFpj1h_j\in\mathbb F_p^{j-1}, the set RdR_d defined by (7) is Kakeya.

Proof. A graph direction (1,a1,,ad1)(1,a_1,\ldots,a_{d-1}) has the line with anchor (0,a12,,ad12)(0,a_1^2,\ldots,a_{d-1}^2) in GdG_d, since

ai2+ait+t2/4=(ai+t/2)2. a_i^2+a_it+t^2/4=(a_i+t/2)^2.

By induction, a witness w+tvw+tv in Rd1R_{d-1} becomes (cd,hd+w)+t(0,v)(c_d,h_d+w)+t(0,v) in the translated slice. These are all remaining directions. Translation changes only the overlap, not this proof. \square

The code reads the complete parameters for every changed cell from the local certificate. Unchanged cells are regenerated by the earlier deterministic transverse recursion. In every case it materializes the point set and checks the displayed witnesses directly.

7.3 The KLSS image-stratum construction

This construction is due to Kopparty et al. (2011).

For f(x)=x2f(x)=x^2, put I(t)={x2+tx:xFp}I(t)=\{x^2+tx:x\in\mathbb F_p\}. Define

Ld=j=0d1{(x1,,xj,t,0,,0):xiI(t)}. L_d=\bigcup_{j=0}^{d-1}\{(x_1,\ldots,x_j,t,0,\ldots,0):x_i\in I(t)\}.

Theorem 7.3. The set LdL_d is Kakeya.

Proof. Let v=(v1,,vd)v=(v_1,\ldots,v_d) be nonzero and let jj be the last index for which vj0v_j\ne0. Along the line with jjth coordinate t=svjt=sv_j, choose, for i<ji<j,

xi(s)=(vivj)2+svi=(vivj)2+tvivjI(t). x_i(s)=\left(\frac{v_i}{v_j}\right)^2+s v_i=\left(\frac{v_i}{v_j}\right)^2+t\frac{v_i}{v_j}\in I(t).

All later coordinates are zero, so this line lies in the jjth stratum. \square

KLSS state the convenient upper bound pj<dSjp\sum_{j<d}S^j. The displayed strata can overlap, so the code also reports the smaller cardinality of their actual union.

7.4 The KLSS missing-digit construction

This construction is also from Kopparty et al. (2011).

Theorem 7.4. For every prime power qq and d1d\ge1,

Md=(Fq×)d{0,1}d M_d=(\mathbb F_q^\times)^d\cup\{0,1\}^d

is Kakeya and has (q1)d+2d1(q-1)^d+2^d-1 points.

Proof. For direction v=(vi)v=(v_i), take the anchor bi=1b_i=1 when vi=0v_i=0 and bi=0b_i=0 otherwise. At parameter zero, b{0,1}db\in\{0,1\}^d. At every nonzero parameter tt, every coordinate of b+tvb+tv is either 11 or the nonzero value tvitv_i, so b+tv(Fq×)db+tv\in(\mathbb F_q^\times)^d. Finally, the two parts intersect only at (1,,1)(1,\ldots,1), proving the size formula. \square

7.5 Sharp planar sets, products, and characteristic three

Proposition 7.5. In dimension two, Theorem 7.1 has size p(p+1)/2+(p1)/2p(p+1)/2+(p-1)/2, the sharp value proved by Blokhuis and Mazzocca (2008). If AFpaA\subseteq\mathbb F_p^a and BFpbB\subseteq\mathbb F_p^b are Kakeya, then A×BA\times B is Kakeya of size AB|A||B|.

Proof. Substitution of d=2d=2 into Theorem 7.1 gives the planar size. For a direction (u,v)(u,v), synchronize witness-line parameters in the two factors; if one component is zero, hold an arbitrary point of that factor constant. This gives a line in A×BA\times B, whose size is exactly the product. \square

The supplied p=3p=3 arrays of sizes 1313, 2727, and 5353 are independently checked below. Only 1313 is used as an exact minimum; 2727 and 5353 are upper bounds.

7.6 Complete 25-cell comparison

The next cell reconstructs every named literature family, uses the 25 Station certificate checks from S2, and checks every AlphaEvolve integer against the commit-pinned provenance map. It then performs the baseline comparison and computes the 14-improvement, 11-tie, 0-loss tally.

Show code

Code cell 29 · In [10]

def assert_explicit_witnesses(points, witnesses, p: int, d: int):
    point_set = {tuple(map(int, point)) for point in points}
    directions = set(canonical_directions(p, d))
    assert set(witnesses) == directions
    for direction, anchor in witnesses.items():
        assert line_points(anchor, direction, p) <= point_set
    return len(directions)


def classical_qr_cone_with_witnesses(p: int, d: int):
    Q = quadratic_residues(p)
    plane = {(0,) + tail for tail in product(range(p), repeat=d - 1)}
    body = {
        (t,) + tail
        for t in range(p)
        for tail in product(range(p), repeat=d - 1)
        if all((y + t * t) % p in Q for y in tail)
    }
    points = plane | body
    inv2 = pow(2, -1, p)
    witnesses = {}
    for direction in canonical_directions(p, d):
        if direction[0] == 0:
            witnesses[direction] = (0,) * d
        else:
            witnesses[direction] = (0,) + tuple(
                (a * inv2) ** 2 % p for a in direction[1:]
            )
    return points, witnesses


def quadratic_graph_block(p: int, d: int):
    slices = {
        t: {(a * a + t * a) % p for a in range(p)}
        for t in range(p)
    }
    points = {
        (t,) + tail
        for t in range(p)
        for tail in product(slices[t], repeat=d - 1)
    }
    witnesses = {
        (1,) + slopes: (0,) + tuple(a * a % p for a in slopes)
        for slopes in product(range(p), repeat=d - 1)
    }
    return points, witnesses


def transverse_extension(p, child, child_witnesses, shift):
    d = len(next(iter(child))) + 1
    graph, witnesses = quadratic_graph_block(p, d)
    embedded = {
        (0,) + tuple((x + h) % p for x, h in zip(point, shift))
        for point in child
    }
    witnesses.update({
        (0,) + direction:
        (0,) + tuple((x + h) % p for x, h in zip(anchor, shift))
        for direction, anchor in child_witnesses.items()
    })
    return graph | embedded, witnesses


def transverse_history_set(p: int, d: int, history):
    points = {(t,) for t in range(p)}
    witnesses = {(1,): (0,)}
    points, witnesses = transverse_extension(p, points, witnesses, (0,))
    assert len(history) == d - 2
    for n, shift in enumerate(history, start=3):
        assert len(shift) == n - 1
        points, witnesses = transverse_extension(
            p, points, witnesses, tuple(shift)
        )
    return points, witnesses


def slice_placement_set(p: int, d: int, placements):
    if d == 1:
        assert len(placements) == 0
        return {(t,) for t in range(p)}, {(1,): (0,)}
    assert len(placements) == d - 1
    c = int(placements[0])
    child, child_witnesses = slice_placement_set(p, d - 1, placements[1:])
    graph, witnesses = quadratic_graph_block(p, d)
    embedded = {(c,) + point for point in child}
    witnesses.update({
        (0,) + direction: (c,) + anchor
        for direction, anchor in child_witnesses.items()
    })
    return graph | embedded, witnesses


def greedy_transverse_set(p: int, d: int):
    if d == 1:
        return {(t,) for t in range(p)}, {(1,): (0,)}, []
    child, child_witnesses, history = greedy_transverse_set(p, d - 1)
    Q_grid = list(product(sorted(quadratic_residues(p)), repeat=d - 1))
    overlap_counts = Counter(
        tuple((s - x) % p for s, x in zip(square_point, child_point))
        for child_point in child
        for square_point in Q_grid
    )
    maximum = max(overlap_counts.values())
    shift = min(h for h, value in overlap_counts.items() if value == maximum)
    points, witnesses = transverse_extension(p, child, child_witnesses, shift)
    return points, witnesses, history + [shift]


def klss_image_strata(p: int, d: int):
    image = {
        t: {(x * x + t * x) % p for x in range(p)}
        for t in range(p)
    }
    points = set()
    for j in range(d):
        for t in range(p):
            for prefix in product(image[t], repeat=j):
                points.add(prefix + (t,) + (0,) * (d - j - 1))
    witnesses = {}
    for direction in canonical_directions(p, d):
        j = max(i for i, value in enumerate(direction) if value)
        inv_last = pow(int(direction[j]), -1, p)
        anchor = tuple(
            (int(direction[i]) * inv_last) ** 2 % p for i in range(j)
        ) + (0,) * (d - j)
        witnesses[direction] = anchor
    return points, witnesses


def missing_digit_set(p: int, d: int):
    nonzero = set(product(range(1, p), repeat=d))
    binary = set(product((0, 1), repeat=d))
    points = nonzero | binary
    witnesses = {
        direction: tuple(1 if value == 0 else 0 for value in direction)
        for direction in canonical_directions(p, d)
    }
    return points, witnesses


CERTIFICATE_HASHES = {
    'benchmark_25_cells.json':
        '286014dfc5d7899dd83b0d7ab102b0e5807236a72e34d1f4b7fe64bf049d78f6',
    'literature_nested_shift_improvements.json':
        '89311a8f76b855878a5a1f88a79c6ff34a1b4e87cb2636afadc5a70df09154c7',
    'station_point_certificates.json':
        'f19d7b03cfdd5956386f2ef4ef2aa41009783ba0b5137d8b54d0ec360b22afad',
    'alphaevolve_provenance.json':
        'eb140c6f3a164e2c3ac6489f8dd932d1757cd8f671f466dbbfef5122bf31dc0a',
}
certificate_data = {}
for filename, expected_hash in CERTIFICATE_HASHES.items():
    path = ARTIFACT_DIR / filename
    if not path.is_file():
        raise FileNotFoundError(f'Required certificate is missing: {path}')
    assert hashlib.sha256(path.read_bytes()).hexdigest() == expected_hash
    certificate_data[filename] = json.loads(path.read_text())

assert station_manifest['benchmark_sha256'] == CERTIFICATE_HASHES['benchmark_25_cells.json']
benchmark_rows = certificate_data['benchmark_25_cells.json']['rows']
assert len(benchmark_rows) == 25
benchmark = {
    (d, p): {'literature': literature, 'ae': ae, 'station': station}
    for d, p, literature, ae, station in benchmark_rows
}
improvement_rows = certificate_data[
    'literature_nested_shift_improvements.json'
]['improvements']
improvements = {(row['d'], row['p']): row for row in improvement_rows}
assert len(improvements) == 12



alphaevolve_provenance = certificate_data['alphaevolve_provenance.json']
alphaevolve_values = {(row['d'], row['p']): row['size']
                      for row in alphaevolve_provenance['records']}
assert len(alphaevolve_values) == 25
assert alphaevolve_values == {(d, p): ae for d, p, _, ae, _ in benchmark_rows}
assert alphaevolve_provenance['official_notebook'].startswith('https://github.com/google-deepmind/')


small_arrays = {}
for filename, p, d, size in [
    ('kakeya_F3_d3_13.npy', 3, 3, 13),
    ('kakeya_F3_d4_27.npy', 3, 4, 27),
    ('kakeya_F3_d5_53.npy', 3, 5, 53),
]:
    path = ARTIFACT_DIR / filename
    if not path.is_file():
        raise FileNotFoundError(f'Required certificate is missing: {path}')
    assert hashlib.sha256(path.read_bytes()).hexdigest() == DATA_HASHES[filename]
    array = np.load(path, allow_pickle=False)
    assert_point_certificate(array, p, d, size)
    small_arrays[d] = {tuple(map(int, row)) for row in array}


def selected_literature_set(d: int, p: int):
    if p == 3 and d in (3, 4):
        points = small_arrays[d]
        return points, find_kakeya_witnesses(points, p, d), 'finite p=3 certificate'
    if p == 3 and d == 5:
        points, witnesses = missing_digit_set(p, d)
        return points, witnesses, 'missing-digit construction'
    if (d, p) in improvements:
        row = improvements[(d, p)]
        if row['method'] == 'transverse_shift_history':
            points, witnesses = transverse_history_set(p, d, row['parameters'])
        elif row['method'] == 'slice_placement_history':
            points, witnesses = slice_placement_set(p, d, row['parameters'])
        else:
            raise AssertionError(f"Unknown method: {row['method']}")
        return points, witnesses, row['method']
    points, witnesses, _ = greedy_transverse_set(p, d)
    return points, witnesses, 'greedy transverse recursion (unchanged cell)'


literature_checks = []
for (d, p), reported in benchmark.items():
    points, witnesses, method = selected_literature_set(d, p)
    direction_count = assert_explicit_witnesses(points, witnesses, p, d)
    assert len(points) == reported['literature'], (
        d, p, len(points), reported['literature']
    )
    literature_checks.append((d, p, len(points), direction_count, method))


family_checks = []
for d, p, _, _, _ in benchmark_rows:
    classical, classical_witnesses = classical_qr_cone_with_witnesses(p, d)
    assert len(classical) == (p - 1) * ((p + 1) // 2) ** (d - 1) + p ** (d - 1)
    assert_explicit_witnesses(classical, classical_witnesses, p, d)

    strata, strata_witnesses = klss_image_strata(p, d)
    assert len(strata) <= p * sum(((p + 1) // 2) ** j for j in range(d))
    assert_explicit_witnesses(strata, strata_witnesses, p, d)

    missing, missing_witnesses = missing_digit_set(p, d)
    assert len(missing) == (p - 1) ** d + 2 ** d - 1
    assert_explicit_witnesses(missing, missing_witnesses, p, d)
    family_checks.append((d, p, len(classical), len(strata), len(missing)))


for p in sorted({p for _, p, _, _, _ in benchmark_rows}):
    planar, planar_witnesses = classical_qr_cone_with_witnesses(p, 2)
    assert len(planar) == p * (p + 1) // 2 + (p - 1) // 2
    assert_explicit_witnesses(planar, planar_witnesses, p, 2)

product_3_5 = {a + b for a in small_arrays[3] for b in small_arrays[5]}
assert len(product_3_5) == 13 * 53
assert len(find_kakeya_witnesses(product_3_5, 3, 8)) == (3**8 - 1) // 2


before_after = [
    (row['d'], row['p'], row['before'], row['after'])
    for row in improvement_rows
]
assert all(benchmark[(d, p)]['literature'] == after for d, p, _, after in before_after)
verdicts = []
for row in benchmark.values():
    prior = min(row['literature'], row['ae'])
    verdicts.append(prior - row['station'])
assert (
    sum(value > 0 for value in verdicts),
    sum(value == 0 for value in verdicts),
    sum(value < 0 for value in verdicts),
) == (14, 11, 0)

print('Verified selected literature sets (d, p, points, directions, method):')
for row in literature_checks:
    print(row)
print('Before/after literature-reference corrections:')
for row in before_after:
    print(row)
print('Verified all named uniform families at all 25 cells.')
print('Verified all 25 Station point sets and every supplied witness line.')
print('Matched all 25 AlphaEvolve integers to the commit-pinned official notebook map.')
print('Verified finite p=3 sizes:', {d: len(points) for d, points in small_arrays.items()})
print('Verified |K_3 x K_5| in dimension 8:', len(product_3_5))
print('Benchmark tally: 14 strictly better, 11 ties, 0 worse.')

Saved output 1

Verified selected literature sets (d, p, points, directions, method):
(3, 3, 13, 13, 'finite p=3 certificate')
(3, 5, 53, 31, 'greedy transverse recursion (unchanged cell)')
(3, 7, 129, 57, 'greedy transverse recursion (unchanged cell)')
(3, 11, 440, 133, 'greedy transverse recursion (unchanged cell)')
(3, 13, 698, 183, 'slice_placement_history')
(3, 19, 2034, 381, 'greedy transverse recursion (unchanged cell)')
(3, 23, 3509, 553, 'greedy transverse recursion (unchanged cell)')
(3, 29, 6837, 871, 'greedy transverse recursion (unchanged cell)')
(3, 31, 8295, 993, 'greedy transverse recursion (unchanged cell)')
(3, 37, 13866, 1407, 'slice_placement_history')
(3, 41, 18708, 1723, 'slice_placement_history')
(3, 43, 21504, 1893, 'greedy transverse recursion (unchanged cell)')
(3, 47, 27899, 2257, 'greedy transverse recursion (unchanged cell)')
(3, 53, 39686, 2863, 'slice_placement_history')
(4, 3, 27, 40, 'finite p=3 certificate')
(4, 5, 163, 156, 'transverse_shift_history')
(4, 7, 528, 400, 'transverse_shift_history')
(4, 11, 2689, 1464, 'greedy transverse recursion (unchanged cell)')
(4, 13, 4972, 2380, 'transverse_shift_history')
(4, 17, 13521, 5220, 'slice_placement_history')
(4, 19, 20586, 7240, 'transverse_shift_history')
(5, 3, 63, 121, 'missing-digit construction')
(5, 5, 497, 781, 'transverse_shift_history')
(5, 7, 2142, 2801, 'transverse_shift_history')
(5, 11, 16307, 16105, 'transverse_shift_history')
Before/after literature-reference corrections:
(3, 13, 699, 698)
(3, 37, 13867, 13866)
(3, 41, 18709, 18708)
(3, 53, 39687, 39686)
(4, 5, 164, 163)
(4, 7, 529, 528)
(4, 13, 4973, 4972)
(4, 17, 13524, 13521)
(4, 19, 20593, 20586)
(5, 5, 503, 497)
(5, 7, 2145, 2142)
(5, 11, 16348, 16307)
Verified all named uniform families at all 25 cells.
Verified all 25 Station point sets and every supplied witness line.
Matched all 25 AlphaEvolve integers to the commit-pinned official notebook map.
Verified finite p=3 sizes: {3: 13, 4: 27, 5: 53}
Verified |K_3 x K_5| in dimension 8: 689
Benchmark tally: 14 strictly better, 11 ties, 0 worse.

Result and table provenance

For the AlphaEvolve column, alphaevolve_provenance.json maps each integer to the displayed constructor in the official commit-pinned notebook. The integers are obtained by executing the indicated constructor, deduplicating its returned points, and counting them; the notebook’s sub-unit random score jitter is not part of the cardinality. All 25 integers were independently replayed from that source.