""" real_support.py -- the carries' real support: every carry of a bounded relaxed path is a linear form in its delta, so its size is at most one L1 norm per d. QUESTION. Over every bounded integer x (t >= 0) with delta_t = x_t - (x_(t+1) + ... + x_(t+d)) in {-1, 0, 1} and first window sum 0 (the carry sequences of pairs of codings, legality dropped), how large can |x_t| and the window sum |x_(t+1) + ... + x_(t+d)| be with integrality dropped too? If the floor of the real sup is C_d = max(1, floor((d - 1)/4)), the triangle wave's carry, the carry bound holds at that d. BY HAND. delta = Delta(F) x, F the shift, Delta(z) = 1 - z - ... - z^d, whose roots are rho_0 = 1/beta inside the unit disc and rho_1, ..., rho_(d-1) outside (their inverses are the d-bonacci polynomial's roots, beta Pisot). 1/Delta = sum_j c_j/(1 - z/rho_j) with c_j = -1/(rho_j Delta'(rho_j)), and its Laurent expansion on |z| = 1 is the bounded inverse h: h_m = sum_(j >= 1) c_j rho_j^(-m) at m >= 0, h_m = -c_0 beta^m at m <= -1, so h_0 = 1 - c_0 (the c_j sum to 1/Delta(0)). On t >= 0, x^p_t = sum_(k >= 0) h_(k - t) delta_k solves every equation (the coefficient of z^(k - t) in Delta(z) h(z) = 1), and is bounded by L(d) = sum_m |h_m|; Delta(F) x = 0 has the bounded solutions c beta^(-t) alone, so every bounded solution is x^p + c beta^(-t), and the first window sum 0 fixes c = -(x^p_0 + ... + x^p_(d-1))/(1 + ... + beta^(-(d-1))). So x_t = g^(t) . delta, g^(t)_k = h_(k - t) + beta^(-t) gam_k, gam = -sum_(s < d) h_(. - s)/sum_(s < d) beta^(-s), and |x_t| <= ||g^(t)||_1; the window sum sigma_(t+1) = x_t - delta_t has the form g^(t) - e_t. H_x(d) and H_s(d) are the sups over t of the two norms. The sups are read at t < T = 81: past T, ||g^(t)||_1 <= L + beta^(-t) ||gam||_1 with ||gam||_1 <= d L and beta^(-T) < 1.84^(-80) < 1e-21, and the window sum's row holds h_0 - 1 at k = t, so its bound is L_s + beta^(-t) d L, L_s = L - |h_0| + |h_0 - 1|. The script prints L <= H_x and L_s <= H_s, so H_x and H_s are the sups over every t, to the float. h two ways: by the roots (mpmath, 40 digits) at d = 3..40, and by the real FFT of 1/Delta on N = 2^lg points at d = 21..100, lg = max(16, ceil(log2(2 d^3)) + 1), the slowest root's modulus being about 1 + 4 pi^2/d^3, so h decays over about d^3 terms. The FFT reads the aliasing, max |h| near N/2, and keeps forward h_m above 1e-18. 1/Delta(z) = (1 - z)/(1 - 2z + z^(d+1)), so h_m = (G_m - G_(m+1))/2, G the Green's function of the walk with steps -1 and +d, each 1/2: L(d) is half G's total variation. PREDICTIONS, CONTROLS AND KILL (fixed before the runs). H1 L(d)/d tends to a constant in [0.25, 0.35] over d = 3..40, and L(d) >= (d - 1)/4 at every d. H2 the route proves the bound (floor H_x, floor H_s <= C_d) at d = 3..8 at most and fails at some d in 9..12. KILL of the route at every d: floor(H_x(d)) > C_d at some d in 3..40. CONTROL, read first: the wave at d = 9, 13, 21 recomputed from g at t < 4(d - 1) equals its x to 1e-9; g's x for a random delta obeys Delta(F) x = delta to 1e-9 on t < 200. H3 (fixed after H's run) the FFT route agrees with the roots to 1e-9 at d = 5, 12, 30, 40 (its control); KILL, floor(H_x) or floor(H_s) > C_d at some d = 21..100. H3b d (L - (d - 1)/4) falls at every step d = 21..100. FINDINGS (copied from the runs' prints). - CONTROL holds: the wave from g within 3.1e-15, 1.2e-14, 1.9e-14 at d = 9, 13, 21; residuals 1.4e-15, 5.6e-15, 1.0e-14 and first window sums below 2.0e-15 at d = 5, 12, 30. - The KILL fires: floor(H_x) > C_d first at d = 6. The route PROVES the bound (floor H_x, floor H_s <= C_d) at d = 3, 4, 9, 10, 13, 14, 15, 17, 18, 19 and every d = 21..100; it leaves d = 5..8, 11, 12, 16, 20. H_x at the open d past 8: 3.038251 (11), 3.243495 (12), 4.113062 (16), 5.030733 (20), against C_d + 1 = 3, 3, 4, 5. - H2 fails the other way: the route passes most d past 12. H1: L(d)/d falls to 0.2469 at d = 40, below 1/4. The tightest d past 20 are d = 0 mod 4, where C_d + 1 - (d - 1)/4 = 1/4; the closest passes are d = 3 (H_s 1.994724 against 2), d = 4 (1.981133) and, past 20, d = 24 (H_x 5.976319, margin 0.0237), then 0.0615 (28), 0.1248 (40), 0.1697 (60), 0.2033 (100) at d = 0 mod 4. - H3b holds: d (L - (d - 1)/4) falls at every step, 5.5016 (21), 5.0014 (40), 4.8169 (60), 4.6717 (100). - H_x = L(d) to 1e-6 at d <= 11 (t past the transient); from d = 12 the sup sits at t = 5..11 to d = 40, above L by 5.0e-4 (12), 4.5e-3 (15, the most), 3.2e-3 (21), 2.0e-4 (40), and moves out with d, t = 26 at d = 97..100, equal to L to 1e-6. L <= H_x and L_s <= H_s at every d = 3..100. - Float: no H within 1e-6 of an integer; the FFT agrees with the roots to 8.9e-13 at d = 5, 12, 30, 40, its aliasing below 3e-17. - Argued, not printed: integrality read one functional at a time pins no carry at the open d: a value C_d + 1 at t forces delta to sign(g^(t)) up to a defect of H - C_d - 1, but delta at the d - 4 places after t (h_1, ..., h_(d-4) small) and more than about four before it costs x_t little and moves each other x_u by whole units, so no one x_u is pinned. RUN RECORD: the bare run 108.1 s, peak commit 247.2 MB. The roots at d = 3..40 first took 6.6 s, 143.5 MB (one matrix of every g^(t) was refused at 541.6 MB at d = 40, so rows are read one at a time); the FFT at d = 21..100 on 2^22 points throughout, with T = 69 + 4d, took 633.5 s, 214.8 MB, the cost per d linear in T and N, hence lg and T as above. """ import math import os import sys os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") import mpmath # noqa: E402 import numpy as np # noqa: E402 mpmath.mp.dps = 40 T = 81 CHECKS = [] def check(name, ok, detail=""): CHECKS.append(bool(ok)) print(f" [{'PASS' if ok else 'FAIL'}] {name}" + (f" -- {detail}" if detail else ""), flush=True) def inverse(d, eps=1e-16): """beta, c_0 and the forward impulse h_m, m = 0..K, by the roots.""" coef = [-1] * d + [1] # Delta, highest power first roots = mpmath.polyroots(coef, maxsteps=400, extraprec=200) dD = lambda z: -sum(k * z ** (k - 1) for k in range(1, d + 1)) c = [-1 / (r * dD(r)) for r in roots] i0 = min(range(d), key=lambda i: abs(roots[i])) beta = float(mpmath.re(1 / roots[i0])) c0 = float(mpmath.re(c[i0])) rest = [(complex(roots[i]), complex(c[i])) for i in range(d) if i != i0] slow = min(abs(r) for r, _ in rest) K = int(math.ceil(math.log(eps) / -math.log(slow))) + d m = np.arange(K) h = np.zeros(K) for r, cj in rest: h += (cj * r ** (-m.astype(float))).real return beta, c0, h def inverse_fft(d): """h by the real FFT of 1/Delta on N = 2^lg points: forward h_m (m >= 0, cut where |h| last exceeds 1e-18), backward h_(-m) (m >= 1, hb[0] unused), and the aliasing, max |h| near N/2.""" lg = max(16, math.ceil(math.log2(2 * d ** 3)) + 1) N = 1 << lg th = -2 * np.pi * np.arange(N // 2 + 1) / N w = np.exp(1j * th) wd = np.exp(1j * th * d) with np.errstate(divide="ignore", invalid="ignore"): D = 1 - w * (1 - wd) / (1 - w) D[0] = 1 - d del w, wd, th h = np.fft.irfft(1 / D, n=N) del D alias = np.abs(h[N // 2 - 1000:N // 2 + 1000]).max() big = np.flatnonzero(np.abs(h[:N // 2]) > 1e-18) hf = h[:big[-1] + 1].copy() hb = np.concatenate([[0.0], h[::-1][:200 * d].copy()]) return hf, hb, alias def norms(d, hf, hb): """L, L_s and, over t < T, ||g^(t)||_1 and ||g^(t) - e_t||_1, from the two halves of h, one row at a time; beta read off hb = -c_0 beta^(-m).""" beta = hb[1] / hb[2] n = len(hf) + T hz = np.concatenate([hf, np.zeros(T)]) k = np.arange(n) def shifted(t): v = np.empty(n) f = k >= t v[f] = hz[k[f] - t] v[~f] = hb[t - k[~f]] return v gam = -sum(shifted(s) for s in range(d)) / sum(beta ** -s for s in range(d)) nx, ns = np.empty(T), np.empty(T) for t in range(T): g = shifted(t) + beta ** -t * gam nx[t] = np.abs(g).sum() g[t] -= 1 ns[t] = np.abs(g).sum() L = np.abs(hf).sum() + np.abs(hb[1:]).sum() return L, L - abs(hf[0]) + abs(hf[0] - 1), nx, ns def by_roots(d): """norms' halves of h from the roots.""" beta, c0, h = inverse(d) hb = np.concatenate([[0.0], -c0 * beta ** -np.arange(1.0, 200 * d)]) return h, hb def gmat(d, n): """g^(t) for t < n as a matrix over k < K + n, by the roots.""" hf, hb = by_roots(d) beta = hb[1] / hb[2] m = len(hf) + n hz = np.concatenate([hf, np.zeros(n)]) k = np.arange(m) def shifted(t): v = np.empty(m) f = k >= t v[f] = hz[k[f] - t] v[~f] = hb[t - k[~f]] return v gam = -sum(shifted(s) for s in range(d)) / sum(beta ** -s for s in range(d)) return np.array([shifted(t) + beta ** -t * gam for t in range(n)]) def wave(d): P, a = d - 1, (d - 1) // 4 x = list(range(a + 1)) + list(range(a - 1, -a - 1, -1)) + \ list(range(-a + 1, 0)) + [0] * (P - 4 * a) return x[:P] def section_control(): print("CONTROL the wave from g; Delta(F) x = delta for a random delta", flush=True) for d in (9, 13, 21): P = d - 1 G = gmat(d, 4 * P) xw = wave(d) dl = np.array([xw[k % P] - xw[(k + 1) % P] for k in range(G.shape[1])]) x = G @ dl err = max(abs(x[t] - xw[t % P]) for t in range(4 * P)) check(f"d = {d}: the wave's x from g", err < 1e-9, f"max |x - x_wave| {err:.2e}") rng = np.random.default_rng(1) for d in (5, 12, 30): G = gmat(d, 200 + d) dl = rng.integers(-1, 2, G.shape[1]) x = G @ dl res = max(abs(x[t] - x[t + 1:t + d + 1].sum() - dl[t]) for t in range(200)) s0 = abs(x[:d].sum()) check(f"d = {d}: a random delta's x solves every equation", res < 1e-9 and s0 < 1e-9, f"max residual {res:.2e}, first window sum {s0:.2e}") def row(d, L, Ls, nx, ns): C = max(1, (d - 1) // 4) hx, hs = nx.max(), ns.max() tail = 1e-21 * d * L proved = (math.floor(max(hx, L + tail)) <= C and math.floor(max(hs, Ls + tail)) <= C) print(f" d = {d:3}: C {C}, L {L:.6f}, L_s {Ls:.6f}, d(L - (d-1)/4) " f"{d * (L - (d - 1) / 4):.4f}, H_x {hx:.6f} (t = " f"{int(nx.argmax())}), H_s {hs:.6f}, margin " f"{C + 1 - max(hx, hs):.4f}, " f"{'PROVED' if proved else 'open'}", flush=True) return proved, L <= hx + 1e-12 and Ls <= hs + 1e-12, hx, hs def section_h(): print("H the roots, d = 3..40", flush=True) proved, past = set(), [] for d in range(3, 41): L, Ls, nx, ns = norms(d, *by_roots(d)) p, ok, *_ = row(d, L, Ls, nx, ns) if p: proved.add(d) if not ok: past.append(d) want = {3, 4, 9, 10, 13, 14, 15, 17, 18, 19} | set(range(21, 41)) check("the sups past T within the sups read: L <= H_x, L_s <= H_s", not past, f"exceptions {past}") check("the bound proved at d = 3, 4, 9, 10, 13..15, 17..19, 21..40, " "left at d = 5..8, 11, 12, 16, 20", proved == want, f"proved {sorted(proved)}") def section_h3(): print("H3 CONTROL the FFT against the roots", flush=True) for d in (5, 12, 30, 40): L0, Ls0, nx0, ns0 = norms(d, *by_roots(d)) hf, hb, al = inverse_fft(d) L, Ls, nx, ns = norms(d, hf, hb) e = max(abs(L - L0), abs(Ls - Ls0), abs(nx.max() - nx0.max()), abs(ns.max() - ns0.max())) check(f"d = {d}: L, L_s, H_x, H_s agree", e < 1e-9 and al < 1e-15, f"largest difference {e:.2e}, aliasing {al:.1e}") print("H3 the FFT, d = 21..100", flush=True) bad, past, scaled = [], [], [] for d in range(21, 101): hf, hb, al = inverse_fft(d) L, Ls, nx, ns = norms(d, hf, hb) p, ok, *_ = row(d, L, Ls, nx, ns) if not p or al >= 1e-15: bad.append(d) if not ok: past.append(d) scaled.append(d * (L - (d - 1) / 4)) check("the sups past T within the sups read: L <= H_x, L_s <= H_s", not past, f"exceptions {past}") check("the bound proved at every d = 21..100, aliasing below 1e-15", not bad, f"exceptions {bad}") check("H3b: d (L - (d - 1)/4) falls at every step d = 21..100", all(b < a for a, b in zip(scaled, scaled[1:])), f"{scaled[0]:.4f} at 21, {scaled[-1]:.4f} at 100") def main(): section_control() section_h() section_h3() print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed") sys.exit(0 if all(CHECKS) else 1) if __name__ == "__main__": main()