majorant.py
Cited on: READING.
To run it: save majorant.py, with edge.py, heat.py, real_support.py, remainder.py beside it; needs mpmath, numpy, scipy.
"""
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()