This notebook is a self-contained mathematical and computational companion to the Discretized Kakeya needle section in the paper.
The task template, including the prompt given to the agents, is available here.
The outline is as follows:
- 1. Problem Definitions
- 2. S1. New upper bounds at .
- 3. S2. Exact optima at and symmetry breaking at .
- 4. Additional findings
The required packages are:
NumPy
1. Problem Definitions
The classical Kakeya needle problem asks how little area is needed to turn a unit line segment through every direction. The finite model below replaces the continuum of directions by equally spaced thin triangles that may slide horizontally.
For real offsets , define
At height , the -th triangle cuts the horizontal line in
We write and .
Lemma 1.1 (slice identity). If are the sorted left endpoints and , then, with ,
Proof. Scan the equal-width intervals from left to right. The first interval contributes . The next left endpoint enlarges the current union by its gap from the previous left endpoint if that gap is below , and by if it is at least . Summing these increments proves the identity.
Lemma 1.2 (exact event sweep). The intersections of the affine endpoint lines, together with , divide the strip into finitely many slabs. On each slab the slice union length is affine in ; hence the trapezoidal rule on consecutive event heights gives exactly.
Proof. Between consecutive endpoint intersections their total order does not change. Consequently the connected components of the interval union have fixed left and right endpoint labels. Each component length, and therefore their sum, is affine in . Integrating an affine function by its endpoint trapezoid is exact. Continuity makes coincident endpoints at slab boundaries harmless.
The construction data are stored in kakeya_needle_offsets.npz. The next cell imports the shared packages, authenticates this local bundle, and loads the seven tested arrays and exact small-case witnesses.
Show code
Code cell 3 · In [1]
from __future__ import annotations
import hashlib
import math
from decimal import Decimal, getcontext
from fractions import Fraction
from pathlib import Path
import numpy as np
getcontext().prec = 50
DATA_FILE = Path("kakeya_needle_offsets.npz")
EXPECTED_SHA256 = "ba68158d6e077aa55d3b5e8132936d2d006cc291a9989b37c05b607a21b4cc8e"
if not DATA_FILE.is_file():
raise FileNotFoundError(f"required artifact is missing: {DATA_FILE.resolve()}")
observed_hash = hashlib.sha256(DATA_FILE.read_bytes()).hexdigest()
if observed_hash != EXPECTED_SHA256:
raise RuntimeError(
f"artifact digest mismatch: expected {EXPECTED_SHA256}, observed {observed_hash}"
)
payload = np.load(DATA_FILE, allow_pickle=False)
expected_keys = {
"tested_n", "offsets_n2", "offsets_n4", "offsets_n8", "offsets_n16",
"offsets_n32", "offsets_n64", "offsets_n128", "n3_numerators",
"n3_denominator", "n4_plateau_c_numerators", "n4_plateau_c_denominators",
"n5_symmetric_numerators", "n5_symmetric_denominator",
"n5_asymmetric_numerators", "n5_asymmetric_denominator",
}
if set(payload.files) != expected_keys:
raise RuntimeError(f"unexpected payload schema: {sorted(payload.files)}")
TESTED_N = tuple(int(v) for v in payload["tested_n"])
assert TESTED_N == (2, 4, 8, 16, 32, 64, 128)
OFFSETS = {n: payload[f"offsets_n{n}"].copy() for n in TESTED_N}
for n, xs in OFFSETS.items():
assert xs.shape == (n,) and xs.dtype == np.float64 and np.all(np.isfinite(xs))
print("payload SHA-256", observed_hash)
print("tested sizes", TESTED_N)
Saved output 1
payload SHA-256 ba68158d6e077aa55d3b5e8132936d2d006cc291a9989b37c05b607a21b4cc8e tested sizes (2, 4, 8, 16, 32, 64, 128)
These shared routines are used by all verification cells below. The next cell defines the exact horizontal slices, union lengths, and event-height area sweep.
Show code
Code cell 5 · In [2]
def interval_at_height(offsets, j, y):
n = len(offsets)
x = offsets[j - 1]
return x + Fraction(j, n) * y, x + Fraction(1, n) + Fraction(j - 1, n) * y
def union_length_exact(offsets, y):
intervals = sorted(interval_at_height(offsets, j, y) for j in range(1, len(offsets) + 1))
total = Fraction(0)
left, right = intervals[0]
for next_left, next_right in intervals[1:]:
if next_left <= right:
right = max(right, next_right)
else:
total += right - left
left, right = next_left, next_right
return total + right - left
def area_exact(offsets):
offsets = tuple(Fraction(v) for v in offsets)
n = len(offsets)
endpoints = []
for j, x in enumerate(offsets, start=1):
endpoints.append((Fraction(j, n), x))
endpoints.append((Fraction(j - 1, n), x + Fraction(1, n)))
breaks = {Fraction(0), Fraction(1)}
for i, (slope_a, intercept_a) in enumerate(endpoints):
for slope_b, intercept_b in endpoints[i + 1:]:
if slope_a == slope_b:
continue
y = (intercept_b - intercept_a) / (slope_a - slope_b)
if 0 < y < 1:
breaks.add(y)
breaks = sorted(breaks)
return sum(
(b - a) * (union_length_exact(offsets, a) + union_length_exact(offsets, b)) / 2
for a, b in zip(breaks, breaks[1:])
)
def decimal_of(q, digits=18):
return Decimal(q.numerator) / Decimal(q.denominator)
assert area_exact([Fraction(2, 9), Fraction(1, 9), 0]) == Fraction(5, 18)
assert area_exact([Fraction(1, 4), Fraction(1, 6), Fraction(1, 12), 0]) == Fraction(1, 4)
print("exact sweep elementary checks passed")
Saved output 1
exact sweep elementary checks passed
2. S1. New upper bounds at .
2.1 Explicit finite upper bounds
Theorem 2.1 (explicit finite upper bounds).
Explicit configurations witness the three bounds; their exact binary-rational areas begin
and
Relative to the separately optimized finite constructions published by AlphaEvolve, the configurations reduce area by , , and , respectively. The configuration also reduces the later HorizonMath value by .
Proof. The bundled offsets encode the three stated configurations. Every stored float is interpreted as its exact binary rational value. Lemma 1.2 then computes the area of that precise configuration with rational arithmetic. Since each explicit configuration is admissible, its area is an upper bound for . The percentage statements are direct arithmetic using the published comparison values.
The next cell interprets every stored offset as an exact binary rational, computes all seven finite areas, and checks the stated reductions.
Show code
Code cell 7 · In [3]
# Published finite comparison values, used only for the arithmetic comparison.
ALPHAEVOLVE = {
2: Decimal("0.333333333333333"),
4: Decimal("0.250000000000000"),
8: Decimal("0.196858619034623"),
16: Decimal("0.162715712988197"),
32: Decimal("0.141647494730520"),
64: Decimal("0.121735177166418"),
128: Decimal("0.114810325818618"),
}
HORIZONMATH_128 = Decimal("0.1091479892")
EXPECTED = {
2: Decimal("0.333333333333333333"),
4: Decimal("0.250000000000000000"),
8: Decimal("0.196858619034622930"),
16: Decimal("0.162715712988192540"),
32: Decimal("0.138596129670274320"),
64: Decimal("0.120900580819172100"),
128: Decimal("0.107066636561634810"),
}
EXACT_FINITE_AREAS = {}
rows = []
for n in TESTED_N:
rational_offsets = [Fraction.from_float(float(v)) for v in OFFSETS[n]]
area_q = area_exact(rational_offsets)
area_d = decimal_of(area_q)
EXACT_FINITE_AREAS[n] = area_q
if abs(area_d - EXPECTED[n]) >= Decimal("2e-16"):
raise AssertionError(f"area mismatch at n={n}: {area_d}")
reduction = (ALPHAEVOLVE[n] - area_d) / ALPHAEVOLVE[n] * 100
if abs(reduction) < Decimal("1e-9"):
reduction = Decimal(0)
rows.append((n, area_d, ALPHAEVOLVE[n], reduction))
print(f"{'n':>4} {'verified area':>22} {'AlphaEvolve':>18} {'reduction':>12}")
for n, area_d, ae, reduction in rows:
print(f"{n:4d} {float(area_d):22.15f} {float(ae):18.15f} {float(reduction):11.6f}%")
assert round(float(rows[4][3]), 2) == 2.15
assert round(float(rows[5][3]), 2) == 0.69
assert round(float(rows[6][3]), 2) == 6.74
horizon_reduction = (HORIZONMATH_128 - decimal_of(EXACT_FINITE_AREAS[128])) / HORIZONMATH_128 * 100
assert round(float(horizon_reduction), 2) == 1.91
assert decimal_of(EXACT_FINITE_AREAS[128]) < Decimal("0.107067")
print(f"HorizonMath reduction at n=128: {float(horizon_reduction):.6f}%")
Saved output 1
n verified area AlphaEvolve reduction 2 0.333333333333333 0.333333333333333 0.000000% 4 0.250000000000000 0.250000000000000 0.000000% 8 0.196858619034623 0.196858619034623 0.000000% 16 0.162715712988193 0.162715712988197 0.000000% 32 0.138596129670274 0.141647494730520 2.154196% 64 0.120900580819172 0.121735177166418 0.685584% 128 0.107066636561635 0.114810325818618 6.744767% HorizonMath reduction at n=128: 1.906909%
Station Reference. Archives #34, #71, and #88; Evaluations #376, #379, #818, #855, and #870. The finite-scope audit and the comparison with the guided AlphaEvolve generalizer are recorded in the corresponding construction reports.
2.2 Scope of the construction benchmark
Benchmark scope. The benchmark tested exactly
The exact results at below are separate mathematical findings.
Station Reference. The Research Center task specification and evaluator define the seven tested sizes, offset-only submission contract, and area-ratio score.
Related Work. Córdoba (1977) and Keich (1999) establish the classical reciprocal-logarithmic framework. Georgiev et al. (2025) reported fixed-size constructions and a guided rule for arbitrary , and Wang et al. (2026) subsequently improved the value. The contribution here is a further improvement at the fixed sizes .
3. S2. Exact optima at and symmetry breaking at .
3.1 Exact global optima at
Translation does not affect area, so fix the gauge . Every endpoint-intersection height is affine in the remaining offsets. The equations saying that such a height equals , equals , or equals another event height form a finite rational hyperplane arrangement in configuration space.
Lemma 3.1 (completeness of the arrangement certificate). Within each full-dimensional cell of , the ordering of all sweep events is fixed and is one rational quadratic polynomial in the offsets. The quadratic need not be convex. The minimum over the closure of a cell is found by solving its full affine stationary system, including singular cases, and recursively minimizing the restricted quadratic on every facet. Therefore the least value over every full-dimensional cell is the global minimum.
Proof. The defining hyperplanes record exactly the values at which an event enters or leaves , or two event heights exchange order. Away from them, the complete sweep schedule is fixed. On each height slab, the slice union length is affine in the offsets and in ; the slab endpoints are affine in the offsets. Exact integration therefore gives a quadratic. A minimum of a differentiable quadratic on a polyhedral cell closure lies in the relative interior of some face on which its projected gradient vanishes. Solving that possibly singular affine system and repeating the argument on every boundary face terminates at edges and vertices, so the recursion is exhaustive. Finally, the closures of the full-dimensional cells cover the entire gauge-fixed space.
Theorem 3.2 (global values).
Proof. The following exact-rational program constructs the complete unbounded arrangements. For it checks all 72 cells. For it finds 9,350 cells and first verifies that reflection pairs them into 4,675 orbits. It then fully minimizes one representative from every orbit, with no stationary-value pruning and with exact treatment of singular stationary affine sets. For each reflection orbit, the representative polynomial is checked at an exact rational interior point against the independent sweep of Lemma 1.2, and reflection transfers the result to its partner. The least values are and , attained at and .
The next cell defines the complete exact-rational hyperplane-arrangement certificate used for and .
Show code
Code cell 12 · In [4]
"""Exact-rational full-space certificates for small Kakeya triangle unions.
The parameter gauge is x_n = 0, so the free coordinates are
``(x_1, ..., x_{n-1})``. For n=3 and n=4 this module enumerates the complete
unbounded rational arrangement cut out by endpoint-crossing heights h=0,
h=1, and h_i=h_j. On each full-dimensional cell the sweep schedule is fixed,
so the union area is an exact quadratic polynomial.
"""
from __future__ import annotations
import json
import math
import os
import time
from fractions import Fraction
from typing import Dict, List, Optional, Sequence, Tuple
Rat = Fraction
Point = Tuple[Rat, ...]
Hyper = Tuple[Rat, ...]
Poly = Dict[Tuple[int, ...], Rat]
def frac_json(x: Rat) -> str:
return f"{x.numerator}/{x.denominator}" if x.denominator != 1 else str(x.numerator)
def parse_frac(s: str) -> Rat:
if "/" in s:
a, b = s.split("/")
return Rat(int(a), int(b))
return Rat(int(s), 1)
def normalize_offsets(offsets: Sequence[Rat]) -> List[Rat]:
m = min(offsets)
return [x - m for x in offsets]
def offsets_from_params(n: int, params: Sequence[Rat]) -> List[Rat]:
if len(params) != n - 1:
raise ValueError("wrong parameter dimension")
return [Rat(v) for v in params] + [Rat(0)]
def endpoint_affines(n: int) -> List[Tuple[List[Rat], Rat, Rat, Tuple[int, str]]]:
d = n - 1
rows: List[Tuple[List[Rat], Rat, Rat, Tuple[int, str]]] = []
for j in range(1, n + 1):
coeff = [Rat(0) for _ in range(d)]
if j < n:
coeff[j - 1] = Rat(1)
rows.append((coeff[:], Rat(0), Rat(j, n), (j, "L")))
rows.append((coeff[:], Rat(1, n), Rat(j - 1, n), (j, "R")))
return rows
def event_affines_frac(n: int) -> List[Tuple[List[Rat], Rat, Tuple[int, int]]]:
lines = endpoint_affines(n)
out: List[Tuple[List[Rat], Rat, Tuple[int, int]]] = []
d = n - 1
for i, (ai, bi, si, _) in enumerate(lines):
for j, (aj, bj, sj, _) in enumerate(lines[i + 1 :], start=i + 1):
if si == sj:
continue
den = si - sj
out.append(([(aj[k] - ai[k]) / den for k in range(d)], (bj - bi) / den, (i, j)))
return out
def canonical_hyper(row: Sequence[Rat]) -> Optional[Hyper]:
if all(v == 0 for v in row[:-1]):
return None
lcm = 1
for v in row:
lcm = lcm * v.denominator // math.gcd(lcm, v.denominator)
vals = [int(v * lcm) for v in row]
g = 0
for v in vals:
g = math.gcd(g, abs(v))
vals = [v // g for v in vals]
first = next(v for v in vals if v != 0)
if first < 0:
vals = [-v for v in vals]
return tuple(Rat(v) for v in vals)
def canonical_hyper_with_sign(row: Sequence[Rat]) -> Optional[Tuple[Hyper, int]]:
h = canonical_hyper(row)
if h is None:
return None
for raw, can in zip(row, h):
if raw != 0:
return h, (1 if raw / can > 0 else -1)
return h, 1
def arrangement_hyperplanes(n: int) -> List[Hyper]:
events = event_affines_frac(n)
d = n - 1
raw: List[List[Rat]] = []
for a, c, _ in events:
raw.append(a + [c])
raw.append(a + [c - 1])
for i, (ai, ci, _) in enumerate(events):
for aj, cj, _ in events[i + 1 :]:
raw.append([ai[k] - aj[k] for k in range(d)] + [ci - cj])
seen = set()
out: List[Hyper] = []
for row in raw:
h = canonical_hyper(row)
if h is not None and h not in seen:
seen.add(h)
out.append(h)
out.sort()
return out
def hyper_value(h: Hyper, p: Point) -> Rat:
return sum(h[i] * p[i] for i in range(len(p))) + h[-1]
def sign_vector(hypers: Sequence[Hyper], p: Point) -> Tuple[int, ...]:
out = []
for h in hypers:
v = hyper_value(h, p)
if v == 0:
raise ValueError("sample on hyperplane")
out.append(1 if v > 0 else -1)
return tuple(out)
def enumerate_cells_1d(points_in: Sequence[Hyper]) -> Dict[Tuple[int, ...], Point]:
points = sorted(-h[1] / h[0] for h in points_in)
samples: List[Rat] = []
if not points:
samples = [Rat(0)]
else:
samples.append(points[0] - 1)
for a, b in zip(points, points[1:]):
samples.append((a + b) / 2)
samples.append(points[-1] + 1)
cells = {}
for x in samples:
cells[tuple(1 if x - p > 0 else -1 for p in points)] = (x,)
return cells
def enumerate_cells_2d(lines_in: Sequence[Hyper]) -> Dict[Tuple[int, ...], Point]:
lines: List[Tuple[Rat, Rat, Rat]] = [(h[0], h[1], h[2]) for h in lines_in]
cells: Dict[Tuple[int, ...], Point] = {}
if not lines:
return {tuple(): (Rat(0), Rat(0))}
for idx, line in enumerate(lines):
a, b, c = line
if b != 0:
p0 = (Rat(0), -c / b)
else:
p0 = (-c / a, Rat(0))
direction = (b, -a)
ts = set()
for j, other in enumerate(lines):
if j == idx:
continue
aa, bb, cc = other
alpha = aa * direction[0] + bb * direction[1]
beta = aa * p0[0] + bb * p0[1] + cc
if alpha != 0:
ts.add(-beta / alpha)
ordered = sorted(ts)
samples: List[Rat] = []
if not ordered:
samples = [Rat(0)]
else:
samples.append(ordered[0] - 1)
for u, v in zip(ordered, ordered[1:]):
samples.append((u + v) / 2)
samples.append(ordered[-1] + 1)
for t in samples:
base = (p0[0] + t * direction[0], p0[1] + t * direction[1])
eps_limit: Optional[Rat] = None
for j, other in enumerate(lines):
if j == idx:
continue
val = other[0] * base[0] + other[1] * base[1] + other[2]
delta = other[0] * a + other[1] * b
if delta != 0:
bound = abs(val / delta) / 2
if bound > 0 and (eps_limit is None or bound < eps_limit):
eps_limit = bound
eps = eps_limit if eps_limit is not None else Rat(1, 2)
if eps == 0:
continue
for side in (-1, 1):
p = (base[0] + side * eps * a, base[1] + side * eps * b)
try:
sv = tuple(1 if (ln[0] * p[0] + ln[1] * p[1] + ln[2]) > 0 else -1 for ln in lines)
except ValueError:
continue
cells.setdefault(sv, p)
return cells
def plane_base_and_basis(plane: Hyper) -> Tuple[Point, Point, Point, Point]:
a, b, c, d = plane
if c != 0:
p0 = (Rat(0), Rat(0), -d / c)
u = (Rat(1), Rat(0), -a / c)
v = (Rat(0), Rat(1), -b / c)
elif b != 0:
p0 = (Rat(0), -d / b, Rat(0))
u = (Rat(1), -a / b, Rat(0))
v = (Rat(0), -c / b, Rat(1))
else:
p0 = (-d / a, Rat(0), Rat(0))
u = (-b / a, Rat(1), Rat(0))
v = (-c / a, Rat(0), Rat(1))
normal = (a, b, c)
return p0, u, v, normal
def map_plane_point(p0: Point, u: Point, v: Point, q: Point) -> Point:
return tuple(p0[i] + q[0] * u[i] + q[1] * v[i] for i in range(3))
def induced_lines_on_plane(hypers: Sequence[Hyper], plane_idx: int) -> Tuple[List[Hyper], List[int]]:
p0, u, v, _ = plane_base_and_basis(hypers[plane_idx])
lines: List[Hyper] = []
owners: List[int] = []
seen = set()
for idx, h in enumerate(hypers):
if idx == plane_idx:
continue
alpha = sum(h[i] * u[i] for i in range(3))
beta = sum(h[i] * v[i] for i in range(3))
gamma = hyper_value(h, p0)
canon = canonical_hyper((alpha, beta, gamma))
if canon is None:
continue
if canon not in seen:
seen.add(canon)
lines.append(canon)
owners.append(idx)
order = sorted(range(len(lines)), key=lambda i: lines[i])
return [lines[i] for i in order], [owners[i] for i in order]
def enumerate_cells_3d(hypers: Sequence[Hyper], deadline: Optional[float] = None) -> Tuple[Dict[Tuple[int, ...], Point], bool]:
cells: Dict[Tuple[int, ...], Point] = {}
for plane_idx, plane in enumerate(hypers):
if deadline is not None and time.time() > deadline:
return cells, False
induced, _owners = induced_lines_on_plane(hypers, plane_idx)
face_cells = enumerate_cells_2d(induced)
p0, u, v, normal = plane_base_and_basis(plane)
for _face_sv, q in face_cells.items():
base = map_plane_point(p0, u, v, q)
eps_limit: Optional[Rat] = None
for idx, h in enumerate(hypers):
if idx == plane_idx:
continue
val = hyper_value(h, base)
delta = sum(h[i] * normal[i] for i in range(3))
if delta != 0:
bound = abs(val / delta) / 2
if bound > 0 and (eps_limit is None or bound < eps_limit):
eps_limit = bound
eps = eps_limit if eps_limit is not None else Rat(1, 2)
if eps == 0:
continue
for side in (-1, 1):
p = tuple(base[i] + side * eps * normal[i] for i in range(3))
try:
sv = sign_vector(hypers, p)
except ValueError:
continue
cells.setdefault(sv, p)
return cells, True
def affine_endpoint_value(n: int, endpoint_index: int, params: Point, y: Rat) -> Rat:
coeff, const, slope, _ = endpoint_affines(n)[endpoint_index]
return sum(coeff[i] * params[i] for i in range(len(params))) + const + slope * y
def interval_at_y(n: int, offsets: Sequence[Rat], j: int, y: Rat) -> Tuple[Rat, Rat]:
x = offsets[j - 1]
return x + Rat(j, n) * y, x + Rat(1, n) + Rat(j - 1, n) * y
def union_length_at_y_frac(n: int, offsets: Sequence[Rat], y: Rat) -> Rat:
intervals = [interval_at_y(n, offsets, j, y) for j in range(1, n + 1)]
intervals.sort()
total = Rat(0)
cur_l: Optional[Rat] = None
cur_r: Optional[Rat] = None
for lft, rgt in intervals:
if cur_l is None:
cur_l, cur_r = lft, rgt
elif lft <= cur_r:
if rgt > cur_r:
cur_r = rgt
else:
total += cur_r - cur_l
cur_l, cur_r = lft, rgt
if cur_l is not None:
total += cur_r - cur_l
return total
def exact_area_offsets(n: int, offsets: Sequence[Rat]) -> Rat:
lines: List[Tuple[Rat, Rat]] = []
for j, x in enumerate(offsets, start=1):
lines.append((Rat(j, n), x))
lines.append((Rat(j - 1, n), x + Rat(1, n)))
breaks = {Rat(0), Rat(1)}
for i, (s1, b1) in enumerate(lines):
for s2, b2 in lines[i + 1 :]:
if s1 == s2:
continue
y = (b2 - b1) / (s1 - s2)
if Rat(0) < y < Rat(1):
breaks.add(y)
bp = sorted(breaks)
area = Rat(0)
for a, b in zip(bp, bp[1:]):
if a != b:
area += (b - a) * (union_length_at_y_frac(n, offsets, a) + union_length_at_y_frac(n, offsets, b)) / 2
return area
def poly_affine(coeff: Sequence[Rat], const: Rat) -> Poly:
d = len(coeff)
out: Poly = {}
if const:
out[(0,) * d] = const
for i, a in enumerate(coeff):
if a:
exp = [0] * d
exp[i] = 1
out[tuple(exp)] = out.get(tuple(exp), Rat(0)) + a
return {k: v for k, v in out.items() if v}
def poly_add(a: Poly, b: Poly, scale: Rat = Rat(1)) -> Poly:
out = dict(a)
for k, v in b.items():
out[k] = out.get(k, Rat(0)) + scale * v
if out[k] == 0:
del out[k]
return out
def poly_mul(a: Poly, b: Poly) -> Poly:
out: Poly = {}
for ea, va in a.items():
for eb, vb in b.items():
exp = tuple(ea[i] + eb[i] for i in range(len(ea)))
if sum(exp) > 2:
raise ValueError("degree exceeded quadratic")
out[exp] = out.get(exp, Rat(0)) + va * vb
return {k: v for k, v in out.items() if v}
def poly_eval(poly: Poly, p: Point) -> Rat:
total = Rat(0)
for exp, coef in poly.items():
term = coef
for i, e in enumerate(exp):
if e == 1:
term *= p[i]
elif e == 2:
term *= p[i] * p[i]
total += term
return total
def union_expression_on_slab(n: int, sample: Point, y_mid: Rat) -> Tuple[List[Rat], Rat, Rat]:
lines = endpoint_affines(n)
intervals = []
for j in range(1, n + 1):
li = 2 * (j - 1)
ri = li + 1
lv = affine_endpoint_value(n, li, sample, y_mid)
rv = affine_endpoint_value(n, ri, sample, y_mid)
intervals.append((lv, rv, li, ri))
intervals.sort(key=lambda row: (row[0], row[1], row[2], row[3]))
clusters: List[Tuple[int, int]] = []
cur_l: Optional[int] = None
cur_r: Optional[int] = None
cur_rv: Optional[Rat] = None
for lv, rv, li, ri in intervals:
if cur_l is None:
cur_l, cur_r, cur_rv = li, ri, rv
elif lv <= cur_rv:
if rv > cur_rv:
cur_r, cur_rv = ri, rv
else:
clusters.append((cur_l, cur_r)) # type: ignore[arg-type]
cur_l, cur_r, cur_rv = li, ri, rv
if cur_l is not None:
clusters.append((cur_l, cur_r)) # type: ignore[arg-type]
d = len(sample)
a = [Rat(0) for _ in range(d)]
b = Rat(0)
cy = Rat(0)
for li, ri in clusters:
ar, br, sr, _ = lines[ri]
al, bl, sl, _ = lines[li]
for k in range(d):
a[k] += ar[k] - al[k]
b += br - bl
cy += sr - sl
return a, b, cy
def schedule_quadratic_exact(n: int, sample: Point) -> Poly:
d = len(sample)
events = event_affines_frac(n)
vals = [(sum(a[i] * sample[i] for i in range(d)) + c, a, c) for a, c, _ in events]
active = [(v, a, c) for v, a, c in vals if Rat(0) < v < Rat(1)]
active.sort(key=lambda row: row[0])
breaks = [(Rat(0), [Rat(0) for _ in range(d)], Rat(0))] + active + [(Rat(1), [Rat(0) for _ in range(d)], Rat(1))]
poly: Poly = {}
for left, right in zip(breaks, breaks[1:]):
v0, a0, c0 = left
v1, a1, c1 = right
if v0 == v1:
continue
y_mid = (v0 + v1) / 2
u_a, u_b, u_cy = union_expression_on_slab(n, sample, y_mid)
h0 = poly_affine(a0, c0)
h1 = poly_affine(a1, c1)
dh = poly_add(h1, h0, scale=-1)
term = poly_mul(poly_affine(u_a, u_b), dh)
h1_sq = poly_mul(h1, h1)
h0_sq = poly_mul(h0, h0)
yterm = poly_add(h1_sq, h0_sq, scale=-1)
yterm = {k: u_cy * v / 2 for k, v in yterm.items()}
poly = poly_add(poly, term)
poly = poly_add(poly, yterm)
return poly
def poly_gradient_linear(poly: Poly, d: int) -> Tuple[List[List[Rat]], List[Rat]]:
mat = [[Rat(0) for _ in range(d)] for _ in range(d)]
vec = [Rat(0) for _ in range(d)]
for exp, coef in poly.items():
for i in range(d):
if exp[i] == 0:
continue
new_exp = list(exp)
new_exp[i] -= 1
deriv_coef = coef * exp[i]
if sum(new_exp) == 0:
vec[i] += deriv_coef
else:
j = new_exp.index(1)
mat[i][j] += deriv_coef
return mat, vec
def solve_linear(mat: List[List[Rat]], rhs: List[Rat]) -> Optional[Point]:
n = len(rhs)
aug = [row[:] + [rhs[i]] for i, row in enumerate(mat)]
pivots = []
r = 0
for c in range(n):
piv = next((i for i in range(r, n) if aug[i][c] != 0), None)
if piv is None:
continue
aug[r], aug[piv] = aug[piv], aug[r]
pv = aug[r][c]
aug[r] = [v / pv for v in aug[r]]
for i in range(n):
if i != r and aug[i][c] != 0:
f = aug[i][c]
aug[i] = [aug[i][j] - f * aug[r][j] for j in range(n + 1)]
pivots.append(c)
r += 1
if r < n:
return None
sol = [Rat(0) for _ in range(n)]
for row, c in enumerate(pivots):
sol[c] = aug[row][-1]
return tuple(sol)
def feasible(hypers: Sequence[Hyper], signs: Tuple[int, ...], p: Point) -> bool:
return all(s * hyper_value(h, p) >= 0 for h, s in zip(hypers, signs))
def q1_coeffs(poly: Poly) -> Tuple[Rat, Rat, Rat]:
return poly.get((2,), Rat(0)), poly.get((1,), Rat(0)), poly.get((0,), Rat(0))
def minimize_cell_1d(poly: Poly, points: Sequence[Hyper], signs: Tuple[int, ...]) -> Tuple[Rat, Point]:
lo: Optional[Rat] = None
hi: Optional[Rat] = None
for h, s in zip(points, signs):
a, c = h
p = -c / a
if s * a > 0:
lo = p if lo is None or p > lo else lo
else:
hi = p if hi is None or p < hi else hi
candidates: List[Point] = []
if lo is not None:
candidates.append((lo,))
if hi is not None:
candidates.append((hi,))
a, b, _ = q1_coeffs(poly)
if a > 0:
candidates.append((-b / (2 * a),))
if (a == 0 and b == 0) or (lo is None and hi is None):
candidates.append((Rat(0),))
best: Optional[Tuple[Rat, Point]] = None
for p in candidates:
if feasible(points, signs, p):
val = poly_eval(poly, p)
if best is None or val < best[0]:
best = (val, p)
if best is None:
raise RuntimeError("no 1D minimizer candidate")
return best
def line_interval(lines: Sequence[Hyper], signs: Tuple[int, ...], boundary_idx: int) -> Tuple[Point, Point, Optional[Rat], Optional[Rat]]:
a, b, c = lines[boundary_idx]
if b != 0:
p0 = (Rat(0), -c / b)
else:
p0 = (-c / a, Rat(0))
direction = (b, -a)
lo: Optional[Rat] = None
hi: Optional[Rat] = None
for idx, (line, s) in enumerate(zip(lines, signs)):
if idx == boundary_idx:
continue
aa, bb, cc = line
alpha = s * (aa * direction[0] + bb * direction[1])
beta = s * (aa * p0[0] + bb * p0[1] + cc)
if alpha == 0:
if beta < 0:
return p0, direction, Rat(1), Rat(0)
continue
bound = -beta / alpha
if alpha > 0:
lo = bound if lo is None or bound > lo else lo
else:
hi = bound if hi is None or bound < hi else hi
return p0, direction, lo, hi
def restrict_poly_to_line(poly: Poly, p0: Point, direction: Point) -> Tuple[Rat, Rat, Rat]:
vals = []
for t in (Rat(0), Rat(1), Rat(2)):
p = tuple(p0[i] + t * direction[i] for i in range(len(p0)))
vals.append(poly_eval(poly, p))
c = vals[0]
a = (vals[2] - 2 * vals[1] + vals[0]) / 2
b = vals[1] - vals[0] - a
return a, b, c
def minimize_cell_2d(
poly: Poly,
lines: Sequence[Hyper],
signs: Tuple[int, ...],
boundary_indices: Optional[Sequence[int]] = None,
) -> Tuple[Rat, Point]:
candidates: List[Point] = []
mat, vec = poly_gradient_linear(poly, 2)
stat = solve_linear(mat, [-vec[0], -vec[1]])
if stat is not None:
candidates.append(stat)
if not lines and stat is None:
candidates.append((Rat(0), Rat(0)))
indices = list(boundary_indices) if boundary_indices is not None else list(range(len(lines)))
for i in indices:
p0, direction, lo, hi = line_interval(lines, signs, i)
if lo is not None and hi is not None and lo > hi:
continue
a, b, _ = restrict_poly_to_line(poly, p0, direction)
ts: List[Rat] = []
if lo is not None:
ts.append(lo)
if hi is not None:
ts.append(hi)
if a > 0:
t0 = -b / (2 * a)
if (lo is None or t0 >= lo) and (hi is None or t0 <= hi):
ts.append(t0)
elif a == 0 and b == 0:
if lo is not None:
ts.append(lo)
elif hi is not None:
ts.append(hi)
else:
ts.append(Rat(0))
for t in ts:
candidates.append((p0[0] + t * direction[0], p0[1] + t * direction[1]))
best: Optional[Tuple[Rat, Point]] = None
for p in candidates:
if feasible(lines, signs, p):
val = poly_eval(poly, p)
if best is None or val < best[0]:
best = (val, p)
if best is None:
raise RuntimeError("no 2D minimizer candidate")
return best
def restrict_poly_affine(poly: Poly, p0: Point, dirs: Sequence[Point]) -> Poly:
k = len(dirs)
out: Poly = {}
d = len(p0)
for exp, coef in poly.items():
terms: Poly = {(0,) * k: coef}
for i in range(d):
power = exp[i]
if power == 0:
continue
aff: Poly = {(0,) * k: p0[i]} if p0[i] else {}
for j, direction in enumerate(dirs):
if direction[i]:
e = [0] * k
e[j] = 1
aff[tuple(e)] = aff.get(tuple(e), Rat(0)) + direction[i]
if power == 1:
terms = poly_mul(terms, aff)
elif power == 2:
terms = poly_mul(poly_mul(terms, aff), aff)
out = poly_add(out, terms)
return out
def hyperplane_param(h: Hyper) -> Tuple[Point, List[Point]]:
d = len(h) - 1
coeff = h[:-1]
const = h[-1]
pivot = next(i for i, a in enumerate(coeff) if a != 0)
p0 = [Rat(0) for _ in range(d)]
p0[pivot] = -const / coeff[pivot]
dirs: List[Point] = []
for free in range(d):
if free == pivot:
continue
v = [Rat(0) for _ in range(d)]
v[free] = Rat(1)
v[pivot] = -coeff[free] / coeff[pivot]
dirs.append(tuple(v))
return tuple(p0), dirs
def minimize_cell(
poly: Poly,
hypers: Sequence[Hyper],
signs: Tuple[int, ...],
dim: Optional[int] = None,
boundary_indices: Optional[Sequence[int]] = None,
) -> Tuple[Rat, Point]:
d = dim if dim is not None else (len(next(iter(poly.keys()))) if poly else len(hypers[0]) - 1)
if d == 1:
return minimize_cell_1d(poly, hypers, signs)
if d == 2:
return minimize_cell_2d(poly, hypers, signs, boundary_indices=boundary_indices)
candidates: List[Point] = []
mat, vec = poly_gradient_linear(poly, d)
stat = solve_linear(mat, [-v for v in vec])
if stat is not None:
candidates.append(stat)
indices = list(boundary_indices) if boundary_indices is not None else list(range(len(hypers)))
for idx in indices:
h = hypers[idx]
p0, dirs = hyperplane_param(h)
sub_poly = restrict_poly_affine(poly, p0, dirs)
sub_hypers: List[Hyper] = []
sub_signs: List[int] = []
sub_boundary: List[int] = []
infeasible_face = False
for j, hh in enumerate(hypers):
if j == idx:
continue
coeff = [sum(hh[i] * direction[i] for i in range(d)) for direction in dirs]
const = hyper_value(hh, p0)
canon_row = canonical_hyper_with_sign(coeff + [const])
if canon_row is None:
if signs[j] * const < 0:
infeasible_face = True
break
else:
canon, orient = canon_row
sub_idx = len(sub_hypers)
sub_hypers.append(canon)
sub_signs.append(signs[j] * orient)
if boundary_indices is not None and j in boundary_indices:
sub_boundary.append(sub_idx)
if infeasible_face:
continue
if len(dirs) == 0:
candidates.append(p0)
else:
try:
sub_val, sub_p = minimize_cell(
sub_poly,
sub_hypers,
tuple(sub_signs),
dim=len(dirs),
boundary_indices=sub_boundary if boundary_indices is not None else None,
)
except RuntimeError:
continue
del sub_val
candidates.append(tuple(p0[i] + sum(sub_p[j] * dirs[j][i] for j in range(len(dirs))) for i in range(d)))
best: Optional[Tuple[Rat, Point]] = None
for p in candidates:
if feasible(hypers, signs, p):
val = poly_eval(poly, p)
if best is None or val < best[0]:
best = (val, p)
if best is None:
raise RuntimeError("no exact minimizer candidate")
return best
def map_affine_point(p0: Point, dirs: Sequence[Point], q: Point) -> Point:
return tuple(p0[i] + sum(q[j] * dirs[j][i] for j in range(len(dirs))) for i in range(len(p0)))
def enumerate_cells_4d(hypers: Sequence[Hyper], deadline: Optional[float] = None) -> Tuple[Dict[Tuple[int, ...], Point], bool]:
cells: Dict[Tuple[int, ...], Point] = {}
for plane_idx, plane in enumerate(hypers):
if deadline is not None and time.time() > deadline:
return cells, False
p0, dirs = hyperplane_param(plane)
induced: List[Hyper] = []
seen = set()
for idx, h in enumerate(hypers):
if idx == plane_idx:
continue
coeff = [sum(h[i] * direction[i] for i in range(4)) for direction in dirs]
const = hyper_value(h, p0)
canon = canonical_hyper(coeff + [const])
if canon is None:
continue
if canon not in seen:
seen.add(canon)
induced.append(canon)
induced.sort()
face_cells, ok = enumerate_cells_3d(induced, deadline=deadline)
if not ok:
return cells, False
normal = plane[:-1]
for q in face_cells.values():
base = map_affine_point(p0, dirs, q)
eps_limit: Optional[Rat] = None
for idx, h in enumerate(hypers):
if idx == plane_idx:
continue
val = hyper_value(h, base)
delta = sum(h[i] * normal[i] for i in range(4))
if delta != 0:
bound = abs(val / delta) / 2
if bound > 0 and (eps_limit is None or bound < eps_limit):
eps_limit = bound
eps = eps_limit if eps_limit is not None else Rat(1, 2)
if eps == 0:
continue
for side in (-1, 1):
p = tuple(base[i] + side * eps * normal[i] for i in range(4))
try:
sv = sign_vector(hypers, p)
except ValueError:
continue
cells.setdefault(sv, p)
return cells, True
def sigma_params(n: int, p: Point) -> Point:
if len(p) != n - 1:
raise ValueError("wrong sigma parameter dimension")
x1 = p[0]
full = list(p) + [Rat(0)]
return tuple(x1 - full[n - j] for j in range(1, n))
def sigma_orbit_report(n: int, hypers: Sequence[Hyper], cells: Dict[Tuple[int, ...], Point]) -> Dict[str, object]:
mapped_missing = 0
non_involutive = 0
fixed = 0
reps = set()
sample_checks = []
for signs, sample in cells.items():
sp = sigma_params(n, sample)
mapped = sign_vector(hypers, sp)
if mapped not in cells:
mapped_missing += 1
if sign_vector(hypers, sigma_params(n, sp)) != signs:
non_involutive += 1
if mapped == signs:
fixed += 1
reps.add(min(signs, mapped))
if len(sample_checks) < 3:
a0 = exact_area_offsets(n, offsets_from_params(n, sample))
a1 = exact_area_offsets(n, offsets_from_params(n, sp))
sample_checks.append({"area": frac_json(a0), "sigma_area": frac_json(a1), "matches": a0 == a1})
return {
"enabled": True,
"orbit_representative_count": len(reps),
"fixed_cell_count": fixed,
"mapped_missing_count": mapped_missing,
"non_involutive_count": non_involutive,
"sample_area_checks": sample_checks,
"bijection_verified": mapped_missing == 0 and non_involutive == 0,
}
def unconstrained_stationary_value(poly: Poly, d: int) -> Optional[Tuple[Rat, Point]]:
mat, vec = poly_gradient_linear(poly, d)
stat = solve_linear(mat, [-v for v in vec])
if stat is None:
return None
return poly_eval(poly, stat), stat
def certify_fullspace_exact(
n: int,
deadline: Optional[float] = None,
incumbent: Optional[Rat] = None,
use_sigma_orbits: bool = False,
) -> Dict[str, object]:
start = time.time()
hypers = arrangement_hyperplanes(n)
d = n - 1
if d == 1:
cells = enumerate_cells_1d(hypers)
complete = True
elif d == 2:
cells = enumerate_cells_2d(hypers)
complete = True
elif d == 3:
cells, complete = enumerate_cells_3d(hypers, deadline=deadline)
elif d == 4:
cells, complete = enumerate_cells_4d(hypers, deadline=deadline)
else:
raise ValueError("fullspace exact certifier currently supports dimensions <= 4")
best_val: Optional[Rat] = None
best_params: Optional[Point] = None
sample_mismatches = []
processed = 0
pruned = 0
degenerate_unconstrained = 0
reps_skipped = 0
cell_sign_set = set(cells)
sigma_report: Dict[str, object] = {"enabled": False}
if use_sigma_orbits:
sigma_report = sigma_orbit_report(n, hypers, cells)
if not sigma_report["bijection_verified"]:
use_sigma_orbits = False
visited_or_represented = 0
for signs, sample in cells.items():
if deadline is not None and time.time() > deadline:
complete = False
break
if use_sigma_orbits:
mapped = sign_vector(hypers, sigma_params(n, sample))
if signs != min(signs, mapped):
reps_skipped += 1
continue
visited_or_represented += 1
poly = schedule_quadratic_exact(n, sample)
sample_area = exact_area_offsets(n, offsets_from_params(n, sample))
formula_area = poly_eval(poly, sample)
if sample_area != formula_area:
sample_mismatches.append(
{"sample": [frac_json(x) for x in sample], "sweep": frac_json(sample_area), "formula": frac_json(formula_area)}
)
continue
stat = unconstrained_stationary_value(poly, d)
if incumbent is not None and stat is not None and stat[0] > incumbent:
pruned += 1
processed += 1
continue
if stat is None:
degenerate_unconstrained += 1
boundary_indices = None
if d == 3 and complete:
local = []
for idx in range(len(hypers)):
flipped = list(signs)
flipped[idx] *= -1
if tuple(flipped) in cell_sign_set:
local.append(idx)
boundary_indices = local
if d == 4 and complete:
local = []
for idx in range(len(hypers)):
flipped = list(signs)
flipped[idx] *= -1
if tuple(flipped) in cell_sign_set:
local.append(idx)
boundary_indices = local
val, params = minimize_cell(poly, hypers, signs, dim=d, boundary_indices=boundary_indices)
if best_val is None or val < best_val:
best_val, best_params = val, params
processed += 1
if sample_mismatches:
raise RuntimeError(f"quadratic/sweep mismatch in {len(sample_mismatches)} cells")
if best_val is None or best_params is None:
if incumbent is not None:
best_val = incumbent
best_params = tuple(Rat(0) for _ in range(d))
else:
raise RuntimeError("no cells processed")
best_offsets = offsets_from_params(n, best_params)
best_sweep = exact_area_offsets(n, best_offsets)
plateau_checks: List[Dict[str, object]] = []
if n == 4:
for c in (Rat(1, 20), Rat(7, 80), Rat(1, 8)):
off = [Rat(1, 4), Rat(1, 4) - c, c, Rat(0)]
plateau_checks.append({"c": frac_json(c), "area": frac_json(exact_area_offsets(4, off))})
return {
"n": n,
"dimension": d,
"area": frac_json(best_val),
"area_float": float(best_val),
"params_xn_zero": [frac_json(x) for x in best_params],
"offsets_normalized": [frac_json(x) for x in normalize_offsets(best_offsets)],
"exact_sweep_at_min": frac_json(best_sweep),
"min_formula_matches_sweep": best_sweep == best_val,
"complete": complete,
"completion_scope": "full_unbounded_parameter_arrangement" if complete else "partial_full_unbounded_parameter_arrangement",
"cell_count": len(cells),
"processed_cell_count": processed,
"orbit_representatives_visited": visited_or_represented,
"sigma_orbit_reduction": sigma_report,
"cell_orbit_representatives_skipped": reps_skipped,
"incumbent_pruned_cell_count": pruned,
"fully_minimized_cell_count": processed - pruned,
"degenerate_unconstrained_cell_count": degenerate_unconstrained,
"incumbent": frac_json(incumbent) if incumbent is not None else None,
"hyperplane_count": len(hypers),
"quadratic_sweep_match_count": processed,
"plateau_checks": plateau_checks,
"wall_clock": time.time() - start,
}
The next cell adds exact minimization on singular affine stationary sets, ensuring that the certificate also handles quadratics with non-isolated stationary points.
Show code
Code cell 14 · In [5]
def affine_solutions(mat, rhs):
"""Return one exact solution and a nullspace basis, or None if inconsistent."""
rows = len(mat)
cols = len(mat[0]) if rows else len(rhs)
aug = [list(mat[i]) + [rhs[i]] for i in range(rows)]
pivots = []
r = 0
for c in range(cols):
pivot = next((i for i in range(r, rows) if aug[i][c] != 0), None)
if pivot is None:
continue
aug[r], aug[pivot] = aug[pivot], aug[r]
q = aug[r][c]
aug[r] = [v / q for v in aug[r]]
for i in range(rows):
if i != r and aug[i][c] != 0:
q = aug[i][c]
aug[i] = [aug[i][j] - q * aug[r][j] for j in range(cols + 1)]
pivots.append(c)
r += 1
for i in range(r, rows):
if all(aug[i][c] == 0 for c in range(cols)) and aug[i][-1] != 0:
return None
free = [c for c in range(cols) if c not in pivots]
p0 = [Rat(0) for _ in range(cols)]
for i, c in enumerate(pivots):
p0[c] = aug[i][-1]
directions = []
for free_col in free:
v = [Rat(0) for _ in range(cols)]
v[free_col] = Rat(1)
for i, c in enumerate(pivots):
v[c] = -aug[i][free_col]
directions.append(tuple(v))
return tuple(p0), directions
def find_feasible_affine(inequalities, dim):
"""Find an exact rational point satisfying a.t+c >= 0, or return None."""
if dim == 0:
return tuple() if all(c >= 0 for _a, c in inequalities) else None
origin = tuple(Rat(0) for _ in range(dim))
if all(sum(a[i] * origin[i] for i in range(dim)) + c >= 0 for a, c in inequalities):
return origin
for a, c in inequalities:
if all(v == 0 for v in a):
continue
p0, directions = hyperplane_param(tuple(a) + (c,))
restricted = []
infeasible = False
for aa, cc in inequalities:
coeff = tuple(sum(aa[i] * v[i] for i in range(dim)) for v in directions)
const = sum(aa[i] * p0[i] for i in range(dim)) + cc
if not directions and const < 0:
infeasible = True
break
restricted.append((coeff, const))
if infeasible:
continue
q = find_feasible_affine(restricted, dim - 1)
if q is not None:
return tuple(
p0[i] + sum(q[j] * directions[j][i] for j in range(dim - 1))
for i in range(dim)
)
return None
def stationary_candidate(poly, hypers, signs, dim):
"""Find a feasible point in the full affine stationary set, including singular cases."""
mat, vec = poly_gradient_linear(poly, dim)
solutions = affine_solutions(mat, [-v for v in vec])
if solutions is None:
return None
p0, directions = solutions
inequalities = []
for h, sign in zip(hypers, signs):
coeff = tuple(
sign * sum(h[i] * v[i] for i in range(dim)) for v in directions
)
const = sign * hyper_value(h, p0)
inequalities.append((coeff, const))
q = find_feasible_affine(inequalities, len(directions))
if q is None:
return None
return tuple(
p0[i] + sum(q[j] * directions[j][i] for j in range(len(directions)))
for i in range(dim)
)
def minimize_cell_robust(poly, hypers, signs, dim=None, boundary_indices=None):
"""Minimize over a polyhedral cell by exact recursion over every relevant face."""
d = dim if dim is not None else len(hypers[0]) - 1
candidates = []
stationary = stationary_candidate(poly, hypers, signs, d)
if stationary is not None:
candidates.append(stationary)
indices = list(boundary_indices) if boundary_indices is not None else list(range(len(hypers)))
for idx in indices:
h = hypers[idx]
p0, directions = hyperplane_param(h)
sub_poly = restrict_poly_affine(poly, p0, directions)
sub_hypers = []
sub_signs = []
sub_boundary = []
infeasible_face = False
for j, hh in enumerate(hypers):
if j == idx:
continue
coeff = [sum(hh[i] * direction[i] for i in range(d)) for direction in directions]
const = hyper_value(hh, p0)
canonical = canonical_hyper_with_sign(coeff + [const])
if canonical is None:
if signs[j] * const < 0:
infeasible_face = True
break
else:
row, orientation = canonical
sub_index = len(sub_hypers)
sub_hypers.append(row)
sub_signs.append(signs[j] * orientation)
if boundary_indices is not None and j in boundary_indices:
sub_boundary.append(sub_index)
if infeasible_face:
continue
if not directions:
candidates.append(p0)
else:
try:
_value, q = minimize_cell_robust(
sub_poly,
sub_hypers,
tuple(sub_signs),
dim=len(directions),
boundary_indices=sub_boundary if boundary_indices is not None else None,
)
except RuntimeError:
continue
candidates.append(tuple(
p0[i] + sum(q[j] * directions[j][i] for j in range(len(directions)))
for i in range(d)
))
best = None
for p in candidates:
if feasible(hypers, signs, p):
value = poly_eval(poly, p)
if best is None or value < best[0]:
best = value, p
if best is None:
raise RuntimeError("no exact minimizer candidate")
return best
# Replace the earlier nonsingular-only routine. Incumbent pruning remains disabled below.
minimize_cell = minimize_cell_robust
The next cell runs the complete and certificates and checks the arrangement sizes, full cell coverage, and exact minima.
Show code
Code cell 16 · In [6]
n3_certificate = certify_fullspace_exact(3)
n4_certificate = certify_fullspace_exact(4, use_sigma_orbits=True)
assert n3_certificate["complete"]
assert n3_certificate["hyperplane_count"] == 15
assert n3_certificate["cell_count"] == 72
assert n3_certificate["processed_cell_count"] == 72
assert n3_certificate["quadratic_sweep_match_count"] == 72
assert n3_certificate["incumbent_pruned_cell_count"] == 0
assert n3_certificate["area"] == "5/18"
assert n3_certificate["min_formula_matches_sweep"]
assert n4_certificate["complete"]
assert n4_certificate["hyperplane_count"] == 59
assert n4_certificate["cell_count"] == 9350
assert n4_certificate["sigma_orbit_reduction"]["bijection_verified"]
assert n4_certificate["orbit_representatives_visited"] == 4675
assert n4_certificate["fully_minimized_cell_count"] == 4675
assert n4_certificate["incumbent_pruned_cell_count"] == 0
assert n4_certificate["quadratic_sweep_match_count"] == 4675
assert n4_certificate["area"] == "1/4"
assert n4_certificate["min_formula_matches_sweep"]
print(f"{'n':>3} {'hyperplanes':>12} {'cells':>8} {'fully minimized':>17} {'minimum':>10}")
for row in (n3_certificate, n4_certificate):
print(
f"{row['n']:3d} {row['hyperplane_count']:12d} {row['cell_count']:8d} "
f"{row['fully_minimized_cell_count']:17d} {row['area']:>10}"
)
Saved output 1
n hyperplanes cells fully minimized minimum 3 15 72 72 5/18 4 59 9350 4675 1/4
Station Reference. Archives #4, #36, #94, and #96; Evaluations #1212, #1220, #1227, and #1238.
3.2 The continuous optimum
Theorem 3.3 (optimal plateau). For
one has . Thus every configuration with is globally optimal by Theorem 3.2.
Proof. Along this line the only event-order changes inside the claimed plateau occur at and . Restricting the exact cell quadratic to the line on each of the three open subintervals gives the identical constant polynomial . Continuity of area includes the two internal walls and both endpoints.
The next cell reconstructs these restricted polynomials symbolically rather than sampling isolated values.
Show code
Code cell 19 · In [7]
plateau_base = (Rat(1, 4), Rat(1, 4), Rat(0))
plateau_direction = (Rat(0), Rat(-1), Rat(1))
def plateau_params(c):
return tuple(plateau_base[i] + c * plateau_direction[i] for i in range(3))
induced_walls = set()
for h in arrangement_hyperplanes(4):
constant = sum(h[i] * plateau_base[i] for i in range(3)) + h[-1]
slope = sum(h[i] * plateau_direction[i] for i in range(3))
if slope:
c_wall = -constant / slope
if Rat(1, 20) <= c_wall <= Rat(1, 8):
induced_walls.add(c_wall)
assert sorted(induced_walls) == [Rat(1, 20), Rat(1, 16), Rat(1, 12), Rat(1, 8)]
# One sample in each open schedule interval. The polynomial itself is then exact on that interval.
plateau_samples = (Rat(9, 160), Rat(7, 96), Rat(5, 48))
plateau_polynomials = []
for c_sample in plateau_samples:
poly = schedule_quadratic_exact(4, plateau_params(c_sample))
coeffs = restrict_poly_to_line(poly, plateau_base, plateau_direction)
plateau_polynomials.append(coeffs)
assert coeffs == (Rat(0), Rat(0), Rat(1, 4))
for c in (Rat(1, 20), Rat(1, 16), Rat(1, 12), Rat(1, 8)):
offsets = [Rat(1, 4), Rat(1, 4) - c, c, Rat(0)]
assert exact_area_offsets(4, offsets) == Rat(1, 4)
print("induced event walls:", sorted(induced_walls))
print("restricted area polynomials:", plateau_polynomials)
print("A_4(x(c)) = 1/4 throughout [1/20, 1/8]")
Saved output 1
induced event walls: [Fraction(1, 20), Fraction(1, 16), Fraction(1, 12), Fraction(1, 8)] restricted area polynomials: [(Fraction(0, 1), Fraction(0, 1), Fraction(1, 4)), (Fraction(0, 1), Fraction(0, 1), Fraction(1, 4)), (Fraction(0, 1), Fraction(0, 1), Fraction(1, 4))] A_4(x(c)) = 1/4 throughout [1/20, 1/8]
Station Reference. Archive #5 and Evaluation #97 give the exact wall-crossing formula; Archives #36 and #96 and Evaluation #1227 connect the plateau to the certified global value.
3.3 Reflection symmetry and its first forced breaking
For a constant , define
Lemma 3.4 (area-preserving reflection). For every and , . Modulo common translation, the fixed locus consists exactly of configurations for which is independent of .
Proof. Put . Reflect horizontally and then translate it by . Its image is
the interval belonging to . Thus every horizontal slice union is carried to the new slice union by an isometry, so their lengths and integrated areas agree. The fixed-locus statement follows by solving , allowing for the irrelevant common translation.
Theorem 3.5 (restricted optimum). On the reflection-fixed locus,
attained at
Proof. In the gauge , the fixed locus is . The exact program below constructs all 37 rational event lines in the -plane, enumerates all 368 full-dimensional cells, builds and independently sweep-checks the area quadratic on each cell, and recursively minimizes it over the cell closure, including singular stationary affine sets. The least value is .
Proposition 3.6 (asymmetric witness).
has and does not lie in the reflection-fixed locus.
Proof. The exact event sweep gives . Its reflected pair sums are not constant: , whereas .
Corollary 3.7 (forced symmetry breaking). Every global minimizer for is asymmetric.
Proof. The exact gap is
Hence an admissible asymmetric configuration lies below every symmetric configuration. The unrestricted value remains open.
The next cell defines the exact reflection-fixed parameterization and full arrangement certificate used at .
Show code
Code cell 22 · In [8]
"""Exact-rational sigma-fixed certificates for small Kakeya triangle unions.
This module avoids the floating LP/QP path used by ``sigma_onset.py``. It
works directly with ``fractions.Fraction`` and enumerates the full affine line
arrangement in sigma parameters for dimensions one and two. The full
unbounded arrangement makes a separate compact-box truncation unnecessary for
the n=3,4,5 runs handled here.
"""
from __future__ import annotations
import json
import math
import os
import time
from dataclasses import dataclass
from fractions import Fraction
from itertools import combinations
from typing import Dict, Iterable, List, Optional, Sequence, Tuple
Rat = Fraction
Point = Tuple[Rat, ...]
Line2 = Tuple[Rat, Rat, Rat]
Poly = Dict[Tuple[int, ...], Rat]
def rat(value: int | Rat) -> Rat:
return value if isinstance(value, Fraction) else Fraction(value, 1)
def frac_json(x: Rat) -> str:
return f"{x.numerator}/{x.denominator}" if x.denominator != 1 else str(x.numerator)
def sigma_dim(n: int) -> int:
return n // 2
def sigma_offsets_frac(n: int, params: Sequence[Rat]) -> List[Rat]:
p = [rat(v) for v in params]
d = sigma_dim(n)
if len(p) != d:
raise ValueError(f"n={n} expects {d} sigma parameters")
c = p[0]
if n % 2 == 0:
left = [c] + p[1:d]
return left + [c - v for v in reversed(p[1:d])] + [Rat(0)]
left = [c] + p[1:d]
return left + [c / 2] + [c - v for v in reversed(p[1:d])] + [Rat(0)]
def normalize_offsets(offsets: Sequence[Rat]) -> List[Rat]:
m = min(offsets)
return [x - m for x in offsets]
def interval_at_y(n: int, offsets: Sequence[Rat], j: int, y: Rat) -> Tuple[Rat, Rat]:
x = offsets[j - 1]
left = x + Rat(j, n) * y
right = x + Rat(1, n) + Rat(j - 1, n) * y
return left, right
def union_length_at_y_frac(n: int, offsets: Sequence[Rat], y: Rat) -> Rat:
intervals = [interval_at_y(n, offsets, j, y) for j in range(1, n + 1)]
intervals.sort()
total = Rat(0)
cur_l: Optional[Rat] = None
cur_r: Optional[Rat] = None
for lft, rgt in intervals:
if cur_l is None:
cur_l, cur_r = lft, rgt
elif lft <= cur_r:
if rgt > cur_r:
cur_r = rgt
else:
total += cur_r - cur_l
cur_l, cur_r = lft, rgt
if cur_l is not None:
total += cur_r - cur_l
return total
def exact_area_offsets(n: int, offsets: Sequence[Rat]) -> Rat:
lines: List[Tuple[Rat, Rat]] = []
for j, x in enumerate(offsets, start=1):
lines.append((Rat(j, n), x))
lines.append((Rat(j - 1, n), x + Rat(1, n)))
breaks = {Rat(0), Rat(1)}
for i, (s1, b1) in enumerate(lines):
for s2, b2 in lines[i + 1 :]:
if s1 == s2:
continue
y = (b2 - b1) / (s1 - s2)
if Rat(0) < y < Rat(1):
breaks.add(y)
bp = sorted(breaks)
area = Rat(0)
for a, b in zip(bp, bp[1:]):
if a == b:
continue
area += (b - a) * (union_length_at_y_frac(n, offsets, a) + union_length_at_y_frac(n, offsets, b)) / 2
return area
def endpoint_affines(n: int) -> List[Tuple[List[Rat], Rat, Rat, Tuple[int, str]]]:
"""Return endpoint lines as intercept affine in params plus slope*y."""
d = sigma_dim(n)
rows: List[Tuple[List[Rat], Rat, Rat, Tuple[int, str]]] = []
zero = [Rat(0) for _ in range(d)]
for j in range(1, n + 1):
# Build x_j affine by probing unit sigma parameters.
coeff = []
for k in range(d):
p = [Rat(0) for _ in range(d)]
p[k] = Rat(1)
coeff.append(sigma_offsets_frac(n, p)[j - 1])
const = sigma_offsets_frac(n, [Rat(0) for _ in range(d)])[j - 1]
rows.append((coeff[:], const, Rat(j, n), (j, "L")))
rows.append((coeff[:], const + Rat(1, n), Rat(j - 1, n), (j, "R")))
return rows
def event_affines_frac(n: int) -> List[Tuple[List[Rat], Rat, Tuple[int, int]]]:
lines = endpoint_affines(n)
out: List[Tuple[List[Rat], Rat, Tuple[int, int]]] = []
for i, (ai, bi, si, _) in enumerate(lines):
for j, (aj, bj, sj, _) in enumerate(lines[i + 1 :], start=i + 1):
if si == sj:
continue
den = si - sj
out.append(([(aj[k] - ai[k]) / den for k in range(len(ai))], (bj - bi) / den, (i, j)))
return out
def canonical_line2(a: Rat, b: Rat, c: Rat) -> Optional[Line2]:
if a == 0 and b == 0:
return None
dens = [a.denominator, b.denominator, c.denominator]
lcm = 1
for den in dens:
lcm = lcm * den // math.gcd(lcm, den)
vals = [int(a * lcm), int(b * lcm), int(c * lcm)]
g = 0
for v in vals:
g = math.gcd(g, abs(v))
vals = [v // g for v in vals]
first = next(v for v in vals if v != 0)
if first < 0:
vals = [-v for v in vals]
return Rat(vals[0]), Rat(vals[1]), Rat(vals[2])
def arrangement_hyperplanes(n: int) -> List[Tuple[Rat, ...]]:
events = event_affines_frac(n)
d = sigma_dim(n)
raw: List[Tuple[Rat, ...]] = []
for a, c, _ in events:
raw.append(tuple(a + [c]))
raw.append(tuple(a + [c - 1]))
for i, (ai, ci, _) in enumerate(events):
for aj, cj, _ in events[i + 1 :]:
raw.append(tuple([ai[k] - aj[k] for k in range(d)] + [ci - cj]))
seen = set()
out: List[Tuple[Rat, ...]] = []
if d == 1:
for a, c in raw:
if a == 0:
continue
point = -c / a
if point not in seen:
seen.add(point)
out.append((point,))
out.sort()
return out
for a, b, c in raw: # d == 2
line = canonical_line2(a, b, c)
if line is not None and line not in seen:
seen.add(line)
out.append(line)
out.sort(key=lambda z: (z[0], z[1], z[2]))
return out
def sign_of_line(line: Line2, p: Tuple[Rat, Rat]) -> int:
v = line[0] * p[0] + line[1] * p[1] + line[2]
if v == 0:
raise ValueError("sample lies on arrangement line")
return 1 if v > 0 else -1
def sign_vector(lines: Sequence[Line2], p: Tuple[Rat, Rat]) -> Tuple[int, ...]:
return tuple(sign_of_line(line, p) for line in lines)
def enumerate_cells_2d(lines_in: Sequence[Tuple[Rat, ...]]) -> Dict[Tuple[int, ...], Tuple[Rat, Rat]]:
lines: List[Line2] = [(a, b, c) for a, b, c in lines_in] # type: ignore[misc]
cells: Dict[Tuple[int, ...], Tuple[Rat, Rat]] = {}
if not lines:
return {tuple(): (Rat(0), Rat(0))}
for idx, line in enumerate(lines):
a, b, c = line
if b != 0:
p0 = (Rat(0), -c / b)
else:
p0 = (-c / a, Rat(0))
direction = (b, -a)
ts = set()
for j, other in enumerate(lines):
if j == idx:
continue
aa, bb, cc = other
alpha = aa * direction[0] + bb * direction[1]
beta = aa * p0[0] + bb * p0[1] + cc
if alpha != 0:
ts.add(-beta / alpha)
ordered = sorted(ts)
samples: List[Rat] = []
if not ordered:
samples = [Rat(0)]
else:
samples.append(ordered[0] - 1)
for u, v in zip(ordered, ordered[1:]):
samples.append((u + v) / 2)
samples.append(ordered[-1] + 1)
for t in samples:
base = (p0[0] + t * direction[0], p0[1] + t * direction[1])
eps_limit: Optional[Rat] = None
for j, other in enumerate(lines):
if j == idx:
continue
val = other[0] * base[0] + other[1] * base[1] + other[2]
delta = other[0] * a + other[1] * b
if delta != 0:
bound = abs(val / delta) / 2
if bound > 0 and (eps_limit is None or bound < eps_limit):
eps_limit = bound
eps = eps_limit if eps_limit is not None else Rat(1, 2)
if eps == 0:
eps = Rat(1, 1024)
for side in (-1, 1):
p = (base[0] + side * eps * a, base[1] + side * eps * b)
try:
sv = sign_vector(lines, p)
except ValueError:
continue
cells.setdefault(sv, p)
return cells
def enumerate_cells_1d(points_in: Sequence[Tuple[Rat, ...]]) -> Dict[Tuple[int, ...], Tuple[Rat]]:
points = [p[0] for p in points_in]
samples: List[Rat] = []
if not points:
samples = [Rat(0)]
else:
samples.append(points[0] - 1)
for a, b in zip(points, points[1:]):
samples.append((a + b) / 2)
samples.append(points[-1] + 1)
cells = {}
for x in samples:
sv = tuple(1 if x - p > 0 else -1 for p in points)
cells[sv] = (x,)
return cells
def affine_endpoint_value(n: int, endpoint_index: int, params: Point, y: Rat) -> Rat:
coeff, const, slope, _ = endpoint_affines(n)[endpoint_index]
return sum(coeff[i] * params[i] for i in range(len(params))) + const + slope * y
def poly_const(d: int, c: Rat) -> Poly:
return {(0,) * d: c} if c else {}
def poly_affine(coeff: Sequence[Rat], const: Rat) -> Poly:
d = len(coeff)
out: Poly = {}
if const:
out[(0,) * d] = const
for i, a in enumerate(coeff):
if a:
exp = [0] * d
exp[i] = 1
out[tuple(exp)] = out.get(tuple(exp), Rat(0)) + a
return {k: v for k, v in out.items() if v}
def poly_add(a: Poly, b: Poly, scale: Rat = Rat(1)) -> Poly:
out = dict(a)
for k, v in b.items():
out[k] = out.get(k, Rat(0)) + scale * v
if out[k] == 0:
del out[k]
return out
def poly_mul(a: Poly, b: Poly) -> Poly:
out: Poly = {}
for ea, va in a.items():
for eb, vb in b.items():
exp = tuple(ea[i] + eb[i] for i in range(len(ea)))
if sum(exp) > 2:
raise ValueError("degree exceeded quadratic")
out[exp] = out.get(exp, Rat(0)) + va * vb
return {k: v for k, v in out.items() if v}
def poly_eval(poly: Poly, p: Point) -> Rat:
total = Rat(0)
for exp, coef in poly.items():
term = coef
for i, e in enumerate(exp):
if e == 1:
term *= p[i]
elif e == 2:
term *= p[i] * p[i]
total += term
return total
def union_expression_on_slab(n: int, sample: Point, y_mid: Rat) -> Tuple[List[Rat], Rat, Rat]:
lines = endpoint_affines(n)
intervals = []
for j in range(1, n + 1):
li = 2 * (j - 1)
ri = li + 1
lv = affine_endpoint_value(n, li, sample, y_mid)
rv = affine_endpoint_value(n, ri, sample, y_mid)
intervals.append((lv, rv, li, ri))
intervals.sort(key=lambda row: (row[0], row[1], row[2], row[3]))
clusters: List[Tuple[int, int]] = []
cur_l: Optional[int] = None
cur_r: Optional[int] = None
cur_rv: Optional[Rat] = None
for lv, rv, li, ri in intervals:
if cur_l is None:
cur_l, cur_r, cur_rv = li, ri, rv
elif lv <= cur_rv:
if rv > cur_rv:
cur_r, cur_rv = ri, rv
else:
clusters.append((cur_l, cur_r)) # type: ignore[arg-type]
cur_l, cur_r, cur_rv = li, ri, rv
if cur_l is not None:
clusters.append((cur_l, cur_r)) # type: ignore[arg-type]
d = len(sample)
a = [Rat(0) for _ in range(d)]
b = Rat(0)
cy = Rat(0)
for li, ri in clusters:
ar, br, sr, _ = lines[ri]
al, bl, sl, _ = lines[li]
for k in range(d):
a[k] += ar[k] - al[k]
b += br - bl
cy += sr - sl
return a, b, cy
def schedule_quadratic_exact(n: int, sample: Point) -> Poly:
d = len(sample)
events = event_affines_frac(n)
vals = [(sum(a[i] * sample[i] for i in range(d)) + c, a, c) for a, c, _ in events]
active = [(v, a, c) for v, a, c in vals if Rat(0) < v < Rat(1)]
active.sort(key=lambda row: row[0])
breaks = [(Rat(0), [Rat(0) for _ in range(d)], Rat(0))] + active + [(Rat(1), [Rat(0) for _ in range(d)], Rat(1))]
poly: Poly = {}
for left, right in zip(breaks, breaks[1:]):
v0, a0, c0 = left
v1, a1, c1 = right
if v0 == v1:
continue
y_mid = (v0 + v1) / 2
u_a, u_b, u_cy = union_expression_on_slab(n, sample, y_mid)
h0 = poly_affine(a0, c0)
h1 = poly_affine(a1, c1)
dh = poly_add(h1, h0, scale=-1)
u_no_y = poly_affine(u_a, u_b)
term = poly_mul(u_no_y, dh)
h1_sq = poly_mul(h1, h1)
h0_sq = poly_mul(h0, h0)
yterm = poly_add(h1_sq, h0_sq, scale=-1)
yterm = {k: u_cy * v / 2 for k, v in yterm.items()}
poly = poly_add(poly, term)
poly = poly_add(poly, yterm)
return poly
def poly_gradient_linear(poly: Poly, d: int) -> Tuple[List[List[Rat]], List[Rat]]:
mat = [[Rat(0) for _ in range(d)] for _ in range(d)]
vec = [Rat(0) for _ in range(d)]
for exp, coef in poly.items():
for i in range(d):
if exp[i] == 0:
continue
new_exp = list(exp)
new_exp[i] -= 1
deriv_coef = coef * exp[i]
if sum(new_exp) == 0:
vec[i] += deriv_coef
else:
j = new_exp.index(1)
mat[i][j] += deriv_coef
return mat, vec
def solve_linear_2(mat: List[List[Rat]], rhs: List[Rat]) -> Optional[Tuple[Rat, Rat]]:
det = mat[0][0] * mat[1][1] - mat[0][1] * mat[1][0]
if det == 0:
return None
x = (rhs[0] * mat[1][1] - mat[0][1] * rhs[1]) / det
y = (mat[0][0] * rhs[1] - rhs[0] * mat[1][0]) / det
return x, y
def feasible_1d(points: Sequence[Tuple[Rat, ...]], signs: Tuple[int, ...], p: Point) -> bool:
x = p[0]
return all((x - pt[0]) * s >= 0 for pt, s in zip(points, signs))
def feasible_2d(lines: Sequence[Tuple[Rat, ...]], signs: Tuple[int, ...], p: Point) -> bool:
x, y = p
return all(s * (line[0] * x + line[1] * y + line[2]) >= 0 for line, s in zip(lines, signs))
def q1_coeffs(poly: Poly) -> Tuple[Rat, Rat, Rat]:
return poly.get((2,), Rat(0)), poly.get((1,), Rat(0)), poly.get((0,), Rat(0))
def minimize_cell_1d(poly: Poly, points: Sequence[Tuple[Rat, ...]], signs: Tuple[int, ...]) -> Tuple[Rat, Point]:
lo: Optional[Rat] = None
hi: Optional[Rat] = None
for pt, s in zip(points, signs):
p = pt[0]
if s > 0:
lo = p if lo is None or p > lo else lo
else:
hi = p if hi is None or p < hi else hi
candidates: List[Point] = []
if lo is not None:
candidates.append((lo,))
if hi is not None:
candidates.append((hi,))
a, b, _ = q1_coeffs(poly)
if a > 0:
x0 = -b / (2 * a)
candidates.append((x0,))
if lo is None and hi is None:
candidates.append((Rat(0),))
best = None
for p in candidates:
if feasible_1d(points, signs, p):
val = poly_eval(poly, p)
if best is None or val < best[0]:
best = (val, p)
if best is None:
raise RuntimeError("no 1D minimizer candidate")
return best
def line_interval(lines: Sequence[Tuple[Rat, ...]], signs: Tuple[int, ...], boundary_idx: int) -> Tuple[Tuple[Rat, Rat], Tuple[Rat, Rat], Optional[Rat], Optional[Rat]]:
a, b, c = lines[boundary_idx]
if b != 0:
p0 = (Rat(0), -c / b)
else:
p0 = (-c / a, Rat(0))
direction = (b, -a)
lo: Optional[Rat] = None
hi: Optional[Rat] = None
for idx, (line, s) in enumerate(zip(lines, signs)):
if idx == boundary_idx:
continue
aa, bb, cc = line
alpha = s * (aa * direction[0] + bb * direction[1])
beta = s * (aa * p0[0] + bb * p0[1] + cc)
if alpha == 0:
if beta < 0:
return p0, direction, Rat(1), Rat(0)
continue
bound = -beta / alpha
if alpha > 0:
lo = bound if lo is None or bound > lo else lo
else:
hi = bound if hi is None or bound < hi else hi
return p0, direction, lo, hi
def restrict_poly_to_line(poly: Poly, p0: Tuple[Rat, Rat], direction: Tuple[Rat, Rat]) -> Tuple[Rat, Rat, Rat]:
vals = []
for t in (Rat(0), Rat(1), Rat(2)):
p = (p0[0] + t * direction[0], p0[1] + t * direction[1])
vals.append(poly_eval(poly, p))
c = vals[0]
a = (vals[2] - 2 * vals[1] + vals[0]) / 2
b = vals[1] - vals[0] - a
return a, b, c
def minimize_cell_2d(poly: Poly, lines: Sequence[Tuple[Rat, ...]], signs: Tuple[int, ...]) -> Tuple[Rat, Point]:
candidates: List[Point] = []
mat, vec = poly_gradient_linear(poly, 2)
stat = solve_linear_2(mat, [-vec[0], -vec[1]])
if stat is not None:
candidates.append(stat)
for i in range(len(lines)):
p0, direction, lo, hi = line_interval(lines, signs, i)
if lo is not None and hi is not None and lo > hi:
continue
a, b, _ = restrict_poly_to_line(poly, p0, direction)
ts: List[Rat] = []
if lo is not None:
ts.append(lo)
if hi is not None:
ts.append(hi)
if a > 0:
t0 = -b / (2 * a)
if (lo is None or t0 >= lo) and (hi is None or t0 <= hi):
ts.append(t0)
elif a == 0 and b == 0:
if lo is not None:
ts.append(lo)
elif hi is not None:
ts.append(hi)
else:
ts.append(Rat(0))
for t in ts:
candidates.append((p0[0] + t * direction[0], p0[1] + t * direction[1]))
best = None
for p in candidates:
if feasible_2d(lines, signs, p):
val = poly_eval(poly, p)
if best is None or val < best[0]:
best = (val, p)
if best is None:
raise RuntimeError("no 2D minimizer candidate")
return best
def certify_sigma_exact(n: int) -> Dict[str, object]:
start = time.time()
d = sigma_dim(n)
hyper = arrangement_hyperplanes(n)
if d == 1:
cells = enumerate_cells_1d(hyper)
elif d == 2:
cells = enumerate_cells_2d(hyper)
else:
raise ValueError("exact certifier currently supports sigma dimension <= 2")
best_val: Optional[Rat] = None
best_params: Optional[Point] = None
sample_mismatches = []
processed = 0
for signs, sample in cells.items():
poly = schedule_quadratic_exact(n, sample)
sample_offsets = sigma_offsets_frac(n, sample)
sample_area = exact_area_offsets(n, sample_offsets)
formula_area = poly_eval(poly, sample)
if sample_area != formula_area:
sample_mismatches.append(
{
"sample": [frac_json(x) for x in sample],
"sweep": frac_json(sample_area),
"formula": frac_json(formula_area),
}
)
continue
if d == 1:
val, params = minimize_cell_1d(poly, hyper, signs)
else:
val, params = minimize_cell_2d(poly, hyper, signs)
if best_val is None or val < best_val:
best_val, best_params = val, params
processed += 1
if sample_mismatches:
raise RuntimeError(f"quadratic/sweep mismatch in {len(sample_mismatches)} cells")
if best_val is None or best_params is None:
raise RuntimeError("no cells processed")
best_offsets = sigma_offsets_frac(n, best_params)
best_sweep = exact_area_offsets(n, best_offsets)
return {
"n": n,
"sigma_dim": d,
"area": frac_json(best_val),
"area_float": float(best_val),
"params": [frac_json(x) for x in best_params],
"offsets": [frac_json(x) for x in normalize_offsets(best_offsets)],
"raw_offsets_xn_zero": [frac_json(x) for x in best_offsets],
"exact_sweep_at_min": frac_json(best_sweep),
"min_formula_matches_sweep": best_sweep == best_val,
"complete": True,
"completion_scope": "full_unbounded_parameter_arrangement",
"cell_count": len(cells),
"processed_cell_count": processed,
"hyperplane_count": len(hyper),
"wall_clock": time.time() - start,
}
def keich_positions(n: int) -> List[float]:
k = int(math.log2(n))
positions = []
for i in range(n):
bits = [int(ch) for ch in bin(i)[2:]]
bits = [0] * (k - len(bits)) + bits
total = 0.0
for j, eps in enumerate(bits, start=1):
total += (1 - j) * eps * (2.0 ** -j) / k
positions.append(total)
m = min(positions)
return [x - m for x in positions]
The next cell adds the singular-safe boundary minimizer used by the reflection-fixed certificate.
Show code
Code cell 24 · In [9]
_boundary_sigma_minimize_cell_2d = minimize_cell_2d
def minimize_cell_2d_robust(poly, lines, signs):
candidates = []
stationary = stationary_candidate(poly, lines, signs, 2)
if stationary is not None:
candidates.append(stationary)
try:
_value, boundary_candidate = _boundary_sigma_minimize_cell_2d(poly, lines, signs)
candidates.append(boundary_candidate)
except RuntimeError:
pass
best = None
for p in candidates:
if feasible_2d(lines, signs, p):
value = poly_eval(poly, p)
if best is None or value < best[0]:
best = value, p
if best is None:
raise RuntimeError("no exact two-dimensional minimizer candidate")
return best
minimize_cell_2d = minimize_cell_2d_robust
The next cell runs the complete reflection-fixed certificate and checks the asymmetric witness and the exact symmetry gap.
Show code
Code cell 26 · In [10]
n5_symmetric_certificate = certify_sigma_exact(5)
assert n5_symmetric_certificate["complete"]
assert n5_symmetric_certificate["hyperplane_count"] == 37
assert n5_symmetric_certificate["cell_count"] == 368
assert n5_symmetric_certificate["processed_cell_count"] == 368
assert n5_symmetric_certificate["area"] == "7/30"
assert n5_symmetric_certificate["exact_sweep_at_min"] == "7/30"
sym = [Fraction(int(v), int(payload["n5_symmetric_denominator"]))
for v in payload["n5_symmetric_numerators"]]
asym = [Fraction(int(v), int(payload["n5_asymmetric_denominator"]))
for v in payload["n5_asymmetric_numerators"]]
assert area_exact(sym) == Fraction(7, 30)
assert area_exact(asym) == Fraction(14, 61)
assert asym[0] + asym[4] != asym[1] + asym[3]
gap = Fraction(7, 30) - Fraction(14, 61)
assert gap == Fraction(7, 1830) and gap > 0
print("fixed-locus lines:", n5_symmetric_certificate["hyperplane_count"])
print("fixed-locus cells:", n5_symmetric_certificate["cell_count"])
print("symmetric minimum:", n5_symmetric_certificate["area"])
print("asymmetric witness:", Fraction(14, 61))
print("exact gap:", gap)
Saved output 1
fixed-locus lines: 37 fixed-locus cells: 368 symmetric minimum: 7/30 asymmetric witness: 14/61 exact gap: 7/1830
Station Reference. Archives #36 and #94; Evaluations #1211, #1212, #1220, #1223, and #1238; Question #36. The fixed-locus certificate is the exact rational enumeration from Evaluation #1220, independently rerun in Evaluation #1223.
Related Work. Falconer (1985) records the classical value , while Schoenberg (1962) gives an construction of area . The results here certify the global optima at , determine the full optimal plateau, and establish forced symmetry breaking at .
4. Additional findings
These are additional results that may be interesting but are not included in the spotlight.
4.1 Exact finite area of the Schoenberg–Keich construction
The dyadic construction of Schoenberg, in the form rederived and analyzed by Keich, gives the classical reciprocal-logarithmic upper bound, but its finite area contains substantial lower-order corrections. We compute that area exactly:
Our proof decomposes the sweep into binary bands and determines the precise height at which each scale changes from merging to separation. The formula proves that the leading constant is exactly one and separates finite-size improvements from genuine improvements in asymptotic behavior.
Let . For , write . The Keich offset assigned to triangle is
Theorem 4.1 (exact finite-area formula).
In particular, .
Proof. At height , remove the common translation and multiply horizontal lengths by ; put . The interval width becomes , while the -th binary displacement is
Process the binary displacements from fine to coarse scale. A bit merges when its translated copy overlaps the current hull, and separates when the two copies are disjoint. On band , where and , the bits merge, the bits separate, and bit changes from merged to separated once, at . These assertions follow from the geometric-sum comparisons
and, writing with ,
Thus band has a constant plateau followed by one linear ramp. On the plateau, the separated bits create congruent components, each of scaled length . If , their total unscaled length is
The plateau lasts units of height. It then decreases with slope for the remaining units. The plateau rectangle minus the triangular ramp correction gives the band contribution
Summing and using and gives the displayed formula. Multiplication by then makes the limit immediate.
The next cell generates the exact Keich offsets for and verifies the formula by the independent event sweep.
Show code
Code cell 30 · In [11]
def keich_offsets_exact(k):
n = 2 ** k
out = []
for i in range(n):
bits = [int(ch) for ch in f"{i:0{k}b}"]
out.append(sum(Fraction(1 - r, k) * bits[r - 1] * Fraction(1, 2 ** r)
for r in range(1, k + 1)))
return out
def keich_formula(k):
return Fraction(1, k * k) * (
k - Fraction(7, 6) + Fraction(1, 2 ** k) + Fraction(1, 6 * 4 ** k)
)
print(f"{'k':>2} {'n':>4} {'exact area':>18} {'decimal':>14}")
for k in range(1, 8):
swept = area_exact(keich_offsets_exact(k))
formula = keich_formula(k)
assert swept == formula
print(f"{k:2d} {2**k:4d} {str(formula):>18} {float(formula):14.10f}")
Saved output 1
k n exact area decimal 1 2 3/8 0.3750000000 2 4 35/128 0.2734375000 3 8 251/1152 0.2178819444 4 16 1483/8192 0.1810302734 5 32 1583/10240 0.1545898438
Saved output 2
6 64 13241/98304 0.1346944173
Saved output 3
7 128 191403/1605632 0.1192072654
Station Reference. Archive #42 and its exact-asymptotic addendum; Evaluations #486, #489, #493, #500, #504, and #519.
Related Work. The construction is due to Schoenberg and was rederived and analyzed asymptotically by Keich (1999). We derive its exact finite area and leading constant.
4.2 Prescribed sweep paths are linear-feasibility problems
A triangle configuration is governed by the order in which its moving endpoints cross and its adjacent gaps reach the triangle width. We prove that any prescribed boundary-labelled sweep path is realizable exactly when a finite system of linear equalities and inequalities is feasible, and that its area is quadratic on this feasible region. This separates the combinatorial choice of a crossing schedule from the continuous optimization of its offsets and event heights, turning each proposed schedule into a finite, generally nonconvex quadratic optimization problem.
A boundary-labelled sweep path specifies finitely many height slabs, the strict order of all affine endpoints in the interior of each slab, which adjacent equal-width gaps are capped by the interval width, and the complete set of endpoint equalities at every slab boundary. Several equalities may occur at the same boundary.
Theorem 4.2 (lifted linear feasibility and quadratic area). The realizability of a prescribed boundary-labelled sweep path is exactly the feasibility of a finite system of linear equalities, weak inequalities, and strict inequalities in the offsets and boundary heights . On the feasible set, the area is a polynomial of degree at most two in .
Proof. Every endpoint has the form , with affine in . An endpoint order or gap-width condition on a slab is an affine inequality in . Because its left side is affine in , it holds throughout the slab exactly when it holds at both boundary heights. A labelled boundary event supplies the full finite collection of affine equalities assigned at that boundary; all remaining labelled orders supply the corresponding weak or strict endpoint inequalities. Adding completes the finite linear system; strict inequalities distinguish open cells and may equivalently use one common positive slack variable.
Once these labels are fixed, the slice union length on slab has the form , where is affine. Its integral is
which has degree at most two. Summing over slabs proves the statement.
Station Reference. Archives #11 and #12; Evaluation #141; accepted solution to Question #2, “Active-Gap Path Realizability for Fixed-Slope Equal-Width Sweeps.”
Related Work. Hyperplane arrangements and their polyhedral cells are classical; see Zaslavsky (1975). The theorem applies this structure to the finite triangle model by lifting each fixed endpoint-and-gap schedule to a system of linear feasibility conditions.
4.3 An exact single-focus family
We analyze the reversal family in which all left endpoints pass through one common focus inside the strip. We derive its exact area as a function of the spacing, prove that the unique optimum occurs at , and obtain the value . This identifies the best possible performance among evenly spaced reversal configurations with an in-strip focus and its limiting barrier . Since the strongest finite constructions already pass below that barrier, their improvement must come from more complicated crossing patterns in which changing adjacency creates nonlocal overlap savings.
For , assign the evenly spaced offsets in reverse order:
The left endpoints are equally spaced at every height, and all meet at the single height .
Theorem 4.3 (reversal-ray area law). For every and ,
The function is strictly convex and has its unique minimum at , where
Proof. Consecutive left endpoints differ by , so their spacing magnitude is . The union of equal-width intervals with equally spaced left endpoints has length
Put . Splitting at and gives
Adding the integral of , which is , proves the formula. The derivative of is , so its unique minimum is at , yielding the stated value.
The next cell checks the area formula at five exact rational spacings for each .
Show code
Code cell 35 · In [12]
def reversal_offsets(n, g):
return [Fraction(n - j) * g for j in range(1, n + 1)]
def reversal_formula(n, g):
a = n * g
return Fraction(1, 2 * n) + Fraction(n - 1, n) * (
a * a / 2 + (1 - a) * (1 - a) / 4
)
for n in range(2, 13):
probes = (Fraction(0), Fraction(1, 6 * n), Fraction(1, 3 * n),
Fraction(1, 2 * n), Fraction(1, n))
for g in probes:
assert area_exact(reversal_offsets(n, g)) == reversal_formula(n, g)
print("reversal-ray formula verified exactly for n=2,...,12 at five rational spacings each")
Saved output 1
reversal-ray formula verified exactly for n=2,...,12 at five rational spacings each
Station Reference. Archives #81, #83, and #86; Evaluations #1019 and #1102, with supporting crossing-family studies in Evaluations #253, #630, #633, and #639.
Related Work. Falconer (1985) describes the finite triangle model and its classical constructions. The theorem derives the exact reversal-ray formula. This family attains the known optima at ; the asymmetric witness and stronger large- constructions lie below it.