norms.py

1import itertools2import math3import sys4import time5from fractions import Fraction67import numpy as np8from mpmath import iv, mp, mpf, mpc, quad, expjpi, cos as mpcos, log as mplog910# DESIGNS1112def cube_minus(q, D, drop):13    return (q, D, [v for v in itertools.product(range(q), repeat=D) if v != drop])1415def base_missing(b, a0):16    return (b, 1, [(a,) for a in range(b) if a != a0])1718GASKET = cube_minus(2, 2, (1, 1))19CARPET = cube_minus(3, 2, (1, 1))20GASKET4 = base_missing(4, 3)21CARPET9 = base_missing(9, 4)2223ROSTER = [24    ("gasket 2D", GASKET),25    ("carpet 2D", CARPET),26    ("gasket base 4", GASKET4),27    ("carpet base 9", CARPET9),28    ("base 10 missing 0", base_missing(10, 0)),29    ("base 10 missing 1", base_missing(10, 1)),30    ("base 10 missing 5", base_missing(10, 5)),31]3233# CHECKS3435FAILS = []3637def chk(value, want, tol):38    if abs(value - want) <= tol:39        return "ok"40    FAILS.append((value, want))41    return "MISMATCH"4243def chk_le(value, cap, name):44    if value <= cap:45        return "ok"46    FAILS.append((name, value, cap))47    return "MISMATCH"4849def close(t0):50    print(f"runtime {time.time() - t0:.1f}s")51    if FAILS:52        raise SystemExit(f"{len(FAILS)} rows off target: {FAILS}")53    print("every asserted row matches")5455# FLOAT TRANSFORM5657def mask(F, D):58    k = len(F)59    def f(t):60        acc = np.zeros(np.broadcast(*t).shape, dtype=complex)61        for v in F:62            acc = acc + np.exp(2j * np.pi * sum(v[i] * t[i] for i in range(D)))63        return np.abs(acc) / k64    return f, k6566def lip_of(F):67    return (2 * math.pi / len(F)) * sum(sum(abs(c) for c in v) for v in F)6869# ONE DIGIT, EXACT7071def digit_sums():72    for name, (q, D, F) in ROSTER:73        k = len(F)74        assert k == q ** D - 175        exact = Fraction(1) + Fraction(q ** D - 1, k)76        f, _ = mask(F, D)77        axes = np.meshgrid(*([np.arange(q)] * D), indexing="ij")78        num = float(f([a / q for a in axes]).sum())79        e = math.log(float(exact)) / math.log(q)80        print(f"{name}: q={q} D={D} k={k} one-digit l1 sum {exact} exact, scan {num:.12f} {chk(num, float(exact), 1e-12)}, log_q sum {e:.6f}, dim {math.log(k) / math.log(q):.6f}")8182def sqrt_split(n):83    a, b = 1, n84    d = 285    while d * d <= b:86        while b % (d * d) == 0:87            b //= d * d88            a *= d89        d += 190    return a, b9192def gasket_anchor():93    q, D, F = GASKET94    c4 = [1, 0, -1, 0]95    terms = {}96    for i in itertools.product(range(4), repeat=2):97        s = []98        for j in (0, 1):99            u = ((2 ** j) * i[0] % 4, (2 ** j) * i[1] % 4)100            s.append(3 + 2 * c4[u[0]] + 2 * c4[u[1]] + 2 * c4[(u[0] - u[1]) % 4])101        a, b = sqrt_split(s[0] * s[1])102        terms[b] = terms.get(b, 0) + a103    body = " + ".join(f"{a} sqrt {b}" if b > 1 else f"{a}" for b, a in sorted(terms.items()))104    val = sum(a * math.sqrt(b) for b, a in terms.items()) / 9105    g = Fraction(sum(terms.values()) if len(terms) == 1 else 0)106    f, k = mask(F, D)107    n = 4108    axes = [a.ravel() for a in np.meshgrid(*([np.arange(n)] * D), indexing="ij")]109    prod = np.ones(axes[0].shape)110    for j in range(2):111        prod *= f([((2 ** j) * a % n) / n for a in axes])112    scan = float(prod.sum())113    print(f"gasket 2D Sigma_2(0) = ({body})/9 = (8 + 2 sqrt 5)/3 exact = {val:.12f}, scan {scan:.12f} {chk(scan, val, 1e-12)}, threshold q^(ND/2) = 4, margin {val - 4:.6f}, exponent {math.log(val) / (2 * math.log(2)):.6f}")114115def entropy_bound():116    mp.dps = 30117    m1 = 2 * quad(lambda u: mplog(2 * mpcos(mp.pi * u)), [0, mpf(1) / 3])118    q, D, F = GASKET119    a = float(m1) - math.log(3)120    print(f"gasket 2D entropy bound: int_(T^2) log|1 + e(t1) + e(t2)| = {float(m1):.12f}, int log|hat F| = {a:.12f}, Jensen gives alpha_1 >= D + that/log q = {2 + a / math.log(2):.6f}, threshold D/2 = 1.0")121122# GRID SUMS123124def grid_sum(q, D, F, L):125    f, k = mask(F, D)126    n = q ** L127    axes = np.meshgrid(*([np.arange(n)] * D), indexing="ij")128    prod = np.ones(axes[0].shape)129    for j in range(L):130        prod *= f([((q ** j) * a % n) / n for a in axes])131    return float(prod.sum())132133GRID = [134    ("gasket 2D", GASKET, [1.0000, 1.0278, 1.0388, 1.0446, 1.0482, 1.0505, 1.0522, 1.0535, 1.0545, 1.0553]),135    ("carpet 2D", CARPET, [0.6309, 0.6931, 0.7174, 0.7301, 0.7378, 0.7429]),136    ("gasket base 4", GASKET4, [0.5000, 0.4903, 0.4876, 0.4862, 0.4853, 0.4848, 0.4844, 0.4841, 0.4838, 0.4836]),137    ("carpet base 9", CARPET9, [0.3155, 0.3298, 0.3344, 0.3367, 0.3380, 0.3389]),138    ("base 10 missing 0", base_missing(10, 0), [None, None, None, None, None, 0.3086]),139    ("base 10 missing 1", base_missing(10, 1), [None, None, None, None, None, 0.3328]),140    ("base 10 missing 5", base_missing(10, 5), [None, None, None, None, None, 0.3420]),141]142143def grids(rows):144    for name, (q, D, F), want in rows:145        for L, w in enumerate(want, start=1):146            S = grid_sum(q, D, F, L)147            e = math.log(S) / (L * math.log(q))148            tag = "" if w is None else " " + chk(e, w, 1e-4)149            print(f"{name} L={L:2d}: sum {S:16.4f} exponent {e:.6f}{tag}")150151# SANDWICH, FLOAT152153def sandwich(name, q, D, F, N, m, want_lo, want_hi, chunk=2048):154    t0 = time.time()155    f, k = mask(F, D)156    lip = lip_of(F)157    n = q ** N158    axes = [a.ravel() for a in np.meshgrid(*([np.arange(n)] * D), indexing="ij")]159    h = 1.0 / (n * m)160    slack = lip * (q ** N - 1) / (q - 1) * (q ** (D * N)) * h / 2161    shifts = [a.ravel() * h for a in np.meshgrid(*([np.arange(m)] * D), indexing="ij")]162    vals = np.empty(len(shifts[0]))163    for s0 in range(0, len(shifts[0]), chunk):164        sh = [s[s0:s0 + chunk] for s in shifts]165        prod = np.ones((len(sh[0]), len(axes[0])))166        for j in range(N):167            prod *= f([((q ** j) * (a[None, :] / n + b[:, None])) % 1.0 for a, b in zip(axes, sh)])168        vals[s0:s0 + chunk] = prod.sum(axis=1)169    lo, hi = float(vals.min()) - slack, float(vals.max()) + slack170    el, eh = math.log(lo) / (N * math.log(q)), math.log(hi) / (N * math.log(q))171    thr = q ** (N * D / 2)172    tl = "" if want_lo is None else " " + chk(el, want_lo, 1e-4)173    th = "" if want_hi is None else " " + chk(eh, want_hi, 1e-4)174    print(f"{name} N={N} m={m}: scan [{vals.min():.6f}, {vals.max():.6f}] slack {slack:.6f} threshold {thr:.4f} -> alpha_1 > {el:.6f}{tl} and < {eh:.6f}{th}  ({time.time() - t0:.1f}s)")175176# INTERVAL KERNEL177178def dn(x):179    return np.nextafter(x, -np.inf)180181def up(x):182    return np.nextafter(x, np.inf)183184PI_UP = up(math.pi)185TABLES = {}186187def costable(Q):188    if Q not in TABLES:189        iv.prec = 96190        lo = np.empty(Q)191        hi = np.empty(Q)192        two = 2 * iv.pi193        for j in range(Q // 2 + 1):194            c = iv.cos(two * j / Q)195            lo[j] = dn(float(c.a))196            hi[j] = up(float(c.b))197        for j in range(Q // 2 + 1, Q):198            lo[j], hi[j] = lo[Q - j], hi[Q - j]199        TABLES[Q] = (lo, hi)200    return TABLES[Q]201202def folded(F):203    mult = {}204    for v in F:205        for w in F:206            d = tuple(a - b for a, b in zip(v, w))207            mult[d] = mult.get(d, 0) + 1208    zero = tuple(0 for _ in F[0])209    out, seen = [], set()210    for d, c in mult.items():211        if d == zero or d in seen:212            continue213        neg = tuple(-a for a in d)214        seen.add(d)215        seen.add(neg)216        out.append((c + mult.get(neg, 0), d))217    return out218219def lip_up(F):220    s = sum(sum(abs(c) for c in v) for v in F)221    return up(up(2 * PI_UP * s) / len(F))222223def hat_iv(kc, Q, dif, k, side):224    lo, hi = costable(Q)225    tab = hi if side else lo226    s = np.full(kc[0].shape, float(k))227    for c, d in dif:228        idx = kc[0] * d[0]229        for j in range(1, len(d)):230            idx = idx + kc[j] * d[j]231        t = tab[idx % Q]232        s = up(s + up(c * t)) if side else dn(s + dn(c * t))233    r = np.sqrt(np.maximum(s, 0.0))234    if side:235        return np.minimum(up(up(r) / k), 1.0)236    return dn(dn(r) / k)237238def kernel_check(pts=400, Q=1944):239    mp.dps = 40240    rng = np.random.default_rng(20250906)241    bad, wide = 0, 0.0242    for name, (q, D, F) in [("gasket 2D", GASKET), ("carpet 2D", CARPET), ("carpet base 9", CARPET9)]:243        dif, k = folded(F), len(F)244        kc = [rng.integers(0, Q, pts) for _ in range(D)]245        lo = hat_iv(kc, Q, dif, k, False)246        hi = hat_iv(kc, Q, dif, k, True)247        for i in range(pts):248            acc = mpc(0)249            for v in F:250                acc += expjpi(2 * mpf(int(sum(int(v[c]) * int(kc[c][i]) for c in range(D)))) / Q)251            t = abs(acc) / k252            if not (mpf(float(lo[i])) <= t <= mpf(float(hi[i]))):253                bad += 1254            wide = max(wide, float(hi[i] - lo[i]))255    print(f"interval kernel against {mp.dps}-digit truth on {3 * pts} arguments mod {Q}: outside enclosure {bad} {chk(bad, 0, 0)}, widest enclosure {wide:.3e}")256257# SANDWICH, CERTIFIED258259def sandwich_iv(name, q, D, F, N, m, chunk=2048):260    t0 = time.time()261    Q = q ** N * m262    costable(Q)263    dif, k = folded(F), len(F)264    n = q ** N265    grid = [a.ravel() * m for a in np.meshgrid(*([np.arange(n)] * D), indexing="ij")]266    shifts = [a.ravel() for a in np.meshgrid(*([np.arange(m)] * D), indexing="ij")]267    slack = up(up(up(lip_up(F) * ((q ** N - 1) // (q - 1))) * q ** (D * N)) / (2 * q ** N * m))268    best_lo, best_hi = np.inf, 0.0269    for s0 in range(0, len(shifts[0]), chunk):270        sh = [s[s0:s0 + chunk] for s in shifts]271        pl = np.ones((len(sh[0]), len(grid[0])))272        ph = np.ones((len(sh[0]), len(grid[0])))273        for j in range(N):274            kc = [((q ** j) * (g[None, :] + b[:, None])) % Q for g, b in zip(grid, sh)]275            pl = dn(pl * hat_iv(kc, Q, dif, k, False))276            ph = up(ph * hat_iv(kc, Q, dif, k, True))277        al = np.zeros(len(sh[0]))278        ah = np.zeros(len(sh[0]))279        for c in range(pl.shape[1]):280            al = dn(al + pl[:, c])281            ah = up(ah + ph[:, c])282        best_lo = min(best_lo, float(al.min()))283        best_hi = max(best_hi, float(ah.max()))284    cl, ch = dn(best_lo - slack), up(best_hi + slack)285    thr = q ** (N * D / 2)286    el = math.floor(math.log(cl) / (N * math.log(q)) * 1e4) / 1e4287    eh = math.ceil(math.log(ch) / (N * math.log(q)) * 1e4) / 1e4288    assert cl > q ** (N * el) and ch < q ** (N * eh)289    print(f"{name} N={N} m={m} certified: min > {cl:.6f}, max < {ch:.6f}, slack {slack:.6f}, threshold q^(ND/2) = {thr:.4f} -> alpha_1 > {el:.4f} and < {eh:.4f}  [{'KILL' if cl > thr else 'no kill'}]  ({time.time() - t0:.1f}s)")290    return el, eh291292# WINDOW STATES293294def window_state(q, D, nd):295    S = q ** D296    dig = np.array(list(itertools.product(range(q), repeat=D)))297    W = S ** nd298    idx = np.arange(W)299    bi = []300    for c in range(D):301        acc = np.zeros(W, dtype=np.int64)302        for i in range(nd):303            acc = acc + dig[:, c][(idx // S ** i) % S] * (q ** (nd - 1 - i))304        bi.append(acc)305    return S, W, idx, bi306307def window_G(q, D, F, nd, m):308    f, k = mask(F, D)309    S, W, idx, bi = window_state(q, D, nd)310    box = 2.0 / q ** nd311    h = box / m312    sub = [a.ravel() for a in np.meshgrid(*([np.arange(m) * h + h / 2] * D), indexing="ij")]313    G = np.empty(W)314    B = max(1, 2_000_000 // len(sub[0]))315    for s0 in range(0, W, B):316        t = [(b[s0:s0 + B, None] / q ** nd + s[None, :]) % 1.0 for b, s in zip(bi, sub)]317        G[s0:s0 + B] = f(t).max(axis=1)318    return np.minimum(G + lip_of(F) * h / 2, 1.0), S, W319320def window_G_iv(q, D, F, nd, m):321    Q = q ** nd * m322    costable(Q)323    dif, k = folded(F), len(F)324    S, W, idx, bi = window_state(q, D, nd)325    sub = [a.ravel() for a in np.meshgrid(*([2 * np.arange(m) + 1] * D), indexing="ij")]326    G = np.empty(W)327    B = max(1, 2_000_000 // len(sub[0]))328    for s0 in range(0, W, B):329        kc = [(b[s0:s0 + B, None] * m + s[None, :]) % Q for b, s in zip(bi, sub)]330        G[s0:s0 + B] = hat_iv(kc, Q, dif, k, True).max(axis=1)331    return np.minimum(up(G + up(lip_up(F) / (q ** nd * m))), 1.0), S, W332333def perron(G, S, nd, t=1.0, iters=500):334    idx = np.arange(len(G))335    succ = [(idx % S ** (nd - 1)) * S + s for s in range(S)]336    Gt = G if t == 1.0 else G ** t337    w = np.ones(len(G))338    lam = 0.0339    for it in range(iters):340        nw = sum(Gt[u] * w[u] for u in succ)341        nl = float(nw.max())342        w = nw / nl343        if abs(nl - lam) < 1e-13:344            lam = nl345            break346        lam = nl347    return lam, np.maximum(w, 1e-12), succ, Gt348349def cw_bound(Gt, w, succ):350    acc = np.zeros(len(w))351    for u in succ:352        acc = up(acc + up(Gt[u] * w[u]))353    plain = np.zeros(len(w))354    for u in succ:355        plain = up(plain + Gt[u])356    return float(up(acc / w).max()), float(plain.max())357358MAY1 = ("published lambda_(1,4) < 2.24190 at 5 digits and box 10^(-5)", 2.24190, 27 / 77)359MAYT = ("published lambda_(235/154,4) < 1.36854", 1.36854, 59 / 433)360KAR1 = ("published 0.3219 at 4 digits", None, 0.3219)361KART = ("published 0.14355 at 4 digits", None, 0.14355)362363WINDOWS = [364    ("base 10 missing 5", base_missing(10, 5), [(4, 64, 2.245878, MAY1), (5, 64, 2.242123, MAY1)]),365    ("base 10 missing 0", base_missing(10, 0), [(4, 64, 2.045493, None), (5, 64, 2.041488, None)]),366    ("carpet 2D", CARPET, [(3, 24, 3.355167, None), (4, 24, 2.670258, None), (5, 24, 2.441254, None)]),367    ("gasket 2D", GASKET, [(4, 24, 2.627210, None), (6, 24, 2.219782, None), (7, 24, 2.153314, None)]),368    ("gasket base 4", GASKET4, [(4, 256, 1.988435, None), (6, 256, 1.952957, None), (8, 256, 1.950717, None)]),369    ("carpet base 9", CARPET9, [(4, 256, 2.134067, None)]),370    ("base 9 missing 0", base_missing(9, 0), [(4, 64, 2.033782, KAR1)]),371]372373def windows(rows):374    for name, (q, D, F), rs in rows:375        for nd, m, want, ref in rs:376            t0 = time.time()377            G, S, W = window_G(q, D, F, nd, m)378            lam, w, succ, Gt = perron(G, S, nd)379            e = math.log(lam) / math.log(q)380            cal = "" if ref is None else f" calibration {ref[0]}, exponent {ref[2]:.6f}, gap {e - ref[2]:+.6f}"381            print(f"{name} windows {nd} digits {W} states sub-scan {m}^{D}: lambda {lam:.6f} {chk(lam, want, 5e-6)} exponent {e:.6f} threshold D/2 = {D / 2}{cal}  ({time.time() - t0:.1f}s)")382383def certify_windows(rows):384    for name, (q, D, F), nd, m in rows:385        t0 = time.time()386        G, S, W = window_G_iv(q, D, F, nd, m)387        lam, w, succ, Gt = perron(G, S, nd)388        mu, plain = cw_bound(Gt, w, succ)389        e = math.ceil(math.log(mu) / math.log(q) * 1e4) / 1e4390        ep = math.ceil(math.log(plain) / math.log(q) * 1e4) / 1e4391        assert mu < q ** e and plain < q ** ep392        print(f"{name} windows {nd} digits {W} states sub-scan {m}^{D} certified: row sum {plain:.6f} -> alpha_1* < {ep:.4f}, weighted {mu:.6f} -> alpha_1* < {e:.4f}, threshold D/2 = {D / 2}  [{'PASS' if mu < q ** (D / 2) else 'no pass'}]  ({time.time() - t0:.1f}s)")393394# MOMENTS395396def crit(b):397    return (1 + math.log(b - 1) / math.log(b) / 2) / 5398399MOMENTS = [400    ("base 10 missing 5", base_missing(10, 5), 4, 64, [(1.0, 0.3514, MAY1), (1.5, 0.1447, None), (235 / 154, 0.1370, MAYT), (1.6, 0.1170, None), (1.7, 0.0937, None)]),401    ("base 9 missing 0", base_missing(9, 0), 4, 64, [(1.5, 0.1446, KART), (235 / 154, 0.1380, None), (1.6, 0.1204, None), (1.7, 0.0995, None)]),402    ("carpet base 9", CARPET9, 4, 64, [(1.0, 0.3451, None), (1.5, 0.1544, None), (235 / 154, 0.1470, None), (1.6, 0.1275, None), (1.7, 0.1042, None)]),403    ("carpet base 9", CARPET9, 5, 64, [(1.0, 0.3437, None), (1.5, 0.1531, None), (235 / 154, 0.1457, None), (1.6, 0.1262, None), (1.7, 0.1031, None), (1.8, 0.0835, None)]),404    ("carpet base 9", CARPET9, 6, 64, [(1.0, 0.3435, None), (1.5, 0.1529, None)]),405    ("gasket base 4", GASKET4, 8, 256, [(1.0, 0.4820, None), (235 / 154, 0.3170, None)]),406]407408def moments(rows):409    for name, (q, D, F), nd, m, ss in rows:410        t0 = time.time()411        G, S, W = window_G(q, D, F, nd, m)412        for s, want, ref in ss:413            lam, w, succ, Gt = perron(G, S, nd, s)414            g = math.log(lam) / math.log(q)415            rhs = crit(q) * (2 - s)416            cal = "" if ref is None else f" calibration {ref[0]}, exponent {ref[2]:.6f}, gap {g - ref[2]:+.6f}"417            print(f"{name} {nd} digits s={s:.6f}: lambda {lam:.6f} g(s) {g:.6f} {chk(g, want, 1e-4)} criterion {rhs:.6f} margin {rhs - g:+.6f} {'PASS' if g < rhs else 'FAIL'}{cal}")418        print(f"  ({time.time() - t0:.1f}s)")419420# LEMMA A PRIME421422def m_of(q, d):423    m, p = 1, q424    while 2 * p <= d:425        p *= q426        m += 1427    return m428429def order(q, d):430    r, x = 1, q % d431    while x != 1:432        x = x * q % d433        r += 1434    return r435436def lemma_a(name, q, D, F, dmax, want):437    f, k = mask(F, D)438    c = 1 - (2 / k) * (1 - math.cos(math.pi / (2 * q)))439    worst = []440    for d in range(3, dmax + 1):441        if math.gcd(d, q) != 1:442            continue443        r, md = order(q, d), m_of(q, d)444        axes = np.meshgrid(*([np.arange(d)] * D), indexing="ij")445        prod = np.ones(axes[0].shape)446        for j in range(r):447            prod *= f([((q ** j % d) * a % d) / d for a in axes])448        prod[(0,) * D] = 0.0449        rate = float(prod.max()) ** (1 / r)450        assert rate <= c ** (1 / md) + 1e-12451        worst.append((rate, d, r, md))452    worst.sort(reverse=True)453    print(f"{name}: c(q,k) = {c:.6f}, moduli 3..{dmax} coprime to q, every per-digit rate at most c^(1/m_d), zero violations")454    for rate, d, r, md in worst[:4]:455        w = want.get(d)456        tag = "" if w is None else " " + chk(rate, w, 5e-7)457        print(f"  d={d:4d} ord={r:3d} m_d={md:2d} rate {rate:.6f}{tag} Lemma A' rate c^(1/m_d) {c ** (1 / md):.6f} Lemma A rate c^(1/ord) {c ** (1 / r):.6f}")458459# COROLLARIES460461def counts_mod(q, D, F, d, n):462    c = np.zeros((d,) * D, dtype=np.int64)463    c[(0,) * D] = 1464    ax = tuple(range(D))465    for j in range(n):466        w = pow(q, j, d)467        nx = np.zeros((d,) * D, dtype=np.int64)468        for v in F:469            nx += np.roll(c, [(w * a) % d for a in v], axis=ax)470        c = nx471    return c472473def corollary(name, q, D, F, dmax, nmax, cells):474    k = len(F)475    assert k ** nmax < 2 ** 62476    c = 1 - (2 / k) * (1 - math.cos(math.pi / (2 * q)))477    worst, rows = (0.0, None), []478    for d in range(3, dmax + 1):479        if math.gcd(d, q) != 1:480            continue481        r, md = order(q, d), m_of(q, d)482        for n in range(2, nmax + 1):483            C = counts_mod(q, D, F, d, n)484            dev = max(abs(Fraction(int(x), k ** n) - Fraction(1, d ** D)) for x in C.ravel())485            bp, ba = c ** (n // md), c ** (n // r)486            assert float(dev) <= bp487            ratio = float(dev) / bp488            if ratio > worst[0]:489                worst = (ratio, (d, n, r, md, float(dev), bp, ba))490            if (d, n) in cells:491                rows.append((d, n, r, md, float(dev), bp, ba, cells[(d, n)]))492    d, n, r, md, dev, bp, ba = worst[1]493    print(f"{name} equidistribution band: moduli 3..{dmax} coprime to q, levels 2..{nmax}, every exact deviation at most c^floor(n/m_d), zero violations")494    print(f"  worst ratio to the window bound {worst[0]:.6f} at d={d} n={n} (ord {r}, m_d {md}): deviation {math.ceil(dev * 1e7) / 1e7:.7f} against A' {math.ceil(bp * 1e6) / 1e6:.6f} and A {math.ceil(ba * 1e6) / 1e6:.6f}")495    for d, n, r, md, dev, bp, ba, want in rows:496        tag = "" if want is None else " " + chk(math.ceil(dev * 1e7) / 1e7, want, 5e-9)497        print(f"  d={d:3d} n={n:3d} ord={r:3d} m_d={md:2d} deviation {math.ceil(dev * 1e7) / 1e7:.7f}{tag} A' bound {math.ceil(bp * 1e6) / 1e6:.6f} A bound {math.ceil(ba * 1e6) / 1e6:.6f}")498499def base_peel(name, q, D, F, ms, nmax):500    k = len(F)501    c = 1 - (2 / k) * (1 - math.cos(math.pi / (2 * q)))502    e = 1503    for p in range(2, q + 1):504        if q % p == 0 and all(p % r for r in range(2, p)):505            e *= p506    ke = sum(1 for v in F if all(a % e == 0 for a in v))507    worst = (0.0, None)508    for m in ms:509        assert math.gcd(m, q) == 1510        mm = m_of(q, m)511        for n in range(2, nmax + 1):512            T = int(counts_mod(q, D, F, e * m, n)[(0,) * D])513            assert T == int(counts_mod(q, D, F, m, n - 1)[(0,) * D])514            err = abs(Fraction(T, k ** n) - Fraction(ke, k) * Fraction(1, m ** D))515            bp = c ** ((n - 1) // mm)516            assert float(err) <= bp517            ratio = float(err) / bp518            if ratio > worst[0]:519                worst = (ratio, (m, n, mm, float(err), bp, c ** ((n - 1) // order(q, m))))520    m, n, mm, err, bp, ba = worst[1]521    print(f"{name} base peel band: rad(q) = {e}, k_e = {ke}, coprime parts {list(ms)}, levels 2..{nmax}, T_(e m)(n) = T_m(n-1) exact and every error at most c^floor((n-1)/m_m), zero violations")522    print(f"  worst ratio {worst[0]:.6f} at m={m} n={n} (m_m {mm}): error {math.ceil(err * 1e7) / 1e7:.7f} against A' {math.ceil(bp * 1e6) / 1e6:.6f} and A {math.ceil(ba * 1e6) / 1e6:.6f}")523524# THE ORDER BAND525526def hat_missing(b, a0, th):527    th = th - np.floor(th)528    s = np.sin(np.pi * th)529    safe = np.where(s == 0.0, 1.0, s)530    K = np.where(s == 0.0, float(b), np.sin(b * np.pi * th) / safe)531    ph = np.cos(2 * np.pi * th * (a0 - (b - 1) / 2))532    return np.sqrt(np.maximum(K * K - 2 * K * ph + 1.0, 0.0)) / (b - 1)533534def hat_missing_check(b, a0, pts=2000):535    f, k = mask(base_missing(b, a0)[2], 1)536    rng = np.random.default_rng(20250906)537    th = rng.random(pts)538    w = float(np.abs(hat_missing(b, a0, th) - f([th])).max())539    print(f"base {b} missing {a0} closed form |sin(b pi th)/sin(pi th) - e((a0 - (b-1)/2) th)|/(b-1) against the character sum on {pts} arguments: worst gap {w:.3e} {chk(w, 0.0, 1e-11)}")540541def order_band(b, a0, N, m, ss, chunk=500):542    lip = (2 * math.pi / (b - 1)) * sum(a for a in range(b) if a != a0)543    n = b ** N544    i = np.arange(n)545    h = 1.0 / (n * m)546    best = {s: np.inf for s in ss}547    for s0 in range(0, m, chunk):548        xg = (np.arange(s0, min(s0 + chunk, m)) + 0.5) * h549        P = np.ones((len(xg), n))550        for j in range(N):551            P *= np.maximum(hat_missing(b, a0, (b ** j) * (xg[:, None] + i[None, :] / n)) - (b ** j) * lip * h / 2, 0.0)552        lg = np.log(np.maximum(P, 1e-300))553        for s in ss:554            best[s] = min(best[s], float(np.exp(s * lg).sum(axis=1).min()))555    return best556557def criterion_band(name, b, a0, N, m, want):558    t0 = time.time()559    ss = sorted(set([round(1.0 + 0.01 * i, 4) for i in range(30)] + [round(1.3 + 0.002 * i, 4) for i in range(151)] + [round(1.6 + 0.01 * i, 4) for i in range(21)] + [1.85, 1.9, 2.0]))560    best = order_band(b, a0, N, m, ss)561    par = (b / (b - 1)) ** N562    print(f"{name} shift sandwich at N={N}, m={m}: Parseval gives Sigma_N^(2)(x) = (b/(b-1))^N = {par:.6f} for every x, scan {best[2.0]:.6f}, slack {par - best[2.0]:.6f}, scan <= anchor {chk_le(best[2.0], par, 'parseval anchor')} ({time.time() - t0:.1f}s)")563    best[2.0] = par564    lo = {s: math.floor(math.log(best[s]) / (N * math.log(b)) * 1e6) / 1e6 for s in ss}565    rhs = {s: math.ceil(crit(b) * (2 - s) * 1e6) / 1e6 for s in ss}566    for s in ss:567        if s in want:568            print(f"  s={s:.4f} min Sigma {best[s]:.6f} -> m_s > {lo[s]:.6f} {chk(lo[s], want[s], 5e-7)} criterion {rhs[s]:.6f} margin {lo[s] - rhs[s]:+.6f}")569    for left in (1.5, 1.0):570        cur, cover = left, []571        while cur < 2.0:572            nxt = max((s for s in ss if s > cur and lo[s] > rhs[cur]), default=None)573            if nxt is None:574                break575            cover.append((cur, nxt))576            cur = nxt577        ok = bool(cover) and cover[-1][1] >= 2.0578        a, z = min(cover, key=lambda p: lo[p[1]] - rhs[p[0]])579        print(f"  monotone cover of [{left}, 2) in {len(cover)} intervals, Sigma_N^(s) falling in s and the criterion falling in s: {'COMPLETE' if ok else 'INCOMPLETE'} {chk(1.0 if ok else 0.0, 1.0, 0.0)}")580        print(f"    tightest cell s={a:.4f} to {z:.4f}: m_s > {lo[z]:.6f} against the criterion {rhs[a]:.6f} at the left end, margin {lo[z] - rhs[a]:+.6f}")581    print(f"  the criterion asks m_s < {crit(b):.6f} (2 - s) for some s in [3/2, 2) and no order in [1, 2) delivers it")582583def order_rows(name, b, a0, N, m, ss, want, chunk):584    t0 = time.time()585    best = order_band(b, a0, N, m, ss, chunk)586    for s in ss:587        lo = math.floor(math.log(best[s]) / (N * math.log(b)) * 1e6) / 1e6588        rhs = math.ceil(crit(b) * (2 - s) * 1e6) / 1e6589        tag = "" if s not in want else " " + chk(lo, want[s], 5e-7)590        print(f"{name} shift sandwich at N={N}, m={m}, s={s:.6f}: min Sigma {best[s]:.6f} -> m_s > {lo:.6f}{tag} criterion {rhs:.6f} margin {lo - rhs:+.6f}  ({time.time() - t0:.1f}s)")591592def two_dim_moments(rows, ss):593    for name, (q, D, F), nd, m, want in rows:594        G, S, W = window_G(q, D, F, nd, m)595        for s in ss:596            lam, w, succ, Gt = perron(G, S, nd, s)597            e = math.log(lam) / (D * math.log(q))598            tag = "" if s not in want else " " + chk(math.ceil(e * 1e6) / 1e6, want[s], 5e-8)599            print(f"{name} {nd} digit-vectors s={s:.6f}: lambda {lam:.6f} exponent in X = q^(D M) units {math.ceil(e * 1e6) / 1e6:.6f}{tag} criterion {math.ceil(crit(q ** D) * (2 - s) * 1e6) / 1e6:.6f}")600601# THE LEAST BASE602603TWIDDLE = {}604605def imul(al, ah, bl, bh):606    p1, p2, p3, p4 = al * bl, al * bh, ah * bl, ah * bh607    return dn(np.minimum(np.minimum(p1, p2), np.minimum(p3, p4))), up(np.maximum(np.maximum(p1, p2), np.maximum(p3, p4)))608609def isq(lo, hi):610    a, c = lo * lo, hi * hi611    return dn(np.where(lo >= 0.0, a, np.where(hi <= 0.0, c, 0.0))), up(np.maximum(a, c))612613def twiddle(Q):614    if Q in TWIDDLE:615        return TWIDDLE[Q]616    B = math.isqrt(Q) + 1617    iv.prec = 96618    def tab(n, step):619        cl, ch, sl, sh = (np.empty(n) for _ in range(4))620        two = 2 * iv.pi621        for j in range(n):622            a = two * ((j * step) % Q) / Q623            c, s = iv.cos(a), iv.sin(a)624            cl[j], ch[j] = dn(float(c.a)), up(float(c.b))625            sl[j], sh[j] = dn(float(s.a)), up(float(s.b))626        return cl, ch, sl, sh627    TWIDDLE[Q] = (B, tab(B, 1), tab(Q // B + 2, B))628    return TWIDDLE[Q]629630def cos_iv(TW, idx):631    B, small, big = TW632    i1 = idx // B633    i0 = idx - i1 * B634    pl, ph = imul(big[0][i1], big[1][i1], small[0][i0], small[1][i0])635    ql, qh = imul(big[2][i1], big[3][i1], small[2][i0], small[3][i0])636    return dn(pl - qh), up(ph - ql)637638def hat_missing_iv(b, a0, j, P, Q, TW):639    jm = j % P640    s1l, s1h = cos_iv(TW, (2 * jm - P) % Q)641    sbl, sbh = cos_iv(TW, (2 * b * jm - P) % Q)642    phl, phh = cos_iv(TW, ((4 * a0 - 2 * b + 2) * jm) % Q)643    if s1l.min() <= 0.0:644        raise SystemExit("denominator enclosure touches zero")645    kl, kh = imul(sbl, sbh, dn(1.0 / s1h), up(1.0 / s1l))646    k2l, k2h = isq(kl, kh)647    ml, mh = imul(kl, kh, phl, phh)648    el = dn(dn(k2l - up(2.0 * mh)) + 1.0)649    eh = up(up(k2h - dn(2.0 * ml)) + 1.0)650    return dn(dn(np.sqrt(np.maximum(el, 0.0))) / (b - 1)), up(up(np.sqrt(np.maximum(eh, 0.0))) / (b - 1))651652def base_windows(b, a0, nd, m, wide=1):653    W = b ** nd654    P = 2 * m * W655    Q = 4 * P656    TW = twiddle(Q)657    lip = up(up(2 * PI_UP * sum(a for a in range(b) if a != a0)) / (b - 1))658    slack = up(up(lip * wide) / (2 * W * m))659    off = (wide * (2 * np.arange(m) + 1)).astype(np.int64)660    Ghi, Glo = np.empty(W), np.empty(W)661    CH = max(1, 400_000 // m)662    for s0 in range(0, W, CH):663        w = np.arange(s0, min(s0 + CH, W), dtype=np.int64)664        j = 2 * m * w[:, None] + off[None, :]665        vl, vh = hat_missing_iv(b, a0, j, P, Q, TW)666        Ghi[s0:s0 + CH] = np.minimum(up(vh.max(axis=1) + slack), 1.0)667        Glo[s0:s0 + CH] = np.maximum(dn(vl.min(axis=1) - slack), 0.0)668    return Ghi, Glo, slack669670def perron_red(G, b, nd, iters=6000, floor_it=300, streak_need=50, tol=1e-13):671    S = b ** (nd - 1)672    U = b ** (nd - 2)673    G3 = G.reshape(b, U, b)674    y = np.ones(S)675    z = np.empty((b, U))676    lam = 0.0677    streak = 0678    for it in range(iters):679        Y = y.reshape(U, b)680        for c in range(b):681            z[c] = (G3[c] * Y).sum(axis=1)682        nl = float(z.max())683        if nl <= 0.0:684            return 0.0, np.ones(S)685        y = (z / nl).reshape(S)686        streak = streak + 1 if abs(nl - lam) <= tol * nl else 0687        lam = nl688        if streak >= streak_need and it >= floor_it:689            break690    return lam, np.maximum(y, 1e-30)691692def cw_red(G, y, b, nd, side):693    U = b ** (nd - 2)694    G3 = G.reshape(b, U, b)695    Y = y.reshape(U, b)696    acc = np.empty((b, U))697    for c1 in range(b):698        r = up(G3[c1] * Y) if side else np.maximum(dn(G3[c1] * Y), 0.0)699        a = np.zeros(U)700        for c0 in range(b):701            a = up(a + r[:, c0]) if side else np.maximum(dn(a + r[:, c0]), 0.0)702        acc[c1] = a703    ratio = acc.reshape(b * U) / y704    return float(up(ratio).max()) if side else float(dn(ratio).min())705706def alpha_band(b, a0, nd, m, wide=1):707    Ghi, Glo, slack = base_windows(b, a0, nd, m, wide)708    lh, yh = perron_red(Ghi, b, nd)709    ll, yl = perron_red(Glo, b, nd)710    mu_hi = cw_red(Ghi, yh, b, nd, True)711    mu_lo = cw_red(Glo, yl, b, nd, False)712    eh = math.ceil(math.log(mu_hi) / math.log(b) * 1e7) / 1e7713    el = None if mu_lo <= 0.0 else math.floor(math.log(mu_lo) / math.log(b) * 1e7) / 1e7714    if not (mu_hi < b ** eh and (el is None or mu_lo > b ** el)):715        raise SystemExit(f"rounding unsafe at base {b} missing {a0}")716    return el, eh, mu_lo, mu_hi, slack717718def shift_check(b, a0, N, m):719    t0 = time.time()720    lip = (2 * math.pi / (b - 1)) * sum(a for a in range(b) if a != a0)721    n = b ** N722    i = np.arange(n)723    h = 1.0 / (n * m)724    lo, hi = np.inf, 0.0725    for s0 in range(0, m, 200):726        xg = (np.arange(s0, min(s0 + 200, m)) + 0.5) * h727        Pl = np.ones((len(xg), n))728        Ph = np.ones((len(xg), n))729        for j in range(N):730            v = hat_missing(b, a0, (b ** j) * (xg[:, None] + i[None, :] / n))731            Pl *= np.maximum(v - (b ** j) * lip * h / 2, 0.0)732            Ph *= np.minimum(v + (b ** j) * lip * h / 2, 1.0)733        lo = min(lo, float(Pl.sum(axis=1).min()))734        hi = max(hi, float(Ph.sum(axis=1).max()))735    el = math.floor(math.log(lo) / (N * math.log(b)) * 1e6) / 1e6736    eh = math.ceil(math.log(hi) / (N * math.log(b)) * 1e6) / 1e6737    print(f"base {b} missing {a0} shift sandwich N={N} m={m}, the grid machine not the window machine: min Sigma {lo:.6f} max Sigma {hi:.6f} -> alpha_1 in [{el:.6f}, {eh:.6f}]  ({time.time() - t0:.1f}s)")738    return el, eh739740def base_cert(b, a0, nd, m, wide=1, tag=""):741    t0 = time.time()742    el, eh, mu_lo, mu_hi, slack = alpha_band(b, a0, nd, m, wide)743    lo = "no positive certificate at this window length" if el is None else f"{el:.7f}"744    verdict = "CLEARS 1/4" if eh < 0.25 else ("FAILS 1/4" if el is not None and el >= 0.25 else "undecided at 1/4")745    print(f"base {b} missing {a0}{tag} windows {nd} digits sub-scan {m} box {wide}/{b}^{nd} certified: lambda in [{mu_lo:.6f}, {mu_hi:.6f}] slack {slack:.3e} -> alpha_1 > {lo} and < {eh:.7f}  [{verdict}]  ({time.time() - t0:.1f}s)")746    return el, eh747748def base_family(b, nd, m, side, want=None):749    t0 = time.time()750    rows = []751    for a0 in range((b + 1) // 2):752        el, eh, mu_lo, mu_hi, slack = alpha_band(b, a0, nd, m)753        rows.append((a0, el, eh))754    worst = max(rows, key=lambda r: r[2])755    best = min(rows, key=lambda r: r[2])756    band = " ".join(f"{a0}:[{-1.0 if el is None else el:.7f},{eh:.7f}]" for a0, el, eh in rows)757    if side:758        ok = all(eh < 0.25 for _, _, eh in rows)759        head = f"every one-missing-digit set in base {b} clears alpha_1 < 1/4" if ok else f"base {b} does not clear at every digit"760    else:761        ok = all(el is not None and el >= 0.25 for _, el, _ in rows)762        head = f"no one-missing-digit set in base {b} clears alpha_1 < 1/4" if ok else f"base {b} clears at some digit"763    print(f"{head}: {nd} digits sub-scan {m}, {len(rows)} digits up to the symmetry a0 <-> {b - 1} - a0, worst a0 = {worst[0]} at [{-1.0 if worst[1] is None else worst[1]:.7f}, {worst[2]:.7f}], best a0 = {best[0]} at [{-1.0 if best[1] is None else best[1]:.7f}, {best[2]:.7f}] {chk(1.0 if ok else 0.0, 1.0, 0.0)}  ({time.time() - t0:.1f}s)")764    print(f"  base {b} certified alpha_1 band by missing digit: {band}")765    return rows766767def family_window(b, a0, nd, m):768    Ghi, Glo, slack = base_windows(b, a0, nd, m)769    lh, yh = perron_red(Ghi, b, nd)770    mu_hi = cw_red(Ghi, yh, b, nd, True)771    eh = math.ceil(math.log(mu_hi) / math.log(b) * 1e7) / 1e7772    if not mu_hi < b ** eh:773        raise SystemExit(f"rounding unsafe at base {b} missing {a0}")774    return eh, mu_hi775776def family_close(lo, hi, m, nds):777    t0 = time.time()778    rows, bad = [], []779    for b in range(lo, hi + 1):780        for nd in nds:781            worst, arg = -1.0, -1782            for a0 in range((b + 1) // 2):783                eh, mu = family_window(b, a0, nd, m)784                if eh > worst:785                    worst, arg = eh, a0786            if worst < 0.25:787                break788        rows.append((b, nd, arg, worst))789        if worst >= 0.25:790            bad.append((b, arg, worst))791    ok = not bad792    head = f"every one-missing-digit set of every base {lo} <= q <= {hi} clears alpha_1 < 1/4" if ok else f"{len(bad)} base(s) in {lo} <= q <= {hi} do not clear at every digit"793    ceiling = max(rows, key=lambda r: r[3])794    print(f"{head}: {hi - lo + 1} bases, {sum((b + 1) // 2 for b in range(lo, hi + 1))} distinct sets, sub-scan {m}, shortest window in {nds} that clears; the band ceiling is q = {ceiling[0]} missing {ceiling[2]} at alpha_1 < {ceiling[3]:.7f} on {ceiling[1]} window digits {chk(1.0 if ok else 0.0, 1.0, 0.0)}  ({time.time() - t0:.1f}s)")795    for b, nd, arg, worst in rows:796        flag = "" if worst < 0.25 else "  FAILS 1/4"797        print(f"  q={b:3d} {(b + 1) // 2:2d} distinct sets {nd} window digits: worst a0 = {arg:2d} at alpha_1 < {worst:.7f}{flag}")798    return rows, bad799800def base_ladder(bases, nd, m, a0s):801    t0 = time.time()802    out = []803    for b in bases:804        a0 = a0s(b)805        el, eh, mu_lo, mu_hi, slack = alpha_band(b, a0, nd, m)806        out.append((b, a0, el, eh))807    print(f"certified alpha_1 ladder at {nd} digits sub-scan {m}, one missing digit per base: " + " ".join(f"q={b} a0={a0} [{-1.0 if el is None else el:.7f},{eh:.7f}]" for b, a0, el, eh in out) + f"  ({time.time() - t0:.1f}s)")808    return out809810# THE PAIR FAMILY811812def pair_missing(b, a, c):813    return (b, 1, [(d,) for d in range(b) if d != a and d != c])814815def hat_pair(b, a, c, th):816    th = th - np.floor(th)817    s = np.sin(np.pi * th)818    safe = np.where(s == 0.0, 1.0, s)819    K = np.where(s == 0.0, float(b), np.sin(b * np.pi * th) / safe)820    p = a - (b - 1) / 2821    r = c - (b - 1) / 2822    v = K * K + 2.0 + 2.0 * np.cos(2 * np.pi * (a - c) * th) - 2.0 * K * (np.cos(2 * np.pi * p * th) + np.cos(2 * np.pi * r * th))823    return np.sqrt(np.maximum(v, 0.0)) / (b - 2)824825def hat_pair_check(b, a, c, pts=2000):826    f, k = mask(pair_missing(b, a, c)[2], 1)827    rng = np.random.default_rng(20250907)828    th = rng.random(pts)829    w = float(np.abs(hat_pair(b, a, c, th) - f([th])).max())830    print(f"base {b} missing {{{a},{c}}} closed form |K - e((a - (b-1)/2) th) - e((c - (b-1)/2) th)|/(b-2) against the character sum on {pts} arguments: worst gap {w:.3e} {chk(w, 0.0, 1e-11)}")831832def pair_count(b):833    return ((b - 2) * (b - 3) // 2 + (b - 2) // 2) // 2 + b // 2834835def pair_sets(b):836    out = [(0, c) for c in range(1, b // 2 + 1)]837    seen = set()838    for a in range(1, b - 1):839        for c in range(a + 1, b - 1):840            if (b - 1 - c, b - 1 - a) in seen:841                continue842            seen.add((a, c))843            out.append((a, c))844    return out845846def hat_pair_iv(b, a, c, j, P, Q, TW):847    jm = j % P848    s1l, s1h = cos_iv(TW, (2 * jm - P) % Q)849    sbl, sbh = cos_iv(TW, (2 * b * jm - P) % Q)850    pal, pah = cos_iv(TW, ((4 * a - 2 * b + 2) * jm) % Q)851    pcl, pch = cos_iv(TW, ((4 * c - 2 * b + 2) * jm) % Q)852    dfl, dfh = cos_iv(TW, (4 * (a - c) * jm) % Q)853    if s1l.min() <= 0.0:854        raise SystemExit("denominator enclosure touches zero")855    kl, kh = imul(sbl, sbh, dn(1.0 / s1h), up(1.0 / s1l))856    k2l, k2h = isq(kl, kh)857    ml, mh = imul(kl, kh, dn(pal + pcl), up(pah + pch))858    el = dn(dn(dn(k2l + 2.0) + dn(2.0 * dfl)) - up(2.0 * mh))859    eh = up(up(up(k2h + 2.0) + up(2.0 * dfh)) - dn(2.0 * ml))860    return dn(dn(np.sqrt(np.maximum(el, 0.0))) / (b - 2)), up(up(np.sqrt(np.maximum(eh, 0.0))) / (b - 2))861862def pair_windows(b, a, c, nd, m, wide=1):863    W = b ** nd864    P = 2 * m * W865    Q = 4 * P866    TW = twiddle(Q)867    lip = up(up(2 * PI_UP * (b * (b - 1) // 2 - a - c)) / (b - 2))868    slack = up(up(lip * wide) / (2 * W * m))869    off = (wide * (2 * np.arange(m) + 1)).astype(np.int64)870    Ghi, Glo = np.empty(W), np.empty(W)871    CH = max(1, 400_000 // m)872    for s0 in range(0, W, CH):873        w = np.arange(s0, min(s0 + CH, W), dtype=np.int64)874        j = 2 * m * w[:, None] + off[None, :]875        vl, vh = hat_pair_iv(b, a, c, j, P, Q, TW)876        Ghi[s0:s0 + CH] = np.minimum(up(vh.max(axis=1) + slack), 1.0)877        Glo[s0:s0 + CH] = np.maximum(dn(vl.min(axis=1) - slack), 0.0)878    return Ghi, Glo, slack879880def pair_band(b, a, c, nd, m, wide=1, side=None):881    Ghi, Glo, slack = pair_windows(b, a, c, nd, m, wide)882    eh = el = None883    mu_hi = mu_lo = 0.0884    if side is None or side:885        lh, yh = perron_red(Ghi, b, nd)886        mu_hi = cw_red(Ghi, yh, b, nd, True)887        eh = math.ceil(math.log(mu_hi) / math.log(b) * 1e7) / 1e7888        if not mu_hi < b ** eh:889            raise SystemExit(f"rounding unsafe at base {b} missing {a},{c}")890    if side is None or not side:891        ll, yl = perron_red(Glo, b, nd)892        mu_lo = cw_red(Glo, yl, b, nd, False)893        el = None if mu_lo <= 0.0 else math.floor(math.log(mu_lo) / math.log(b) * 1e7) / 1e7894        if el is not None and not mu_lo > b ** el:895            raise SystemExit(f"rounding unsafe at base {b} missing {a},{c}")896    return el, eh, mu_lo, mu_hi, slack897898def pair_cert(b, a, c, nd, m, wide=1, tag="", thr=0.25, name="1/4"):899    t0 = time.time()900    el, eh, mu_lo, mu_hi, slack = pair_band(b, a, c, nd, m, wide)901    lo = "no positive certificate at this window length" if el is None else f"{el:.7f}"902    verdict = f"CLEARS {name}" if eh < thr else (f"FAILS {name}" if el is not None and el >= thr else f"undecided at {name}")903    print(f"base {b} missing {{{a},{c}}}{tag} windows {nd} digits sub-scan {m} box {wide}/{b}^{nd} certified: lambda in [{mu_lo:.6f}, {mu_hi:.6f}] slack {slack:.3e} -> alpha_1 > {lo} and < {eh:.7f}  [{verdict}]  ({time.time() - t0:.1f}s)")904    return el, eh905906def pair_census(bs):907    th = (np.arange(1001) + 0.5) / 1001908    for b in bs:909        fp = {}910        for a in range(b):911            for c in range(a + 1, b):912                fp.setdefault(tuple(np.round(hat_pair(b, a, c, th), 11)), []).append((a, c))913        sets = pair_sets(b)914        reps = {tuple(np.round(hat_pair(b, a, c, th), 11)) for (a, c) in sets}915        print(f"base {b}: {b * (b - 1) // 2} excluded pairs fall into {len(fp)} distinct transforms against (C(q-2,2) + floor((q-2)/2))/2 + floor(q/2) = {pair_count(b)} {chk(float(len(fp)), float(pair_count(b)), 0.0)}, the scan list holds {len(sets)} {chk(float(len(sets)), float(pair_count(b)), 0.0)} and meets every class {chk(float(len(reps)), float(len(fp)), 0.0)}; the edge class {{0,c}} is the one-missing-digit set of the {b - 1}-digit interval and collapses {{0,c}} with {{0,{b}-c}}")916917def pair_shortest(b, nds, m, side, thr=0.25, name="1/4"):918    t0 = time.time()919    todo = pair_sets(b)920    n0 = len(todo)921    done, wnd = [], {}922    for nd in nds:923        nxt = []924        for (a, c) in todo:925            el, eh, mu_lo, mu_hi, slack = pair_band(b, a, c, nd, m, side=side)926            v = eh if side else (-1.0 if el is None else el)927            if (v < thr) if side else (v >= thr):928                done.append((a, c, v))929                wnd[nd] = wnd.get(nd, 0) + 1930            else:931                nxt.append((a, c))932        todo = nxt933        if not todo:934            break935    ok = not todo936    key = (max(done, key=lambda r: r[2]) if side else min(done, key=lambda r: r[2])) if done else (-1, -1, -1.0)937    spread = " ".join(f"{nd}:{wnd[nd]}" for nd in sorted(wnd))938    verb = f"clears alpha_1 < {name} at every pair" if side else f"fails {name} at every pair"939    tail = ""940    if not ok:941        ups = []942        for (a, c) in todo:943            el, eh, mu_lo, mu_hi, slack = pair_band(b, a, c, nds[-1], m, side=not side)944            ups.append(f"{{{a},{c}}} S={a + c - b + 1} D={c - a} alpha_1 {'>' if side else '<'} {(-1.0 if el is None else el) if side else eh:.7f}")945        note = "not shown to clear" if side else "no positive lower certificate"946        tail = f"; {len(todo)} UNCLEAR at {nds[-1]} digits, {note}: " + ", ".join(ups)947    print(f"  q={b:3d} {n0:5d} distinct sets {verb}: window digits used {spread}, {'worst' if side else 'closest'} {{{key[0]},{key[1]}}} at alpha_1 {'<' if side else '>'} {key[2]:.7f}{tail} {chk(1.0 if ok else 0.0, 1.0, 0.0)}  ({time.time() - t0:.1f}s)", flush=True)948    return b, ok, key, wnd949950def pair_first(b, nds, m, thr=0.25, name="1/4", want=None):951    t0 = time.time()952    best, arg, wnd = None, None, None953    for (a, c) in pair_sets(b):954        for nd in nds:955            el, eh, mu_lo, mu_hi, slack = pair_band(b, a, c, nd, m, side=True)956            if eh < thr:957                break958        if best is None or eh < best:959            best, arg, wnd = eh, (a, c), nd960    ok = best < thr961    fl = pair_band(b, arg[0], arg[1], wnd, m, side=False)[0]962    floor = "no positive lower certificate" if fl is None else f"{fl:.7f}"963    pin = "" if want is None else " " + chk(1.0 if ok == want[0] else 0.0, 1.0, 0.0) + chk(best, want[1], 5e-7)964    print(f"  q={b:3d} {len(pair_sets(b)):5d} distinct sets, best pair {{{arg[0]},{arg[1]}}} at {wnd} window digits: alpha_1 < {best:.7f} [{'CLEARS' if ok else 'does not clear'} {name}], the infimum matrix of that set brackets it from below at alpha_1 > {floor} {chk(1.0 if fl is None or fl <= best else 0.0, 1.0, 0.0)}{pin}  ({time.time() - t0:.1f}s)", flush=True)965    return b, ok, arg, best966967def pair_some(b, nds, m, thr=0.25, name="1/4"):968    t0 = time.time()969    best, arg, wnd = -1.0, None, None970    for nd in nds:971        for (a, c) in pair_sets(b):972            el, eh, mu_lo, mu_hi, slack = pair_band(b, a, c, nd, m, side=False)973            v = -1.0 if el is None else el974            if v > best:975                best, arg, wnd = v, (a, c), nd976            if v >= thr:977                break978        if best >= thr:979            break980    print(f"  q={b:3d} {len(pair_sets(b)):5d} distinct sets, witness pair {{{arg[0]},{arg[1]}}} at {wnd} window digits: alpha_1 > {best:.7f} [{'FAILS' if best >= thr else 'no witness at these windows'} {name}] {chk(1.0 if best >= thr else 0.0, 1.0, 0.0)}  ({time.time() - t0:.1f}s)", flush=True)981    return b, best >= thr, arg, best982983def pair_ladder(bases, a, c, nd, m):984    t0 = time.time()985    out = []986    for b in bases:987        el, eh, mu_lo, mu_hi, slack = pair_band(b, a, c, nd, m)988        out.append((b, el, eh))989    print(f"certified alpha_1 ladder at {nd} digits sub-scan {m}, the interval class {{{a},{c}}} in each base: " + " ".join(f"q={b} [{-1.0 if el is None else el:.7f},{eh:.7f}]" for b, el, eh in out) + f"  ({time.time() - t0:.1f}s)")990    return out991992# VERBS993994def main():995    verb = sys.argv[1] if len(sys.argv) > 1 else "check"996    t0 = time.time()997    if verb == "check":998        digit_sums()999        gasket_anchor()1000        entropy_bound()1001        grids(GRID[:2])1002        lemma_a("gasket 2D", *GASKET, 129, {127: 0.808166, 129: 0.809637})1003        windows(WINDOWS[:1])1004        kernel_check()1005        sandwich_iv("gasket 2D", *GASKET, 2, 256)1006        certify_windows([("carpet 2D", CARPET, 4, 24)])1007        corollary("gasket 2D", *GASKET, 9, 10, {(5, 10): 0.0014741})1008    elif verb == "grid":1009        grids(GRID)1010    elif verb == "sandwich":1011        sandwich("gasket 2D", *GASKET, 2, 256, 1.0105, 1.1029)1012        sandwich("gasket 2D", *GASKET, 3, 256, 1.0126, 1.1022)1013        sandwich("gasket 2D", *GASKET, 4, 300, 1.0096, 1.1046)1014        sandwich("carpet 2D", *CARPET, 2, 256, 0.5828, 0.9507)1015        sandwich("carpet 2D", *CARPET, 3, 400, None, 0.9421)1016        sandwich("gasket base 4", *GASKET4, 4, 20000, None, 0.4864)1017        sandwich("carpet base 9", *CARPET9, 3, 20000, None, 0.3704)1018    elif verb == "certify":1019        gasket_anchor()1020        entropy_bound()1021        kernel_check()1022        sandwich_iv("gasket 2D", *GASKET, 2, 256)1023        sandwich_iv("gasket 2D", *GASKET, 3, 256)1024        certify_windows([("carpet 2D", CARPET, 4, 24), ("carpet 2D", CARPET, 5, 24)])1025    elif verb == "windows":1026        windows(WINDOWS)1027    elif verb == "moments":1028        moments(MOMENTS)1029    elif verb == "lemma":1030        lemma_a("gasket 2D", *GASKET, 301, {257: 0.830915, 255: 0.830253, 129: 0.809637, 127: 0.808166})1031    elif verb == "corollary":1032        corollary("gasket 2D", *GASKET, 15, 12, {(3, 2): 0.2222223, (5, 10): 0.0014741})1033        corollary("carpet 2D", *CARPET, 11, 8, {(4, 2): 0.0625, (5, 8): 0.0002842})1034        base_peel("gasket 2D", *GASKET, [3, 5, 7], 11)1035        base_peel("carpet 2D", *CARPET, [2, 5, 7], 7)1036    elif verb == "criterion":1037        hat_missing_check(9, 4)1038        hat_missing_check(10, 5)1039        criterion_band("carpet base 9", 9, 4, 4, 4000, {1.0: 0.334604, 1.46: 0.159988, 1.5: 0.148588, 1.6: 0.122780, 1.8: 0.081955, 2.0: 0.053605})1040        order_rows("carpet base 9", 9, 4, 5, 3000, [1.5, 235 / 154], {1.5: 0.149397, 235 / 154: 0.142274}, 60)1041        two_dim_moments([("carpet 2D", CARPET, 5, 24, {1.0: 0.406200, 1.5: 0.195631}), ("carpet base 9", CARPET9, 5, 64, {1.0: 0.343674, 1.5: 0.153069})], [1.0, 1.5, 235 / 154, 1.6, 1.7, 1.8])1042    elif verb == "least":1043        hat_missing_check(21, 0)1044        hat_missing_check(34, 16)1045        base_cert(10, 5, 4, 16, 2, " calibration against the study float row lambda 2.245878 and Maynard lambda_(1,4) < 2.24190")1046        base_cert(10, 5, 5, 16, 2, " calibration against the study float row lambda 2.242123 and 27/77 = 0.3506494")1047        base_cert(10, 5, 6, 16, 2, " calibration: 27/77 = 0.3506494 is a finite-window upper bound, not the exponent")1048        base_cert(9, 0, 4, 16, 2, " calibration against the study float row lambda 2.033782 and Karwatowski 0.3219")1049        base_cert(21, 0, 5, 8, 1, " THE LEAST BASE CARRYING ONE SUCH SET")1050        base_family(20, 4, 8, False)1051        base_family(34, 4, 8, True)1052        base_cert(33, 15, 4, 8, 1, " THE WITNESS THAT 34 IS LEAST FOR THE WHOLE FAMILY")1053        base_ladder([10, 14, 18, 20, 21, 22, 26, 30, 33, 34], 4, 8, lambda b: 0)1054        base_ladder([10, 14, 18, 20, 21, 22, 26, 30, 33, 34], 4, 8, lambda b: b // 2)1055        shift_check(21, 0, 4, 400)1056        shift_check(34, 16, 3, 600)1057        grids([("base 21 missing 0", base_missing(21, 0), [None] * 4)])1058    elif verb == "six":1059        base_cert(21, 0, 6, 8, 1, " THE LEAST BASE CARRYING ONE SUCH SET, THE HEADLINE WINDOW")1060        base_cert(10, 5, 7, 16, 2, " calibration: the exponent is strictly under 27/77 = 0.3506494")1061    elif verb == "pairs":1062        hat_pair_check(21, 0, 1)1063        hat_pair_check(32, 7, 19)1064        pair_census([6, 9, 10, 12])1065        pair_cert(32, 0, 1, 4, 8, 1, " THE LEAST BASE CARRYING ONE CERTIFIED SUCH SET")1066        pair_cert(31, 0, 1, 4, 8, 1, " THE WITNESS THAT 32 IS LEAST FOR THE INTERVAL CLASS")1067        pair_cert(20, 3, 11, 5, 8, 1, " THE WITNESS THAT 21 IS LEAST FOR THE WHOLE FAMILY AGAINST 1/3", 1 / 3, "1/3")1068        pair_first(13, [3], 8, 1 / 3, "1/3", (True, 0.3318819))1069        pair_first(12, [5], 8, 1 / 3, "1/3", (False, 0.3371162))1070        pair_ladder([20, 24, 28, 30, 31, 32, 33, 36, 40], 0, 1, 4, 8)1071    elif verb == "pairone":1072        v = sys.argv[2:]1073        thr, name = (1 / 3, "1/3") if v and v[-1] == "third" else (0.25, "1/4")1074        n = [int(x) for x in v if x != "third"]1075        pair_cert(n[0], n[1], n[2], n[3], n[4], n[5] if len(n) > 5 else 1, "", thr, name)1076    elif verb == "pairfail":1077        v = sys.argv[2:]1078        thr, name = (1 / 3, "1/3") if len(v) > 2 and v[2] == "third" else (0.25, "1/4")1079        for b in range(int(v[0]), int(v[1]) + 1):1080            pair_shortest(b, [2, 3, 4, 5], 8, False, thr, name)1081    elif verb == "pairclear":1082        v = sys.argv[2:]1083        thr, name = (1 / 3, "1/3") if len(v) > 2 and v[2] == "third" else (0.25, "1/4")1084        for b in range(int(v[0]), int(v[1]) + 1):1085            pair_shortest(b, [2, 3, 4], 8, True, thr, name)1086    elif verb == "pairsome":1087        v = sys.argv[2:]1088        thr, name = (1 / 3, "1/3") if len(v) > 2 and v[2] == "third" else (0.25, "1/4")1089        for b in range(int(v[0]), int(v[1]) + 1):1090            pair_some(b, [3, 4, 5], 8, thr, name)1091    elif verb == "pairfirst":1092        v = sys.argv[2:]1093        thr, name = (1 / 3, "1/3") if len(v) > 2 and v[2] == "third" else (0.25, "1/4")1094        for b in range(int(v[0]), int(v[1]) + 1):1095            pair_first(b, [2, 3, 4, 5], 8, thr, name)1096    elif verb == "family":1097        family_close(35, 125, 8, [2, 3, 4])1098    else:1099        raise SystemExit("verbs: check grid sandwich certify windows moments lemma corollary criterion least six family pairs pairfail pairclear")1100    close(t0)11011102main()