uniform.py
Cited on: READING.
To run it: save uniform.py, with edge.py, heat.py, majorant.py, real_support.py, remainder.py beside it; needs mpmath, numpy, scipy.
"""
uniform.py -- the remainder R_maj that majorant.py's per-term majorants
leave, read at every P >= 399 at once: its limit C_R = lim P R_maj, the
side it is approached from, and sup_P R_maj + F against 1/4.
QUESTION X. Does P R_maj(P) tend to a constant C_R, from which side, and
does sup_P R_maj + F stay below 1/4 with c = 400 while F <= 0.01 (P below
1.8e10)? Notation as majorant.py; A and A_k are read at P0 = 399.
BY HAND.
(X1) The tail past M. R_maj stops at M = P^3 and adds 2 sum_(m >= M)
|h_m|, every root by the Q-bound with |c| <= 2/P: ln Q_k ~ 8 pi^2
k^2/d^2, so the tail ~ (8/P)(P^3/(4 pi^2)) e^(-4 pi^2) = 1.46e-18
P^2, 1/4 near P = 4.1e8. The cut moves: M = T P^3, T(P) = 1 +
ln(P/P0)/(2 pi^2), and [1, T) is read by the modes' majorant.
(X2) The limit. Near a sign change x0 in {0, 1/2}, r = x0 P + rho, P^3
phi = -K_xx(x0, tau) rho + O(rho^3/P^2), and the sawtooth's sign
gives s phi = Gamma |rho|/P^3, Gamma = |K_xx(x0, tau)|, on the
lattice Z at x0 = 0, and at x0 = 1/2 at even P; at odd P the lattice
at 1/2 is Z + 1/2 and the class rho = -1/2 is wrong-signed, s phi =
-Gamma/(2 P^3). P^3 B tends to b(x0, tau), P-free: the images sum_k
|G(x0 + k)| A(|x0 + k|/(2 tau), tau), G = (4 pi tau)^(-1/2)
e^(-xi^2/(4 tau)) xi/(2 tau), the modes sum_(k <= 22) 4 pi k
e^(-4 pi^2 k^2 tau) A_k(tau). Away from x0, |phi| ~ P^-2 against B
~ P^-3. With P^2 values of q per unit tau,
C_R = int_(1/c)^T sum_(x0) S dtau, S = 2b #{rho : Gamma|rho| <= b},
plus Gamma at odd P, x0 = 1/2, which integrates to K(1/2, T) -
K(1/2, 1/c).
(X3) At finite P the image term is min(|phi_k| A/P, crude), the crude
branch binding where A/P > 1, near x0 = 1/2 at tau < 0.0063, a
weight below e^(-1/(16 tau)); it leaves as P grows. TRANSPLANT from
the heat flow's C + a/P: P R_maj - C_R = O(1/P).
(X4) A theorem needs one side: R_maj(P) <= (C_R + sup D)/399 with D =
P R_maj - C_R, and 1/4 - F leaves C_R + sup D < 99.
(X5) Windows: at x0 = 1/2 and small tau, b/Gamma ~ 1/(16 tau), 25
classes at 1/400; past 1/4, b/Gamma ~ 4 pi^2 tau. Windows of 160
classes per side hold both to T(1.8e10) = 1.89.
INSTRUMENT. B at m = qP + rho with q a float and rho signed (m overflows
int64 past P near 2e6), images as majorant.images_maj, modes as
modes_maj; the walk's tail past image K is e^(-g) with g about
K^2/(4 tau), P^3-scaled, so K grows: K = ceil(sqrt(4 tau_hi (64 + 3
ln(P/P0)))) + 2, majorant.py's at P0. The roots' tail past K4 = 22 by
exact Q_k (1 - cos as 2 sin^2) for k <= 2000, beyond by ln Q >=
min(1.196 phi^2, 0.846) (phi <= 0.841: sin(phi/2) >= (phi/2)(0.97053),
ln(1 + u) >= 0.63548 u on u <= 4/3; Q rising in phi past it). R by
classes in the windows, the q sum by a midpoint rule in tau (geometric
nodes per band, P^2 per unit tau).
PREDICTIONS, CONTROL AND KILL (frozen before code; X0e, X0f after the
fourth run).
X0 CONTROL, read first, at d = 400: (a) B in (q, rho) equals
majorant.py's images_maj and modes_maj (tails apart) at 2000 sampled
m to 1e-12 relative; (b) no Bad class in the outer 10 of each window
at any tau read, at every P read (at P = 399 the windows are every
class); (c) the quadrature at P = 399, T = 1, within 0.5% of
majorant.py's R_maj by band (0.04496, 0.01836, 1.970e-4), and n_tau
doubled moves it by under 0.1%; (d) PLANT: the window cut to 20
classes fails (b) at d = 800; and C_R's odd term read alone equals
K(1/2, 1) - K(1/2, 1/400) to 1e-6. (e) On every block of 256 nodes
the run reads below tau = 0.004, at P = 1e4, 1e6, 1e10 (odd, M =
T(P) P^3), a dense scan of every class within 5000 of each sign
change finds no Bad class outside the adaptive windows' inner part
(the outer 10 excluded); (f) PLANT: a coarse scan of 6 points
leaves one there.
X1 past_M/P^2 at P = 399, 1e4, 1e6 within 10% of 1.46e-18.
X2 |C_R(odd) - P R_maj(399)| <= 3, C_R(odd) above it; C_R(even)
printed, not predicted.
X3 P |P R_maj(P) - C_R| over P = 799, 1599, 3199 (odd P, c = 400, T
= 1) within a factor 3 of each other (O(1/P)).
KILL: max over P in {399, 799, 1599, 3199, 1e4, 1e5, 1e6, 1e8, 1e10,
1.8e10} (odd) of R_maj(P) + F(400, P), M = T(P) P^3, at or above
1/4; or the X0 control failing. Then the per-term route needs P0
raised or c moved before 1.8e10.
AFTER THE FIRST RUN: X0b failed at P = 399 (852 Bad (class, node)
pairs in the outer 10), where the windows are every class, so R there is
whole (X0c); X0b no longer stops the run, the instrument's checks (X0a,
X0c, X0d) do, and the largest tau of an outer Bad class is printed. By
hand: where A(xi, tau) >= P the image term's crude branch puts B above
|phi_k|, so every class with A >= P is Bad, at tau below about
1/(8 sqrt P) near x0 = 1/2, a weight near P e^(-sqrt(P)/2).
AFTER THE SECOND RUN: X0b failed at every P read, to 1.8e10 (outer Bad
at tau <= 0.0035), as the positive-l image's A = P0(e^(l/P0) - 1), about
1.4e5 at xi = -1/2, tau = 1/400, widens the window at 1/2 past 160
classes at every P read; and C_R with A at P0 diverges (8.84e98 in [1/c, 1/16)),
G A growing like e^(xi^3/(8 tau^2 P0) - xi^2/(4 tau)) over the far
images, which at finite P the crude branch caps. So the windows adapt
(R_quad's coarse scan, the plant X0d run unadapted), the kill and X3
read n_tau = 12000 (24000 to 48000 moved R by 5e-6 at P = 399), and the
limit of the majorants read at P0 = P (A = |l|, A_k = |beta1 - beta2|)
is printed beside, not predicted.
FINDINGS (copied from the fourth run's prints, X0e and X0f the fifth's).
- The instrument holds: B in (q, rho) equals majorant.py's to 1.7e-15
relative; at P = 399 the quadrature reads R_maj by band 0.04496,
0.01836, 1.9696e-4 (R 0.06351, P R 25.342), within 0.0002 relative
of majorant.py, n_tau doubled moving R by 4.8e-6 relative; the plant (windows of
20, unadapted, d = 800) puts 52283 Bad pairs in the outer 10; the odd
term integrates to 1.00000000. X0b fails at P = 399 alone (852 outer
Bad pairs, tau <= 0.00269), where every class is read.
- X1 fails at P = 399: past_M/P^2 over 1.46e-18 is 1.359, 1.006, 0.994
at P = 399, 1e4, 1e6, by hand about 1.35 from X1's 4 pi^2 dropping (P/d)^3. The tail past
P^3 is 0.2438 at P = 4.1e8 and 469.9 at 1.8e10; past T(P) P^3 it is
2.3e-13 to 3.2e-13 at every P read.
- X2 and X3 fail on their premise: with A at P0 the far images' G A
grows like e^(xi^3/(8 tau^2 P0) - xi^2/(4 tau)), and C_R diverges
(8.84e98 in [1/c, 1/16)); its bands past 1/16 are finite, 7.3229
and 0.0783, which P R_maj's reach (7.3234, 0.0783 at P = 3199). The
first orders, the limit of majorants read at P0 = P, give C_R
17.9415 at odd P (11.1199, 6.7480, 0.0736 by band), 16.8206 at even.
- P R_maj at c = 400, M = T(P) P^3: 25.3408, 22.0587, 21.1002, 20.8247,
20.7245, 20.7086, 20.7064, 20.7068, 20.7077, 20.7136 at P = 399, 799,
1599, 3199, 1e4, 1e5, 1e6, 1e8, 1e10, 1.8e10 (odd): within 20.706 to
20.725 past 1e4 (n_tau checked at 399 alone), about 2.77 above the
first orders' C_R, the price of one A for
every P; P R_maj - C_R(first) is 4.12, 3.16, 2.88 at T = 1, P = 799,
1599, 3199.
- The KILL, its X0b clause demoted after the first run, does not fire
over the ten odd P read: R_maj + F is largest at P = 399, 0.063511,
and at 1.8e10 is 0.009999, F's. Even P is unread past 399.
- The Bad set at 1/2 below tau near 0.0035 reaches the edge of a
160-class window at every P read (the second run's X0b); the
adaptive windows read 718 to 3120 classes, and no Bad class sits at
a window's edge from P = 799 to 1.8e10, which the scan's own widening
makes nearly automatic. X0e reads it against a dense scan: on the 12
blocks below tau = 0.004 at P = 1e4, 1e6, 1e10 the farthest Bad class
sits at least 23 classes inside the window's outer 10 (853, 783, 749
classes from 1/2 at the first block, windows 950, 879, 837); the
6-point plant leaves one 615 classes past (X0f).
RUN RECORD: 17/21 checks, the misses X0b at P = 399, X1, X2 and X3, as
recorded; 728.1 s under the 512 MB guard, peak commit 200.9 MB. The
fourth run, before X0e and X0f, read 15/19 in 584.6 s. The third run was
stopped at P = 1e8, its evenly spaced scan having grown a window without
bound, and had failed X0b at P = 1e6 (6732 outer Bad pairs), its 600
evenly spaced classes passing over the Bad set; the scan is now geometric and
capped.
"""
import math
import os
import sys
import time
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import numpy as np # noqa: E402
from heat import kernel # noqa: E402
from edge import kx, signs # noqa: E402
from remainder import F_edge, g_and_dg # noqa: E402
from majorant import (image_A, mode_A, mode_consts, images_maj, # noqa: E402
modes_maj, mode_tails)
CHECKS, NAMES = [], []
MISSED = ("X0b no Bad class in the outer 10 at P = 399",
"X1 past_M/P^2 within 10% of 1.46e-18 at P = 399, 1e4, 1e6",
"X2 |C_R(odd) - P R_maj(399)| <= 3, C_R above it",
"X3 P |P R_maj - C_R| at P = 799, 1599, 3199 within a factor 3")
C, P0, K4 = 400, 399, 22
W, CAP = 160, 5000
BAND_399 = (0.04496, 0.01836, 1.970e-4) # majorant.py's R_maj by band
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 kimg(P, tau_hi):
return int(math.ceil(math.sqrt(4 * tau_hi * (64 + 3 * math.log(P / P0))))
) + 2
def img_B(q, rho, P, tau_hi):
"""images_maj at m = qP + rho, q float, rho signed, K grown with P."""
x = rho / P
tau = (q + x) / P ** 2
K = kimg(P, tau_hi)
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(tau.shape, 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 = rho + 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()
A = image_A(a, tau, np.sign(xi), P0)
nf = np.where(ok, n, 1.0)
y = np.abs(j) / 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(-j * j / (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 lq_exact(k, d):
s = np.sin(np.pi * k / d)
return np.log1p(np.minimum(8 * s * s, 2.0))
def lq_beyond(d):
"""ln Q's lower bound past k = 2000 (or None when d/2 <= 2000)."""
if d // 2 <= 2000:
return None
ph = 2 * np.pi * 2001 / d
return min(1.196 * ph * ph, 0.846)
def root_tail(tau, P, k_lo, factor):
"""factor/P * sum_(k_lo <= k <= d/2) exp(-tau P^3 ln Q_k/(2(P + 1)))."""
d = P + 1
k = np.arange(k_lo, min(d // 2, 2000) + 1)
out = np.zeros(len(tau))
if len(k):
e = -np.outer(tau, lq_exact(k, d)) * (P ** 3 / (2 * (P + 1)))
out += np.exp(e).sum(axis=1)
lb = lq_beyond(d)
if lb is not None:
out += (d // 2 - 2000) * np.exp(-tau * P ** 3 * lb / (2 * (P + 1)))
return factor / P * out
def mode_B(tau, P, cst, tails=True):
"""modes_maj's per-mode sum at each tau, its tails at tau itself."""
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
if tails:
ks = np.arange(K4 + 1, K4 + 400)
B += (4 * np.pi * ks * np.exp(-4 * np.pi ** 2 * np.outer(tau, ks * ks))
).sum(axis=1) / P ** 2
B += root_tail(tau, P, K4 + 1, 4.0)
return B
def past_T(P, T):
"""2 sum_(m >= T P^3) |h_m|, every root by the Q-bound, |c| <= 2/P."""
d = P + 1
k = np.arange(1, min(d // 2, 2000) + 1)
lr = -lq_exact(k, d) / (2 * (P + 1))
out = (8 / P * np.exp(T * P ** 3 * lr) / -np.expm1(lr)).sum()
lb = lq_beyond(d)
if lb is not None:
lr = -lb / (2 * (P + 1))
out += (d // 2 - 2000) * 8 / P * math.exp(T * P ** 3 * lr) / \
-math.expm1(lr)
return float(out)
def nodes(c, T, n_tau):
"""Midpoint nodes and widths in tau, geometric per band."""
cuts = [1 / c] + [v for v in (1 / 64, 1 / 16, 1 / 4, 1.0) if 1 / c < v < T] \
+ [T]
span = math.log(T * c)
ts, ws = [], []
for lo, hi in zip(cuts[:-1], cuts[1:]):
n = max(8, int(n_tau * math.log(hi / lo) / span))
e = np.geomspace(lo, hi, n + 1)
ts.append(np.sqrt(e[:-1] * e[1:]))
ws.append(np.diff(e))
return np.concatenate(ts), np.concatenate(ws)
def sign_of(r, P):
return np.where(r < (P - 1) / 2, 1.0, -1.0)
def dists(r, P):
"""Each class's distance to the sign change at 0 and at 1/2."""
return np.minimum(r, P - r), np.abs(r - (P - 1) / 2)
def summand(t, rho, P, cst):
"""At nodes t and signed classes rho: Bad, and (|phi| - s phi + 2B)
where Bad."""
r = np.mod(rho, P)
s = sign_of(r, P)
q = np.floor(t * P * P)[:, None] + 0 * rho[None, :]
x = rho[None, :] / P
tau = (q + x) / P ** 2
B = np.full(q.shape, np.inf)
lo = t < 1 / 4
if lo.any():
B[lo] = img_B(q[lo], np.broadcast_to(rho, q.shape)[lo], P,
float(tau[lo].max()))
hi = t >= 1 / 64
if hi.any(): # x-free, read at the node's tau
Bf = mode_B(t[hi], P, cst)
B[hi] = np.minimum(B[hi], Bf[:, None])
assert np.isfinite(B).all()
ph = -kx(np.broadcast_to(x, q.shape).ravel(),
tau.ravel()).reshape(q.shape) / P ** 2
bad = s[None, :] * ph <= B
return bad, np.where(bad, np.abs(ph) - s[None, :] * ph + 2 * B, 0.0)
def classes(P, w0, wh):
"""Signed classes within w0 of 0 and wh of (P - 1)/2, every class when
the two windows meet; with each class's nearer window's slack."""
h = (P - 1) // 2
if w0 + wh >= h - 1:
r = np.arange(P)
rho = np.where(r > h + (P - h) // 2, r - P, r)
else:
rho = np.concatenate([np.arange(-w0, w0 + 1),
h + np.arange(-wh, wh + 1)])
d0, dh = dists(np.mod(rho, P), P)
slack = np.where(d0 <= dh, w0 - d0, wh - dh)
if w0 + wh >= h - 1:
slack = np.full(len(rho), np.inf)
far = np.minimum(d0, dh)
slack = far.max() - far
return rho.astype(float), slack
def scan_classes(P, coarse):
"""The coarse scan's signed classes, at `coarse` distances from each
sign change geometric to P/4 (ratio under 1.07 past 20 at 600), with
their distances to 0 and to 1/2."""
dg = np.unique(np.geomspace(1, max(2, P / 4), coarse).round())
h = (P - 1) // 2
rhoc = np.unique(np.concatenate([dg, -dg, h + dg, h - dg, [0, h]]))
return (rhoc,) + dists(np.mod(rhoc, P).astype(int), P)
def adapt(t, P, cst, w, scan, cap):
"""The windows at nodes t: each widened to 1.1 times the farthest Bad
class of the coarse scan near its sign change, plus 20, capped at
`cap` a side; with the coarse Bad classes past the cap (count,
farthest distance)."""
rhoc, d0c, dhc = scan
w0 = wh = w
badc = summand(t, rhoc, P, cst)[0].any(axis=0)
near0 = badc & (d0c <= dhc)
nearh = badc & (dhc < d0c)
if near0.any():
w0 = max(w, int(1.1 * d0c[near0].max()) + 20)
if nearh.any():
wh = max(w, int(1.1 * dhc[nearh].max()) + 20)
past = (near0 & (d0c > cap)) | (nearh & (dhc > cap))
far = float(np.minimum(d0c, dhc)[past].max()) if past.any() else 0.0
return min(w0, cap), min(wh, cap), int(past.sum()), far
def dense_far(t, P, cst, cap):
"""The farthest Bad class from each sign change at nodes t, every class
within `cap` of each read (every class when the two meet)."""
rho, _ = classes(P, cap, cap)
d0, dh = dists(np.mod(rho, P).astype(int), P)
bad = np.zeros(len(rho), bool)
sub = max(1, int(160000 // len(rho)))
for j in range(0, len(t), sub):
bad |= summand(t[j:j + sub], rho, P, cst)[0].any(axis=0)
f0 = d0[bad & (d0 <= dh)]
fh = dh[bad & (dh < d0)]
return (float(f0.max()) if len(f0) else -1.0,
float(fh.max()) if len(fh) else -1.0)
def R_quad(P, c, T, n_tau, cst, w=W, block=256, adaptive=True, coarse=600,
cap=CAP):
"""R_maj over [1/c, T) by the midpoint rule over windowed classes, by
band, with the Bad count in each window's outer 10 classes. Adaptive
(after the first runs): below tau = 1/16 the coarse scan (adapt)
widens each block's windows, capped at `cap` classes a side (after
the third run, stopped at P = 1e8 on a window grown past it): the
coarse Bad classes past the cap are counted into the outer Bad, the
farthest printed."""
ts, ws = nodes(c, T, n_tau)
edges = np.array([1 / 16, 1 / 4, 1.0])
bands = np.zeros(4)
bad_outer, t_outer, widest = 0, 0.0, 0
R_quad.far = (0.0, 0.0)
scan = scan_classes(P, coarse)
for i in range(0, len(ts), block):
t, wt = ts[i:i + block], ws[i:i + block]
w0 = wh = w
if adaptive and t.min() < 1 / 16:
w0, wh, n_past, far = adapt(t, P, cst, w, scan, cap)
bad_outer += n_past
if far > R_quad.far[0]:
R_quad.far = (far, float(t.min()))
rho, slack = classes(P, w0, wh)
widest = max(widest, len(rho))
outer = slack < 10
sub = max(1, int(160000 // len(rho)))
for j in range(0, len(t), sub):
bad, rr = summand(t[j:j + sub], rho, P, cst)
bo = bad[:, outer].any(axis=1)
bad_outer += int(bad[:, outer].sum())
if bo.any():
t_outer = max(t_outer, float(t[j:j + sub][bo].max()))
band = np.searchsorted(edges, t[j:j + sub], side="right")
bands += np.bincount(band, weights=rr.sum(axis=1) *
wt[j:j + sub] * P * P, minlength=4)
R_quad.t_outer = t_outer
R_quad.widest = widest
return bands, bad_outer
def b_limit(x0, tau, cst, arm="P0"):
"""b(x0, tau) = lim P^3 B near x0: images below 1/4, modes from 1/64,
the smaller kept; arm "first" (after the first runs) puts the first
orders |l| and |beta1 - beta2| for A and A_k, the limit of the
majorants read at P0 = P."""
b = np.full(len(tau), np.inf)
lo = tau < 1 / 4
t = tau[lo]
bi = np.zeros(len(t))
for k in range(-12, 13):
xi = x0 + k
if xi == 0:
continue
G = (4 * np.pi * t) ** -0.5 * np.exp(-xi * xi / (4 * t)) * abs(xi) \
/ (2 * t)
a = abs(xi) / (2 * t)
if arm == "first":
A = np.abs(a * (1.5 - a * a * t))
else:
A = np.minimum(image_A(a, t, np.sign(xi), P0), 1e300)
bi += G * A
b[lo] = bi
hi = tau >= 1 / 64
if arm == "first":
t = tau[hi]
bf = np.zeros(len(t))
for k in range(1, K4 + 1):
bf += 4 * np.pi * k * np.exp(-4 * np.pi ** 2 * k * k * t) * \
np.abs(3 * np.pi * k - (2 * np.pi * k) ** 3 * t)
else:
bf = mode_B(tau[hi], 1.0, cst, tails=False)
b[hi] = np.minimum(b[hi], bf)
return b
def C_limit(odd, c, T, n_tau, cst, parity_only=False, arm="P0"):
"""C_R = int sum_(x0) S dtau, by band; and the odd term alone."""
ts, ws = nodes(c, T, n_tau)
edges = np.array([1 / 16, 1 / 4, 1.0])
band = np.searchsorted(edges, ts, side="right")
S = np.zeros(len(ts))
par = np.zeros(len(ts))
for x0 in (0.0, 0.5):
G = np.abs(kernel(x0, ts)[0])
if parity_only:
if x0 == 0.5 and odd:
par = G
continue
b = b_limit(x0, ts, cst, arm)
ratio = b / G
if x0 == 0.5 and odd:
nh = np.floor(ratio + 0.5)
S += G + 2 * b * (1 + nh + np.maximum(nh - 1, 0))
else:
S += 2 * b * (2 * np.floor(ratio) + 1)
if parity_only:
return float((par * ws).sum())
return np.bincount(band, weights=S * ws, minlength=4)
def K_half(tau):
k = np.arange(1, 400)
return 1 + 2 * float(((-1.0) ** k * np.exp(-4 * np.pi ** 2 * k * k * tau)
).sum())
def T_of(P):
return 1 + math.log(P / P0) / (2 * math.pi ** 2)
def section_control(cst):
print("X0 control at d = 400", flush=True)
P = 399
rng = np.random.default_rng(5)
m = np.unique(rng.integers(P ** 3 // C, P ** 3, 2000))
worst_i, worst_f = 0.0, 0.0
for lo, hi in ((P ** 3 // C, P ** 3 // 4), (P ** 3 // 64, P ** 3)):
mm = m[(m >= lo) & (m < hi)]
thi = float(mm.max() + 1) / P ** 3
if hi == P ** 3 // 4:
ref = images_maj(mm, P, P0, thi)
q, r = (mm // P).astype(float), (mm % P).astype(float)
got = img_B(q, r, P, thi)
worst_i = float(np.abs(got / ref - 1).max())
else:
tau = mm / P ** 3
m0 = int(mm.min())
ref = modes_maj(tau, P + 1, P, P0, cst, m0) - \
mode_tails(P + 1, P, K4, m0)
got = mode_B(tau, P, cst, tails=False)
worst_f = float(np.abs(got / ref - 1).max())
check("X0a B in (q, rho) equals majorant.py's (1e-12)",
max(worst_i, worst_f) <= 1e-12,
f"images {worst_i:.1e}, modes {worst_f:.1e} relative")
t0 = time.time()
b1, bo1 = R_quad(P, C, 1.0, 24000, cst)
b2, bo2 = R_quad(P, C, 1.0, 48000, cst)
tail = past_T(P, 1.0)
R1, R2 = b1.sum() + tail, b2.sum() + tail
ref = sum(BAND_399)
print(f" P = 399: by band [1/c, 1/16) [1/16, 1/4) [1/4, 1): "
f"{b1[0]:.5f} {b1[1]:.5f} {b1[2]:.4e}; R {R1:.5f} (P R {P * R1:.3f}),"
f" n_tau doubled {R2:.5f}; past P^3 {tail:.2e} "
f"({time.time() - t0:.0f} s)", flush=True)
dev = max(abs(b1[i] / BAND_399[i] - 1) for i in range(3))
check("X0c quadrature at P = 399 within 0.5% of majorant.py by band",
dev <= 0.005, f"largest band deviation {dev:.4f}; total "
f"{R1 / ref - 1:+.4f}")
check("X0c n_tau doubled moves R by under 0.1%",
abs(R2 / R1 - 1) < 1e-3, f"{R2 / R1 - 1:+.2e}")
check("X0b no Bad class in the outer 10 at P = 399", bo1 == 0,
f"{bo1} Bad (class, node) pairs, the largest such tau "
f"{R_quad.t_outer:.5f}")
_, bo = R_quad(801 - 2, C, 1.0, 6000, cst, w=20, adaptive=False)
check("X0d PLANT: the window cut to 20 fails X0b at d = 800", bo > 0,
f"{bo} Bad (class, node) pairs in the outer 10")
par = C_limit(True, C, 1.0, 24000, cst, parity_only=True)
want = K_half(1.0) - K_half(1 / C)
check("X0d the odd term integrates to K(1/2, 1) - K(1/2, 1/400) (1e-6)",
abs(par - want) <= 1e-6, f"{par:.8f} against {want:.8f}")
return R1
def section_scan(cst):
print("X0e the coarse scan against a dense scan, the blocks below "
"tau = 0.004", flush=True)
miss = {600: [], 6: []}
for P in (10001, 1000001, 10 ** 10 + 1):
t0 = time.time()
ts, _ = nodes(C, T_of(P), 12000)
scans = {k: scan_classes(P, k) for k in miss}
for i in range(0, len(ts), 256):
t = ts[i:i + 256]
if t.min() >= 0.004:
break
f0, fh = dense_far(t, P, cst, CAP)
row = [f" P = {P:.4g}, tau from {t.min():.5f}: dense farthest "
f"Bad {f0:.0f} from 0, {fh:.0f} from 1/2"]
for k, sc in scans.items():
w0, wh, _, _ = adapt(t, P, cst, W, sc, CAP)
miss[k].append(max(f0 - (w0 - 10), fh - (wh - 10)))
row.append(f"coarse {k}: windows {w0}, {wh}")
print("; ".join(row) + f" ({time.time() - t0:.0f} s)", flush=True)
check("X0e the adaptive windows hold every dense Bad class inside their "
"outer 10, P = 1e4, 1e6, 1e10", max(miss[600]) <= 0,
f"largest overreach {max(miss[600]):.0f} classes")
check("X0f PLANT: a 6-point coarse scan leaves a dense Bad class in or "
"past the outer 10", max(miss[6]) > 0,
f"largest overreach {max(miss[6]):.0f} classes")
def section_tail():
print("X1 the tail past M = P^3", flush=True)
rows = []
for P in (399, 10001, 1000001):
t = past_T(P, 1.0)
rows.append(t / P ** 2)
print(f" P = {P}: past P^3 {t:.3e}, over P^2 {t / P ** 2:.4e}; "
f"with T(P) = {T_of(P):.4f}: {past_T(P, T_of(P)):.3e}",
flush=True)
check("X1 past_M/P^2 within 10% of 1.46e-18 at P = 399, 1e4, 1e6",
all(abs(v / 1.46e-18 - 1) <= 0.10 for v in rows),
" ".join(f"{v / 1.46e-18:.3f}" for v in rows))
for P in (4.1e8 + 1, 1.8e10 + 1):
print(f" P = {P:.2e}: past P^3 {past_T(int(P), 1.0):.3e}", flush=True)
def section_limit(cst, R399):
print("X2, X3 the limit C_R against P R_maj(P), c = 400, T = 1",
flush=True)
Co = C_limit(True, C, 1.0, 24000, cst)
Ce = C_limit(False, C, 1.0, 24000, cst)
print(f" C_R odd {Co.sum():.4g} (by band {Co[0]:.4g} {Co[1]:.4f} "
f"{Co[2]:.4f}), even {Ce.sum():.4g}", flush=True)
C1 = C_limit(True, C, 1.0, 24000, cst, arm="first")
E1 = C_limit(False, C, 1.0, 24000, cst, arm="first")
print(f" first orders (A = |l|, A_k = |beta1 - beta2|): C_R odd "
f"{C1.sum():.4f} (by band {C1[0]:.4f} {C1[1]:.4f} {C1[2]:.4f}), "
f"even {E1.sum():.4f}", flush=True)
check("X2 |C_R(odd) - P R_maj(399)| <= 3, C_R above it",
0 <= Co.sum() - 399 * R399 <= 3,
f"C_R {Co.sum():.4f}, P R_maj(399) {399 * R399:.4f}")
gaps = []
for P in (799, 1599, 3199):
t0 = time.time()
b, bo = R_quad(P, C, 1.0, 12000, cst)
PR = P * (b.sum() + past_T(P, 1.0))
gaps.append(P * abs(PR - Co.sum()))
print(f" P = {P}: P R_maj {PR:.4f} (by band {P * b[0]:.4f} "
f"{P * b[1]:.4f} {P * b[2]:.4f}), P (P R - C_R) "
f"{P * (PR - Co.sum()):+.4g}, P R - C_R(first) "
f"{PR - C1.sum():+.4f}; outer Bad {bo}, widest read "
f"{R_quad.widest} classes ({time.time() - t0:.0f} s)",
flush=True)
check(f"X0b no Bad class in the outer 10 at P = {P}", bo == 0,
f"{bo}, the largest such tau {R_quad.t_outer:.5f}")
check("X3 P |P R_maj - C_R| at P = 799, 1599, 3199 within a factor 3",
max(gaps) <= 3 * min(gaps), " ".join(f"{g:.1f}" for g in gaps))
return Co
def section_kill(cst):
print("KILL: R_maj(P) + F(400, P), M = T(P) P^3", flush=True)
worst = 0.0
for P in (399, 799, 1599, 3199, 10001, 100001, 1000001, 10 ** 8 + 1,
10 ** 10 + 1, 18 * 10 ** 9 + 1):
t0 = time.time()
T = T_of(P)
b, bo = R_quad(P, C, T, 12000, cst)
tail = past_T(P, T)
R = b.sum() + tail
F = F_edge(C, P)[0]
Cl = C_limit(True, C, T, 24000, cst, arm="first").sum()
worst = max(worst, R + F)
print(f" P = {P:.4g}: T {T:.4f}, R_maj {R:.6f} (P R {P * R:.4f}, "
f"C_R(T, first) {Cl:.4f}), past T {tail:.1e}, F {F:.2e}, R + F "
f"{R + F:.6f}; outer Bad {bo}, widest read {R_quad.widest}, "
f"farthest Bad past the cap {R_quad.far[0]:.0f} classes at tau "
f"{R_quad.far[1]:.5f} ({time.time() - t0:.0f} s)",
flush=True)
if P > 3199:
check(f"X0b no Bad class in the outer 10 at P = {P:.4g}", bo == 0,
f"{bo}, the largest such tau {R_quad.t_outer:.5f}")
check("KILL not fired: max R_maj + F below 1/4", worst < 0.25,
f"{worst:.6f}")
def main():
cst = mode_consts(P0)
R399 = section_control(cst)
section_scan(cst)
if not all(ok for n, ok in zip(NAMES, CHECKS) if n not in MISSED):
print("control failed; stopping")
sys.exit(1)
section_tail()
section_limit(cst, R399)
section_kill(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()