""" excess.py -- the larger sup's excess over L past 400, bounded through D_1 = sum_m |h_m - h_(m+1)|, h's total variation, by the walk's bumps and the roots. QUESTION Z. Does an explicit bound D_1^+ >= D_1 leave the excess inside the carry bound's budget at d = 400, and how does it move with P? BY HAND. (Z1) The excess. H_s - L <= (L_s - L) + sup_t eps_t (D_(t+1) + 1 - L_s)^+ (the real support). D_s <= s D_1, so the term is 0 while (t + 1) D_1^+ <= L_s - 1, i.e. t < t1 = floor(t0), t0 = (L_s - 1)/D_1^+. At t = t1 + k, D_(t+1) + 1 - L_s <= (k + 1) D_1^+ and <= 2L + 1 - L_s, so the excess is at most X = (L_s - L) + eps_t1 min(2L + 1 - L_s, (2/beta) D_1^+), max_k (k + 1) beta^-k = 2/beta for 3/2 <= beta < 2. The coarser X_R = eps_t1 (2L + 1) is printed beside. H_x's excess vanishes up to L/D_1^+ >= t0, the smaller. (Z2) The bumps. On m >= 0, h_m = [m = 0]/2 + sum_(u>=1) T_u(m), T_u = b(u, n) j/(2n), j = m - uP, n = 2u - j, on j <= u, 0 beyond (the walk). As a function of j, T_u does not involve P: its total variation V_u is one sequence. By the local limit theorem T_u ~ -phi_n'(j), phi_n the N(0, n) density, n ~ 2u, so u V_u -> kappa = 2/sqrt(2 pi e) = 0.48394, a TRANSPLANT the script reads. (Z3) The bound. With U = floor(tau_c P^2), M = U P: D_1 <= |h_-1| + |h_0 - h_-1| + 1/2 + sum_(u<=U) V_u + 2 sum_(u>U) sum_(m<=M) |T_u(m)| + sum_k wt_k |c_k| |1 - x_k| |x_k|^M/(1 - |x_k|), h on m <= 0 being -c_0 beta^m, monotone; the 1/2 the jump at m = 0; the bumps by the triangle inequality, those past U by their L1 at m <= M (the pairs with m < M); the roots past M, about 2 sum_k e^(-4 pi^2 k^2 tau_c), 4.64 at 1/400, 0.17 at 1/16. Each tau_c gives a bound; D_1^+ is the least over the grid. (Z4) Sizes: D_1 ~ 2 kappa ln P + b, 2 kappa = 0.968, about 6.5 at d = 400; the root tail's triangle inequality loses 1 to 3; t0 ~ 98.75/D_1^+, near 12 to 15. PREDICTIONS, CONTROL AND KILL (frozen before code). Z0 CONTROL, read first. (a) h from the bumps equals the FFT's (real_support.inverse_fft) on [0, M] at d = 40, 100, tau_c = 1/16, to 1e-12; (b) D_1 rebuilt from those pieces (h on m <= 0 by c_0 and beta, the walk to M, the FFT past M) equals the FFT's D_1 to 1e-10, and route.py's streamed D_1 equals the FFT's to 1e-9; (c) PLANT: j shifted by one in T_u fails (a). Z1 u V_u within 2% of kappa at u = 1000, |u V_u - kappa| falling over u = 10, 100, 1000, 10^4. Z2 D_1 - 0.968 ln P within a band of width 0.6 over d = 101 and 110..400 by 10; D_1 in [5.5, 7.5] at d = 400. Z3 D_1^+ - D_1 in [0.5, 3.5] at d = 400; the best tau_c printed. Z4 at d = 400: X_R in [3e-3, 0.1] (the planned bar, X_R below 1e-3, predicted to miss: a transplant of route.py's X1a, itself missed; the budget is what the theorem needs), X in [1e-4, 5e-3]. Z5 X falls with d over the d = 0 mod 4 read. KILL: (K1) D_1^+ below a printed D_1 at any d read (those of Z2 and 800); or (K2) R_maj + F + X at d = 400 at or above 1/4 + 1/(4P) - 1/P^2, R_maj + F = 0.063511 (uniform.py). Past 400, X is read at P = 799, 1599, 3199, 9999 with V_u past its computed range as kappa'/u, kappa' the largest u V_u on its last decade: printed. FINDINGS (copied from the second run's prints). - Z0 holds: the bumps' h equals the FFT's on [0, M] to 4.6e-16 (d = 40) and 1.1e-15 (d = 100); D_1 rebuilt to 1.2e-14 and 5.1e-14, route.py's streamed D_1 to 1.5e-13 and 3.3e-13; the shifted plant misses by 0.25. - Z1 holds: u V_u is 0.625, 0.500314, 0.486121, 0.484145, 0.483959, 0.483947 at u = 1, 10, 100, 1000, 10^4, 4e4 (0.65625 at u = 2), falling onto kappa = 0.483941 from above past u = 2; the largest on 4000..4e4 is 0.483996. - Z2 holds, tighter than frozen: D_1 - 0.968 ln P lies in [0.7780, 0.7790] at every d read (0.7778 at d = 800), D_1 = 6.575283 at d = 400, 7.247336 at 800: D_1 ~ 0.968 ln P + 0.778 (2 kappa = 0.96788), an observation. - Z3 fails on the tight side: D_1^+ - D_1 is 0.1382 at d = 400, not 0.5 to 3.5; the least tau_c is 1/16 at every d, the root tail there 0.1698 (by hand 0.17), the far bumps 1.1e-4. The by-hand loss of 1 to 3 was the root tail at small tau_c, which the grid does not pick. - Z4 holds: at d = 400, t0 = 14.711, X_R = 6.119e-3 (the planned 1e-3 bar missed, as predicted), X = 2.049e-4 (1.003e-4 with D_1 itself). - Z5 fails on its premise: X steps down as t1 = floor(t0) steps and rises slightly between steps (eps_t1 fixed, D_1^+ growing), 8.66e-2 at d = 120, 2.05e-4 at 400; it falls only in envelope. - The KILL does not fire: D_1^+ >= D_1 at every d read; R_maj + F + X = 0.063716 at d = 400 against 0.250620. - Past 400, V_u past 4e4 as 0.483996/u: D_1^+ - 0.968 ln P = 0.9159 at d = 800 (U = 39900 inside the table, so the exact read), 1600, 3200 and 0.9158 at 10^4; X = 5.5e-8, 7.2e-15, 1.8e-27, 1.7e-76. RUN RECORD: 10/12 checks, the misses Z3 and Z5, as recorded; 116.2 s under the 512 MB guard, peak commit 232.4 MB. The first run stopped in its last section: heat.roots' Newton step for beta overflowed b^d at d = 1600; heat.roots now takes beta by its fixed point. Each bump is read on a fixed window, its ends asserted below 1e-30 (the top only where the window stops short of j = u), monotone beyond; the lower end below j_lo is unpriced in far(), 1.2e-38 at u = 1 and 2.7e-83 at u = 1000 by a call of bump() after the run. """ import math import os import sys import time os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") import numpy as np # noqa: E402 from heat import roots, coef # noqa: E402 from real_support import inverse_fft # noqa: E402 from route import stream # noqa: E402 CHECKS, NAMES = [], [] MISSED = ("Z3 D_1^+ - D_1 in [0.5, 3.5] at d = 400", "Z5 X falls with d over the d = 0 mod 4 read") KAPPA = 2 / math.sqrt(2 * math.pi * math.e) TAUS = (1 / 400, 1 / 200, 1 / 100, 1 / 64, 1 / 48, 1 / 32, 1 / 24, 1 / 16, 1 / 12, 1 / 8) UMAX = 40000 RMAJ_F = 0.063511 # uniform.py, P = 399 def check(name, ok, detail=""): CHECKS.append(bool(ok)) NAMES.append(name) print(f" [{'PASS' if ok else 'FAIL'}] {name}" + (f" -- {detail}" if detail else ""), flush=True) def bump(u): """(j, T_u(j)) over j in [j_lo, j_hi], j_hi = min(u, 30 sqrt u + 100), j_lo = -(30 sqrt u + 100): log C(n, u) by lgamma at the top, then C(n + 1, u) = C(n, u)(n + 1)/(n + 1 - u) as j falls.""" w = int(math.ceil(30 * math.sqrt(u))) + 100 j_hi = min(u, w) j = np.arange(j_hi, -w - 1, -1) n = 2 * u - j n0 = int(n[0]) lb0 = math.lgamma(n0 + 1) - math.lgamma(u + 1) - math.lgamma(n0 - u + 1) \ - n0 * math.log(2) step = np.log((n[1:]) / (2.0 * (n[1:] - u))) lb = lb0 + np.concatenate([[0.0], np.cumsum(step)]) return j, np.exp(lb) * j / (2.0 * n), j_hi def V_table(umax): """V_u, u = 1..umax: the total variation of T_u over j, its end jumps included (both below 1e-30 but at j = u, asserted).""" V = np.zeros(umax + 1) for u in range(1, umax + 1): j, T, j_hi = bump(u) end_lo = abs(T[-1]) assert end_lo < 1e-30, (u, end_lo) if j_hi < u: assert abs(T[0]) < 1e-30, (u, T[0]) V[u] = abs(T[0]) + np.abs(np.diff(T)).sum() + end_lo return V def walk_h(P, M, shift=0): """h_m on m = 0..M by the bumps (shift moves each bump's j: the plant).""" h = np.zeros(M + 1) h[0] += 0.5 u = 1 while True: j, T, _ = bump(u) m = u * P + j + shift if m.min() > M: break keep = (m >= 0) & (m <= M) np.add.at(h, m[keep], T[keep]) u += 1 return h def far(P, U): """2 sum_(u > U) sum_(m <= M) |T_u(m)|, M = U P.""" M, s, u = U * P, 0.0, U + 1 while True: j, T, _ = bump(u) m = u * P + j if m.min() > M: break s += np.abs(T[m <= M]).sum() u += 1 return 2 * s def root_tail(d, M, rt=None): x, wt, b, _ = rt or roots(d) c = coef(x, d) ax = np.abs(x) return float((wt * np.abs(c) * np.abs(1 - x) * np.exp(M * np.log(ax)) / (1 - ax)).sum()) def edge(d, rt=None): """|h_-1| + |h_0 - h_-1| + 1/2, and beta, c_0.""" x, wt, b, _ = rt or roots(d) c0 = float(coef(b, d)) hm1, h0 = -c0 / b, 1 - c0 return abs(hm1) + abs(h0 - hm1) + 0.5, b, c0 def D1_bound(d, Vc, kappa_x=None): """min over tau_c of the bound (Z3), with its pieces; V_u past the table by kappa_x/u when kappa_x is given.""" P = d - 1 rt = roots(d) A, b, c0 = edge(d, rt) best = None for tc in TAUS: U = int(tc * P * P) if U > UMAX: if kappa_x is None: continue bumps = Vc[UMAX] + kappa_x * math.log(U / UMAX) # >= sum 1/u else: bumps = Vc[U] fa = far(P, U) rtl = root_tail(d, U * P, rt) tot = A + bumps + fa + rtl if best is None or tot < best[0]: best = (tot, tc, A, bumps, fa, rtl) return best def excess(d, L, h0, D1p, b): Ls = L - abs(h0) + abs(h0 - 1) B = sum(b ** -s for s in range(d)) t0 = (Ls - 1) / D1p t1 = math.floor(t0) eps = b ** -t1 / B XR = eps * (2 * L + 1) X = (Ls - L) + eps * min(2 * L + 1 - Ls, (2 / b) * D1p) return X, XR, t0, t1 def section_control(): print("Z0 CONTROL", flush=True) for d in (40, 100): P = d - 1 U = int(P * P / 16) M = U * P hf, hb, _ = inverse_fft(d) hw = walk_h(P, M) err = np.abs(hw - hf[:M + 1]).max() check(f"Z0a the bumps' h equals the FFT's on [0, M] at d = {d}", err < 1e-12, f"max {err:.1e}, M = {M}") full = np.concatenate([hb[:0:-1], hf, [0.0]]) D1_fft = np.abs(np.diff(full)).sum() A, b, c0 = edge(d) rest = np.abs(np.diff(np.concatenate([hf[M:], [0.0]]))).sum() D1_re = A - 0.5 + np.abs(np.diff(hw)).sum() + rest L, E, h0s, D, tail, bs, integ = stream(d) check(f"Z0b D_1 rebuilt and route.py's streamed D_1 equal the FFT's " f"at d = {d}", abs(D1_re - D1_fft) < 1e-10 and abs(D[1] - D1_fft) < 1e-9, f"FFT {D1_fft:.12f}, rebuilt {D1_re - D1_fft:+.1e}, " f"streamed {D[1] - D1_fft:+.1e}") if d == 40: hp = walk_h(P, M, shift=1) errp = np.abs(hp - hf[:M + 1]).max() check("Z0c PLANT: j shifted by one fails Z0a at d = 40", errp >= 1e-12, f"max {errp:.1e}") def section_V(): print("Z1 V_u, the bumps' total variation", flush=True) t0 = time.time() V = V_table(UMAX) Vc = np.cumsum(V) devs = [] for u in (1, 2, 5, 10, 100, 1000, 10000, 40000): print(f" u = {u}: V_u {V[u]:.6e}, u V_u {u * V[u]:.6f}, " f"sum to u {Vc[u]:.6f}", flush=True) if u in (10, 100, 1000, 10000): devs.append(abs(u * V[u] - KAPPA)) kx = float((np.arange(4000, UMAX + 1) * V[4000:]).max()) print(f" kappa {KAPPA:.6f}; largest u V_u on u = 4000..{UMAX}: {kx:.6f}" f" ({time.time() - t0:.0f} s)", flush=True) check("Z1 u V_u within 2% of kappa at u = 1000, the gap falling", abs(1000 * V[1000] / KAPPA - 1) <= 0.02 and all(b_ < a_ for a_, b_ in zip(devs, devs[1:])), " ".join(f"{v:.2e}" for v in devs)) return Vc, kx def section_read(Vc): print("Z2-Z5 and KILL: D_1, D_1^+ and the excess at d <= 400 and 800", flush=True) ds = [101] + list(range(110, 401, 10)) + [800] rows = {} for d in ds: t0 = time.time() P = d - 1 L, E, h0, D, tail, b, integ = stream(d) L = L + tail D1 = D[1] + 2 * tail tot, tc, A, bumps, fa, rtl = D1_bound(d, Vc) X, XR, t0_, t1 = excess(d, L, h0, tot, b) Xa, _, _, _ = excess(d, L, h0, D1, b) rows[d] = (D1, tot, X, XR, Xa, L, b) print(f" d = {d}: D_1 {D1:.6f} (- 0.968 ln P {D1 - 0.968 * math.log(P):.4f})," f" D_1^+ {tot:.6f} at tau_c 1/{1 / tc:.0f} (edge {A:.4f}, " f"bumps {bumps:.4f}, far {fa:.1e}, roots {rtl:.4f}); t0 {t0_:.3f}," f" X {X:.3e}, X_R {XR:.3e}, X by D_1 itself {Xa:.3e}; beta " f"{b:.12f} ({time.time() - t0:.1f} s)", flush=True) low = [d for d in ds if rows[d][1] < rows[d][0]] check("KILL K1 not fired: D_1^+ >= D_1 at every d read", not low, f"below at {low}") sp = [rows[d][0] - 0.968 * math.log(d - 1) for d in ds if d <= 400] check("Z2 D_1 - 0.968 ln P in a band of width 0.6 over d <= 400, D_1 in " "[5.5, 7.5] at 400", max(sp) - min(sp) <= 0.6 and 5.5 <= rows[400][0] <= 7.5, f"band [{min(sp):.4f}, {max(sp):.4f}], D_1(400) {rows[400][0]:.4f}") g = rows[400][1] - rows[400][0] check("Z3 D_1^+ - D_1 in [0.5, 3.5] at d = 400", 0.5 <= g <= 3.5, f"{g:.4f}") X, XR = rows[400][2], rows[400][3] check("Z4 at d = 400 X_R in [3e-3, 0.1] and X in [1e-4, 5e-3]", 3e-3 <= XR <= 0.1 and 1e-4 <= X <= 5e-3, f"X_R {XR:.3e} (1e-3 bar: {'met' if XR < 1e-3 else 'missed'}), " f"X {X:.3e}") q = [rows[d][2] for d in ds if d <= 400 and d % 4 == 0] check("Z5 X falls with d over the d = 0 mod 4 read", all(b_ < a_ for a_, b_ in zip(q, q[1:])), " ".join(f"{v:.2e}" for v in q)) P = 399 budget = 1 / 4 + 1 / (4 * P) - 1 / P ** 2 check("KILL K2 not fired: R_maj + F + X below the budget at d = 400", RMAJ_F + X < budget, f"{RMAJ_F + X:.6f} against {budget:.6f}") return rows def section_past(Vc, kx): print("PAST 400 (V_u past the table as kappa'/u): an observation", flush=True) for d in (800, 1600, 3200, 10000): P = d - 1 tot, tc, A, bumps, fa, rtl = D1_bound(d, Vc, kx) L = (P * P - 1) / (4 * P) + 1 / P ** 2 # W at odd P; E > 0 left out rt = roots(d) h0 = 1 - float(coef(rt[2], d)) X, XR, t0_, t1 = excess(d, L, h0, tot, rt[2]) print(f" d = {d}: D_1^+ {tot:.4f} at tau_c 1/{1 / tc:.0f} (- 0.968 ln" f" P {tot - 0.968 * math.log(P):.4f}); t0 {t0_:.2f}, X {X:.3e}, " f"X_R {XR:.3e}", flush=True) def main(): t0 = time.time() section_control() if not all(CHECKS): print("control failed; stopping") sys.exit(1) Vc, kx = section_V() section_read(Vc) section_past(Vc, kx) print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed, {time.time() - t0:.0f} s") sys.exit(as_recorded(MISSED)) def as_recorded(missed): """0 iff the checks that failed are exactly the record's misses, by name.""" got = {n for n, ok in zip(NAMES, CHECKS) if not ok} print("misses %s; the record's %s: %s" % ( "; ".join(sorted(got)) or "none", "; ".join(sorted(missed)) or "none", "as recorded" if got == set(missed) else "NOT AS RECORDED")) return 0 if got == set(missed) else 1 if __name__ == "__main__": main()