""" majorant.py -- majorants of |h - phi| per image and per mode that hold for every P >= P0 at once, and the remainder R they leave at d = 400. QUESTION W. Do majorants of the per-term errors that hold for every P >= P0 at once (images: |T_k - phi_k| <= |phi_k| A(xi, tau)/P; modes: |root_k - mode_k| <= env_k(tau) A_k(tau)/P), put in place of the exact per-term errors in remainder.py's route, keep R_maj + F below 1/4 at d = 400, P0 = 399? Notation as remainder.py; a = |xi|/(2 tau), w = xi^2/(4 tau) = a^2 tau, s = sgn xi. BY HAND. (W1) Images. Stirling with Robbins' bracket (needs 1 <= u <= n - 1): T/phi_k = (n0/n)^(3/2) (1 - y^2)^(-1/2) exp(w - n D(y)) e^eps, eps in (-1/(3n(1 - y^2)), 0]. With delta = xi(P + 2)/n0 = s a (1 + 2/P)/P, n = n0(1 - delta), y = j/n = s a/(P(1 - delta)), n y^4 = 2 a^4 tau/(P^2 (1 - delta)^3), exactly L = ln(T/phi_k) = delta(3/2 - w) + E2, E2 = (3/2)(-ln(1 - delta) - delta) - (1/2) ln(1 - y^2) - w delta^2/(1 - delta) - R_D + eps, R_D = n sum_(i>=2) y^(2i)/(2i(2i - 1)) in [n y^4/12, n y^4/(12(1 - y^2))], the first two terms in [0, (3/4)delta^2/(1 - |delta|)] and [0, y^2/(2(1 - y^2))]. First order: P L -> l = s a (3/2 - w) = 3xi/(4 tau) - xi^3/(8 tau^2), re-derived as [(3/2)d_x^2 + tau d_x^4]G / d_x G on the Gaussian image (He_2, He_4: the 1/xi cancels), heat.py's correction, so not a TRANSPLANT. Core: a <= P/4, so |delta| <= dbar = (1 + 2/P0) min(a/P0, 1/4), |y| <= ybar = min(a/P0, 1/4)/(1 - dbar), and u, n - u >= n(1 - ybar)/2 >= 1 at tau >= 1/c. For every P >= P0, P L in [V, U]: Hbar = (3/4)a^2(1 + 2/P0)^2/(1 - dbar) + a^2/(2(1 - dbar)^2(1 - ybar^2)) - a^4 tau/(1 + dbar) - a^4 tau/(6(1 + dbar)^3), Lbar = -a^4 tau (1 + 2/P0)^2/(1 - dbar) - a^4 tau/(6(1 - dbar)^3 (1 - ybar^2)) - 1/(6 tau (1 - dbar)(1 - ybar^2)), U = l + max(Hbar + max(2l, 0), 0)/P0, V = l + min(Lbar + min(2l, 0), 0)/P0, and A = max(f(U), f(V)), f(t) = P0(e^(t/P0) - 1) at t >= 0 (P(e^(t/P) - 1) falls with P), |t| at t < 0 (P(1 - e^(t/P)) rises to |t|). Outside the core, |T| <= (1/2) exp(-w/(1 + |delta|)) (Hoeffding, |y| <= 1), w >= P^2 tau/16 there; the images past K and phi's image tail as in remainder.py, both uniform (g(KP) rises with P). (W2) Modes. 2Re(c x^m) - mode_k = 2Re(c0 e^(mu0 m) e^(2 pi i k x)((c/c0) e^((mu - mu0) m) - 1)), c0 = -i theta/P, mu0 = -theta^2/P, theta = 2 pi k/P, |2 c0| e^(mu0 m) = env_k. P mu = g(lambda), g = -lambda - Log(2 - e^lambda) = lambda^2 + lambda^3 + G, G = sum_(n>=4) g_n lambda^n, g_n exact to n = 40 (Fractions), the rest by Cauchy on |lambda| = 0.6 with M = 0.6 - ln(2 - e^0.6) = 2.3266: Ghat(r). The disc |mu| <= mubar = sigma theta^2/P maps into itself when (1 + sigma theta/P)^2 (1 + theta + sigma theta^2/P) + Ghat(theta + mubar)/ theta^2 <= sigma, |g'| <= 11.2 < P0 makes it a contraction; sigma_k read at P0, theta0 = 2 pi k/P0, the left side falling in P. Then (mu - mu0) m = i alpha1 + z2, alpha1 = -(2 pi k)^3 tau/P, |z2| <= zeta = rhobar tau P^2, rhobar = 2 theta mubar + mubar^2 + 3 theta^2 mubar + 3 theta mubar^2 + mubar^3 + Ghat(theta + mubar). c/c0 = 1 + 1.5 i theta + e, |e| <= theta^2/2 + |E_N|(1 + theta + |E_D'|) + (1 + theta/2) |E_D'|, |E_v| <= theta^3/6 + e^mubar - 1, |E_N| <= |E_v|/theta, |v| <= e^(theta + mubar) - 1, q = (1 + 2/P)|v|, |E_D'| <= (2/P)|v| + theta^2/2 + |E_v| + q^2/(1 - q). With w1 = (1 + 1.5 i theta)e^(i alpha1), beta1 = 3 pi k, beta2 = (2 pi k)^3 tau, P^2 |w1 - 1|^2 <= (beta1 - beta2)^2 + beta1 beta2^3/(3 P0^2) (1 - cos x <= x^2/2, sin x >= x - x^3/6), so A_k = sqrt((beta1 - beta2)^2 + beta1 beta2^3/(3 P0^2)) + E_c e^zeta0 + sqrt(1 + 2.25 theta0^2) P0 (e^zeta0 - 1), E_c = P0 |e| and zeta0 = zeta at P0 (each P-scaled term falls with P). First order beta1 - beta2 matches heat.py's (3/(2P)) d_x^2 + (tau/P) d_x^4 on mode k; the triangle inequality on c's and the drift's errors separately would lose |beta1| + |beta2|. Modes k <= K4, theta0 <= 0.35; past K4 every root x (Im x >= 0, arg phi, d phi = 2 pi k + arg(1/(2 - x)), that arg in [0, pi/6]) has |x|^(P+1)|2 - x| = 1 and |2 - x|^2 >= Q(phi) = min(5 - 4 cos phi, 3), so |x|^m <= Q^(-m/(2(P + 1))), phi >= 2 pi k/d, and |c| <= 2/P. Owed, not proved here: that root index k <= K4 is the disc's fixed point (one root per sector); read as the solver's roots lying in the disc. Past tau = 1 the same Q-bound summed geometrically in m. (W3) R_maj = remainder.py's sum with B_maj = min(B_img_maj (tau < 1/4), B_four_maj (tau >= 1/64)); B_maj >= B_prov term by term wherever the majorants dominate, so R_maj >= R_prov. (False as written: B_prov carries floors, the float |sum phi_k - phi| and the dropped roots' 1e-20 sum, that the majorants' analytic tails undercut; FINDINGS.) PREDICTIONS, CONTROL AND KILL. W0 CONTROL, read first, at d = 400, P0 = 399, on 4000 m sampled in [P^3/400, P^3/4) and 2000 in [P^3/64, P^3): (a) every core image's |T - phi_k| <= |phi_k| A/P (1e-12 relative slack) and every mode k <= K4's |root - mode| <= env A_k/P; (b) the PLANT: with A and A_k halved, (a) fails on some sample of each kind; (c) tightness: the largest P|T/phi - 1|/A over images with w <= 10, and P|root - mode|/(env A_k) over k <= 3, are both at least 0.9; (d) the solver's roots k <= K4 lie in the disc |mu - 0| <= mubar. W1 Uniformity: A and A_k read at P0 = 399 dominate the exact relative errors at d = 800 and 1600 on the same tau samples. W2 R_maj/R_prov in [1, 1.10] at d = 400, P0 = 399; B_maj >= B_prov (1e-12) at every m. At d = 200 (P0 = 199) printed, not predicted. AMENDED BEFORE CODE (W1, W2): the core a <= P/4 holds the images' scale sqrt(30/tau) only when tau >> 30/P^2, so at tau0 = 1/400 not below P near 110: d = 100 leaves W2 and d = 200 is read only. Outside the core |T| <= (2 pi n (1 - |y|))^(-1/2) |y| exp(-j^2/(2n)) (Robbins, eps <= 0, D(y) >= y^2/2, at 1 <= u <= n - 1; u = n gives 2^-n/2), P^-2-scaled, not Hoeffding's (1/2)e^(-w/(1 + |delta|)), which is not on phi_k's P^-2 scale. The frozen band "P R_maj in [20, 30] at each d" left with d = 100 and 200, unsaid in the amendment (at d = 400 it holds). AFTER THE FIRST RUN: W1 draws its own tau samples (seed 11, 1500 and 600), not W0's; the mode samples are log-uniform in tau over [1/64, 1), uniform before, which reached tau near 0.02 too rarely to catch a dropped beta1 (W0e, the planted break, added then). KILL: R_maj + F(400, 399) >= 1/4 at d = 400; then a P-uniform route needs the cancellation between images or modes, and the per-term route is closed. FINDINGS (copied from the run's prints). - The control holds. g = lambda^2 + lambda^3 + (13/12)lambda^4 + ..., M = 2.3266, ghat(0.35) = 2.775e-2. Over 65424 core image terms (4000 m) and 25075 mode terms (2000 m, K4 = 22) at d = 400, the exact (30 digit) errors stay within the majorants, largest ratios 0.999074 and 0.976082; the largest ratio at w <= 10 is 0.9991 and at k <= 3 0.9761, so both are tight to first order. W0b reads the same ratios past 0.5 (44084 and 4233 terms), W0c's corollary, not a plant. The plant is W0e's break: with the bracket's first order l's sign flipped the largest image ratio is 7.4e17, and with beta1 = 0 in A_k the largest mode ratio 1.2899, W0a failing on both. The solver's roots k <= 22 lie in their discs (|mu|/mubar at most 0.9840). sigma runs 1.0161 (k = 1) to 1.5723 (k = 22), E_c 0.269 to 210.5. - W1 holds: the constants read at P0 = 399 dominate at d = 800 and 1600, largest ratios 0.998047, 0.998425 (images, about 24800 terms each) and 0.973942, 0.974154 (modes, about 7550 each). - The KILL does not fire: R_maj + F = 0.06351 at d = 400 (P R 25.34), a quarter of the 1/4 margin, against R_prov 0.04988 and the oracle's 0.04979; by band [1/c, 1/16), [1/16, 1/4), [1/4, 1): 0.04496, 0.01836, 1.970e-4 against 0.03264, 0.01706, 1.836e-4. B_maj sits above |h - phi| at every m but by float, 2.3e-16 at most. - W2 fails twice, both on its own premises. The ratio is 1.2734, the excess in the images' band, and W4 reads it as the price of one A for every P >= P0: there the exact per-term R is 0.03264, the uniform A's 0.04496, its clamps kept with f at P0 0.03906, the bracket at P0 alone 0.03351, the first order |l| (no bound) 0.03034. So, in that order, f(t) = |t| at t < 0, the sup over P, costs 0.00590 (the clamps kept) and the clamps of U and V at 0 cost 0.00555 (f at P0); together 0.01145 of the band's 0.01232. A is fixed at P0, so neither price leaves as P grows: the exact relative error tends to |l|/P while A/P keeps the bracket's 1/P0 terms, and P R_maj need not fall toward the oracle's. B_maj >= B_prov fails (0.94 relative) because B_prov's floors are not the majorants' tails: at m = 63521018 (tau 0.999997) B_prov is 5.874e-21, of it 5.538e-21 the dropped roots' crude floor, the roots' per-term sum 3.362e-22, and B_maj 3.472e-22 sits above that sum. - At d = 200 (P0 = 199) R_maj is 0.22683 against R_prov 0.12158, the excess in [1/c, 1/16) (0.18628 against 0.08681); the core a <= P/4 cutting the images' scale and the bracket's 1/P0 terms at P0 = 199 are not separated there. RUN RECORD: 15/17 checks, the misses both W2's, as recorded; 320.1 s under the 512 MB guard, peak commit 146.9 MB. W4 (read only, no check), W0e and W1's own samples were added after the first run. Three faults were fixed before any verdict: a core term with phi_k underflowed and A overflowed gave a NaN that max() hid (now the smaller of the majorant and the Stirling bound, and B asserted finite); heat.roots overflows on beta at d = 1600 (the roots k <= K4 now start from its angle formula); the modes' e^zeta0 overflowed where env underflows (now in logs). """ import math import os import sys import time from fractions import Fraction os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") import numpy as np # noqa: E402 import mpmath as mp # noqa: E402 from heat import roots, coef # noqa: E402 from edge import signs, kx # noqa: E402 from remainder import images, modes, F_edge # noqa: E402 CHECKS, NAMES = [], [] MISSED = ("W2 R_maj/R_prov in [1, 1.10] at d = 400", "W2 B_maj >= B_prov at every m (1e-12)") # the record's two misses C = 400 R0, NG = 0.6, 40 # Cauchy disc for g, exact coefficients M_G = R0 - math.log(2 - math.exp(R0)) mp.mp.dps = 30 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_coeffs(N): """Taylor coefficients of g = -lambda - Log(2 - e^lambda), exact.""" v = [Fraction(0)] + [Fraction(1, math.factorial(n)) for n in range(1, N + 1)] g = [Fraction(0)] * (N + 1) g[1] -= 1 pw = [Fraction(1)] + [Fraction(0)] * N for j in range(1, N + 1): new = [Fraction(0)] * (N + 1) for a in range(N + 1): if pw[a]: for b in range(1, N + 1 - a): new[a + b] += pw[a] * v[b] pw = new for n in range(N + 1): g[n] += pw[n] / j return g GC = g_coeffs(NG) GABS = [abs(float(c)) for c in GC] def ghat(r): """W2: |G(lambda)| <= ghat(|lambda|), G = g - lambda^2 - lambda^3.""" assert r < R0 s = sum(GABS[n] * r ** n for n in range(4, NG + 1)) return s + M_G * (r / R0) ** (NG + 1) / (1 - r / R0) def image_A(a, tau, s, P0, arm="uniform"): """W1: P |T/phi_k - 1| <= A for every P >= P0 with a <= P/4; the arm "at P0" reads the same bracket at P = P0 alone, unclamped, and "clamped" keeps the clamps of U and V with f at P0 on both signs.""" e0 = 1 + 2 / P0 rr = np.minimum(a / P0, 0.25) db = e0 * rr yb = rr / (1 - db) w = a * a * tau lo = s * a * (1.5 - w) a4t = a ** 4 * tau Hb = (0.75 * a * a * e0 ** 2 / (1 - db) + a * a / (2 * (1 - db) ** 2 * (1 - yb ** 2)) - a4t / (1 + db) - a4t / (6 * (1 + db) ** 3)) Lb = (-a4t * e0 ** 2 / (1 - db) - a4t / (6 * (1 - db) ** 3 * (1 - yb ** 2)) - 1 / (6 * tau * (1 - db) * (1 - yb ** 2))) if arm == "at P0": U = lo + (Hb + np.maximum(2 * lo, 0)) / P0 V = lo + (Lb + np.minimum(2 * lo, 0)) / P0 if arm == "clamped": U = lo + np.maximum(Hb + np.maximum(2 * lo, 0), 0) / P0 V = lo + np.minimum(Lb + np.minimum(2 * lo, 0), 0) / P0 if arm != "uniform": return P0 * np.maximum(np.abs(np.expm1(np.minimum(U / P0, 700))), np.abs(np.expm1(np.minimum(V / P0, 700)))) U = lo + np.maximum(Hb + np.maximum(2 * lo, 0), 0) / P0 V = lo + np.minimum(Lb + np.minimum(2 * lo, 0), 0) / P0 def f(t): return np.where(t >= 0, P0 * np.expm1(np.minimum(np.maximum(t, 0) / P0, 700)), -t) return np.maximum(f(U), f(V)) def images_maj(m, P, P0, tau_hi, arm="uniform"): """W1: B_img_maj at each m: core images |phi_k| A/P, the rest |T|'s Stirling bound + |phi_k|, the walk past K and phi's images past K.""" 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 Bw = 1 + 2 / P J = K * P g = J * J / (2 * (n0max + Bw * J)) dg = J * (2 * n0max + Bw * J) / (2 * (n0max + Bw * J) ** 2) assert g >= 60, g pref = (4 * np.pi * tau) ** -0.5 / P ** 2 Bm = np.full(len(m), math.exp(-g) / (-math.expm1(-P * dg))) for N in range(K, K + 40): Bm += 2 * pref * (N + 1) / (2 * tau) * np.exp(-N * N / (4 * tau)) for k in range(-K, K + 1): xi = x + k ph = pref * np.exp(-xi * xi / (4 * tau)) * xi / (2 * tau) a = np.abs(xi) / (2 * tau) core = a <= P / 4 j = r + k * P u = q - k n = 2 * u - j ok = (u >= 0) & (n >= np.maximum(u, 1)) assert (ok[core] & (u[core] >= 1) & (u[core] <= n[core] - 1)).all() if arm == "first": A = np.abs(np.sign(xi) * a * (1.5 - a * a * tau)) else: A = image_A(a, tau, np.sign(xi), P0, arm) nf = np.where(ok, n, 1).astype(float) jf = j.astype(float) y = np.abs(jf) / nf inner = ok & (u >= 1) & (u < n) Tb = np.where(inner, (2 * np.pi * nf * np.maximum(1 - y, 1e-300)) ** -0.5 * y * np.exp(-jf * jf / (2 * nf)), 0.0) Tb = np.where(ok & ((u == n) | (u == 0)), 0.5 * 2.0 ** (-nf), Tb) crude = Tb + np.abs(ph) maj = np.abs(ph) * np.minimum(A, 1e300) / P Bm += np.where(core, np.minimum(maj, crude), crude) return Bm def mode_consts(P0): """W2: per mode k <= K4, (k, theta0, sigma, mubar0, E_c, zeta0/tau).""" K4 = int(0.35 * P0 / (2 * math.pi)) out = [] for k in range(1, K4 + 1): th = 2 * math.pi * k / P0 def lhs(sg): return ((1 + sg * th / P0) ** 2 * (1 + th + sg * th * th / P0) + ghat(th + sg * th * th / P0) / th ** 2) sg = 1.0 for _ in range(50): sg = lhs(sg) sg *= 1 + 1e-9 assert lhs(sg) <= sg mb = sg * th * th / P0 rho = (2 * th * mb + mb * mb + 3 * th * th * mb + 3 * th * mb * mb + mb ** 3 + ghat(th + mb)) Ev = th ** 3 / 6 + math.expm1(mb) EN = Ev / th v = math.expm1(th + mb) qq = (1 + 2 / P0) * v assert qq < 1 ED = 2 / P0 * v + th * th / 2 + Ev + qq * qq / (1 - qq) ec = th * th / 2 + EN * (1 + th + ED) + (1 + th / 2) * ED out.append((k, th, sg, mb, P0 * ec, rho * P0 * P0)) return out def mode_A(k, th0, Ec, zc, tau, P0): b1 = 3 * math.pi * k b2 = (2 * math.pi * k) ** 3 * tau z0 = zc * tau return (np.sqrt((b1 - b2) ** 2 + b1 * b2 ** 3 / (3 * P0 * P0)) + Ec * np.exp(z0) + math.sqrt(1 + 2.25 * th0 * th0) * P0 * np.expm1(z0)) def mode_tails(d, P, K4, m0): """Phi's modes past K4 and every root past K4, at the chunk's least m.""" tau0 = m0 / P ** 3 ks = np.arange(K4 + 1, K4 + 400) four = (4 * np.pi * ks * np.exp(-4 * np.pi ** 2 * ks * ks * tau0)).sum() \ / P ** 2 kr = np.arange(K4 + 1, d // 2 + 1) lq = np.log(np.minimum(5 - 4 * np.cos(2 * np.pi * kr / d), 3.0)) root = (4 / P * np.exp(-m0 * lq / (2 * (P + 1)))).sum() return four + root def modes_maj(tau, d, P, P0, cst, m0): """env_k A_k/P summed, in logs: past zeta0 = 600, ln A_k <= ln 2 + max(ln of the root term, ln(E_c + sqrt(1 + 2.25 theta0^2) P0) + zeta0).""" B = np.zeros(len(tau)) for k, th0, sg, mb, Ec, zc in cst: lnenv = math.log(4 * math.pi * k / P ** 2) - 4 * np.pi ** 2 * k * k * tau z0 = zc * tau b1, b2 = 3 * math.pi * k, (2 * math.pi * k) ** 3 * tau S = np.sqrt((b1 - b2) ** 2 + b1 * b2 ** 3 / (3 * P0 * P0)) cz = Ec + math.sqrt(1 + 2.25 * th0 * th0) * P0 small = z0 < 600 lnA = np.log(2) + np.maximum(np.log(S), math.log(cz) + z0) if small.any(): lnA[small] = np.log(mode_A(k, th0, Ec, zc, tau[small], P0)) B += np.exp(lnenv + lnA) / P return B + mode_tails(d, P, len(cst), m0) def past_M(d, P): """The sum over m >= P^3 of 2|h_m|, every root by the Q-bound.""" kr = np.arange(1, d // 2 + 1) lq = np.log(np.minimum(5 - 4 * np.cos(2 * np.pi * kr / d), 3.0)) lr = -lq / (2 * (P + 1)) return float((2 * 4 / P * np.exp(P ** 3 * lr) / -np.expm1(lr)).sum()) def stream(d, P0, cst): """R_prov (remainder.py's B), the oracle and R_maj, by tau band.""" 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", "maj")} worst_pm, worst_hm = 0.0, -np.inf 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] Bi = np.full(n, np.inf) Bf = np.full(n, np.inf) Bim = np.full(n, np.inf) Bfm = np.full(n, np.inf) if m0 / P ** 3 < 1 / 4: Ts, Ps, D, tailT = images(m, P, thi) Bi = D + np.abs(Ps - ph) + tailT Bim = images_maj(m, P, P0, thi) if thi >= 1 / 64: D4, S4 = modes(m, P, xs, cw, act + 1) Bf = D4 + np.abs(S4 - ph) + drop Bfm = modes_maj(tau, d, P, P0, cst, m0) Bp = np.minimum(Bi, Bf) Bm = np.minimum(Bim, Bfm) assert np.isfinite(Bp).all() and np.isfinite(Bm).all() err = np.abs(h - ph) worst_pm = max(worst_pm, float(((Bp - Bm) / Bp).max())) worst_hm = max(worst_hm, float((err - Bm).max())) band = np.searchsorted(edges, tau, side="right") for key, B in (("prov", Bp), ("orac", err), ("maj", Bm)): bad = sm * ph <= B rr = np.where(bad, np.abs(ph) - sm * ph + 2 * B, 0.0) out[key] += np.bincount(band, weights=rr, minlength=3) m0 += n out["tail"] = 2 * (np.abs(cw) * absx ** M / (1 - absx)).sum() out["tail_maj"] = past_M(d, P) out["worst_pm"] = worst_pm out["worst_hm"] = worst_hm return out def exact_T(m, k, P): """T_(q-k)(m) and phi_k in mpmath, or None outside 1 <= u <= n - 1.""" q, r = divmod(int(m), P) j, u = r + k * P, q - k n = 2 * u - j if not (1 <= u <= n - 1): return None T = mp.binomial(n, u) * mp.mpf(2) ** (-n) * j / (2 * n) tau = mp.mpf(int(m)) / P ** 3 xi = mp.mpf(j) / P ph = (4 * mp.pi * tau) ** -0.5 / P ** 2 * mp.exp(-xi * xi / (4 * tau)) \ * xi / (2 * tau) return T, ph def exact_roots(d, K4): """The roots k <= K4 refined to 30 digits from heat.py's start, with their c (heat.roots itself overflows on beta past d near 1000).""" out = [] for k in range(1, K4 + 1): ph = 2 * math.pi * k / d for _ in range(60): ph = (2 * math.pi * k + math.atan2(math.sin(ph), 2 - math.cos(ph))) / d x = mp.mpc(complex((5 - 4 * math.cos(ph)) ** (-0.5 / d) * np.exp(1j * ph))) for _ in range(40): f = x ** (d + 1) - 2 * x ** d + 1 fp = (d + 1) * x ** d - 2 * d * x ** (d - 1) x = x - f / fp out.append((x, (x - 1) / ((d + 1) * x - 2 * d))) return out def image_ratios(d, P0, ms): """Over the core images of each m: P |T - phi_k| / (|phi_k| A), and w.""" P = d - 1 rat, ws = [], [] for mm in ms: tau = mm / P ** 3 K = int(math.ceil(math.sqrt(4 * tau * 64))) + 2 x = (mm % P) / P for k in range(-K, K + 1): xi = x + k a = abs(xi) / (2 * tau) if a > P / 4 or xi == 0: continue e = exact_T(mm, k, P) if e is None: raise AssertionError("core image outside 1 <= u <= n - 1") T, ph = e if abs(ph) < mp.mpf(10) ** -280: continue A = float(image_A(np.array([a]), tau, np.sign(xi), P0)[0]) rat.append(float(P * abs(T - ph) / (abs(ph) * A))) ws.append(xi * xi / (4 * tau)) return np.array(rat), np.array(ws) def mode_ratios(d, P0, cst, ms): """Over k <= K4 and each m: P |root_k - mode_k| / (env_k A_k), and k.""" P = d - 1 rts = exact_roots(d, len(cst)) rat, kk = [], [] for mm in ms: tau = mm / P ** 3 x = mp.mpf(int(mm % P)) / P for (k, th0, sg, mb, Ec, zc), (xr, c) in zip(cst, rts): env = 4 * mp.pi * k * mp.exp(-4 * mp.pi ** 2 * k * k * mp.mpf(tau)) / P ** 2 if env < mp.mpf(10) ** -280: continue root = 2 * mp.re(c * xr ** int(mm)) mode = env * mp.sin(2 * mp.pi * k * x) A = float(mode_A(k, th0, Ec, zc, tau, P0)) rat.append(float(P * abs(root - mode) / (env * A))) kk.append(k) return np.array(rat), np.array(kk) def samples(P, n_img, n_mode, seed): rng = np.random.default_rng(seed) t1 = rng.uniform(1 / 400, 1 / 4, n_img) t2 = np.exp(rng.uniform(math.log(1 / 64), 0, n_mode)) mi = np.unique(np.round(t1 * P ** 3).astype(np.int64)) mo = np.unique(np.round(t2 * P ** 3).astype(np.int64)) return mi, mo def section_control(): print("W0 control at d = 400, P0 = 399", flush=True) d, P0 = 400, 399 print(f" g_2 = {GC[2]}, g_3 = {GC[3]}, g_4 = {GC[4]}, M = {M_G:.4f}, " f"ghat(0.35) = {ghat(0.35):.3e}", flush=True) check("W0 g's series starts lambda^2 + lambda^3", GC[1] == 0 and GC[2] == 1 and GC[3] == 1, f"g_4 = {GC[4]}") cst = mode_consts(P0) mi, mo = samples(d - 1, 4000, 2000, 7) t0 = time.time() ri, wi = image_ratios(d, P0, mi) rm, km = mode_ratios(d, P0, cst, mo) print(f" {len(ri)} core image terms over {len(mi)} m, {len(rm)} mode " f"terms over {len(mo)} m, K4 = {len(cst)} ({time.time() - t0:.0f} s)", flush=True) check("W0a every core image within |phi_k| A/P (1e-12)", ri.max() <= 1 + 1e-12, f"largest ratio {ri.max():.6f}") check("W0a every mode k <= K4 within env A_k/P (1e-12)", rm.max() <= 1 + 1e-12, f"largest ratio {rm.max():.6f}") hi = (ri > 0.5).sum() hm = (rm > 0.5).sum() check("W0b (W0c's corollary) some image ratio past 1/2", hi > 0, f"{hi} terms") check("W0b (W0c's corollary) some mode ratio past 1/2", hm > 0, f"{hm} terms") ti = ri[wi <= 10].max() tm = rm[km <= 3].max() check("W0c tight: largest image ratio at w <= 10 at least 0.9", ti >= 0.9, f"{ti:.4f}") check("W0c tight: largest mode ratio at k <= 3 at least 0.9", tm >= 0.9, f"{tm:.4f}") xs = roots(d)[0] worst = 0.0 for k, th0, sg, mb, Ec, zc in cst: mu = complex(np.log(xs[k - 1])) - 1j * 2 * math.pi * k / (d - 1) worst = max(worst, abs(mu) / mb) check("W0d the solver's roots k <= K4 lie in the discs |mu| <= mubar", worst <= 1, f"largest |mu|/mubar {worst:.4f}") for k, th0, sg, mb, Ec, zc in cst[:3] + cst[-1:]: print(f" k = {k}: sigma {sg:.4f}, E_c {Ec:.3f}, zeta0/tau {zc:.3e}", flush=True) return cst def section_plant(cst): """W0e (added after the first run): the PLANT, W0a's check run on broken majorants, l's sign flipped in A and beta1 = 0 in A_k; both must fail.""" global image_A, mode_A print("W0e PLANT: broken majorants at d = 400", flush=True) good_A, good_Ak = image_A, mode_A def bad_A(a, tau, s, P0, arm="uniform"): return good_A(a, tau, -s, P0, arm) def bad_Ak(k, th0, Ec, zc, tau, P0): b2 = (2 * math.pi * k) ** 3 * tau z0 = zc * tau return (np.sqrt(b2 ** 2) + Ec * np.exp(z0) + math.sqrt(1 + 2.25 * th0 * th0) * P0 * np.expm1(z0)) image_A, mode_A = bad_A, bad_Ak try: mi, mo = samples(399, 4000, 2000, 7) ri, wi = image_ratios(400, 399, mi) rm, km = mode_ratios(400, 399, cst, mo) finally: image_A, mode_A = good_A, good_Ak check("W0e PLANT: l's sign flipped, W0a fails on the images", ri.max() > 1 + 1e-12, f"largest ratio {ri.max():.4g}") check("W0e PLANT: beta1 = 0, W0a fails on the modes", rm.max() > 1 + 1e-12, f"largest ratio {rm.max():.4f}") def section_uniform(cst): print("W1 constants at P0 = 399 against d = 800, 1600", flush=True) for d in (800, 1600): mi, mo = samples(d - 1, 1500, 600, 11) t0 = time.time() ri, wi = image_ratios(d, 399, mi) rm, km = mode_ratios(d, 399, cst, mo) check(f"W1 images within A/P at d = {d}", ri.max() <= 1 + 1e-12, f"largest ratio {ri.max():.6f} over {len(ri)} terms") check(f"W1 modes within A_k/P at d = {d}", rm.max() <= 1 + 1e-12, f"largest ratio {rm.max():.6f} over {len(rm)} terms " f"({time.time() - t0:.0f} s)") def section_R(cst): print("W2, KILL: R_maj against R_prov", flush=True) for d, P0 in ((200, 199), (400, 399)): t0 = time.time() P = d - 1 c = cst if P0 == 399 else mode_consts(P0) o = stream(d, P0, c) Rp = o["prov"].sum() + o["tail"] Ro = o["orac"].sum() + o["tail"] Rm = o["maj"].sum() + o["tail_maj"] sh = lambda v: " ".join(f"{t:.4e}" for t in v) # noqa: E731 print(f" d = {d} (P0 = {P0}): R_maj {Rm:.5f} (P R {P * Rm:.2f}), " f"R_prov {Rp:.5f}, oracle {Ro:.5f}; by band [1/c, 1/16) [1/16, " f"1/4) [1/4, 1): maj {sh(o['maj'])}, prov {sh(o['prov'])}; " f"tails {o['tail_maj']:.1e}, {o['tail']:.1e}; B_maj below B_prov " f"by at most {o['worst_pm']:.1e} relative; |h - phi| - B_maj at " f"most {o['worst_hm']:.1e} ({time.time() - t0:.0f} s)", flush=True) if d == 400: check("W2 R_maj/R_prov in [1, 1.10] at d = 400", 1 <= Rm / Rp <= 1.10, f"{Rm / Rp:.4f}") check("W2 B_maj >= B_prov at every m (1e-12)", o["worst_pm"] <= 1e-12, f"{o['worst_pm']:.1e}") F = F_edge(C, 399)[0] check("KILL not fired: R_maj + F < 1/4 at d = 400", Rm + F < 0.25, f"{Rm + F:.5f}") def section_arms(cst): """W4 (added after the run, read only): the images' band [1/c, 1/16) at d = 400 under the exact per-term B, the uniform A, its clamps kept with f at P0, the bracket at P0 alone, and the first order |l| (no bound); and the m where B_maj falls furthest below B_prov on the modes past tau = 1/4.""" print("W4 arms at d = 400", flush=True) d, P0 = 400, 399 P = d - 1 xs, wt, b, _ = roots(d) cw = coef(xs, d) * wt absx = np.abs(xs) s = signs(P) M = P ** 3 acc = {key: 0.0 for key in ("exact", "uniform", "clamped", "at P0", "first")} m0 = M // C while m0 < M // 16: n = int(min(2 ** 17, max(4096, m0 // 16), M // 16 - m0)) m = np.arange(m0, m0 + n, dtype=np.int64) tau = m / P ** 3 thi = (m0 + n) / P ** 3 ph = -kx((m % P) / P, tau) / P ** 2 sm = s[m % P] Ts, Ps, D, tailT = images(m, P, thi) Bs = {"exact": D + np.abs(Ps - ph) + tailT, "uniform": images_maj(m, P, P0, thi), "clamped": images_maj(m, P, P0, thi, "clamped"), "at P0": images_maj(m, P, P0, thi, "at P0"), "first": images_maj(m, P, P0, thi, "first")} for key, B in Bs.items(): bad = sm * ph <= B acc[key] += float(np.where(bad, np.abs(ph) - sm * ph + 2 * B, 0).sum()) m0 += n print(" images over [1/c, 1/16): " + ", ".join( f"{k} {v:.5f}" for k, v in acc.items()), flush=True) best = (0.0, None) m0 = M // C 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)) if m0 / P ** 3 >= 1 / 4: m = np.arange(m0, m0 + n, dtype=np.int64) tau = m / P ** 3 ph = -kx((m % P) / P, tau) / P ** 2 drop = (np.abs(cw) * 1e-20).sum() if len(act) < len(xs) else 0.0 D4, S4 = modes(m, P, xs, cw, act + 1) Bf = D4 + np.abs(S4 - ph) + drop Bfm = modes_maj(tau, d, P, P0, cst, m0) r = (Bf - Bfm) / Bf i = int(r.argmax()) if r[i] > best[0]: best = (float(r[i]), int(m[i]), float(tau[i]), float(Bf[i]), float(Bfm[i]), float(D4[i]), drop) m0 += n print(f" modes past 1/4: B_maj below B_prov by {best[0]:.4f} at m = " f"{best[1]} (tau {best[2]:.6f}): B_prov {best[3]:.3e} of which the " f"dropped roots' floor {best[6]:.3e} and the roots' per-term sum " f"{best[5]:.3e}; B_maj {best[4]:.3e}", flush=True) def main(): cst = section_control() section_plant(cst) if not all(CHECKS): print("control failed; stopping") sys.exit(1) section_uniform(cst) section_R(cst) section_arms(cst) 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()