route.py
Checks: READING: The carry bound rule.
Cited on: READING.
To run it: save route.py, with heat.py, real_support.py beside it; needs mpmath, numpy.
"""
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<T) eps_t (D'_(t+1) - L)^+, eps_T L),
B_s = L_s + max(max_(t<T) eps_t (D'_(t+1) + 1 - L_s)^+, eps_T (2L + 1)),
T = 64, D'_s = D_s read at s <= S = 16 and min(2L, q D_S + D_r) at s =
qS + r > 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()