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