""" variation.py -- h's total variation D_1 bounded at every P >= 399 at once, each piece a closed form or a P-free constant, and the larger sup's excess X it leaves. QUESTION Q. Is there a bound D1U(P) >= D_1 holding at every P >= 399, each piece a closed form or a P-free constant, whose X stays below 1/4 - 0.0636 at every P >= 399? Notation as excess.py: T_u(j) = C(n, u) 2^-n j/(2n), n = 2u - j, j <= u; U = floor(P^2/16), M = UP. BY HAND. (Q1) The bumps' variation in closed form. For 1 <= j <= u - 1, T(j+1)/T(j) = 2(u - j)(j + 1)/(j(2u - j - 1)) > 1 iff j^2 + j < 2u; for j = -i, i >= 2, |T(-i+1)|/|T(-i)| > 1 iff i^2 - i > 2u. So T is 0 at j = 0, rises to its max at j* = min{j : j^2 + j >= 2u}, falls to T(u) = 2^-(u+1) and jumps to 0; on j < 0 it is negative, |T| largest at i* = max{i : i^2 - i <= 2u}, falling to 0. Exactly V_u = 2 T(j*) + 2 |T(-i*)|. Robbins with its correction negative (n >= m + 1) and nH(m/n) <= n ln 2 - 2(m - n/2)^2/n give, for 1 <= u <= n - 1, C(n, u) 2^-n <= sqrt(n/(2 pi u (n - u))) e^(-j^2/(2n)), and sup_x x e^(-x^2/(2n)) = sqrt(n) e^(-1/2), so u V_u <= K(u) = (2 pi e)^(-1/2) [(1 - j*/u)^(-1/2) + (1 + i*/u)^(-1/2)], -> kappa. j*/u <= (1 + sqrt(1 + 8u))/(2u), falling in u, and the second bracket term is below 1, so past a = 4e4 sum_(a < u <= U) V_u <= K_a ln(U/a), K_a = (2 pi e)^(-1/2) [(1 - (1 + sqrt(1 + 8a))/(2a))^(-1/2) + 1] = 0.484804 by hand; below a the table V_u (excess.py). (Q2) The far bumps, sharper than excess.py's 2 sum |T|. Bump u = U + s meets m <= M only at j <= -sP, and there i = -j has i^2 - i > 2u (s^2 P^2 - sP > 2(U + s) as U <= P^2/16), so |T_u| is monotone on that window: its variation over m in [0, M] is at most |T_(U+s)(-sP)|. far' = sum_(s >= 1) |T_(U+s)(-sP)|, by (Q1)'s bound with u(n - u) = u(u + i): |T_u(-i)| <= i/(2 sqrt(n)) (2 pi u (u + i))^(-1/2) e^(-i^2/(2n)). By hand about 16 e^-4/(sqrt(pi) P^2), 1.0e-6 at P = 399, falling like P^-2; s past 50 under Hoeffding (|T| <= e^(-i^2/(2n))/2) below e^(-2000). (Q3) The edge, exactly: 1.5 + c_0 eps/beta, eps = 2 - beta = beta^-d, c_0 = (1 - eps)/(2 - (d + 1) eps) in (0, 1); at most 1.5 + (2/3)^d. L_s - L = 2 c_0 - 1 = (d - 1) eps/(2 - (d + 1) eps). (Q4) The root tail, a P-free constant. Modes with theta = 2 pi k/P <= 0.35: the disc |mu| <= sigma theta^2/P (majorant.py W2, sigma read at theta0 = 2 pi k/P0 for k <= 22 and at theta = 0.35 past, the self-map's left side rising in theta, falling in P) holds the root, and |P mu + theta^2 + i theta^3| <= rhobar, so Re(P mu) <= -theta^2 (1 - r), r = rhobar/theta^2 read at (theta0, P0) (rhobar/theta^2 rises in theta, falls in P). With |c| <= (theta/P)C1, C1 = sqrt(1 + 2.25 theta0^2) + E_c/P0; |1 - x| <= theta C2, C2 = expm1(theta0 + mubar0)/theta0; 1/(1 - |x|) <= 1 + P/(theta^2 (1 - r)) (1 - e^-a >= a/(1 + a)); U theta^2 >= (pi^2 k^2/4)(1 - 16/P^2): mode k <= 2 C1 C2 [theta0^2/P0 + 1/(1 - r)] exp(-(pi^2 k^2/4)(1 - 16/P0^2)(1 - r)), whatever P >= P0. Modes with theta > 0.35: phi >= 2 pi k/d > 0.349125, Q >= Q0 = 5 - 4 cos 0.349125, |c| <= 2/P (|(d+1)x - 2d| >= P), |1 - x| <= 2, at most d/2 of them: (4d/P)(1 + 2d/ln Q0) e^(-UP ln Q0/(2d)), falling in P (e^-1072 at 399). Sum: R_r = 2 sum_k e^(-pi^2 k^2/4) (1 + 1%) = 0.1715 by hand (k = 1: C1 1.00095, C2 1.0079, r 3.5e-4), the exact tail 0.1698 at d = 400. That the roots are one per sector, so the fixed points are all of them, is (Q6). (Q5) D1U(P) = edge + (table to U, K_a ln(U/a) past) + far' + R_r, so past P = 800, D1U <= 2K_a ln P + 0.9069 + far' (2K_a = 0.969606, U = floor(P^2/16) <= P^2/16). X(P) as excess.py (Z1) with L >= W >= (P^2 - 1)/(4P) (E >= 0): X falls as L rises ((2/beta)D1U <= 2L + 1 - L_s), rises with D1U, so both substitutions are safe. Past the read range, t1 >= (P/4 - 2)/D1U - 1 grows like P/ln P and X <= P eps + 1.99^-t1 2.02 D1U falls. (Q6) The roots, one per sector (written after the first two runs, read by Q6's check in the third). On the closed unit disc, Re(2 - x) >= 1 and |2 - x| >= 1, so g_k(x) = e^(2 pi i k/d) (2 - x)^(-1/d) (principal branch) maps the disc into itself with |g_k'| <= 1/d: one fixed point x_k each, k = 0..d - 1 (Banach). A root in the disc has x^d = 1/(2 - x), so it is x_k for exactly one k; x_0 = 1, and |x_k| = 1 forces |2 - x_k| = 1, so x_k = 1: the P roots inside are x_1..x_P, one per sector, dphi - 2 pi k = arg(1/(2 - x)) in [-pi/6, pi/6] (the largest |Im x|/(2 - Re x) on the disc is 1/sqrt 3, at cos phi = 1/2). A root in majorant.py's disc about i theta_k, theta_k <= 0.35, has dphi - 2 pi k = theta_k + d Im mu, within 0.35 + sigma theta^2 (d/P) < 0.6 of 0, so in no other sector's window (2 pi - pi/6 away): it is x_k. The owed clause of (Q4) and of majorant.py's W2 is this. PREDICTIONS, CONTROL AND KILL (frozen before code). Q0 CONTROL, read first. (a) (Q1)'s identity: 2T(j*) + 2|T(-i*)| equals excess.py's V_table at every u = 1..4e4, 1e-12 relative; PLANT: j* + 1 in its place fails at some u. (b) far''s monotone window: the windowed variation of bump U + s over m in [0, M] equals |T(-sP)| - |T(-uP)| at d = 40, 100, s = 1, 2, 3, to 1e-15 absolute. (c) X recomputed from excess.py's D_1^+ = 6.713490 and its L at d = 400 gives its 2.049e-4. Q1 K(u) >= u V_u at every u = 4..4e4; K(4e4)/(4e4 V) - 1 in [0, 1e-4]; K_a in [0.4847, 0.4849]. Q2 far' at P = 399 in [5e-7, 2e-6], below excess.py's far 1.1e-4. Q3 R_r in [0.1698, 0.175], every mode bound above the exact root term at d = 400, 800, 1600 (k <= 22, 30 digits not needed: doubles); PLANT: the exponent's 16/P^2 read as 16/P fails at d = 400. Q4 D1U(399) - D_1^+(399) in [0, 0.01] (D_1^+ = 6.713490). Q5 X(P) at every P = 399..10^6: largest within t1 = 14's plateau, at most 2.1e-4; D1U(P) <= 0.9696 ln P + 0.91 past P = 800. KILL, refrozen from ROAD's (its "below the exact D_1^+" assumed the far piece unchanged; (Q2) moves it, so the wrong-bound test is per piece): (K1) a piece's P-uniform bound below its exact counterpart at d = 400 or 800 (Sum V_u, far' against exact binomials, the root tail, the edge), or D1U below the exact D_1 there; (K2) X(399) from D1U at or above 1/4 - 0.0636, or X(P) at or above it at any P read. AMENDED BEFORE THE RUN: Q3's plant as frozen moves the bound up (1 - 16/P widens the exponent's slack), so it could not fail; the plant reads the cut as U = P^2/16 + P, the factor 1 + 16/P0. AFTER THE SECOND RUN, frozen before its own: Q6 g_k's fixed points, iterated 40 times from 0, equal heat.roots' x_k, k = 1..d/2, at d = 400 and 1600 to 1e-12; |d arg x_k - 2 pi k| <= pi/6 at both; the window 0.35 + sigma(0.35) 0.35^2 (d/P) at P = 399 in [0.5, 0.6), by hand 0.544; PLANT: g_(k+1)'s fixed points in x_k's place miss at the float arm's 1e-9 bar. AFTER THE FIRST RUN: Q0a's 1e-12 bar missed, 3.2e-10 at u = 36242; against 40 digits the table and the closed form err there by -1.3e-10 and +1.8e-10, and by 1e-11 to 1e-10 both from u = 10^4 on (a 40-digit diagnostic run beside the script, not kept): lgamma's float wall near 8e5 (ulp 1.2e-10), not the identity. The float arm's bar is now 1e-9, and the identity is read against a 40-digit brute variation at u = 10, 100, 1000. FINDINGS (copied from the second run's prints, Q6's from the fourth). - Q0 holds: the identity to 3.2e-10 against the float table (the wall) and to 2.3e-40 against 40 digits; j* + 1 misses by 0.23; the far bumps' windowed variation equals the end value to 2.1e-25 at d = 40, 100, s = 1, 2, 3; excess.py's X at d = 400 recomputes, 2.0488e-4. - Q1 holds: K(u) >= u V_u at every u = 4..4e4, 0.666854 against 0.550781 at u = 4, 0.483951 against 0.483947 at 4e4 (7.6e-6); K_a = 0.484803, kappa 0.483941. - Q2 holds: far' is 1.0910e-6 at P = 399 (exact 1.0907e-6), 2.6550e-7 at 799, P^2 far' 0.1737 and 0.1695: a hundred times below excess.py's far. - Q3 holds: R_r = 0.171483, mode 1's bound 0.171376 against the exact 0.169726, 0.169639, 0.169617 at d = 400, 800, 1600 (C1 1.000953, C2 1.007956, r 3.56e-4); the exact tails 0.169831, 0.169743, 0.169721; past k = 22 and past theta = 0.35 below e^-1000 (e^-1063). The plant (cut read as P^2/16 + P) falls below at modes 1 to 4. - Q4 holds: D1U(399) = 6.715030, 0.00154 above excess.py's D_1^+ and 0.1397 above the exact D_1 (6.575283); at d = 800 7.387126 against 7.247336. - Q5 holds: X is largest at P = 408, 2.0559e-4 (t1 = 14), 2.0493e-4 at 399, 5.50e-8 at 799, 1.70e-76 at 9999, 10^-622.7 at 10^5; D1U - 0.9696 ln P lies in [0.9069, 0.9081] at the P printed, at most 0.9070 past 800. - The KILL does not fire: D1U, far' and the root tail sit above the exact at d = 400, 800 (Sum V_u is the table there, the same numbers; the edge's bound 1.5 + (2/3)^d is read as 1.5, below the doubles' resolution), and X stays below 1/4 - 0.0636 at every P = 399..10^6. Past 10^6 by (Q5) by hand, not run. Unproved: the table and lgamma in doubles; terms below the doubles' resolution left out of the sums, the edge's (2/3)^d, L_s - L's P (2/3)^d, and beta read as 2 in eps_t1 = beta^-(t1+1) (B = beta exactly) and in 2/beta, which lower X by a factor near 1 + (t1 + 2) 2^-(d+1). - Q6 holds: g_k's fixed points equal heat.roots' x_k to 7.3e-16 at d = 400 and 7.5e-16 at 1600; |d arg x_k - 2 pi k| is at most 0.5228 and 0.5234 against pi/6 = 0.5236; sigma(0.35) = 1.5822, the window 0.5443; g_(k+1)'s fixed points miss x_k by at least 1.6e-2 and 3.9e-3, the plant caught at each d (its nearest, 3.9e-3). The plant reads that the fixed points are indexed by k alone: g_(k+1)'s is the root x_(k+1), caught by construction. RUN RECORD: 15/15 checks; 60.0 s under the 512 MB guard, peak commit 162.1 MB. The first run stopped at its control, Q0a's 1e-12 bar (above). The second read 13/13 in 59.5 s, 164.0 MB; K1's check was renamed after it to the pieces it reads, its condition unchanged, and Q6's two checks were added after it, frozen before the third run. The third read 15/15 with the plant's nearest taken as the larger of the two d's (1.6e-2), so a plant missed at one d passed; the fourth takes the smaller. """ import math import os import sys import time 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 majorant import mode_const # noqa: E402 from excess import V_table, edge, excess # noqa: E402 from route import stream # noqa: E402 CHECKS, NAMES = [], [] MISSED = () P0, A_TAB, TH_CAP = 399, 40000, 0.35 G = (2 * math.pi * math.e) ** -0.5 BUDGET_X = 0.25 - 0.0636 D1P_399 = 6.713490 # excess.py, d = 400 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 logT(u, j): """ln |T_u(j)| exactly (lgamma), j != 0, j <= u.""" n = 2 * u - j return (math.lgamma(n + 1) - math.lgamma(u + 1) - math.lgamma(n - u + 1) - n * math.log(2) + math.log(abs(j) / (2 * n))) def logT_rob(u, j): """(Q1)'s Robbins bound on ln |T_u(j)|, 1 <= u <= n - 1.""" n = 2 * u - j return (math.log(abs(j) / (2 * n)) + 0.5 * math.log(n / (2 * math.pi * u * (n - u))) - j * j / (2 * n)) def jstar(u): j = max(1, math.ceil((-1 + math.sqrt(1 + 8 * u)) / 2) - 1) while j * j + j < 2 * u: j += 1 return j def istar(u): i = math.floor((1 + math.sqrt(1 + 8 * u)) / 2) + 1 while i * i - i > 2 * u: i -= 1 return i def V_closed(u, shift=0): j = min(jstar(u) + shift, u) return 2 * math.exp(logT(u, j)) + 2 * math.exp(logT(u, -istar(u))) def K_of(u): return G * ((1 - jstar(u) / u) ** -0.5 + (1 + istar(u) / u) ** -0.5) def K_a(a): return G * ((1 - (1 + math.sqrt(1 + 8 * a)) / (2 * a)) ** -0.5 + 1) def far_prime(P, exact=False): """sum_(s <= 50) |T_(U+s)(-sP)|, by the Robbins bound or exactly.""" U = P * P // 16 tot = 0.0 for s in range(1, 51): u, i = U + s, s * P assert i * i - i > 2 * u tot += math.exp((logT if exact else logT_rob)(u, -i)) return tot def mode_bound(k, th0, plant=False): sg, mb, ec, rho = mode_const(th0, P0) r = rho / th0 ** 2 C1 = math.sqrt(1 + 2.25 * th0 * th0) + ec C2 = math.expm1(th0 + mb) / th0 cut = (1 + 16 / P0) if plant else (1 - 16 / P0 ** 2) return (2 * C1 * C2 * (th0 * th0 / P0 + 1 / (1 - r)) * math.exp(-(math.pi * k) ** 2 / 4 * cut * (1 - r))), r, C1, C2 def q_part(P): """ln of the modes past theta = 0.35: (4d/P)(1 + 2d/ln Q0) e^(-UP lnQ0/(2d)).""" d, U = P + 1, P * P // 16 lq = math.log(5 - 4 * math.cos(TH_CAP * P0 / (P0 + 1))) return (math.log(4 * d / P) + math.log1p(2 * d / lq) - U * P * lq / (2 * d)) def root_const(plant=False): K4 = int(TH_CAP * P0 / (2 * math.pi)) per = [mode_bound(k, 2 * math.pi * k / P0, plant)[0] for k in range(1, K4 + 1)] sg, mb, ec, rho = mode_const(TH_CAP, P0) rm = rho / TH_CAP ** 2 C1 = math.sqrt(1 + 2.25 * TH_CAP ** 2) + ec C2 = math.expm1(TH_CAP + mb) / TH_CAP k = K4 + 1 ex = (math.pi * k) ** 2 / 4 * (1 - 16 / P0 ** 2) * (1 - rm) ratio = math.exp(-(math.pi ** 2 / 4) * (2 * k + 1) * (1 - 16 / P0 ** 2) * (1 - rm)) past = 2 * C1 * C2 * (TH_CAP ** 2 / P0 + 1 / (1 - rm)) * math.exp(-ex) \ / (1 - ratio) return per, past, rm def exact_modes(d): P = d - 1 U = P * P // 16 x, wt, b, _ = roots(d) c = coef(x, d) ax = np.abs(x) return wt * np.abs(c) * np.abs(1 - x) * np.exp(U * P * np.log(ax)) / (1 - ax) def section_control(V): print("Q0 CONTROL", flush=True) t0 = time.time() err = max(abs(V_closed(u) / V[u] - 1) for u in range(1, A_TAB + 1)) check("Q0a V_u = 2T(j*) + 2|T(-i*)| at every u = 1..4e4 (1e-9 rel)", err < 1e-9, f"largest {err:.1e} ({time.time() - t0:.0f} s)") mp.mp.dps = 40 worst = 0 for u in (10, 100, 1000): lo = -int(30 * math.sqrt(u)) - 100 T = [mp.binomial(2 * u - j, u) * mp.mpf(2) ** (j - 2 * u) * j / (2 * (2 * u - j)) for j in range(lo, u + 1)] var = abs(T[0]) + sum(abs(T[k + 1] - T[k]) for k in range(len(T) - 1)) + abs(T[-1]) cl = 2 * T[jstar(u) - lo] + 2 * abs(T[-istar(u) - lo]) worst = max(worst, abs(var / cl - 1)) check("Q0a the identity by 40-digit brute variation at u = 10, 100, 1000", worst < 1e-30, f"largest {mp.nstr(worst, 3)}") errp = max(abs(V_closed(u, 1) / V[u] - 1) for u in range(1, A_TAB + 1)) check("Q0a PLANT: j* + 1 in its place fails", errp >= 1e-12, f"largest {errp:.1e}") worst = 0.0 for d in (40, 100): P = d - 1 U = P * P // 16 for s in (1, 2, 3): u = U + s i = np.arange(u * P, s * P - 1, -1) n = 2 * u + i lg = np.array([math.lgamma(v + 1) for v in n]) - math.lgamma(u + 1) \ - np.array([math.lgamma(v - u + 1) for v in n]) \ - n * math.log(2) + np.log(i / (2.0 * n)) T = np.exp(lg) var = np.abs(np.diff(T)).sum() gap = abs(var - (T[-1] - T[0])) worst = max(worst, gap) print(f" d = {d}, s = {s}: windowed variation {var:.6e}, " f"|T(-sP)| - |T(-uP)| {T[-1] - T[0]:.6e}", flush=True) check("Q0b the far bumps' variation is the end value at d = 40, 100", worst < 1e-15, f"largest gap {worst:.1e}") L, E, h0, D, tail, b, integ = stream(400) L = L + tail X, XR, t0_, t1 = excess(400, L, h0, D1P_399, b) check("Q0c X from excess.py's D_1^+ at d = 400 is its 2.049e-4", abs(X - 2.049e-4) < 5e-8, f"{X:.4e}, t0 {t0_:.3f}") return L, h0, b, D[1] + 2 * tail def section_bumps(V): print("Q1 THE BUMPS: K(u) >= u V_u, K_a past the table", flush=True) bad = [u for u in range(4, A_TAB + 1) if K_of(u) < u * V[u] * (1 - 1e-13)] u = A_TAB rel = K_of(u) / (u * V[u]) - 1 ka = K_a(A_TAB) for v in (4, 10, 100, 1000, 10000, 40000): print(f" u = {v}: j* {jstar(v)}, i* {istar(v)}, u V_u " f"{v * V[v]:.6f}, K(u) {K_of(v):.6f}", flush=True) print(f" K_a at a = {A_TAB}: {ka:.6f}; kappa {2 * G:.6f}", flush=True) check("Q1 K(u) >= u V_u at every u = 4..4e4, K(4e4) within 1e-4, K_a " "in [0.4847, 0.4849]", not bad and 0 <= rel <= 1e-4 and 0.4847 <= ka <= 0.4849, f"below at {bad[:5]}, K(4e4)/(u V) - 1 {rel:.2e}, K_a {ka:.6f}") return ka def section_far(): print("Q2 THE FAR BUMPS by their end values", flush=True) rows = {} for P in (399, 799): fb, fe = far_prime(P), far_prime(P, exact=True) rows[P] = (fb, fe) print(f" P = {P}: far' bound {fb:.4e}, exact {fe:.4e}, P^2 far' " f"{P * P * fb:.4f}", flush=True) check("Q2 far' at P = 399 in [5e-7, 2e-6]", 5e-7 <= rows[399][0] <= 2e-6, f"{rows[399][0]:.3e}") return rows def section_roots(): print("Q3 THE ROOT TAIL, a P-free constant", flush=True) per, past, rm = root_const() qp = q_part(P0) R = sum(per) + past + math.exp(qp) print(f" per mode k = 1, 2, 3: {per[0]:.6f} {per[1]:.3e} {per[2]:.3e};" f" past k = 22 at theta 0.35 (r {rm:.4f}): {past:.1e}; theta >" f" 0.35: e^{qp:.0f}", flush=True) b1, r1, C1, C2 = mode_bound(1, 2 * math.pi / P0) print(f" k = 1: C1 {C1:.6f}, C2 {C2:.6f}, r {r1:.2e}; R_r {R:.6f}, " f"2 sum e^(-pi^2 k^2/4) " f"{2 * sum(math.exp(-(math.pi * k) ** 2 / 4) for k in range(1, 9)):.6f}", flush=True) low, plant_low, tails = [], [], {} pp, _, _ = root_const(plant=True) for d in (400, 800, 1600): ex = exact_modes(d) tails[d] = float(ex.sum()) K4 = len(per) low += [(d, k + 1) for k in range(K4) if per[k] < ex[k]] plant_low += [(d, k + 1) for k in range(K4) if pp[k] < ex[k]] print(f" d = {d}: exact tail {tails[d]:.6f}, mode 1 {ex[0]:.6f} " f"(bound {per[0]:.6f}), modes past 22 {ex[K4:].sum():.1e}", flush=True) check("Q3 R_r in [0.1698, 0.175], every mode bound above the exact " "term at d = 400, 800, 1600", 0.1698 <= R <= 0.175 and not low and all(tails[d] <= R for d in tails), f"R_r {R:.6f}, below at {low}") check("Q3 PLANT: the cut read as P^2/16 + P fails at d = 400", any(d == 400 for d, k in plant_low), f"below at {plant_low[:4]}") return R, tails def fixed_points(d, shift=0): """x_k = g_(k+shift)'s fixed point, iterated from 0, k = 1..d/2.""" w = np.exp(2j * np.pi * (np.arange(1, d // 2 + 1) + shift) / d) x = np.zeros(d // 2, dtype=complex) for _ in range(40): x = w * (2 - x) ** (-1.0 / d) return x def section_sectors(): print("Q6 THE ROOTS, one per sector, by g_k's contraction", flush=True) worst, plant = 0.0, math.inf arg = 0.0 for d in (400, 1600): x, _, _, _ = roots(d) k = np.arange(1, d // 2 + 1) err = float(np.abs(fixed_points(d) - x).max()) mis = float(np.abs(fixed_points(d, 1) - x).min()) dev = float(np.abs(d * np.angle(x) - 2 * np.pi * k).max()) worst, plant, arg = max(worst, err), min(plant, mis), max(arg, dev) print(f" d = {d}: |fixed point - heat.roots| at most {err:.1e}; " f"g_(k+1)'s at least {mis:.1e}; |d arg x_k - 2 pi k| at most " f"{dev:.4f}", flush=True) sg = mode_const(TH_CAP, P0)[0] win = TH_CAP + sg * TH_CAP ** 2 * (P0 + 1) / P0 print(f" sigma(0.35) {sg:.4f}; the window 0.35 + sigma theta^2 (d/P) " f"{win:.4f}; pi/6 {math.pi / 6:.4f}", flush=True) check("Q6 g_k's fixed points are heat.roots' x_k at d = 400, 1600 " "(1e-12), each within pi/6 of its sector, the window in [0.5, 0.6)", worst < 1e-12 and arg <= math.pi / 6 and 0.5 <= win < 0.6, f"largest {worst:.1e}, {arg:.4f}, window {win:.4f}") check("Q6 PLANT: g_(k+1)'s fixed points miss x_k at the 1e-9 bar", plant >= 1e-9, f"nearest {plant:.1e}") def D1U_all(Ps, Vc, ka, R): U = Ps * Ps // 16 bumps = np.where(U <= A_TAB, Vc[np.minimum(U, A_TAB)], Vc[A_TAB] + ka * np.log(np.maximum(U, 1) / A_TAB)) far = np.zeros(len(Ps)) for s in range(1, 51): u = (U + s).astype(float) i = (s * Ps).astype(float) n = 2 * u + i far += np.exp(np.log(i / (2 * n)) + 0.5 * np.log(n / (2 * np.pi * u * (u + i))) - i * i / (2 * n)) return 1.5 + bumps + far + R, bumps, far def X_all(Ps, D): L = (Ps * Ps - 1) / (4.0 * Ps) Ls = L # + 2 c_0 - 1 <= P (2/3)^d t1 = np.floor((Ls - 1) / D) B = 2.0 # (1 - 2^-d)/(1/2), beta = 2 in doubles lg = -t1 * math.log(2) - math.log(B) + np.log(np.minimum(2 * L + 1 - Ls, D)) return np.exp(lg), t1, lg / math.log(10) def section_read(Vc, ka, R, far_rows, ctrl, tails): print("Q4, Q5 and KILL: D1U and X at every P = 399..10^6", flush=True) L400, h0, b, D1_400 = ctrl _, _, _, D, tail, _, _ = stream(800) D1_800 = D[1] + 2 * tail Ps = np.arange(399, 10 ** 6 + 1, dtype=np.int64) D, bumps, far = D1U_all(Ps, Vc, ka, R) X, t1, lx = X_all(Ps, D) for P in (399, 400, 799, 800, 1599, 3199, 9999, 10 ** 5, 10 ** 6): k = P - 399 print(f" P = {P}: D1U {D[k]:.6f} (bumps {bumps[k]:.6f}, far' " f"{far[k]:.2e}), - 0.9696 ln P {D[k] - 0.9696 * math.log(P):.4f};" f" t1 {t1[k]:.0f}, X {X[k]:.3e} (log10 {lx[k]:.1f})", flush=True) for d, Dex in ((400, D1_400), (800, D1_800)): k = d - 1 - 399 fb, fe = far_rows[d - 1] print(f" d = {d}: D1U {D[k]:.6f} against the exact D_1 {Dex:.6f};" f" far' {fb:.3e} against {fe:.3e}; R_r {R:.6f} against " f"{tails[d]:.6f}", flush=True) wrong = [d for d, Dex in ((400, D1_400), (800, D1_800)) if D[d - 400] < Dex or far_rows[d - 1][0] < far_rows[d - 1][1] or R < tails[d]] check("KILL K1 not fired: D1U, far' and the root tail above the exact " "at d = 400, 800", not wrong, f"wrong at {wrong}") g = D[0] - D1P_399 check("Q4 D1U(399) - D_1^+(399) in [0, 0.01]", 0 <= g <= 0.01, f"{g:.5f}") km = int(np.argmax(X)) past = Ps > 800 sl = float((D[past] - 0.9696 * np.log(Ps[past])).max()) check("Q5 X largest within t1 = 14's plateau, at most 2.1e-4; D1U <= " "0.9696 ln P + 0.91 past 800", t1[km] == 14 and X[km] <= 2.1e-4 and sl <= 0.91, f"largest X {X[km]:.4e} at P = {Ps[km]} (t1 " f"{t1[km]:.0f}), D1U - 0.9696 ln P at most {sl:.4f}") over = Ps[X >= BUDGET_X] check("KILL K2 not fired: X below 1/4 - 0.0636 at every P read", X[0] < BUDGET_X and len(over) == 0, f"X(399) {X[0]:.4e}, at or above at {over[:5]}") def main(): t0 = time.time() V = V_table(A_TAB) Vc = np.cumsum(V) ctrl = section_control(V) if not all(ok for n, ok in zip(NAMES, CHECKS) if n not in MISSED): print("control failed; stopping") sys.exit(1) ka = section_bumps(V) far_rows = section_far() R, tails = section_roots() section_sectors() section_read(Vc, ka, R, far_rows, ctrl, tails) print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed, {time.time() - t0:.0f} s") 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()