""" remainder.py -- a bound on |h - phi| past the edge built from per-term triangle inequalities, images below and modes above, and the edge's per-term bound Ebar in closed form. QUESTION V. Does a bound B_m >= |h_m - phi_m| built from per-term triangle inequalities (no h, no cancellation between terms) keep R_prov = sum_(m >= m1, s phi <= B) (|phi| - s phi) + 2B below 1/4 - Ebar at d = 400; and does Ebar have a closed form, and how does it scale with P at fixed c? Notation as edge.py: P = d - 1, m = qP + r, x = r/P, tau = m/P^3, phi_m = -P^-2 K_x(x, tau), s the sawtooth's sign, m1 = P^3/c, c = 400. BY HAND. (V1) Images. j = m - uP = kP + r with k = q - u, xi = j/P = x + k. The walk's u-term T_u(m) = b(u, n) j/(2n), n = 2u - j = n0 - xi(P + 2), n0 = 2 tau P^2, and phi's image k, P^-2 (4 pi tau)^(-1/2) e^(-xi^2/(4 tau)) xi/(2 tau), are the same function at n = n0. So B_img = sum_k |T_(q-k) - phi_k| + |sum_k phi_k - phi| + tail, T = 0 where u < 0 or n < max(u, 1). Past |j| = J, b <= e^(-j^2/(2n)) (Hoeffding), |j|/(2n) <= 1/2 and n <= n0 + B|j|, B = 1 + 2/P, so with g(a) = a^2/(2(n0 + B a)), convex, the walk's terms beyond the images read cost at most e^(-g(J))/(1 - e^(-P g'(J))) per side. TRANSPLANT from Stirling: T_k/phi_k - 1 ~ (3 xi/(4 tau) - xi^3/(8 tau^2))/P, so B_img ~ P^-3 C/tau with no e^(-4 pi^2 tau) decay, while |phi| has it: past tau* ~ ln P/(4 pi^2) every class is Bad. (V2) Modes. h_m = sum over the root pairs of wt Re(c_k x_k^m), x_k near e^(2 pi i k/P), and phi_m = P^-2 sum_(k>=1) 4 pi k sin(2 pi k x) e^(-4 pi^2 k^2 tau); pair root k with mode k: B_four = sum_(k <= K) |wt Re(c_k x_k^m) - mode_k| + |sum_(k <= K) mode_k - phi| + sum_(dropped) wt |c| |x|^m, a root dropped once |x|^m < 1e-20. Per mode the relative miss is first order, k^3 tau/P (the drift i theta^3 m/P). (V3) B = min(B_img, B_four), B_img read at tau < 1/4, B_four at tau >= 1/64. Past M = P^3, |h| - sh <= 2|h| <= 2 sum wt |c| |x|^M/(1 - |x|). Roots and the root sum are floating, as in route.py. (V4) The edge. b(u, n)|j|/n = Pr(tau_j = n), the simple walk's first passage to j (the hitting time theorem, S_n = 2u - n = j). For fixed j the m < m1 run over n = 2u - j <= N_j = 2(m1 - j)/P - j, so the wrong-signed sum at j is at most Pr(tau_|j| <= N) <= 2 Pr(S_N >= |j|) <= 2 e^(-j^2/(2N)) (reflection, Hoeffding). With N <= A + B a, A = 2 m1/P, a = |j| >= J0 = ceil((P - 1)/2), and g as above with A for n0, Ebar <= F = 4 e^(-g(J0))/(1 - e^(-g'(J0))), g'(a) = a(2A + B a)/(2(A + B a)^2). At fixed c, g(J0) -> c/16 and g'(J0) -> c/(4P), so F ~ (16P/c) e^(-c/16). F(400, 399) = 5.6e-8. PREDICTIONS, CONTROL AND KILL (fixed before the code). V0 CONTROL, read first. (a) The images' T sum equals h (root sum) to 1e-14 and their phi sum equals the Fourier phi to 1e-14, at d = 100 on m sampled over tau in [1/400, 1/4). (b) With B = |h - phi| the R machinery returns edge.py's R: 0.27568, 0.12140, 0.04979 at d = 100, 200, 400 (4 digits). (c) F >= Ebar at d = 100, 200, 400 (edge.py's 7.070e-7, 1.504e-8, 9.931e-10). V1 Images alone (B = B_img on [1/c, 1/4), nothing past): its R over [1/16, 1/4) exceeds the oracle's there by a factor above 3 at d = 400. V2 The min route: P R_prov in [20, 60] at d = 100, 200, 400. V3 F(400, 399) <= 1e-7; F(400, 10^6) / ((16P/c) e^(-c/16)) in [0.8, 1.25]. V0d (a sanity read added with the code, not frozen): B >= |h - phi| at every m read, to 1e-15, the root sum's float error (edge.py Y0a). V0e (added after the run, not frozen): h to 40 digits by the binomial series on 300 sampled m at tau < 1/4, against the root sum over all roots and the float binomial sum. KILL of the per-term route: R_prov + F >= 1/4 at d = 400; then the comparison function must carry the first order (phi1, heat.py's correction) before any majorant is worth deriving. FINDINGS (copied from the run's prints). - V0 holds but V0d at d = 400. (a) The images' T sum meets the root sum to 9.1e-16 and their phi sum the Fourier phi to 3.0e-18 at d = 100. (b) The oracle R is 0.27568, 0.12140, 0.04979, edge.py's. (c) F = 7.000e-5, 1.195e-6, 5.651e-8 >= Ebar at d = 100, 200, 400. (d) B falls below |h - phi| by at most 2.9e-16, 9.0e-16, 1.5e-15 at d = 100, 200, 400, the last past the bar set. On the image branch the float error is of that size (V0e): against h to 40 digits on 300 m in [P^3/400, P^3/4), the root sum over all roots misses by up to 8.84e-16, 1.20e-15 and the float binomial sum by up to 1.30e-15, 1.06e-15 at d = 200, 400. V0d's worst m, the streamed h and B_four's float error are not read, so the shortfall is unattributed, though of the float's size; an exact B owes outward rounding. - The per-term route is tight: R_prov = 0.27603, 0.12158, 0.04988 (P R 27.33, 24.19, 19.90) against the oracle's 0.27568, 0.12140, 0.04979 at d = 100, 200, 400, within 0.23% in every band. V2 fails at d = 400 by being tighter than its band (19.90 < 20). The KILL does not fire: R_prov + F = 0.04988 at d = 400. - V1 holds: images alone over [1/16, 1/4) are 0.659, 0.505, 0.373 at d = 100, 200, 400 against the oracle's 0.0693, 0.0343, 0.0170, a factor 9.5, 14.7, 21.9, rising with P: the modes carry that band. - V3 holds: F(400, 399) = 5.65e-8; F/((16P/c) e^(-c/16)) = 1.307, 1.0027, 1.0000 at P = 1e4, 1e6, 1e8. The per-term edge grows like P at fixed c, past 0.01 only near P = 2e10 at c = 400. RUN RECORD: 15/17 checks, the misses V0d and V2 at d = 400; 52.2 s under the 512 MB guard, peak commit 115.7 MB. The images-alone column is masked at tau >= 1/4, where a chunk straddling 1/4 would carry its values. V0e's tail reads 4.9e-41, 4.8e-43 at d = 200, 400. The exit is 0 iff the misses are the record's, V0d and V2. """ 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 edge import signs, logb, kx # noqa: E402 CHECKS, NAMES = [], [] MISSED = ("V0d B >= |h - phi| at every m, d = 400, to 1e-15", "V2 P R_prov in [20, 60] at d = 400") # the record's two misses C = 400 EDGE_R = {100: 0.27568, 200: 0.12140, 400: 0.04979} EDGE_EBAR = {100: 7.070e-7, 200: 1.504e-8, 400: 9.931e-10} 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 g_and_dg(a, A, B): g = a * a / (2 * (A + B * a)) dg = a * (2 * A + B * a) / (2 * (A + B * a) ** 2) return g, dg def F_edge(c, P): """V4: Ebar over m < P^3/c is at most 4 e^-g(J0)/(1 - e^-g'(J0)).""" A, B = 2 * P * P / c, 1 + 2 / P J0 = math.ceil((P - 1) / 2) g, dg = g_and_dg(J0, A, B) return 4 * math.exp(-g) / (-math.expm1(-dg)), g, dg def images(m, P, tau_hi): """V1: the images' T sum, phi sum and sum |T_k - phi_k|, and the walk's tail past the images read, at each m (tau_hi bounds the chunk's tau).""" q, r = m // P, m % P x = r / P tau = m / P ** 3 K = int(math.ceil(math.sqrt(4 * tau_hi * 64))) + 2 n0max = 2 * tau_hi * P * P B = 1 + 2 / P J = K * P # every image past K has |j| >= K P g, dg = g_and_dg(J, n0max, B) assert g >= 60, g tailT = 2 * math.exp(-g) / (-math.expm1(-P * dg)) / 2 Ts = np.zeros(len(m)) Ps = np.zeros(len(m)) D = np.zeros(len(m)) pref = (4 * np.pi * tau) ** -0.5 / P ** 2 for k in range(-K, K + 1): xi = x + k ph = pref * np.exp(-xi * xi / (4 * tau)) * xi / (2 * tau) j = r + k * P u = q - k n = 2 * u - j ok = (u >= 0) & (n >= np.maximum(u, 1)) T = np.zeros(len(m)) if ok.any(): uu, nn, jj = u[ok].astype(float), n[ok].astype(float), j[ok] T[ok] = np.exp(logb(uu, nn)) * jj / (2 * nn) Ts += T Ps += ph D += np.abs(T - ph) return Ts, Ps, D, tailT def modes(m, P, xs, cw, act_k): """V2: sum_k |wt Re(c_k x_k^m) - mode_k| over the active pairs, and the modes' sum.""" x = (m % P) / P tau = m / P ** 3 D = np.zeros(len(m)) S = np.zeros(len(m)) for k in act_k: root = (cw[k - 1] * np.exp(m * np.log(xs[k - 1]))).real mode = 4 * np.pi * k * np.sin(2 * np.pi * k * x) * \ np.exp(-4 * np.pi ** 2 * k * k * tau) / P ** 2 D += np.abs(root - mode) S += mode return D, S def stream(d): """R_prov (B = min of the two), R with B = |h - phi| (the oracle) and with B = B_img alone on [1/c, 1/4), each by tau band; and the checks that B >= |h - phi|.""" P = d - 1 xs, wt, b, resid = roots(d) cw = coef(xs, d) * wt absx = np.abs(xs) lx = np.log(xs) s = signs(P) M = P ** 3 m1 = M // C edges = np.array([1 / 16, 1 / 4]) out = {key: np.zeros(3) for key in ("prov", "orac", "img")} worst = 0.0 nB = {"img": 0, "four": 0} m0 = m1 while m0 < M: act = np.nonzero(absx ** m0 > 1e-20)[0] n = int(min(2 ** 17, max(4096, m0 // 16), 2 ** 21 // max(1, min(len(act), 64)), M - m0)) m = np.arange(m0, m0 + n, dtype=np.int64) tau = m / P ** 3 thi = (m0 + n) / P ** 3 h = np.zeros(n) for a in range(0, len(act), 64): sl = act[a:a + 64] h += (np.exp(np.outer(m - m0, lx[sl])) @ (cw[sl] * np.exp(m0 * lx[sl]))).real drop = (np.abs(cw) * 1e-20).sum() if len(act) < len(xs) else 0.0 ph = -kx((m % P) / P, tau) / P ** 2 sm = s[m % P] Bimg = np.full(n, np.inf) Bfour = np.full(n, np.inf) if m0 / P ** 3 < 1 / 4: Ts, Ps, D, tailT = images(m, P, thi) Bimg = D + np.abs(Ps - ph) + tailT nB["img"] += n if thi >= 1 / 64: D4, S4 = modes(m, P, xs, cw, act + 1) Bfour = D4 + np.abs(S4 - ph) + drop nB["four"] += n Bmin = np.minimum(Bimg, Bfour) err = np.abs(h - ph) worst = max(worst, float((err - Bmin * (1 + 1e-9) - 1e-300).max())) band = np.searchsorted(edges, tau, side="right") Bimg_alone = np.where(tau < 1 / 4, Bimg, np.inf) for key, B in (("prov", Bmin), ("orac", err), ("img", Bimg_alone)): if key == "img" and not np.isfinite(B).any(): continue fin = np.isfinite(B) bad = fin & (sm * ph <= B) rr = np.where(bad, np.abs(ph) - sm * ph + 2 * np.where(fin, B, 0), 0.0) out[key] += np.bincount(band, weights=rr, minlength=3) m0 += n tail = 2 * (np.abs(cw) * absx ** M / (1 - absx)).sum() out["tail"] = tail out["worst"] = worst out["nB"] = nB return out def section_control(): print("V0a images against the root sum and the Fourier phi, d = 100", flush=True) d = 100 P = d - 1 xs, wt, b, _ = roots(d) cw = coef(xs, d) * wt rng = np.random.default_rng(1) m = np.sort(rng.integers(P ** 3 // 400, P ** 3 // 4, 2000)) h = (np.exp(np.outer(m, np.log(xs))) @ cw).real Ts, Ps, D, tailT = images(m, P, m.max() / P ** 3) tau = m / P ** 3 x = (m % P) / P ph4 = np.zeros(len(m)) for k in range(1, 80): ph4 += 4 * np.pi * k * np.sin(2 * np.pi * k * x) * \ np.exp(-4 * np.pi ** 2 * k * k * tau) / P ** 2 e1, e2 = np.abs(Ts - h).max(), np.abs(Ps - ph4).max() check("V0a images' T sum = h and phi sum = Fourier phi to 1e-14", e1 < 1e-14 and e2 < 1e-14, f"{e1:.1e}, {e2:.1e}; walk tail " f"{tailT:.1e}") def section_edge(): print("V0c, V3 the edge in closed form", flush=True) for d in (100, 200, 400): P = d - 1 F, g, dg = F_edge(C, P) check(f"V0c F >= Ebar at d = {d}", F >= EDGE_EBAR[d], f"F {F:.3e} (g {g:.3f}, g' {dg:.4f}), Ebar {EDGE_EBAR[d]:.3e}") F, g, dg = F_edge(C, 399) check("V3a F(400, 399) <= 1e-7", F <= 1e-7, f"{F:.3e}") for P in (10 ** 4, 10 ** 6, 10 ** 8): F, g, dg = F_edge(C, P) ratio = F / (16 * P / C * math.exp(-C / 16)) print(f" P = {P:.0e}: F {F:.3e}, F/((16P/c) e^(-c/16)) " f"{ratio:.4f}", flush=True) if P == 10 ** 6: check("V3b F(400, 1e6)/((16P/c) e^(-c/16)) in [0.8, 1.25]", 0.8 <= ratio <= 1.25, f"{ratio:.4f}") Pc = 0.01 * C / 16 * math.exp(C / 16) print(f" (16P/c) e^(-c/16) = 0.01 at P = {Pc:.2e}, c = {C}", flush=True) def section_R(): print("V0b, V1, V2, KILL at d = 100, 200, 400", flush=True) res = {} for d in (100, 200, 400): t0 = time.time() P = d - 1 o = stream(d) Rp = o["prov"].sum() + o["tail"] Ro = o["orac"].sum() + o["tail"] res[d] = (o, Rp, Ro) sh = lambda v: " ".join(f"{t:.4e}" for t in v) # noqa: E731 print(f" d = {d}: R_prov {Rp:.5f} (P R {P * Rp:.2f}), oracle " f"{Ro:.5f} (P R {P * Ro:.2f}); by band [1/c, 1/16) [1/16, " f"1/4) [1/4, 1): prov {sh(o['prov'])}, oracle {sh(o['orac'])}," f" images alone {sh(o['img'])}; B below |h - phi| by at most " f"{o['worst']:.1e}; tail {o['tail']:.1e} ({time.time() - t0:.0f}" f" s)", flush=True) check(f"V0b oracle R = edge.py's {EDGE_R[d]} at d = {d} (4 digits)", f"{Ro:.4g}" == f"{EDGE_R[d]:.4g}", f"{Ro:.5f}") check(f"V0d B >= |h - phi| at every m, d = {d}, to 1e-15", o["worst"] <= 1e-15, f"{o['worst']:.1e}") o, Rp, Ro = res[400] f = o["img"][1] / o["orac"][1] check("V1 images alone over [1/16, 1/4) exceed the oracle by a factor " "above 3 at d = 400", f > 3, f"factor {f:.2f}") for d in (100, 200, 400): o, Rp, Ro = res[d] check(f"V2 P R_prov in [20, 60] at d = {d}", 20 <= (d - 1) * Rp <= 60, f"{(d - 1) * Rp:.2f}") F = F_edge(C, 399)[0] o, Rp, Ro = res[400] check("KILL not fired: R_prov + F < 1/4 at d = 400", Rp + F < 0.25, f"{Rp + F:.5f}") def section_float(): """V0e: h to 40 digits by the binomial series (mpmath) on sampled m, against the root sum and the float binomial terms B_img is built of.""" import mpmath as mp mp.mp.dps = 40 print("V0e float error against h to 40 digits", flush=True) for d in (200, 400): P = d - 1 xs, wt, b, _ = roots(d) cw = coef(xs, d) * wt m = np.sort(np.random.default_rng(3).integers(P ** 3 // 400, P ** 3 // 4, 300)) hr = np.concatenate([(np.exp(np.outer(m[a:a + 50], np.log(xs))) @ cw).real for a in range(0, len(m), 50)]) thi = m.max() / P ** 3 Ts, _, _, tl = images(m, P, thi) K = int(math.ceil(math.sqrt(4 * thi * 64))) + 2 ex = [] for mm in m.tolist(): q, r = divmod(mm, P) t = mp.mpf(0) for k in range(-K, K + 1): u, j = q - k, r + k * P n = 2 * u - j if u >= 0 and n >= max(u, 1): t += mp.binomial(n, u) * mp.mpf(2) ** (-n) * j / (2 * n) ex.append(t) er = max(abs(float(mp.mpf(a) - e)) for a, e in zip(hr, ex)) eb = max(abs(float(mp.mpf(a) - e)) for a, e in zip(Ts, ex)) print(f" d = {d}: over 300 m in [P^3/400, P^3/4), the root sum " f"misses h by up to {er:.2e}, the float binomial sum by up to " f"{eb:.2e} (the walk's tail past the images {tl:.1e})", flush=True) def main(): section_control() section_edge() if not all(CHECKS): print("control failed; stopping") sys.exit(1) section_R() section_float() print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed") 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()