heat.py
Checks: READING: The real support theorem; READING: The heat flow pattern.
Cited on: READING.
To run it: save heat.py, with real_support.py beside it; needs mpmath, numpy.
"""
heat.py -- P E's limit from the first-order heat flow: the bulk of the
real support h, read by classes mod P = d - 1, is the derivative of the
periodic heat kernel, and its first-order correction moves the two
zeros by z(tau) classes; it predicts P E tends to one constant per
parity of P.
QUESTION Z. Does the first-order term of the heat-flow reading predict
P E(d) at large d, one constant per parity of P?
BY HAND.
(Z1) m >= 0: h_m = sum_k c_k x_k^m over the P roots of x^d (x - 2) = -1
inside the disc, c = (x - 1)/((d + 1) x - 2d); m < 0: h_m = -c_0
beta^m, c_0 the same formula at x = beta (partial fractions of
(1 - z)/(1 - 2z + z^(d+1)) at rho = 1/x, rho^d = 2 - x).
(Z2) x = e^lam, lam = i theta + mu, theta = 2 pi k/P: the root equation
is exactly e^(P mu) = 1/(1 - (1 - e^lam)^2), so P mu = lam^2 +
lam^3 + O(theta^4) = -theta^2 - i theta^3 + O(theta^4, theta^3/P),
and at m = qP + r, x^m = e^(i theta r) e^(m mu).
(Z3) c = -(i theta/P)(1 + (3/2) i theta + O(theta^2)), the denominator
being -P(1 - i theta + O(P^-2)). With i theta <-> d_x/P these give
(Z4); the expansion is a TRANSPLANT from k << P, the small-tau end
(k ~ P) where it may fail.
(Z4) In the bulk, x = r/P, tau = m/P^3,
h_m = -(1/P^2) [d_x + (3/(2P)) d_x^2 + (tau/P) d_x^4] K(x, tau),
K the periodic heat kernel, K_tau = K_xx.
(Z5) Each zero x0 (0, and 1/2) moves by z(tau) = -3/2 - tau K_xxxx/K_xx
classes, measured from r = 0 and r = P/2.
(Z6) P E -> C = 2 int_0^inf sum_(x0) |K_xx(x0, tau)| Phi(z(tau)) d tau,
Phi(z) the sum of |j - z| over the lattice points j between the
class sums' sign change and z: at x0 = 0 lattice Z, switch -1/2;
at x0 = 1/2 the same at P even, and at P odd lattice Z + 1/2 with
switch -1 (the class at -1/2 has H = -1/P^2 < 0).
PREDICTIONS AND KILL (fixed before the runs).
Z0 CONTROL, read first: the root sum equals real_support.inverse_fft's
h to 1e-12 on m in [-200, 4 d^3) at d = 21, 60, and reproduces
classes.py's E at d = 99, 100 to 1e-9.
Z1 LOCAL: at d = 100, 101 the measured zero of r -> h_(qP + r) near r =
0 lies within 0.15 classes of z(tau) at tau = 0.02, 0.05, 0.1, 0.2.
KILL of the reading: off by more than 0.5 class at two of the four
tau at both d.
Z2 THE CONSTANT: P E at d = 199..202, 399..402; per parity, C + a/P
fitted to two d. KILL: Z6's C off the fitted C by more than 0.05,
or the parity gap's sign wrong. Not blind: before these were
fixed, a fit C + a/P to classes.py's P E at d = 95..100 already put
the limits near 4.46 and 4.71 (a ~ 16), the gap near 0.25.
FINDINGS (copied from the runs' prints).
- Z0 holds: the root sum equals the FFT's h to 4.8e-16 at d = 21 and
9.4e-16 at d = 60 (m from -200 to where the FFT's h is cut at 1e-18,
32768 and 524286, short of 4 d^3 at both d, the cut where the FFT's
noise floor meets N/2); streamed E 0.047210793845
and 0.049140307484 against the FFT's ...852 and ...488 at d = 99, 100.
- Z1 holds at a window of +-12 classes (at +-6 the crossing at tau =
0.2, where z = 6.39, falls outside the window and reads nan): the
zero near r = 0 is off z(tau) by 0.001, 0.005, 0.033, 0.091 at tau =
0.02, 0.05, 0.1, 0.2, at d = 100 and 101 alike. Near r = P/2
(not frozen) measured -1.9709, 0.4140, 2.4818, 6.4882 against
-1.9349, 0.4095, 2.4471, 6.3949 at d = 100. The error grows with tau,
1.3% and 1.4% of z near r = 0 at tau = 0.1, 0.2.
- Z2 holds, far inside the kill: C = 4.4570 at P even (1.5211 from x0
= 0, 2.9359 from 1/2) and 4.7070 at P odd (1.5211 + 3.1859); C + a/P
fitted to d = 199/399 and 201/401 gives 4.4571 (a 16.68), to 200/400
and 202/402 4.7071 (a 15.68, 15.69). P E = 4.49904, 4.74645, 4.49883,
4.74625 at d = 399..402. E's shares by tau over [0, 1/16), [1/16,
1/8), [1/8, 1/4), [1/4, inf) are 0.3979 0.4317 0.1658 0.0045 at d =
400 against C's 0.3975 0.4322 0.1657 0.0045, and 0.4679 0.3690 0.1586
0.0046 at d = 399 against 0.4674 0.3694 0.1586 0.0046.
- The parity gap C_odd - C_even = 0.2500150. BY HAND: at x0 = 1/2,
g(tau) = K_xx > 0 with int g d tau = 1 and int tau g'(tau) d tau =
-1 (g' = K_xxxx), so int g (z - j) d tau = -(j + 1/2),
and Phi(z) = sum_(j in L) (z - j)^+ - sum_(j in L, j < s) (z - j)
makes C_(1/2)/2 = sum_(j in L) Psi(j), Psi(j) = F(j) - (-j - 1/2)^+,
F(j) = int g (z - j)^+ d tau, one Psi for both parities (the second
term vanishes at j = -1/2). Poisson: sum over Z minus sum over Z +
1/2 is 2 sum_(n odd) Psi^(n); the kink of slope 1 at -1/2 gives
-2 sum_(n odd) 1/(4 pi^2 n^2) = -1/8, F'' (the density of z under g)
the rest. So the gap is 2/8 = 1/4 - 4 sum_(n odd) F^(n), the second
term 1.5e-5 by the quadrature's total.
- S (after Z2, no kill): g^(t)_k = h_(k-t) - (beta^-t/B) h_(k+1) on k
>= 0, since sum_(s=1..d) h_(n-s) = h_n - [n = 0]; read to 1.6e-15 at
d = 13, 40, 100. So ||g^(t)||_1 <= (1 - eps) L + eps D_(t+1), eps =
beta^-t/B and D_s = sum_m |h_m - h_(m+s)|, and H_x <= L + sup_t eps
(D_(t+1) - L)^+. D_s/L at s = 1, 2, 4, 8, 16: 0.9448 1.3286 1.7557
1.7772 1.7530 (d = 13), 0.4386 0.6814 1.0473 1.5275 1.9514 (40),
0.2108 0.3436 0.5616 0.8974 1.3570 (100).
RUN RECORD: the bare run 22.4 s, peak commit 241.4 MB; d = 399..402
4.1 s each by the root sum with roots dropped below
|x|^m = 1e-20.
"""
import os
import sys
import time
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import numpy as np # noqa: E402
from real_support import inverse_fft # noqa: E402
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 roots(d):
"""The roots of x^d (x - 2) = -1 inside the disc, one per conjugate
pair (upper half plane) with weight 2, a real root weight 1; and beta."""
ks = np.arange(1, d // 2 + 1)
phi = 2 * np.pi * ks / d
for _ in range(60):
phi = (2 * np.pi * ks + np.arctan2(np.sin(phi), 2 - np.cos(phi))) / d
x = (5 - 4 * np.cos(phi)) ** (-0.5 / d) * np.exp(1j * phi)
for _ in range(30):
xd = x ** (d - 1)
f = xd * x * (x - 2) + 1
fp = xd * ((d + 1) * x - 2 * d)
x = x - f / fp
wt = np.full(len(x), 2.0)
if d % 2 == 0:
x[-1] = x[-1].real
wt[-1] = 1.0
b = 2.0
for _ in range(100):
b = b - (b ** d * (b - 2) + 1) / (b ** (d - 1) * ((d + 1) * b - 2 * d))
resid = np.abs(x ** d * (x - 2) + 1).max()
return x, wt, b, resid
def coef(x, d):
return (x - 1) / ((d + 1) * x - 2 * d)
def h_at(d, m, rt=None):
"""h at an array of m >= 0 by the root sum (all roots)."""
x, wt, b, _ = rt or roots(d)
c = coef(x, d) * wt
lx = np.log(x)
return (np.exp(np.outer(m, lx)) @ c).real
def stream_E(d, mmax=None):
"""L, E and E's share by tau band, streaming m over [-80 d, P^3),
roots dropped once |x|^m < 1e-20."""
P = d - 1
x, wt, b, resid = roots(d)
c = coef(x, d) * wt
lx = np.log(x)
rr = np.arange(P)
s = np.sign(((P - 1) / 2 - rr) / P - 1 / P ** 2)
mmax = mmax or P ** 3
L = E = 0.0
bands = np.zeros(4)
edges = np.array([1 / 16, 1 / 8, 1 / 4]) * P ** 3
# m < 0
mneg = -np.arange(1, 80 * d)
hneg = -coef(b, d) * b ** mneg.astype(float)
eneg = np.abs(hneg) - s[mneg % P] * hneg
L += np.abs(hneg).sum()
E += eneg.sum()
bands[0] += eneg.sum()
m0 = 0
order = np.argsort(-np.abs(x)) # slowest decay first
absx = np.abs(x)
while m0 < mmax:
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),
mmax - 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
m = m0 + j
e = np.abs(h) - s[m % P] * h
L += np.abs(h).sum()
E += e.sum()
idx = np.searchsorted(edges, m, side="right")
bands += np.bincount(idx, weights=e, minlength=4)
m0 += n
return L, E, bands / E, resid
# ---- the continuum: K_xx and K_xxxx at x0 = 0, 1/2 ----
def kernel(x0, tau):
"""K_xx, K_xxxx of the periodic heat kernel at x0, by images below
tau = 0.05 and by Fourier above."""
tau = np.asarray(tau, float)
k2 = np.zeros_like(tau)
k4 = np.zeros_like(tau)
lo = tau < 0.05
t = tau[lo]
for n in range(-6, 7):
y = x0 + n
G = (4 * np.pi * t) ** -0.5 * np.exp(-y * y / (4 * t))
k2[lo] += G * (y * y / (4 * t * t) - 1 / (2 * t))
k4[lo] += G * (y ** 4 / (16 * t ** 4) - 3 * y * y / (4 * t ** 3)
+ 3 / (4 * t * t))
t = tau[~lo]
for k in range(1, 40):
w = 2 * np.pi * k
g = 2 * np.cos(w * x0) * np.exp(-w * w * t)
k2[~lo] += -w ** 2 * g
k4[~lo] += w ** 4 * g
return k2, k4
def phi(z, a, s):
"""Sum of |j - z| over j in a + Z strictly between s and z."""
z = np.asarray(z, float)
out = np.zeros_like(z)
up = z > s
n1 = np.floor(s - a) + 1
n2 = np.ceil(z - a) - 1
N = np.maximum(n2 - n1 + 1, 0)
out = np.where(up, N * (z - a) - (n1 + n2) * N / 2, out)
n1d = np.floor(z - a) + 1
n2d = np.ceil(s - a) - 1
Nd = np.maximum(n2d - n1d + 1, 0)
out = np.where(~up, (n1d + n2d) * Nd / 2 - Nd * (z - a), out)
return out
def zshift(x0, tau):
k2, k4 = kernel(x0, tau)
return -1.5 - tau * k4 / k2, k2
def constant(Podd):
tau = np.exp(np.linspace(np.log(1e-5), np.log(2.0), 400001))
with np.errstate(divide="ignore", invalid="ignore"):
z0, k20 = zshift(0.0, tau)
zh, k2h = zshift(0.5, tau)
f0 = np.abs(k20) * phi(z0, 0.0, -0.5)
fh = np.abs(k2h) * phi(zh, 0.5 if Podd else 0.0,
-1.0 if Podd else -0.5)
f0 = np.where(k20 != 0, f0, 0.0) # K_xx underflows: no mass there
fh = np.where(k2h != 0, fh, 0.0)
i0 = 2 * np.trapezoid(f0, tau)
ih = 2 * np.trapezoid(fh, tau)
prof = []
for lo_, hi_ in ((0, 1 / 16), (1 / 16, 1 / 8), (1 / 8, 1 / 4),
(1 / 4, 3)):
sel = (tau >= lo_) & (tau < hi_)
prof.append(2 * np.trapezoid((f0 + fh)[sel], tau[sel]))
return i0 + ih, i0, ih, np.array(prof) / (i0 + ih)
def section_control():
print("Z0 CONTROL: the root sum against the FFT")
for d in (21, 60):
hf, hb, alias = inverse_fft(d)
n = min(len(hf), 4 * d ** 3)
m = np.arange(n)
rt = roots(d)
hr = np.concatenate([h_at(d, m[a:a + 20000], rt)
for a in range(0, n, 20000)])
b = rt[2]
mneg = np.arange(1, 201)
hnr = -coef(b, d) * b ** (-mneg.astype(float))
err = max(np.abs(hr - hf[:n]).max(), np.abs(hnr - hb[1:201]).max())
check(f"Z0 root sum = FFT h on m in [-200, {n}) at d = {d}",
err < 1e-12, f"max error {err:.1e}, root residual "
f"{rt[3]:.1e}")
for d, Eref in ((99, 4.721e-02), (100, 4.914e-02)):
L, E, bands, resid = stream_E(d)
hf, hb, alias = inverse_fft(d)
P = d - 1
m = np.concatenate([-np.arange(len(hb) - 1, 0, -1),
np.arange(len(hf))])
h = np.concatenate([hb[:0:-1], hf])
H = np.bincount(m % P, weights=h, minlength=P)
Ef = (np.abs(h) - np.sign(H)[m % P] * h).sum()
del hf, hb, h, m
check(f"Z0 streamed E = the FFT's E at d = {d}", abs(E - Ef) < 1e-9,
f"E {E:.12f} vs {Ef:.12f}, L {L:.9f}")
def section_local():
print("Z1 LOCAL: the zero near r = 0 (and r = P/2) against z(tau)")
taus = (0.02, 0.05, 0.1, 0.2)
off = {}
for d in (100, 101):
P = d - 1
rt = roots(d)
off[d] = []
for tau in taus:
q = round(tau * P ** 2)
row = []
for centre in (0, P / 2):
js = np.arange(-12, 13) + (0 if centre == 0 else
(0.5 if P % 2 else 0))
r = np.round(centre + js).astype(int) % P
m = q * P + r
h = h_at(d, m, rt)
sgn = 1 if centre == 0 else -1 # slope sign at the zero
k = np.flatnonzero(np.diff(np.sign(sgn * h)) > 0)
zm = (js[k] - h[k] * (js[k + 1] - js[k]) /
(h[k + 1] - h[k])) if len(k) == 1 else [np.nan]
zp, _ = zshift(0.0 if centre == 0 else 0.5,
np.array([q * P / P ** 3]))
row.append((float(np.ravel(zm)[0]), float(zp[0])))
off[d].append(abs(row[0][0] - row[0][1]))
print(f" d = {d}, tau = {tau}: near 0 measured {row[0][0]:+.4f}"
f" predicted {row[0][1]:+.4f}; near P/2 measured "
f"{row[1][0]:+.4f} predicted {row[1][1]:+.4f}", flush=True)
check("Z1 the zero near r = 0 within 0.15 classes, 4 tau, d = 100, 101",
all(o < 0.15 for d in off for o in off[d]),
" ".join(f"{o:.3f}" for d in off for o in off[d]))
killed = all(sum(o > 0.5 for o in off[d]) >= 2 for d in off)
check("Z1 KILL not fired (off by 0.5 at two tau at both d)", not killed)
def section_constant():
print("Z2 THE CONSTANT")
Ce, i0e, ihe, prof_e = constant(False)
Co, i0o, iho, prof_o = constant(True)
print(f" C (P even, d odd) = {Ce:.4f} = {i0e:.4f} (x0 = 0) + {ihe:.4f}"
f" (x0 = 1/2); profile {' '.join(f'{v:.4f}' for v in prof_e)}")
print(f" C (P odd, d even) = {Co:.4f} = {i0o:.4f} (x0 = 0) + {iho:.4f}"
f" (x0 = 1/2); profile {' '.join(f'{v:.4f}' for v in prof_o)}")
PE = {}
for d in (199, 200, 201, 202, 399, 400, 401, 402):
t0 = time.time()
L, E, bands, resid = stream_E(d)
PE[d] = (d - 1) * E
print(f" d = {d}: L {L:.9f} E {E:.6e} P E {PE[d]:.5f} "
f"profile {' '.join(f'{v:.4f}' for v in bands)} "
f"({time.time() - t0:.1f} s)", flush=True)
fits = {}
for par, pairs in (("even", ((199, 399), (201, 401))),
("odd", ((200, 400), (202, 402)))):
fits[par] = []
for d1, d2 in pairs:
P1, P2 = d1 - 1, d2 - 1
a = (PE[d1] - PE[d2]) / (1 / P1 - 1 / P2)
fits[par].append((PE[d2] - a / P2, a))
print(f" P {par}: C + a/P fits " + "; ".join(
f"C {C:.4f}, a {a:.2f}" for C, a in fits[par]))
ok_e = all(abs(C - Ce) <= 0.05 for C, _ in fits["even"])
ok_o = all(abs(C - Co) <= 0.05 for C, _ in fits["odd"])
check("Z2 C (P even) within 0.05 of the fits", ok_e,
f"{Ce:.4f} vs {[round(C, 4) for C, _ in fits['even']]}")
check("Z2 C (P odd) within 0.05 of the fits", ok_o,
f"{Co:.4f} vs {[round(C, 4) for C, _ in fits['odd']]}")
check("Z2 the parity gap's sign", (Co - Ce) * (PE[400] - PE[399]) > 0,
f"C gap {Co - Ce:+.7f}, P E gap at 399/400 "
f"{PE[400] - PE[399]:+.4f}")
def section_sup():
"""S (after Z2, no kill): g^(t)_k = h_(k-t) - (beta^-t/B) h_(k+1) on k
>= 0, the window sum of h being h_(k+1) (Delta(F) h = delta); and the
shift norms D_s = sum_m |h_m - h_(m+s)| against L."""
print("S THE SUP OVER t: the window sum of h is one term")
for d in (13, 40, 100):
hf, hb, alias = inverse_fft(d)
h = np.concatenate([hb[:0:-1], hf])
off, n = len(hb) - 1, len(hf)
x, wt, b, _ = roots(d)
B = sum(b ** -s for s in range(d))
k = np.arange(n - 2)
win = sum(h[off + k - s] for s in range(d))
err = np.abs(win - h[off + k + 1]).max()
L = np.abs(h).sum()
D = [np.abs(h[:-s] - h[s:]).sum() / L for s in (1, 2, 4, 8, 16)]
print(f" d = {d}: identity error {err:.1e}, L {L:.4f}, D_s/L at "
f"s = 1, 2, 4, 8, 16: " + " ".join(f"{v:.4f}" for v in D))
del hf, hb, h, win
check(f"S the window sum of h is h_(k+1) at d = {d}", err < 1e-13,
f"{err:.1e}")
def main():
section_control()
if not all(CHECKS):
print("control failed; stopping")
sys.exit(1)
section_local()
section_constant()
section_sup()
print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed")
sys.exit(0 if all(CHECKS) else 1)
if __name__ == "__main__":
main()