remainder.py
Checks: READING: The walk theorem.
Cited on: READING.
To run it: save remainder.py, with edge.py, heat.py, real_support.py beside it; needs mpmath, numpy, scipy.
"""
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()