real_support.py
Checks: READING: The real support theorem; READING: The carry bound rule.
Cited on: READING.
To run it: save real_support.py; needs mpmath, numpy.
"""
real_support.py -- the carries' real support: every carry of a bounded
relaxed path is a linear form in its delta, so its size is at most one
L1 norm per d.
QUESTION. Over every bounded integer x (t >= 0) with delta_t = x_t -
(x_(t+1) + ... + x_(t+d)) in {-1, 0, 1} and first window sum 0 (the
carry sequences of pairs of codings, legality dropped), how large can
|x_t| and the window sum |x_(t+1) + ... + x_(t+d)| be with integrality
dropped too? If the floor of the real sup is C_d = max(1, floor((d -
1)/4)), the triangle wave's carry, the carry bound holds at that d.
BY HAND. delta = Delta(F) x, F the shift, Delta(z) = 1 - z - ... - z^d,
whose roots are rho_0 = 1/beta inside the unit disc and rho_1, ...,
rho_(d-1) outside (their inverses are the d-bonacci polynomial's
roots, beta Pisot). 1/Delta = sum_j c_j/(1 - z/rho_j) with c_j =
-1/(rho_j Delta'(rho_j)), and its Laurent expansion on |z| = 1 is the
bounded inverse h: h_m = sum_(j >= 1) c_j rho_j^(-m) at m >= 0, h_m =
-c_0 beta^m at m <= -1, so h_0 = 1 - c_0 (the c_j sum to 1/Delta(0)).
On t >= 0, x^p_t = sum_(k >= 0) h_(k - t) delta_k solves every equation
(the coefficient of z^(k - t) in Delta(z) h(z) = 1), and is bounded by
L(d) = sum_m |h_m|; Delta(F) x = 0 has the bounded solutions c beta^(-t)
alone, so every bounded solution is x^p + c beta^(-t), and the first
window sum 0 fixes c = -(x^p_0 + ... + x^p_(d-1))/(1 + ... +
beta^(-(d-1))). So x_t = g^(t) . delta, g^(t)_k = h_(k - t) + beta^(-t)
gam_k, gam = -sum_(s < d) h_(. - s)/sum_(s < d) beta^(-s), and |x_t| <=
||g^(t)||_1; the window sum sigma_(t+1) = x_t - delta_t has the form
g^(t) - e_t. H_x(d) and H_s(d) are the sups over t of the two norms.
The sups are read at t < T = 81: past T, ||g^(t)||_1 <= L +
beta^(-t) ||gam||_1 with ||gam||_1 <= d L and beta^(-T) < 1.84^(-80) <
1e-21, and the window sum's row holds h_0 - 1 at k = t, so its bound is
L_s + beta^(-t) d L, L_s = L - |h_0| + |h_0 - 1|. The script prints L
<= H_x and L_s <= H_s, so H_x and H_s are the sups over every t, to the
float.
h two ways: by the roots (mpmath, 40 digits) at d = 3..40, and by the
real FFT of 1/Delta on N = 2^lg points at d = 21..100, lg =
max(16, ceil(log2(2 d^3)) + 1), the slowest root's modulus being about 1 +
4 pi^2/d^3, so h decays over about d^3 terms. The FFT reads the
aliasing, max |h| near N/2, and keeps forward h_m above 1e-18.
1/Delta(z) = (1 - z)/(1 - 2z + z^(d+1)), so h_m = (G_m - G_(m+1))/2,
G the Green's function of the walk with steps -1 and +d, each 1/2:
L(d) is half G's total variation.
PREDICTIONS, CONTROLS AND KILL (fixed before the runs).
H1 L(d)/d tends to a constant in [0.25, 0.35] over d = 3..40, and
L(d) >= (d - 1)/4 at every d.
H2 the route proves the bound (floor H_x, floor H_s <= C_d) at
d = 3..8 at most and fails at some d in 9..12.
KILL of the route at every d: floor(H_x(d)) > C_d at some d in 3..40.
CONTROL, read first: the wave at d = 9, 13, 21 recomputed from g at
t < 4(d - 1) equals its x to 1e-9; g's x for a random delta obeys
Delta(F) x = delta to 1e-9 on t < 200.
H3 (fixed after H's run) the FFT route agrees with the roots to 1e-9
at d = 5, 12, 30, 40 (its control); KILL, floor(H_x) or floor(H_s)
> C_d at some d = 21..100. H3b d (L - (d - 1)/4) falls at every
step d = 21..100.
FINDINGS (copied from the runs' prints).
- CONTROL holds: the wave from g within 3.1e-15, 1.2e-14, 1.9e-14 at
d = 9, 13, 21; residuals 1.4e-15, 5.6e-15, 1.0e-14 and first window
sums below 2.0e-15 at d = 5, 12, 30.
- The KILL fires: floor(H_x) > C_d first at d = 6. The route PROVES
the bound (floor H_x, floor H_s <= C_d) at d = 3, 4, 9, 10, 13, 14,
15, 17, 18, 19 and every d = 21..100; it leaves d = 5..8, 11, 12,
16, 20. H_x at the open d past 8: 3.038251 (11), 3.243495 (12),
4.113062 (16), 5.030733 (20), against C_d + 1 = 3, 3, 4, 5.
- H2 fails the other way: the route passes most d past 12. H1: L(d)/d
falls to 0.2469 at d = 40, below 1/4. The tightest d past 20 are d =
0 mod 4, where C_d + 1 - (d - 1)/4 = 1/4; the closest passes are
d = 3 (H_s 1.994724 against 2), d = 4 (1.981133) and, past 20,
d = 24 (H_x 5.976319, margin 0.0237), then 0.0615 (28), 0.1248
(40), 0.1697 (60), 0.2033 (100) at d = 0 mod 4.
- H3b holds: d (L - (d - 1)/4) falls at every step, 5.5016 (21),
5.0014 (40), 4.8169 (60), 4.6717 (100).
- H_x = L(d) to 1e-6 at d <= 11 (t past the transient); from d = 12
the sup sits at t = 5..11 to d = 40, above L by 5.0e-4 (12), 4.5e-3
(15, the most), 3.2e-3 (21), 2.0e-4 (40), and moves out with d, t =
26 at d = 97..100, equal to L to 1e-6. L <= H_x and L_s <= H_s at
every d = 3..100.
- Float: no H within 1e-6 of an integer; the FFT agrees with the
roots to 8.9e-13 at d = 5, 12, 30, 40, its aliasing below 3e-17.
- Argued, not printed: integrality read one functional at a time pins
no carry at the open d: a value C_d + 1 at t forces delta to sign(g^(t)) up to a defect
of H - C_d - 1, but delta at the d - 4 places after t (h_1, ...,
h_(d-4) small) and more than about four before it costs x_t little
and moves each other x_u by whole units, so no one x_u is pinned.
RUN RECORD: the bare run 108.1 s, peak commit 247.2 MB. The roots at
d = 3..40 first took 6.6 s, 143.5 MB (one matrix of every g^(t) was
refused at 541.6 MB at d = 40, so rows are read one at a time); the FFT
at d = 21..100 on 2^22 points throughout, with T = 69 + 4d, took
633.5 s, 214.8 MB, the cost per d linear in T and N, hence lg and T as
above.
"""
import math
import os
import sys
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import mpmath # noqa: E402
import numpy as np # noqa: E402
mpmath.mp.dps = 40
T = 81
CHECKS = []
def check(name, ok, detail=""):
CHECKS.append(bool(ok))
print(f" [{'PASS' if ok else 'FAIL'}] {name}" +
(f" -- {detail}" if detail else ""), flush=True)
def inverse(d, eps=1e-16):
"""beta, c_0 and the forward impulse h_m, m = 0..K, by the roots."""
coef = [-1] * d + [1] # Delta, highest power first
roots = mpmath.polyroots(coef, maxsteps=400, extraprec=200)
dD = lambda z: -sum(k * z ** (k - 1) for k in range(1, d + 1))
c = [-1 / (r * dD(r)) for r in roots]
i0 = min(range(d), key=lambda i: abs(roots[i]))
beta = float(mpmath.re(1 / roots[i0]))
c0 = float(mpmath.re(c[i0]))
rest = [(complex(roots[i]), complex(c[i])) for i in range(d) if i != i0]
slow = min(abs(r) for r, _ in rest)
K = int(math.ceil(math.log(eps) / -math.log(slow))) + d
m = np.arange(K)
h = np.zeros(K)
for r, cj in rest:
h += (cj * r ** (-m.astype(float))).real
return beta, c0, h
def inverse_fft(d):
"""h by the real FFT of 1/Delta on N = 2^lg points: forward h_m
(m >= 0, cut where |h| last exceeds 1e-18), backward h_(-m) (m >= 1,
hb[0] unused), and the aliasing, max |h| near N/2."""
lg = max(16, math.ceil(math.log2(2 * d ** 3)) + 1)
N = 1 << lg
th = -2 * np.pi * np.arange(N // 2 + 1) / N
w = np.exp(1j * th)
wd = np.exp(1j * th * d)
with np.errstate(divide="ignore", invalid="ignore"):
D = 1 - w * (1 - wd) / (1 - w)
D[0] = 1 - d
del w, wd, th
h = np.fft.irfft(1 / D, n=N)
del D
alias = np.abs(h[N // 2 - 1000:N // 2 + 1000]).max()
big = np.flatnonzero(np.abs(h[:N // 2]) > 1e-18)
hf = h[:big[-1] + 1].copy()
hb = np.concatenate([[0.0], h[::-1][:200 * d].copy()])
return hf, hb, alias
def norms(d, hf, hb):
"""L, L_s and, over t < T, ||g^(t)||_1 and ||g^(t) - e_t||_1, from the
two halves of h, one row at a time; beta read off hb = -c_0 beta^(-m)."""
beta = hb[1] / hb[2]
n = len(hf) + T
hz = np.concatenate([hf, np.zeros(T)])
k = np.arange(n)
def shifted(t):
v = np.empty(n)
f = k >= t
v[f] = hz[k[f] - t]
v[~f] = hb[t - k[~f]]
return v
gam = -sum(shifted(s) for s in range(d)) / sum(beta ** -s
for s in range(d))
nx, ns = np.empty(T), np.empty(T)
for t in range(T):
g = shifted(t) + beta ** -t * gam
nx[t] = np.abs(g).sum()
g[t] -= 1
ns[t] = np.abs(g).sum()
L = np.abs(hf).sum() + np.abs(hb[1:]).sum()
return L, L - abs(hf[0]) + abs(hf[0] - 1), nx, ns
def by_roots(d):
"""norms' halves of h from the roots."""
beta, c0, h = inverse(d)
hb = np.concatenate([[0.0], -c0 * beta ** -np.arange(1.0, 200 * d)])
return h, hb
def gmat(d, n):
"""g^(t) for t < n as a matrix over k < K + n, by the roots."""
hf, hb = by_roots(d)
beta = hb[1] / hb[2]
m = len(hf) + n
hz = np.concatenate([hf, np.zeros(n)])
k = np.arange(m)
def shifted(t):
v = np.empty(m)
f = k >= t
v[f] = hz[k[f] - t]
v[~f] = hb[t - k[~f]]
return v
gam = -sum(shifted(s) for s in range(d)) / sum(beta ** -s
for s in range(d))
return np.array([shifted(t) + beta ** -t * gam for t in range(n)])
def wave(d):
P, a = d - 1, (d - 1) // 4
x = list(range(a + 1)) + list(range(a - 1, -a - 1, -1)) + \
list(range(-a + 1, 0)) + [0] * (P - 4 * a)
return x[:P]
def section_control():
print("CONTROL the wave from g; Delta(F) x = delta for a random delta",
flush=True)
for d in (9, 13, 21):
P = d - 1
G = gmat(d, 4 * P)
xw = wave(d)
dl = np.array([xw[k % P] - xw[(k + 1) % P] for k in range(G.shape[1])])
x = G @ dl
err = max(abs(x[t] - xw[t % P]) for t in range(4 * P))
check(f"d = {d}: the wave's x from g", err < 1e-9,
f"max |x - x_wave| {err:.2e}")
rng = np.random.default_rng(1)
for d in (5, 12, 30):
G = gmat(d, 200 + d)
dl = rng.integers(-1, 2, G.shape[1])
x = G @ dl
res = max(abs(x[t] - x[t + 1:t + d + 1].sum() - dl[t])
for t in range(200))
s0 = abs(x[:d].sum())
check(f"d = {d}: a random delta's x solves every equation",
res < 1e-9 and s0 < 1e-9,
f"max residual {res:.2e}, first window sum {s0:.2e}")
def row(d, L, Ls, nx, ns):
C = max(1, (d - 1) // 4)
hx, hs = nx.max(), ns.max()
tail = 1e-21 * d * L
proved = (math.floor(max(hx, L + tail)) <= C and
math.floor(max(hs, Ls + tail)) <= C)
print(f" d = {d:3}: C {C}, L {L:.6f}, L_s {Ls:.6f}, d(L - (d-1)/4) "
f"{d * (L - (d - 1) / 4):.4f}, H_x {hx:.6f} (t = "
f"{int(nx.argmax())}), H_s {hs:.6f}, margin "
f"{C + 1 - max(hx, hs):.4f}, "
f"{'PROVED' if proved else 'open'}", flush=True)
return proved, L <= hx + 1e-12 and Ls <= hs + 1e-12, hx, hs
def section_h():
print("H the roots, d = 3..40", flush=True)
proved, past = set(), []
for d in range(3, 41):
L, Ls, nx, ns = norms(d, *by_roots(d))
p, ok, *_ = row(d, L, Ls, nx, ns)
if p:
proved.add(d)
if not ok:
past.append(d)
want = {3, 4, 9, 10, 13, 14, 15, 17, 18, 19} | set(range(21, 41))
check("the sups past T within the sups read: L <= H_x, L_s <= H_s",
not past, f"exceptions {past}")
check("the bound proved at d = 3, 4, 9, 10, 13..15, 17..19, 21..40, "
"left at d = 5..8, 11, 12, 16, 20", proved == want,
f"proved {sorted(proved)}")
def section_h3():
print("H3 CONTROL the FFT against the roots", flush=True)
for d in (5, 12, 30, 40):
L0, Ls0, nx0, ns0 = norms(d, *by_roots(d))
hf, hb, al = inverse_fft(d)
L, Ls, nx, ns = norms(d, hf, hb)
e = max(abs(L - L0), abs(Ls - Ls0), abs(nx.max() - nx0.max()),
abs(ns.max() - ns0.max()))
check(f"d = {d}: L, L_s, H_x, H_s agree", e < 1e-9 and al < 1e-15,
f"largest difference {e:.2e}, aliasing {al:.1e}")
print("H3 the FFT, d = 21..100", flush=True)
bad, past, scaled = [], [], []
for d in range(21, 101):
hf, hb, al = inverse_fft(d)
L, Ls, nx, ns = norms(d, hf, hb)
p, ok, *_ = row(d, L, Ls, nx, ns)
if not p or al >= 1e-15:
bad.append(d)
if not ok:
past.append(d)
scaled.append(d * (L - (d - 1) / 4))
check("the sups past T within the sups read: L <= H_x, L_s <= H_s",
not past, f"exceptions {past}")
check("the bound proved at every d = 21..100, aliasing below 1e-15",
not bad, f"exceptions {bad}")
check("H3b: d (L - (d - 1)/4) falls at every step d = 21..100",
all(b < a for a, b in zip(scaled, scaled[1:])),
f"{scaled[0]:.4f} at 21, {scaled[-1]:.4f} at 100")
def main():
section_control()
section_h()
section_h3()
print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed")
sys.exit(0 if all(CHECKS) else 1)
if __name__ == "__main__":
main()