""" 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()