""" edge.py -- E at small tau by the walk's binomial series, and the remainder the leading heat term leaves past it. QUESTION Y. Is E's share at small tau, m < m1 = P^3/c, bounded rigorously and exponentially small; and does the leading heat term phi = -P^-2 K_x(r/P, m/P^3), as an oracle (h known), leave a remainder below 1/4 on m >= m1 at d = 400? BY HAND. (Y1) G_m = sum_n Pr(S_n = m), S the walk with steps -1, +d (each 1/2). u up-steps in n steps reach u(d + 1) - n, so n = u(P + 2) - m, down-steps u(P + 1) - m >= 0, and G_m = sum_u b(u, n_u), b(u, n) = C(n, u) 2^-n. As C(n - 1, u) = C(n, u)(n - u)/n, b(u, n) - b(u, n - 1) = b(u, n)(2u - n)/n, and 2u - n = m - uP, so with h_m = (G_m - G_(m+1))/2, h_m = [m = 0]/2 + sum_(u >= 0, n_u >= max(u, 1)) T_u(m), T_u(m) = b(u, n_u) (m - uP)/(2 n_u). Read: h_0 = 1/2 - (P/(2(P + 2))) sum_(u>=1) b(u, u(P + 2)), so 1 - 2 h_0 ~ P 2^-(P+2), the u = 1 term; m < 0 has the u = 0 term -2^(m-1), against -c_0 beta^m. (Y2) T_u(m) has the sign of j = m - uP, j = m (mod P). The sawtooth's sign s is + on r < (P - 1)/2, - on r > (P - 1)/2 and on r = (P - 1)/2 at odd P, so T_u is right-signed at j in [-(P + 1)/2, (P - 1)/2), and wrong only where |j| >= (P - 1)/2. Per m, |h| - s h <= sum_u (|T_u| - s T_u), so E over any m-range is at most Ebar = sum over that range and the u with T_u wrong of b |j|/n. Hoeffding: b(u, n) <= exp(-j^2/(2n)) (one-sided tail at |u - n/2| = |j|/2), so Ebar <= Ebar_H = the same with b replaced. At m < m1 = P^3/c, u <= P^2/c, n ~ 2u: the exponent is at least ((P - 1)/2)^2/(4P^2/c) ~ c/16, over ~P^2/c bumps; c = 400 puts Ebar near 400 e^-25 ~ 1e-8 at P = 399, uniformly c/16 in P. (Y3) Oracle remainder on m >= m1: |h| - s h <= (|phi| - s phi) + 2|h - phi| at every m, zero where s h > 0, so E_(>=m1) <= R = sum over Bad of (|phi| - s phi) + 2|h - phi|, Bad = {m >= m1: s phi <= |h - phi|}. TRANSPLANT (heat.py Z4): h - phi ~ P^-3 K_xx z(tau) near each zero x0, phi ~ -P^-3 K_xx rho, rho classes from x0, so Bad is |rho| <= |z| and R ~ (2/P) sum_(x0) int |K_xx| |z| (2|z| + 1) d tau: P R of order 10 to 50, flat in P; z -> 0 at small tau, so tau0 = 1/c adds little. PREDICTIONS, CONTROL AND KILL (fixed before the code). Y0 CONTROL, read first. (a) The binomial series equals the root sum's h to 1e-12 on m in [-80, 40P] at d = 13, 40, 100. (b) At c = 4 (tau0 = 1/4) E_(= E_(= E_(>=m1) at every d read. Y1 THE EDGE: at c = 400, Ebar <= 1e-6 and Ebar_H <= 1e-6 at d = 100, 200, 400. Y2 THE REMAINDER: P R between 10 and 50 at d = 100, 200, 400. KILL of the leading phi: R >= 1/4 at d = 400 (P = 399, d = 0 mod 4, margin 1/4 + 1/(4P) - 1/P^2 less E and the excess); phi must then carry the first order. FINDINGS (copied from the run's prints). - Y0 holds. (a) The binomial series equals the root sum on m in [-80, 40P] to 3.3e-16, 4.4e-16, 9.2e-16 at d = 13, 40, 100; 1 - 2 h_0 = 7.336e-4 against P 2^-(P+2) = 7.324e-4 at d = 13, 1.774e-11 against 1.774e-11 at d = 40. (b) At c = 4, E's share below m1 at d = 400 is 0.99545. (c) Ebar >= E below m1 and R >= E past m1 at every d read; only Ebar's half at c = 400 can fail, R's holding term by term and Ebar at c = 4 being 125 against E's 0.0118. The window W drops terms with j^2/(2n) above about 64 (by hand), so truncation cannot move Ebar at the printed digits. - Ebar is tight: E below m1 against Ebar at c = 400 is 7.070e-7 / 7.070e-7 (d = 100), 1.503e-8 / 1.504e-8 (200), 9.835e-10 / 9.931e-10 (400): the edge's wrong-signed mass is single bumps' tails. - Y1 fails as frozen at d = 100 only: Ebar_H is 1.68e-5 there (Ebar 7.07e-7 holds); 4.61e-7 and 4.48e-8 at d = 200, 400. Hoeffding drops the local factor of b, about (pi n/2)^(-1/2), and the bump count P^2/c multiplies e^(-c/16); the bar is met from d = 200. - Y2 holds, and the KILL does not fire: P R = 27.29, 24.16, 19.87 at d = 100, 200, 400, falling, not flat; R = 0.0498 at d = 400, a fifth of 1/4. R's shares over [tau0, 1/16), [1/16, 1/4), [1/4, 1): 0.746 0.252 0.002 (d = 100), 0.714 0.283 0.003 (200), 0.655 0.342 0.004 (400). RUN RECORD: 14/15 checks, the miss Y1 at d = 100, which the exit answers for (0 iff the misses are the record's); 32.3 s, peak 218.3 MB. """ import math import os import sys import time os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") import numpy as np # noqa: E402 from scipy.special import gammaln # noqa: E402 from heat import roots, coef # noqa: E402 CHECKS, NAMES = [], [] MISSED = ("Y1",) # the record's one miss LN2 = math.log(2) 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 signs(P): r = np.arange(P) return np.sign(((P - 1) / 2 - r) / P - 1 / P ** 2) def logb(u, n): """log of C(n, u) 2^-n.""" return gammaln(n + 1) - gammaln(u + 1) - gammaln(n - u + 1) - n * LN2 def h_binom(d, ms): """h_m = [m = 0]/2 + sum_u b(u, n)(m - uP)/(2n), n = u(P + 2) - m (Y1), at each m of ms.""" P = d - 1 out = [] for m in ms: u0 = max(0, -(-m // (P + 1))) u = np.arange(u0, max(u0, m // P) + 12, dtype=float) n = u * (P + 2) - m ok = (n >= np.maximum(u, 1)) u, n = u[ok], n[ok] t = np.exp(logb(u, n)) * (m - u * P) / (2 * n) out.append(t.sum() + (0.5 if m == 0 else 0.0)) return np.array(out) def ebar(d, m1, mlo): """Ebar and Ebar_H (Y2) over m in [mlo, m1): sum over the wrong-signed T_u of b |j|/n, and the same with b <= exp(-j^2/(2n)).""" P = d - 1 s = signs(P) tot = totH = 0.0 u = 0 while True: W = int(130 + 16 * math.sqrt(u + 1)) jlo = max(mlo - u * P, -W) jhi = min(u, m1 - 1 - u * P, W) if u * P - W >= m1: break if jlo <= jhi: j = np.arange(jlo, jhi + 1) j = j[j != 0] n = (2 * u - j).astype(float) j, n = j[n >= max(u, 1)], n[n >= max(u, 1)] wrong = s[j % P] != np.sign(j) j, n = j[wrong], n[wrong] if len(j): aj = np.abs(j) tot += (np.exp(logb(float(u), n)) * aj / n).sum() totH += (np.exp(-j.astype(float) ** 2 / (2 * n)) * aj / n).sum() u += 1 return tot, totH def kx(x, tau): """K_x of the periodic heat kernel (K_tau = K_xx, period 1), by images below tau = 0.05 and by Fourier above.""" out = np.zeros_like(tau) lo = tau < 0.05 if lo.any(): t, xx = tau[lo], x[lo] acc = np.zeros_like(t) for k in range(-3, 4): y = xx - k acc += (4 * np.pi * t) ** -0.5 * np.exp(-y * y / (4 * t)) \ * (-y / (2 * t)) out[lo] = acc hi = ~lo if hi.any(): t, xx = tau[hi], x[hi] q = np.exp(-4 * np.pi ** 2 * t) qk, q2k = q.copy(), q ** 3 # q^(k^2), q^(2k + 1) acc = np.zeros_like(t) for k in range(1, 11): acc += -4 * np.pi * k * np.sin(2 * np.pi * k * xx) * qk qk = qk * q2k q2k = q2k * q * q out[hi] = acc return out def stream(d, m1): """E below m1 and at or past it, and the oracle remainder R past m1 (Y3) with its tau bands, h by the root sum over [-80 d, P^3).""" P = d - 1 x, wt, b, resid = roots(d) c = coef(x, d) * wt lx = np.log(x) absx = np.abs(x) s = signs(P) out = dict(Elt=0.0, Ege=0.0, R=0.0, Rb=np.zeros(3)) edges = np.array([1 / 16, 1 / 4]) def take(m, h): e = np.abs(h) - s[m % P] * h lt = m < m1 out["Elt"] += e[lt].sum() out["Ege"] += e[~lt].sum() if (~lt).any(): mm, hh = m[~lt], h[~lt] ph = -kx((mm % P) / P, mm / P ** 3) / P ** 2 sm = s[mm % P] bad = sm * ph <= np.abs(hh - ph) rr = (np.abs(ph) - sm * ph + 2 * np.abs(hh - ph))[bad] out["R"] += rr.sum() band = np.searchsorted(edges, mm[bad] / P ** 3, side="right") out["Rb"] += np.bincount(band, weights=rr, minlength=3) mneg = np.arange(-80 * d, 0) take(mneg, -coef(b, d) * b ** mneg.astype(float)) M = P ** 3 order = np.argsort(-absx) m0 = 0 while m0 < M: act = order[absx[order] ** m0 > 1e-20] R = max(len(act), 1) n = int(min(max(m0 // 8, 4096), 2 ** 20, 2 ** 21 // min(R, 64), M - m0)) j = np.arange(n) base = c[act] * np.exp(m0 * lx[act]) h = np.zeros(n) for a in range(0, len(act), 64): sl = slice(a, a + 64) h += (np.exp(np.outer(j, lx[act][sl])) @ base[sl]).real take(m0 + j, h) m0 += n ac = np.abs(c) tail = (ac * absx ** M / (1 - absx)).sum() + (ac * 1e-20 / (1 - absx)).sum() out["tail"] = tail return out def section_control(): print("Y0 CONTROL", flush=True) for d in (13, 40, 100): P = d - 1 ms = np.arange(-80, 40 * P + 1) hb = h_binom(d, ms) x, wt, b, _ = roots(d) c = coef(x, d) * wt pos = ms[ms >= 0] hr = np.concatenate([-coef(b, d) * b ** ms[ms < 0].astype(float), (np.exp(np.outer(pos, np.log(x))) @ c).real]) err = np.abs(hb - hr).max() check(f"Y0a binomial series = root sum on m in [-80, 40P] at d = {d}", err < 1e-12, f"max difference {err:.1e}; 1 - 2 h_0 = " f"{1 - 2 * hb[80]:.3e} against P 2^-(P+2) = " f"{P * 2.0 ** -(P + 2):.3e}") def section_edge(): print("Y0b, Y0c, Y1, Y2 at d = 100, 200, 400", flush=True) res = {} for d in (100, 200, 400): t0 = time.time() P = d - 1 for cc in ((4, 400) if d == 400 else (400,)): m1 = P ** 3 // cc o = stream(d, m1) eb, ebH = ebar(d, m1, -80 * d) E = o["Elt"] + o["Ege"] R = o["R"] + 2 * o["tail"] res[(d, cc)] = (o, eb, ebH, R) print(f" d = {d}, c = {cc}: P E {P * E:.5f}, E below m1 " f"{o['Elt']:.3e} (share {o['Elt'] / E:.5f}), Ebar {eb:.3e}," f" Ebar_H {ebH:.3e}; E past m1 {o['Ege']:.4e}, R {R:.4e}," f" P R {P * R:.3f}, R's bands [tau0, 1/16) [1/16, 1/4) " f"[1/4, 1) {' '.join(f'{v / o['R']:.4f}' for v in o['Rb'])}" f"; tail {o['tail']:.1e} ({time.time() - t0:.1f} s)", flush=True) check(f"Y0c Ebar >= E below m1 and R >= E past m1, d = {d}, " f"c = {cc}", eb >= o["Elt"] - 1e-15 and R >= o["Ege"]) o, eb, ebH, R = res[(400, 4)] sh = o["Elt"] / (o["Elt"] + o["Ege"]) check("Y0b at c = 4, E's share below m1 at d = 400 is 1 - 0.0045 " "within 0.001", abs(sh - (1 - 0.0045)) <= 0.001, f"{sh:.5f}") for d in (100, 200, 400): o, eb, ebH, R = res[(d, 400)] check(f"Y1 Ebar, Ebar_H <= 1e-6 at c = 400, d = {d}", eb <= 1e-6 and ebH <= 1e-6, f"{eb:.2e}, {ebH:.2e}") for d in (100, 200, 400): o, eb, ebH, R = res[(d, 400)] check(f"Y2 P R in [10, 50] at d = {d}", 10 <= (d - 1) * R <= 50, f"P R {(d - 1) * R:.3f}") o, eb, ebH, R = res[(400, 400)] check("KILL not fired: R < 1/4 at d = 400", R < 0.25, f"R {R:.4f}") def main(): section_control() if not all(CHECKS): print("control failed; stopping") sys.exit(1) section_edge() 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.""" got = {n.split()[0] 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()