""" route.py -- the carry bound past 100 by the shift norms: the real support's sup over t bounded by L and D_s alone, h streamed by the root sum and never held whole. QUESTION X. At each d = 101..400, is floor(B_x) <= C_d and floor(B_s) <= C_d, where B_x = L + max(max_(t S (D_0 = 0), eps_t = beta^-t/B, B = 1 + beta^-1 + ... + beta^-(d-1), L_s = L - |h_0| + |h_0 - 1| and D_s = sum over all m of |h_m - h_(m+s)|? BY HAND. (X1) g^(t)_k = h_(k-t) - eps_t h_(k+1) on k >= 0 (heat.py, S), so g^(t)_k - [k = t] = (1 - eps)(h_(k-t) - [k = t]) + eps (h_(k-t) - [k = t] - h_(k+1)), and over k >= 0, a part of the sum over m, ||g^(t) - e_t||_1 <= (1 - eps) L_s + eps (D_(t+1) + 1); without e_t, (1 - eps) L + eps D_(t+1). 0 < eps_t <= 1. (X2) D is subadditive, D_(a+b) <= D_a + D_b (the triangle inequality through h_(m+a)), and D <= 2L, so D_s <= D'_s; t >= T: eps_t (D - L)^+ <= eps_T L and eps_t (D + 1 - L_s)^+ <= eps_T (2L + 1). T is 64, not S: eps_S (2L + 1) is 1.5e-3 at d = 400, above X1's bar and rising with L. (X3) past m = M = P^3, |h_m| <= sum_k |c_k| |x_k|^m, so the tail is at most tau = sum_k |c_k| |x_k|^M/(1 - |x_k|); a root dropped once |x|^m < 1e-20 owes at most |c| 1e-20/(1 - |x|); below m = -80 d, sum |c_0| beta^m <= |c_0| beta^(-80 d)/(1 - 1/beta). All are added to L and L_s, twice to each D_s. (X4) C_d + 1 - W at d = 0 mod 4 is 1/4 + 1/(4P) - 1/P^2, the least of the four classes of d mod 4; with E ~ (4.7071 + 15.68/P)/P there, the least margin is at d = 104, 0.2052 less the excess. (X5) streamed D_s: with ext = h on [m0 - S, m0 + n), D_s += sum |ext[S - s : S - s + n] - ext[S : S + n]|, each pair (m, m + s) counted once, as m + s enters; a closing chunk of S zeros flushes the pairs past the end. Integrity at every d: the class sums of the streamed h against the sawtooth ((P - 1)/2 - r)/P - 1/P^2, and h_0 against 1 - c_0 (the residues sum to 1/Delta(0) = 1), c_0 the coefficient at beta. PREDICTIONS, CONTROL AND KILL (fixed before the runs, X0d after). X0 CONTROL, read first. (a) The streamed L, E, L_s and D_1..D_16 at d = 40, 100 equal the FFT h's (real_support.inverse_fft) to 1e-9. (b) The route is never below the true sups: B_x >= H_x and B_s >= H_s (real_support.norms on the FFT h) at d = 13, 24, 40, 60, 100. (c) The route prints OPEN at d = 11, 12, 16, 20, where H_x > C_d + 1, and (d) closed at d = 100, where H_x and H_s close the bound. X1 every d = 101..400 closes; the excess max(B_x - L, B_s - L_s) is below 1e-3 at every d (X1a) and falls (X1b). X2 the least margin C_d + 1 - max(B_x, B_s) is at d = 104, within 0.002 of 0.2052. KILL of the route: a d in 101..400 printed OPEN. FINDINGS (copied from the runs' prints). - X0 holds. (a) The streamed L, E, L_s, D_1..D_16 equal the FFT's to 8.6e-13 (d = 40) and 9.3e-12 (d = 100). (b) B_x - H_x = 0.281, 0.150, 0.0593, 0.0162, 0.00115 and B_s - H_s = 0.531, 0.275, 0.0915, 0.0247, 0.00178 at d = 13, 24, 40, 60, 100: never below, loose at small d. (c) OPEN at d = 11, 12, 16, 20 (B_x 3.355, 3.545, 4.350, 5.232). (d) closed at d = 100, margin 0.2015. - The KILL does not fire: every d = 101..400 closes. Integrity at every d: the class sums on the sawtooth to 2.6e-11, h_0 = 1 - c_0 to 1.1e-16, root residuals at most 3.6e-13, roots at least 0.016 apart, the tail tau at most 9.3e-16. - X1a fails as frozen: the excess max(B_x - L, B_s - L_s) is 1.63e-3 at d = 101 and at least 1e-3 at d = 101..108. The frozen estimate, eps_15 0.357 L, about 1.4e-4 at d = 100 from heat.py's D_16/L, put the excess at t = 15; it sits at t = 9 at d = 101, where eps is larger, and t rises with d, never falling, to 28 at d = 400. B_s's excess is the larger at every d. X1b holds: the excess falls at every step, to 7.36e-9 at d = 400. - X2 holds: the least margin is 0.2038 at d = 104. By d mod 4 (0, 1, 2, 3) the least margins are 0.2038 (104), 0.9521 (101), 0.7027 (102), 0.4533 (103), rising to 0.2387, 0.9886, 0.7387, 0.4887 at d = 400, 397, 398, 399: the offsets C_d + 1 - W (1/4, 1, 3/4, 1/2 up to O(1/P)) less E and the excess. - L_s = L to the printed 9 decimals at every d: 1 - 2 h_0 = 2 c_0 - 1 is of order d 2^-d. P E falls at every step of 2 within each parity of d, 4.6233 to 4.49904 (odd d), 4.8618 to 4.74645 (even d). RUN RECORD: the bare run, d = 101..400, 624.4 s, peak commit 233.1 MB; 0.1 s per d at d = 101, 6.1 s at d = 400; 16 of 17 checks, the miss X1a as recorded. """ 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, norms # noqa: E402 S, T = 16, 64 CHECKS, NAMES = [], [] MISSED = ("X1a",) # the record's one miss 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) class Acc: """L, E, the class sums and D_1..D_S over a stream of h in m order.""" def __init__(self, P): self.P = P r = np.arange(P) self.sign = np.sign(((P - 1) / 2 - r) / P - 1 / P ** 2) self.prev = np.zeros(S) self.L = self.E = 0.0 self.H = np.zeros(P) self.D = np.zeros(S + 1) def add(self, m0, h): n = len(h) cls = (m0 + np.arange(n)) % self.P self.L += np.abs(h).sum() self.E += (np.abs(h) - self.sign[cls] * h).sum() self.H += np.bincount(cls, weights=h, minlength=self.P) ext = np.concatenate([self.prev, h]) for s in range(1, S + 1): self.D[s] += np.abs(ext[S - s:S - s + n] - ext[S:S + n]).sum() self.prev = ext[-S:].copy() def close(self): """Flush the pairs (m, m + s) with m + s past the end, h = 0 there.""" ext = np.concatenate([self.prev, np.zeros(S)]) for s in range(1, S + 1): self.D[s] += np.abs(ext[S - s:S - s + S] - ext[S:2 * S]).sum() def stream(d): """h by the root sum over m in [-80 d, P^3), roots dropped once |x|^m < 1e-20: L, E, h_0, D_1..D_S, the tail bound, beta, the integrity.""" P = d - 1 x, wt, b, resid = roots(d) c = coef(x, d) * wt lx = np.log(x) absx = np.abs(x) acc = Acc(P) mneg = np.arange(-80 * d, 0) acc.add(int(mneg[0]), -coef(b, d) * b ** mneg.astype(float)) M = P ** 3 order = np.argsort(-absx) m0, h0 = 0, None 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 if h0 is None: h0 = h[0] acc.add(m0, h) m0 += n acc.close() ac = np.abs(c) tail = (ac * absx ** M / (1 - absx)).sum() + (ac * 1e-20 / (1 - absx)).sum() tail += abs(coef(b, d)) * b ** (-80.0 * d) / (1 - 1 / b) r = np.arange(P) saw = np.abs(acc.H - (((P - 1) / 2 - r) / P - 1 / P ** 2)).max() herr = abs(h0 - (1 - coef(b, d))) dist = np.abs(x[:, None] - x[None, :]) + np.eye(len(x)) integ = (saw, herr, resid, dist.min()) return acc.L, acc.E, h0, acc.D, tail, b, integ def route(d, L, h0, D, b, tail=0.0): """B_x, B_s and the two excesses, the tail folded in (X1-X3).""" B = sum(b ** -s for s in range(d)) eps = b ** -np.arange(T + 1.0) / B L = L + tail Ls = L - abs(h0) + abs(h0 - 1) Dt = D[1:S + 1] + 2 * tail Dq = np.empty(T + 1) Dq[:S + 1] = np.concatenate([[0.0], Dt]) for s in range(S + 1, T + 1): Dq[s] = min(2 * L, (s // S) * Dq[S] + Dq[s % S]) vx = np.append(eps[:T] * np.maximum(Dq[1:] - L, 0), eps[T] * L) vs = np.append(eps[:T] * np.maximum(Dq[1:] + 1 - Ls, 0), eps[T] * (2 * L + 1)) ex, es = float(vx.max()), float(vs.max()) tm = int(vx.argmax()) if ex >= es else int(vs.argmax()) return L + ex, Ls + es, ex, es, L, Ls, tm def fft_stats(d): """L, E, h_0 and D_1..D_S from the FFT h, held whole (d <= 100).""" hf, hb, _ = inverse_fft(d) P = d - 1 h = np.concatenate([hb[:0:-1], hf, np.zeros(S)]) off = len(hb) - 1 m = np.arange(len(h)) - off H = np.bincount(m % P, weights=h, minlength=P) E = (np.abs(h) - np.sign(H)[m % P] * h).sum() D = np.zeros(S + 1) hz = np.concatenate([np.zeros(S), h]) for s in range(1, S + 1): D[s] = np.abs(hz[S - s:len(hz) - s] - hz[S:]).sum() b = hb[1] / hb[2] return np.abs(h).sum(), E, hf[0], D, b, hf, hb def row(d, Bx, Bs): C = max(1, (d - 1) // 4) closed = math.floor(Bx) <= C and math.floor(Bs) <= C return C, closed, C + 1 - max(Bx, Bs) def section_control(): print("X0 CONTROL", flush=True) for d in (40, 100): L, E, h0, D, tail, b, integ = stream(d) Lf, Ef, h0f, Df, bf, _, _ = fft_stats(d) Ls = L - abs(h0) + abs(h0 - 1) Lsf = Lf - abs(h0f) + abs(h0f - 1) err = max(abs(L - Lf), abs(E - Ef), abs(Ls - Lsf), np.abs(D[1:] - Df[1:]).max()) check(f"X0a streamed L, E, L_s, D_1..D_16 = the FFT's at d = {d}", err < 1e-9, f"max difference {err:.1e}; tail {tail:.1e}; " f"sawtooth {integ[0]:.1e}, h_0 {integ[1]:.1e}") print(" d H_x B_x H_s B_s", flush=True) for d in (13, 24, 40, 60, 100): L, E, h0, D, tail, b, integ = stream(d) Bx, Bs, ex, es, *_ = route(d, L, h0, D, b, tail) _, _, _, _, _, hf, hb = fft_stats(d) _, _, nx, ns = norms(d, hf, hb) Hx, Hs = nx.max(), ns.max() print(f" {d:3} {Hx:.8f} {Bx:.8f} {Hs:.8f} {Bs:.8f}", flush=True) check(f"X0b the route bounds the sups at d = {d}", Bx >= Hx - 1e-12 and Bs >= Hs - 1e-12, f"B_x - H_x {Bx - Hx:.2e}, B_s - H_s {Bs - Hs:.2e}") for d in (11, 12, 16, 20): L, E, h0, D, tail, b, integ = stream(d) Bx, Bs, *_ = route(d, L, h0, D, b, tail) C, closed, mg = row(d, Bx, Bs) check(f"X0c the route prints OPEN at d = {d}", not closed, f"B_x {Bx:.6f}, B_s {Bs:.6f}, C_d + 1 = {C + 1}") L, E, h0, D, tail, b, integ = stream(100) Bx, Bs, *_ = route(100, L, h0, D, b, tail) C, closed, mg = row(100, Bx, Bs) check("X0d the route prints closed at d = 100, closed by the FFT's sups", closed, f"B_x {Bx:.6f}, B_s {Bs:.6f}, margin {mg:.4f}") def section_past(lo, hi): print(f"X THE ROUTE, d = {lo}..{hi}", flush=True) opened, exs, margins, bad = [], [], {}, [] for d in range(lo, hi + 1): t0 = time.time() L, E, h0, D, tail, b, integ = stream(d) Bx, Bs, ex, es, Lt, Ls, tm = route(d, L, h0, D, b, tail) C, closed, mg = row(d, Bx, Bs) margins[d] = mg exs.append(max(ex, es)) if not closed: opened.append(d) if integ[0] > 1e-9 or integ[1] > 1e-9 or integ[2] > 1e-9 \ or integ[3] < 1e-9: bad.append(d) print(f" d = {d}: C {C}, L {Lt:.9f}, L_s {Ls:.9f}, P E " f"{(d - 1) * E:.5f}, excess {ex:.2e} / {es:.2e} at t = {tm}, " f"B_x {Bx:.6f}," f" B_s {Bs:.6f}, margin {mg:.4f}, " f"{'closed' if closed else 'OPEN'}; tail {tail:.1e}, sawtooth " f"{integ[0]:.1e}, h_0 {integ[1]:.1e}, roots {integ[2]:.1e} " f"apart {integ[3]:.1e} ({time.time() - t0:.1f} s)", flush=True) check("integrity at every d: sawtooth, h_0, residual to 1e-9, roots " "distinct", not bad, f"exceptions {bad}") check(f"KILL not fired: every d = {lo}..{hi} closed", not opened, f"open {opened}") check("X1a the excess below 1e-3 at every d", max(exs) < 1e-3, f"largest {max(exs):.2e}") check("X1b the excess falls at every step", all(b_ < a_ for a_, b_ in zip(exs, exs[1:])), f"{exs[0]:.2e} at {lo}, {exs[-1]:.2e} at {hi}") dmin = min(margins, key=margins.get) check("X2 the least margin at d = 104 within 0.002 of 0.2052", dmin == 104 and abs(margins[104] - 0.2052) <= 0.002, f"least {margins[dmin]:.4f} at d = {dmin}") def main(): lo = int(sys.argv[1]) if len(sys.argv) > 1 else 101 hi = int(sys.argv[2]) if len(sys.argv) > 2 else 400 section_control() if not all(CHECKS): print("control failed; stopping") sys.exit(1) section_past(lo, hi) 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()