files / problems / lonely-runner-fourteen / compute / q2 / check_cells.py
#!/usr/bin/env python3
"""check_cells.py -- independent check of the cell/window reduction in cells.c.
Written from the definitions, sharing no code with the C. Four checks:
cells B_i(t) = floor(m*frac(i*t)) really is constant on each half-open cell
[c_r, c_{r+1}) of the arrangement generated by the breakpoints
union_i (1/(m*i))Z, sampled at random interior points.
window for every q, every r and every integer p >= q, the half-open window
[c_r, c_r + 1/q) contains a point of (1/p)Z, and that point's cell
lies in W_r = { cells meeting the window }. This is the whole
content of the reduction: it is why an unsaved v must satisfy the
window constraint.
schema schema.py's interval arithmetic for the family v = e_a + e_b agrees,
pair by pair and p by p, with exhausting all (s,j) directly.
brute T2(k,p) decided by enumerating all of Z_m^k directly from ST26's
definitions -- no covering, no search -- for small k.
scan the certified claim "T2(k,p) holds for every p >= q0" is not
contradicted by brute force over a range of p.
usage: python3 check_cells.py --selftest
python3 check_cells.py --k 5 --p 91 # one brute-force decision
"""
import argparse, random
from math import gcd
def lcm_upto(k):
L = 1
for i in range(1, k + 1):
L = L * i // gcd(L, i)
return L
def cells(k):
"""half-open cells [c_r, c_{r+1}) as integer numerators over D."""
m, L = k + 1, lcm_upto(k)
D = m * L
cut = set()
for i in range(1, k + 1):
step = L // i # D/(m*i)
cut.update(range(0, D, step))
return D, sorted(cut)
def Bvec(k, num, den):
"""B_i at the exact rational time num/den."""
m = k + 1
return tuple(m * ((i * num) % den) // den for i in range(1, k + 1))
def check_cells_constant(k, trials=200, seed=1):
"""B is constant on each half-open cell."""
rnd = random.Random(seed)
D, cut = cells(k)
n = len(cut)
bad = 0
for r in range(n):
lo, hi = cut[r], (cut[r + 1] if r + 1 < n else D)
ref = Bvec(k, lo, D)
for _ in range(trials):
# a random rational strictly inside [lo, hi), fine denominator
f = rnd.randrange(0, 10**6)
num = lo * 10**6 + (hi - lo) * f
if Bvec(k, num, D * 10**6) != ref:
bad += 1
return n, bad
def check_window(k, q, plist):
"""every window [c_r, c_r+1/q) meets (1/p)Z, in a cell of W_r."""
D, cut = cells(k)
n = len(cut)
bad = 0
for r in range(n):
lo = cut[r]
# W_r = cells r..hi with q*(c_{r'} - c_r) < D (wrapping)
hi = r
while q * (cut[(hi + 1) % n] + (D if hi + 1 >= n else 0) - lo) < D:
hi += 1
W = set((x % n) for x in range(r, hi + 1))
for p in plist:
# smallest j with j/p >= lo/D
j = -(-lo * p // D)
if q * (j * D - lo * p) >= D * p: # j/p - lo/D >= 1/q
bad += 1 # window missed the grid
continue
c = j * D // p # numerator of j/p over D
# locate its cell
import bisect
cr = bisect.bisect_right(cut, c % D) - 1
if cr % n not in W:
bad += 1
return n, bad
def eligible(k, v):
"""ST26 Def 2.1: needs a witness iff some v_i = 0 and not gcd-proper."""
m = k + 1
if all(x != 0 for x in v):
return False
q = 2
x = m
ps = []
while q * q <= x:
if x % q == 0:
ps.append(q)
while x % q == 0:
x //= q
q += 1
if x > 1:
ps.append(x)
for pr in ps:
if sum(1 for a in v if a % pr) <= 1:
return False # gcd-proper, no witness needed
return True
def saved(k, p, v):
"""some (s,j) puts s*v + r_k(j/p) inside {1,...,k-1}^k."""
m = k + 1
for j in range(p):
B = Bvec(k, j, p)
for s in range(m):
if all((s * v[i] + B[i]) % m not in (0, m - 1) for i in range(k)):
return True
return False
def brute(k, p, report=False):
"""all obstructions in Z_m^k; T2(k,p) holds iff there are none."""
m = k + 1
bad = []
v = [0] * k
def rec(i):
if i == k:
t = tuple(v)
if eligible(k, t) and not saved(k, p, t):
bad.append(t)
return
for a in range(m):
v[i] = a
rec(i + 1)
rec(0)
if report and bad:
print(" first obstruction:", bad[0])
return len(bad)
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--k", type=int, default=5)
ap.add_argument("--p", type=int)
ap.add_argument("--selftest", action="store_true")
a = ap.parse_args()
if a.p:
n = brute(a.k, a.p, report=True)
print("k=%d p=%d obstructions=%d -> T2 %s"
% (a.k, a.p, n, "FAILS" if n else "HOLDS"))
return
if not a.selftest:
print("nothing to do; pass --selftest or --k/--p")
return
print("[cells] B constant on every half-open cell")
for k in (5, 7, 13):
n, bad = check_cells_constant(k)
print(" k=%-3d cells=%-5d mismatches=%d %s" % (k, n, bad, "ok" if not bad else "FAIL"))
print("[window] every window [c_r, c_r+1/q) meets (1/p)Z inside W_r")
for k, q in ((5, 90), (7, 84), (13, 200), (13, 300)):
pl = list(range(q, q + 40)) + [q * 3, q * 7 + 1, 10007]
n, bad = check_window(k, q, pl)
print(" k=%-3d q=%-4d cells=%-5d p tested=%-3d violations=%d %s"
% (k, q, n, len(pl), bad, "ok" if not bad else "FAIL"))
print("[schema] the e_a+e_b family: interval arithmetic vs all (s,j)")
import schema
from itertools import combinations
for k in (5, 13):
pl = [p for p in range(6, 130)]
mism = n = 0
for T in combinations(range(1, k + 1), 2):
Z = [i for i in range(1, k + 1) if i not in T]
D, good = schema.good_set(k, Z)
v = tuple(1 if i in T else 0 for i in range(1, k + 1))
if not eligible(k, v):
continue
for p in pl:
n += 1
if (not schema.has_grid_point(D, good, p)) != (not saved(k, p, v)):
mism += 1
print(" k=%-3d pairs x p checked=%-6d mismatches=%d %s"
% (k, n, mism, "ok" if not mism else "FAIL"))
print("[brute] T2(5,p) by direct enumeration of all of Z_6^5")
fails = [p for p in range(6, 121) if brute(5, p)]
print(" p in 6..120 that FAIL:", fails)
print(" last failure p=%d, so T2(5,p) holds for every p in %d..120"
% (max(fails), max(fails) + 1))
print(" cells.c certifies every p >= 90:", "consistent" if max(fails) < 90 else "CONTRADICTED")
print("[scan] brute force above the certified threshold, k=5, q0=90")
bad = [p for p in range(90, 400) if brute(5, p)]
print(" failures in 90..399:", bad if bad else "none -- consistent")
if __name__ == "__main__":
main()