""" 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/399)))) + 2, majorant.py's at 399. 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 # noqa: E402 from remainder import F_edge # 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): """K grown by ln(P/399), the module's P0, whatever P0 the majorants are read at: at P0 = P a growth of ln(P/P0) = 0 would leave the walk's tail past K, about P^3 e^-64, growing like P^3.""" return int(math.ceil(math.sqrt(4 * tau_hi * (64 + 3 * math.log(P / P0)))) ) + 2 def img_B(q, rho, P, tau_hi, P0=P0): """images_maj at m = qP + rho, q float, rho signed, K grown with P, A read at P0.""" 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, P0=P0): """modes_maj's per-mode sum at each tau, its tails at tau itself; cst and A_k read at P0.""" 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, P0=P0): """At nodes t and signed classes rho: Bad, and (|phi| - s phi + 2B) where Bad, the majorants read at P0.""" 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()), P0) hi = t >= 1 / 64 if hi.any(): # x-free, read at the node's tau Bf = mode_B(t[hi], P, cst, P0=P0) 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, P0=P0): """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, P0)[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, P0=P0): """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, P0)[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, P0=P0): """R_maj over [1/c, T) by the midpoint rule over windowed classes, by band, the majorants read at P0 (cst read there too), 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, P0) 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, P0) 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()