box.py
Checks: READING: The carry bound rule.
Cited on: READING.
To run it: save box.py, with edge.py, every.py, heat.py, majorant.py, real_support.py, remainder.py, uniform.py beside it; needs mpmath, numpy, scipy.
"""
box.py -- R_maj + F + X against the real support's budget on boxes
[Pa, Pb] of P from 399 to 1.3388e30: per box an upper bound U on P R_maj
at EVERY P of the box, from sups of monotone pieces over (offset, tau
cell), the majorants read at P0 = Pa.
QUESTION. At d = P + 1 = 0 mod 4 the real support's margin is positive
when R_maj(P) + F(c, P) + X(P) < budget(P) = 1/4 + 1/(4P) - 1/P^2
(uniform.py's R_maj, remainder.py's F_edge, variation.py's X; every.py's
budget). P R_maj read at grid P (every.py) bounds it at those P alone.
So P is cut in boxes, and per box the margin
m(Pa, Pb) = budget(Pb) - F_box(c, Pa, Pb) - X - U(Pa, Pb)/Pa
is read: budget falls past P = 8, U/Pa >= R_maj(P) at every P of the
box, F_box >= F(c, P) there (D10), X <= 2.06e-4 at every P >= 399
(variation.py: run to 1e6, past it by its (Q5) by hand). Is every box's
margin above 0? Scope: P in [399, 399 2^25 10^20) = [399, 1.3388e30),
under the majorants' own conditions (majorant.py's W1, W2), D6's and
D8's monotonicities (not proved), and doubles.
BY HAND.
(D1) The q sum. P R_maj(P) = sum_m [Bad](|phi| - s phi + 2B) P over m
in [P^3/c, T(P) P^3), m = qP + r, plus P past_T. Per class rho the
tau = (q + rho/P)/P^2 run on a lattice of step P^-2, so a tau cell
[t1, t2) holds at most (t2 - t1)P^2 + 1 of them: P R_maj <= sum over
cells (t2 - t1 + Pa^-2) sum_rho sup_cell P^3 term + P past_T, the cells
geometric on [1/c, T(Pb)], cut at 1/64, 1/4 and 1 (B's arms switch
at 1/64 and 1/4: images below 1/4, modes from 1/64, the smaller
kept).
(D2) The classes. rho = x0 P + o, x0 in {0, 1/2}: o in Z, |o| < P/4 at
0; at 1/2, o in Z at even P, Z + 1/2 at odd, |o| <= P/4. A box holds
both parities: each lattice summed over o, the larger lattice at 1/2
kept (every.py's combine). u = o/P lies in [o/Pb, o/Pa]: the shift
is the only way P enters an image's xi = x0 + k + u. Each offset's
u range is cut in sub-boxes at most DU = 0.005 wide.
(D3) phi in P^3 units: P^3 phi = -P K_x(x0 + u, tau). A right-signed
class (s phi > 0) is Bad only if P |K_x| <= P^3 B, and is cut on
a (sub-box, tau cell) when a lower bound on P |K_x| exceeds the sup
UB of P^3 B there; else it adds 2 UB. Two lower bounds, the larger
kept. Direct: Plo inf |K_x|, Plo = max(Pa, 4|o|, |o|/max|u|) (a
class at offset o lives at P >= 4|o|, and at P = |o|/|u| for its u
in the sub-box), inf |K_x| by a Taylor at the sub-box's centre,
|K_x| - |K_xx| du/2 - |K_xxx| dt/2, less the second order by local
sups of K's derivatives (K's Gaussian images n = -3..3 below tau =
0.05, Hermite coefficients bounding each derivative; the Fourier sum
sum_k 2 (2 pi k)^n e^(-4 pi^2 k^2 t1) above), and |u|min times mv
beside it. Per class: K_x(x0) = 0, so where mv > 0 bounds |K_xx|
below on the class's segment [x0, x0 + u], P |K_x(x0 + o/P)| >= P
|u| mv = |o| mv, at every P. The class o = 0 (s phi = 0) adds 2 UB.
The one wrong-signed class (x0 = 1/2, odd P, o = -1/2, s = -1, K_x <
0) adds 2 UB + 2 P |K_x(1/2 - 1/(2P))| <= 2 UB + sup |K_xx| over [1/2
- 1/(2Pa), 1/2] by the mean value (K_x(1/2) = 0).
(D4) Images, core (a = |xi|/(2 tau) <= P/4): P^3 |phi_k| A/P = G(xi, tau)
A(a, tau, s, Pa), G = (4 pi tau)^(-1/2) (|xi|/(2 tau)) e^(-xi^2/(4
tau)). Sup of G over |xi| in [l, h], tau in [t1, t2]: (4 pi t1)^(-1/2)
(2 t1)^-1 times max |xi| e^(-xi^2/(4 t2)), at |xi| = clip(sqrt(2 t2),
l, h), read in logs. A by its pieces (majorant.py image_A): a in
[l/(2 t2), h/(2 t1)]; lo = s a(3/2 - a^2 tau) exactly ranged (falls
in tau, one maximum in a at a^2 = 1/(2 tau)); Hbar's two added terms
at a_hi, its two subtracted at a_lo, t1, db at a_hi; Lbar's terms at
a_hi, t2 and 1/(6 tau) at t1; U rises in lo and Hbar, V in lo and
Lbar; A <= max(f(U_hi+), f(V_lo-)), f falling on t < 0 and rising on
t >= 0, in logs: ln f(U) <= ln P0 + U/P0 past U/P0 = 700, so nothing
is clamped. An image core at every P of the box (a_hi <= Pa/4)
takes G A; one non-core at every P (a_lo > Pb/4) the crude arm; one
mixed, the larger.
(D5) Images, non-core (|xi| > tau P/2). With j = P xi, n = 2 tau P^2 -
(P + 2) xi: n(1 - y) = n - |j| >= 2, so Tb <= (4 pi)^(-1/2) e^(-E),
E = min(xi^2/(4 tau), P^2 |xi|/(10P + 4)) (xi >= 0: n <= 2 tau P^2;
xi < 0: n < (5P + 2)|xi| by |xi| > tau P/2); the end terms 0.5 2^(-P
|xi|). So P^3 crude <= P^3 ((4 pi)^(-1/2) e^(-E) + 0.5 2^(-P |xi|)) +
P G, P^3 and P G at P_eff = min(Pb, 2h/t1) (non-core needs P < 2|xi|
/tau), E at max(l, t1 Pa/2), t2 and Pa. Past P = 2(K + 1)c every
image is core (|xi| <= K + 1 < tau P/2).
(D6) The walk and the images past K. Past |j| = J the walk's terms
cost at most e^(-g(J))/(1 - e^(-P g'(J))) a side, g(a) = a^2/(2(n0 +
B a)), n0 = 2 tau P^2, B = 1 + 2/P (remainder.py, Hoeffding): g is
convex and falls in n0, so the bound holds at any J and with tau
read at any tau_hi >= tau, which is how uniform.py's img_B is used
here: K = kimg(Pb, 1/4) at every cell (its g >= 60 a guard, not a
condition), g(K P) rising in P (W1), P^3 e^(-g) at Pb^3 and g at Pa,
the denominator at the smaller of P g' at Pa, Pb (P g' not proved
monotone); the N sum, P (4 pi tau)^(-1/2) ..., at Pb, its prefactor
at t1 and e^(-N^2/(4 tau)) at t2. Both below 1e-15 a cell by g >=
64 + 3 ln(P/399), so P^3 e^(-g) <= 399^3 e^-64 at every P.
(D7) Modes k <= 22: P^3 env_k A_k/P = 4 pi k e^(-4 pi^2 k^2 tau) A_k(tau;
Pa), P-free: env at t1; A_k's sqrt((b1 - b2)^2 + b1 b2^3/(3 Pa^2))
with (b1 - b2)^2 the larger at t1, t2 and b2^3 at t2, the rest rising
in tau, at t2. Modes 23..421: P sum 4 pi k e^(-4 pi^2 k^2 t1) at Pb.
The roots past 22: sin(pi k/d) >= 2k/d and ln(1 + y) >= y ln3/2 on
[0, 2] give the exponent >= 8.72 tau k^2 (k <= d/4) and >= 0.546 tau
P^2 (past), so P^3 root_tail <= 4 Pb^2 (sum_(k >= 23) e^(-8.72 t1
k^2) + (Pb/4 + 1) e^(-0.546 t1 Pa^2)): 4 Pb^2 e^-72.1 at t1 = 1/64,
2e-5 at 1e13, passing 1 near 2.2e15 (by hand), where the images'
smaller arm is kept below 1/4, and about 1e-440 at t1 >= 1/4 to P =
1.34e30 (it underflows to 0 in doubles), passing 1 near P = e^576
at t1 = 1/4, k = 23 (by hand). The range stops at 1.3388e30 by
choice.
(D8) Past T: P past_T(P, T(P)) grows about like P (past_T near 2.3e-13,
uniform.py): read at Pb, its rise in P not proved.
(D9) Boxes and the window. The layout (REGIONS): [399, 798) in 16
boxes of ratio 2^(1/16), [798, 1596) in 8 of 2^(1/8), [1596, 3192)
in 4 of 2^(1/4), [3192, 6384) in 2 of 2^(1/2), [6384, 399 2^25) in
3 of 128, each at 500 tau cells, and [399 2^25, 1.3388e30) in 5 of
1e4 at 1000, floors integer: 38 boxes. A box's cost grows with its
sub-boxes, o(1/Pa - 1/Pb)/DU an offset, so below 6384 they are
narrow; past it, to 399 2^25, W/Pa is small and the ratio nearly
free. Offsets |o| <= W = 200; past W, per u-cell [u1, u2] of [W/Pb,
1/4] (geometric to 0.02, DU wide above) on each of the four sides
(x0, sign of u), at most (u2 - u1)Pb + 1 classes at any P of the
box (a lattice of step 1 in [u1 P, u2 P]), each right-signed, adding
2 UB a class where its cut fails at that tau cell, the cut read at P
>= W/max|u| and with W mv. At fixed uncut cells this past-W add
over Pa grows like the ratio Pb/Pa at one Pb; at one ratio it
grows with Pb, by a split this run does not print (the uncut cells,
and UB's walk and tail terms read at Pb), and it falls as the tau
cells narrow: so the open range takes ratio 1e4 at 1000 cells.
(D10) c per box. F(c, P) = 4 e^(-g(a))/(1 - e^(-g'(a))), a = J0 >= (P -
1)/2, A = 2P^2/c, B = 1 + 2/P (remainder.py's F_edge); F falls in a.
Along a = (P - 1)/2, g = (P - 1)^2/(8D), D = 2P^2/c + (P - 1)(P +
2)/(2P), and 2D - (P - 1)D' = 4P/c + (P - 1)(1/2 + 2/P - 1/P^2) > 0,
so g rises in P: read at Pa. g' = (a/A)(2 + x)/(2(1 + x)^2) >=
(a/A)/(1 + x)^2, x = Ba/A, a/A = c(P - 1)/(4P^2) falling in P (at
Pb), x <= (1 + 2/Pa)c/(4Pa). So F(c, P) <= F_box(c, Pa, Pb) on the
box; c is the least integer >= 400 with F_box <= 0.01, and the cells
start at 1/c.
PREDICTIONS (frozen before the run).
(P1) U(399, 399), the degenerate box, over every.py's 25.3379 in
[1.15, 1.30].
(P2) The largest U below 399 2^25 at range 0's first box [399, 416],
in [38, 55], and the least margin there, in [0.11, 0.16].
(P3) The past-W add at most 5 a box below 399 2^25.
(P4) The open boxes' margins at least 0.237, their least in [0.2398,
0.2400], their past-W add over Pa at most 1e-6.
(P5) The wall at most 600 s.
CONTROLS.
(D0) The dominance: at P in a box, at the centre of every third tau
cell wide enough (2.5/P^2) to hold its classes' tau, the bound's
cell sum over the classes against the true sum (uniform.py's
summand), both with the majorants at P0 = Pa: the largest ratio at
most 1. Read at P = 399 (Pa = Pb) and at five P inside four boxes,
830, 869, 1777, 600001 and 3e10 + 1, the last two over the classes
within 400 of each sign change (every summand is nonnegative, so
that can only lower the truth). Sharing P0, the control cannot see
the majorants' validity at P > Pa, which rests on W1, W2.
PLANTS. The wrong-signed class's Gamma dropped: a ratio past 1 at
every odd control P. |o| read as 3|o| in the per-class cut (D3):
past 1 at 399 or 869. The last u-cell dropped (D9): the coverage
(the u-cells gapless from W/Pb to 1/4, each lattice's count at
Pa, the middle and Pb - 1 within (u2 - u1)Pb + 1, in rationals,
since past 2^53 a double's spacing passes 1) leaves a gap.
(D10) F_box at or above F_edge at Pa, the middle and Pb of every box,
to 1e-13.
The past-W cut (W mv) has no control: the truth past W is 0 at every
control P.
KILL. A box's margin at or below 0, or a control ratio above 1. (K2)
D0's U below 25.3379: a bug.
FINDINGS (copied from the run's prints).
- The kill misses: every box's margin is above 0, the least 0.142075
at [399, 416]. Per region, U max, least margin, past-W add max:
range 0, 43.2173 at [399, 416], 0.142075, 0; range 1, 26.5044 at
[798, 870], 0.216867, 0.1357; range 2, 22.7682, 0.235660, 0.0985;
range 3, 22.0612, 0.242938, 0.00272; ranges 4..24, 21.9922 at
[6384, 817152], 0.242356 at [104595456, 13388218368], 0.0001115;
open, c 538..1116, U max 4.28979e17 at the last box, U/Pa at most
3.83e-9, least margin 0.239865 at the first, past-W add 1.199e17,
over Pa at most 8.952e-10. So the real support's margin at d = P +
1 = 0 mod 4 is positive at every P of [399, 1.3388e30), under the
QUESTION's conditions.
- The controls hold: the dominance 0.9577 at 399, 0.9608 and 0.9618
at 830, 869 in [798, 870], 0.9655 at 1777 in [1596, 1897], 0.9674
at 600001 in [6384, 817152], 0.9821 at 3e10 + 1 in [13388218368,
133882183680000]; the Gamma plant 1.9798 at 399 and 2.0932,
2.1368, 2.1548, 2.3572 at the odd P; the |o| plant 3.5494, 3.3740;
the coverage over 22 boxes, least slack 0.013, its plant leaving a
gap; F_box over F_edge least 1 to fifteen digits.
- P1 misses: U(399, 399) = 33.1196, ratio 1.3071 (31.2886, 1.2349 at
1000 cells: the +Pa^-2 a cell and the cells' sups). P2..P5 hold.
- How the layout was found, from reads at 1000 cells by this design.
Below 6384 a box's cost grows with its sub-boxes: range 0's first
box at ratio 2^(1/16), 2^(1/4), 2^(1/2) reads U 40.8616, 55.0896,
74.3924, and at 2 its arrays passed 512 MB (before a pass's rows
were scaled to its columns). Past it the ratio is
nearly free: [6384, Pb] at ratio 2, 8, 64 reads 20.4102, 20.4864,
20.4857, and [1634304, Pb] at 2, 64, 1024 reads 20.1789, 20.2094,
20.2194. The per-class cut (D3's |o| mv) took the ratio-2 boxes from
6384 and 408576 from 38.69 and 38.08 to 20.4102 and 20.1842: the
mean-value arm read on a sub-box loses Pa/Pb. W = 200 against
every class to 1600 at [798, 870]: U 24.2051 both. Past 399 2^25
the past-W add over Pa grows with Pb at one ratio and like the ratio
at one Pb (D9): at ratio 1e10 and 500 cells the box [1.3388e20,
1.3388e30] reads margin -7.732848, the add 1.067e21 all of U (add/Pa
0.0133 at the first open box, 7.97 there); at 1000 cells the last
box at ratio 1e4 and at 10 reads add/Pa 8.952e-10 and 8.952e-13,
margin 0.239996 at both, and at ratio 1e4 the first box reads
4.793e-40.
MISSED: P1, by name in MISSED; the run exits 0 iff the checks that fail
are exactly it.
RUN RECORD: 19 of 20 checks, the miss P1 as recorded; 453.4 s, peak
commit 135.2 MB under the 512 MB guard (estimate 460 s): the controls
and coverage 130 s, ranges 0..24 211 s, the open range 112 s. A
rehearsal, the first box of each region, ran first: 193.5 s, peak 135.2
MB, every check as in the run.
"""
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 uniform as UN # noqa: E402
from every import at_P, T_of, K4, budget # noqa: E402
from heat import kernel # noqa: E402
from remainder import F_edge # noqa: E402
C, W = 400, 200
DU = 0.005 # (D2) a sub-box's width in u, at most
CHUNK, ELEMS = 250, 10 ** 6 # tau cells a pass, and cells x columns
P_OPEN = 399 * 2 ** 25 # 13388218368: the open range's start
REGIONS = [("range 0", 399, 2, 1 / 16, 16, 500),
("range 1", 798, 2, 1 / 8, 8, 500),
("range 2", 1596, 2, 1 / 4, 4, 500),
("range 3", 3192, 2, 1 / 2, 2, 500),
("ranges 4..24", 6384, 2, 7, 3, 500),
("open", P_OPEN, 10, 4, 5, 1000)] # (D9) start, ratio base^step,
# boxes, tau cells
INSIDE = [(798, 870, 830, None), (798, 870, 869, None),
(1596, 1897, 1777, None), (6384, 817152, 600001, 400),
(P_OPEN, P_OPEN * 10 ** 4, 30000000001, 400)] # the control at P inside boxes
R399 = 25.3379 # every.py's P R_maj at P = 399, P0 = P
CHECKS, NAMES = [], []
MISSED = ("P1 U(399, 399) over 25.3379 in [1.15, 1.30]",) # the record's misses
REHEARSE = "--rehearse" in sys.argv
PLANT = {"gamma": True, "edge": True}
PW = {} # (D9) per box: the past-W add, uncut count
OF = [1.0] # (D3) |o|'s factor in the cut; its PLANT 3
F_MAX, X_MAX = 0.01, 2.06e-4 # (D10) F_box's cap; X (variation.py)
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 cells(T, n, c=C):
"""Geometric tau cells on [1/c, T], cut at 1/64, 1/4 and 1: (t1, t2)."""
cuts = [1 / c] + [v for v in (1 / 64, 1 / 4, 1.0) if 1 / c < v < T] + [T]
span = math.log(T * c)
e1, e2 = [], []
for lo, hi in zip(cuts[:-1], cuts[1:]):
m = max(8, int(n * math.log(hi / lo) / span))
e = np.geomspace(lo, hi, m + 1)
e1.append(e[:-1])
e2.append(e[1:])
return np.concatenate(e1), np.concatenate(e2)
def dsup(t, n):
"""sup over x of |d_x^n K(x, t)| <= sum_k 2 (2 pi k)^n e^(-4 pi^2 k^2 t)."""
k = np.arange(1, 200)
w = 2 * np.pi * k
return (2 * w ** n * np.exp(-np.outer(t, w * w))).sum(axis=1)
def G_sup(l, h, t1, t2):
"""(D4) sup of (4 pi tau)^(-1/2) (|xi|/(2 tau)) e^(-xi^2/(4 tau)) over
|xi| in [l, h], tau in [t1, t2]."""
z = np.clip(np.sqrt(2 * t2), l, h)
return (4 * np.pi * t1) ** -0.5 / (2 * t1) * z * np.exp(-z * z / (4 * t2))
def A_sup(al, ah, s, t1, t2, P0):
"""(D4) ln of the sup of majorant.py's image_A (uniform arm) over a in
[al, ah], tau in [t1, t2], sign s, at P0: ln fU <= ln P0 + U/P0 past
U/P0 = 700, so no clamp understates it."""
e0 = 1 + 2 / P0
rr = np.minimum(ah / P0, 0.25)
db = e0 * rr
yb = rr / (1 - db)
a_pk = np.clip(1 / np.sqrt(2 * t1), al, ah)
mx = a_pk * (1.5 - a_pk ** 2 * t1) # max of a(3/2 - a^2 tau)
mn = np.minimum(al * (1.5 - al ** 2 * t2), ah * (1.5 - ah ** 2 * t2))
lo_hi = np.where(s > 0, mx, -mn)
lo_lo = np.where(s > 0, mn, -mx)
Hb = (0.75 * ah ** 2 * e0 ** 2 / (1 - db)
+ ah ** 2 / (2 * (1 - db) ** 2 * (1 - yb ** 2))
- al ** 4 * t1 / (1 + db) - al ** 4 * t1 / (6 * (1 + db) ** 3))
Lb = (-ah ** 4 * t2 * e0 ** 2 / (1 - db)
- ah ** 4 * t2 / (6 * (1 - db) ** 3 * (1 - yb ** 2))
- 1 / (6 * t1 * (1 - db) * (1 - yb ** 2)))
U = lo_hi + np.maximum(Hb + np.maximum(2 * lo_hi, 0), 0) / P0
V = lo_lo + np.minimum(Lb + np.minimum(2 * lo_lo, 0), 0) / P0
x = np.maximum(U, 0) / P0
with np.errstate(divide="ignore"):
lnU = math.log(P0) + np.where(x <= 700, np.log(np.expm1(
np.minimum(x, 700))), x)
lnV = np.log(-np.minimum(V, 0))
return np.maximum(lnU, lnV)
def lnG_sup(l, h, t1, t2):
"""ln of G_sup, in closed form (no underflow)."""
z = np.clip(np.sqrt(2 * t2), l, h)
with np.errstate(divide="ignore"):
return (-0.5 * np.log(4 * np.pi * t1) - np.log(2 * t1) + np.log(z)
- z * z / (4 * t2))
def img_sup(x0, ulo, uhi, Pa, Pb, t1, t2):
"""(D4)-(D6) sup of P^3 img_B over P in [Pa, Pb], u in [ulo, uhi] (per
column), tau in [t1, t2] (per row)."""
K = UN.kimg(Pb, 0.25)
T1, T2 = t1[:, None], t2[:, None]
out = np.zeros((len(t1), len(ulo)))
for k in range(-K, K + 1):
xa, xb = x0 + k + ulo, x0 + k + uhi
s = np.sign(xa + xb)
l = np.where(xa * xb <= 0, 0.0, np.minimum(np.abs(xa), np.abs(xb)))
h = np.maximum(np.abs(xa), np.abs(xb))
l2, h2 = l[None, :], h[None, :]
al, ah = l2 / (2 * T2), h2 / (2 * T1)
lnA = A_sup(al, np.minimum(ah, Pb / 4), s[None, :], T1, T2, Pa)
lnc = lnG_sup(l2, h2, T1, T2) + lnA
with np.errstate(over="ignore"):
core = np.exp(lnc)
Pe = np.minimum(Pb, 2 * h2 / T1)
ln = np.maximum(l2, T1 * Pa / 2) # non-core: |xi| > tau P/2
E = np.minimum(ln * ln / (4 * T2), Pa ** 2 * ln / (10 * Pa + 4))
crude = np.where(ln <= h2, Pe ** 3 * (
(4 * np.pi) ** -0.5 * np.exp(-E) + 0.5 * np.exp(-Pa * ln * math.log(2)))
+ Pe * G_sup(ln, h2, T1, T2), 0.0)
out += np.where(ah <= Pa / 4, core,
np.where(al > Pb / 4, crude, np.maximum(core, crude)))
# (D6) the walk past K and the N sum
walk = np.zeros(len(t1))
for i, t in enumerate(t2):
gs, xs = [], []
for P in (Pa, Pb):
J = K * P
n0 = 2 * t * P * P
Bw = 1 + 2 / P
gs.append(J * J / (2 * (n0 + Bw * J)))
xs.append(P * J * (2 * n0 + Bw * J) / (2 * (n0 + Bw * J) ** 2))
assert min(gs) >= 60, gs
walk[i] = Pb ** 3 * math.exp(-gs[0]) / (-math.expm1(-min(xs)))
N = np.arange(K, K + 40)
ntail = Pb * 2 * (4 * np.pi * t1) ** -0.5 * (
(N + 1)[None, :] / (2 * T1) * np.exp(-(N * N)[None, :] / (4 * T2))
).sum(axis=1)
return out + (walk + ntail)[:, None]
def mode_sup(Pa, Pb, t1, t2):
"""(D7) sup of P^3 mode_B over P in [Pa, Pb], tau in [t1, t2]."""
out = np.zeros(len(t1))
for k, th0, sg, mb, Ec, zc in at_P(Pa):
lnenv = math.log(4 * math.pi * k) - 4 * np.pi ** 2 * k * k * t1
b1 = 3 * np.pi * k
b2a, b2b = (2 * np.pi * k) ** 3 * t1, (2 * np.pi * k) ** 3 * t2
S = np.sqrt(np.maximum((b1 - b2a) ** 2, (b1 - b2b) ** 2)
+ b1 * b2b ** 3 / (3 * Pa * Pa))
z = zc * t2
cz = math.sqrt(1 + 2.25 * th0 * th0) * Pa
lnA = np.where(z < 600, np.log(S + Ec * np.exp(np.minimum(z, 600))
+ cz * np.expm1(np.minimum(z, 600))),
np.log(2) + np.maximum(np.log(S), math.log(Ec + cz) + z))
out += np.exp(lnenv + lnA)
ks = np.arange(K4 + 1, K4 + 400)
out += Pb * (4 * np.pi * ks * np.exp(-4 * np.pi ** 2 * np.outer(t1, ks * ks))
).sum(axis=1)
kr = np.arange(K4 + 1, 4000)
out += 4 * Pb ** 2 * (np.exp(-8.72 * np.outer(t1, kr * kr)).sum(axis=1)
+ (Pb / 4 + 1) * np.exp(-0.546 * t1 * Pa ** 2))
return out
def kder(x, t, n):
"""d_x^n K(x, t) by its Fourier series, t >= 1/1200 (k < 60)."""
out = np.zeros(np.broadcast(x, t).shape)
for k in range(1, 60):
w = 2 * np.pi * k
ph = w * x + n * np.pi / 2
out += 2 * w ** n * np.cos(ph) * np.exp(-w * w * t)
return out
HERM = [np.polynomial.hermite.herm2poly([0] * m + [1]) for m in range(8)]
def loc_sup(xl, xh, t1, t2, m):
"""sup of |d_x^m K| over x in [xl, xh], tau in [t1, t2]: below tau =
0.05 by K's Gaussian images n = -3..3, each d_y^m of (4 pi tau)^(-1/2)
e^(-y^2/(4 tau)) = (4 pi tau)^(-1/2) (-1)^m (4 tau)^(-m/2) H_m(y/(2
sqrt tau)) e^(-y^2/(4 tau)), |H_m| <= its coefficients' absolute sum
at |z|; at and above 0.05 by dsup at t1."""
shape = np.broadcast(xl, xh, t1).shape
T1 = np.broadcast_to(t1, shape)
T2 = np.broadcast_to(t2, shape)
lo = T2 < 0.05
c = np.abs(HERM[m])
acc = np.zeros(shape)
for n in range(-3, 4):
a, b = xl - n, xh - n
l = np.where(a * b <= 0, 0.0, np.minimum(np.abs(a), np.abs(b)))
h = np.maximum(np.abs(a), np.abs(b))
z = h / (2 * np.sqrt(T1))
poly = sum(ck * z ** k for k, ck in enumerate(c))
acc = acc + ((4 * np.pi * T1) ** -0.5 * (4 * T1) ** (-m / 2) * poly
* np.exp(-l * l / (4 * T2)))
hi = np.broadcast_to(dsup(np.ravel(T1[..., :1] if T1.ndim > 1 else T1),
m).reshape(T1.shape[:-1] + (1,) if T1.ndim > 1
else T1.shape), shape)
return np.where(lo, acc, hi)
def kx_inf(x0, ulo, uhi, t1, t2):
"""(D3) inf of |K_x(x0 + u, tau)| over the box, the larger of two
bounds, and mv, the mean value's inf of |K_xx| (P |K_x| >= |o| mv).
Direct: Taylor at the box's centre, |K_x| - |K_xx| du/2 - |K_xxx|
dt/2 there, less the second order by loc_sup over the box.
Mean value from x0 (K_x(x0) = 0): |u|min times inf |K_xx| over x
between x0 and x0 + u, by Taylor at that segment's centre alike."""
T1, T2 = t1[:, None], t2[:, None]
xl, xh = x0 + ulo[None, :], x0 + uhi[None, :]
xm = (xl + xh) / 2
tm = (T1 + T2) / 2
du = xh - xl
dt = T2 - T1
first = (np.abs(kder(xm, tm, 1)) - np.abs(kder(xm, tm, 2)) * du / 2
- np.abs(kder(xm, tm, 3)) * dt / 2)
sec = (du ** 2 / 8 * loc_sup(xl, xh, T1, T2, 3)
+ du * dt / 4 * loc_sup(xl, xh, T1, T2, 4)
+ dt ** 2 / 8 * loc_sup(xl, xh, T1, T2, 5))
um = np.minimum(np.abs(ulo), np.abs(uhi))[None, :]
um = np.where((ulo * uhi)[None, :] <= 0, 0.0, um)
uf = np.maximum(np.abs(ulo), np.abs(uhi))[None, :]
sg = np.sign(ulo + uhi)[None, :]
sl, sh = np.minimum(x0, x0 + sg * uf), np.maximum(x0, x0 + sg * uf)
xc = (sl + sh) / 2
mv = (np.abs(kder(xc, tm, 2)) - np.abs(kder(xc, tm, 3)) * uf / 2
- np.abs(kder(xc, tm, 4)) * dt / 2
- (uf ** 2 / 8 * loc_sup(sl, sh, T1, T2, 4)
+ uf * dt / 4 * loc_sup(sl, sh, T1, T2, 5)
+ dt ** 2 / 8 * loc_sup(sl, sh, T1, T2, 6)))
return np.maximum(np.maximum(first - sec, um * mv), 0.0), np.maximum(mv, 0.0)
def gamma_sup(Pa, t1, t2):
"""(D3) sup of |K_xx| over x in [1/2 - 1/(2Pa), 1/2], tau in the cell."""
x1 = 0.5 - 0.5 / Pa
cs = [np.abs(kernel(x, t)[0]) for x in (x1, 0.5) for t in (t1, t2)]
return (np.max(cs, axis=0) + 0.25 / Pa * dsup(t1, 3)
+ (t2 - t1) / 2 * dsup(t1, 4))
def lattice(x0, half, Pa, Pb):
"""Offsets |o| <= W of one lattice live at some P of the box, and
their u sub-boxes."""
lim = min(W, Pb / 4)
if half:
o = np.arange(-math.floor(lim + 0.5), math.floor(lim + 0.5)) + 0.5
else:
o = np.arange(-math.floor(lim), math.floor(lim) + 1, dtype=float)
o = o[np.abs(o) <= Pb / 4]
ulo = np.where(o >= 0, o / Pb, o / Pa)
uhi = np.where(o >= 0, o / Pa, o / Pb)
n = np.maximum(1, np.ceil((uhi - ulo) / DU)).astype(int)
start = np.concatenate([[0], np.cumsum(n)[:-1]])
i = np.repeat(np.arange(len(o)), n)
f = np.arange(n.sum()) - np.repeat(start, n)
w = (uhi - ulo)[i] / n[i]
return o, ulo[i] + f * w, ulo[i] + (f + 1) * w, i, start
def box_bound(Pa, Pb, n_cells, parities=(0, 1), per_cell=False):
"""U(Pa, Pb) and its parts; parities: the P parities the box holds;
per_cell: also the cells and each cell's sum over the classes."""
t1, t2 = cells(T_of(Pb), n_cells, c_box(Pa, Pb))
h = t2 - t1 + Pa ** -2.0
img = t2 <= 0.25 + 1e-15
mod = t1 >= 1 / 64 - 1e-15
UM = np.full(len(t1), np.inf)
UM[mod] = mode_sup(Pa, Pb, t1[mod], t2[mod])
lats = [("0", 0.0, False)]
if 0 in parities:
lats.append(("hZ", 0.5, False))
if 1 in parities:
lats.append(("hZ2", 0.5, True))
sums, live = {}, 0
for key, x0, half in lats:
o, ulo, uhi, gi, start = lattice(x0, half, Pa, Pb)
always = (o == 0) | (half & (o == -0.5))
tsum = np.zeros(len(t1))
rows = max(8, min(CHUNK, ELEMS // len(ulo)))
for c in range(0, len(t1), rows):
sl = slice(c, c + rows)
a1, a2, im = t1[sl], t2[sl], img[sl]
UB = np.broadcast_to(UM[sl, None], (len(a1), len(ulo))).copy()
if im.any():
UB[im] = np.minimum(UB[im], img_sup(x0, ulo, uhi, Pa, Pb,
a1[im], a2[im]))
um = np.maximum(np.abs(ulo), np.abs(uhi))
Plo = np.maximum(np.maximum(Pa, 4 * np.abs(o[gi])), np.divide(
np.abs(o[gi]), um, out=np.zeros(len(um)), where=um > 0))[None, :]
kxi, mv = kx_inf(x0, ulo, uhi, a1, a2)
cut = (Plo * kxi > UB) | (OF[0] * np.abs(o[gi])[None, :] * mv > UB)
cut &= ~always[gi][None, :]
term = np.maximum.reduceat(np.where(cut, 0.0, 2 * UB), start,
axis=1)
cut = np.logical_and.reduceat(cut, start, axis=1)
if half and PLANT["gamma"]:
k = int(np.flatnonzero(o == -0.5)[0])
term[:, k] += gamma_sup(Pa, a1, a2)
live += int((~cut).sum())
tsum[sl] = term.sum(axis=1)
sums[key] = tsum
half_sum = np.maximum(sums.get("hZ", 0), sums.get("hZ2", 0))
pw, n_pw = past_W(Pa, Pb, t1, t2, UM, img)
cell = sums["0"] + half_sum + pw
tot = (cell * h).sum()
T = T_of(Pb)
tot += Pb * UN.past_T(Pb, T)
PW[(Pa, Pb)] = (float((pw * h).sum()), n_pw)
if per_cell:
return float(tot), live, (t1, t2, cell)
return float(tot), live
def F_box(c, Pa, Pb):
"""(D10) F(c, P) at every P of [Pa, Pb]: g at a = (P - 1)/2, rising in
P, read at Pa; g' >= (a/A)/(1 + x)^2, a/A = c(P - 1)/(4P^2) least at
Pb, x = Ba/A <= (1 + 2/Pa) c/(4Pa)."""
a, A, B = (Pa - 1) / 2, 2 * Pa * Pa / c, 1 + 2 / Pa
g = a * a / (2 * (A + B * a))
dg = c * (Pb - 1) / (4 * Pb * Pb) / (1 + B * c / (4 * Pa)) ** 2
return 4 * math.exp(-g) / -math.expm1(-dg)
def c_box(Pa, Pb):
"""(D10) the least integer c >= C with F_box <= F_MAX."""
c = C
while F_box(c, Pa, Pb) > F_MAX:
c += 1
return c
def u_edges(Pb):
"""(D9) the u-cells' edges on [W/Pb, 1/4], one side."""
lo = W / Pb
e = (np.geomspace(lo, 0.02, max(2, int(math.log(0.02 / lo) / 0.08)) + 1)
if lo < 0.02 else np.array([lo]))
e = np.concatenate([e, np.arange(max(lo, 0.02) + DU, 0.25, DU), [0.25]])
return e[:-1] if PLANT["edge"] is False else e
def coverage(Pa, Pb):
"""(D9) the u-cells run without a gap from at most W/Pb to 1/4, and at P
= Pa, the middle and Pb - 1, each lattice's classes in a u-cell (o in
[u1 P, u2 P] on Z or Z + 1/2, counted in closed form) number at most
(u2 - u1)Pb + 1: the least slack (negative is a breach) and the
gap's verdict. Counted in rationals (each edge's double read exactly):
past P = 2^53 a double's spacing passes 1."""
e = u_edges(Pb)
whole = bool(e[0] <= W / Pb and e[-1] >= 0.25 and np.all(np.diff(e) > 0))
q = [Fraction(float(v)) for v in e]
worst = math.inf
for P in (Pa, (Pa + Pb) // 2, Pb - 1):
for c in (Fraction(0), Fraction(1, 2)):
for lo, hi in zip(q[:-1], q[1:]):
n = math.floor(hi * P - c) - math.ceil(lo * P - c) + 1
worst = min(worst, float((hi - lo) * Pb + 1 - n))
return worst, whole
def past_W(Pa, Pb, t1, t2, UM, img):
"""(D9) the classes past W, per tau cell: on each of the four sides (x0
in {0, 1/2}, sign of u) a u-cell [u1, u2] of [W/Pb, 1/4] holds at most
(u2 - u1)Pb + 1 classes at any P of the box, each right-signed, adding
2 sup P^3 B where the cut fails at that tau cell (P >= W/|u|). The
per-cell sum, and the count of (side, u-cell, tau cell) uncut."""
add = np.zeros(len(t1))
if W / Pb >= 0.25:
return add, 0
e = u_edges(Pb)
n_u = len(e) - 1
cnt = (e[1:] - e[:-1]) * Pb + 1
n_fail = 0
for x0 in (0.0, 0.5):
for sg in (1, -1):
a, b = sg * e[:-1], sg * e[1:]
ulo, uhi = np.minimum(a, b), np.maximum(a, b)
Plo = np.maximum(Pa, W / np.maximum(np.abs(ulo), np.abs(uhi)))
for c in range(0, len(t1), CHUNK):
sl = slice(c, c + CHUNK)
a1, a2, im = t1[sl], t2[sl], img[sl]
UB = np.broadcast_to(UM[sl, None], (len(a1), n_u)).copy()
if im.any():
UB[im] = np.minimum(UB[im], img_sup(x0, ulo, uhi, Pa, Pb,
a1[im], a2[im]))
kxi, mv = kx_inf(x0, ulo, uhi, a1, a2)
ok = (Plo[None, :] * kxi > UB) | (OF[0] * W * mv > UB)
n_fail += int((~ok).sum())
add[sl] += np.where(ok, 0.0, 2 * UB * cnt[None, :]).sum(axis=1)
return add, n_fail
def geo(lo, base, step, n):
"""n boxes of ratio base^step from lo, floors integer."""
fl = sorted({int(math.floor(lo * base ** (i * step))) for i in range(n + 1)})
return list(zip(fl[:-1], fl[1:]))
def layout():
"""(D9) the boxes, (region, Pa, Pb, tau cells), gapless from 399 to
1.3388e30."""
out = []
for reg, lo, base, step, n, nc in REGIONS:
out += [(reg, Pa, Pb, nc) for Pa, Pb in geo(lo, base, step, n)]
assert all(a[2] == b[1] for a, b in zip(out, out[1:])), "a gap"
return out
def cells_at(P):
"""(D9) the tau cells of the region holding P."""
return [r[5] for r in REGIONS if r[1] <= P][-1]
def dominance(n_cells, Pa=399, Pb=399, P=399, w=None):
"""At P in [Pa, Pb], the true sum over the classes (within w of each
sign change, every class when w is None) of P^3 [Bad](|phi| - s phi +
2B), uniform.py's summand with the majorants at P0 = Pa, at the centre
of every third cell wide enough (2.5/P^2) to hold its classes' tau,
over the box bound's cell sum: the largest ratio."""
par = (P % 2,) if Pa == Pb else (0, 1)
U, live, (t1, t2, S) = box_bound(Pa, Pb, n_cells, par, per_cell=True)
idx = np.flatnonzero((t2 - t1) >= 2.5 / P ** 2)[::3]
tm = (t1[idx] + t2[idx]) / 2
rho = UN.classes(P, P, P)[0] if w is None else UN.classes(P, w, w)[0]
cst = at_P(Pa)
far = np.minimum(np.abs(rho), np.abs(rho - P / 2)) > W
tru, tfar = np.zeros(len(tm)), np.zeros(len(tm))
for a in range(0, len(tm), 16):
v = P ** 3 * UN.summand(tm[a:a + 16], rho, P, cst, Pa)[1]
tru[a:a + 16] = v.sum(axis=1)
tfar[a:a + 16] = v[:, far].sum(axis=1)
r = tru / S[idx]
share = float(np.max(np.divide(tfar, tru, out=np.zeros(len(tm)),
where=tru > 0)))
return U, live, float(r.max()), float(tm[r.argmax()]), len(idx), share
def section_control(n_cells):
print("D0 CONTROL: the degenerate box Pa = Pb = 399", flush=True)
t0 = time.time()
U, live, r, at, n, _ = dominance(n_cells)
check("D0 the bound's cell sums dominate the true class sums at P = 399",
r <= 1, f"largest ratio {r:.4f} at tau {at:.4f} over {n} cells "
f"({time.time() - t0:.0f} s)")
PLANT["gamma"] = False
try:
Up, _, rp, atp, _, _ = dominance(n_cells)
finally:
PLANT["gamma"] = True
check("D0 PLANT: the wrong-signed class's Gamma dropped, a ratio past 1",
rp > 1, f"largest ratio {rp:.4f} at tau {atp:.4f}; U {Up:.4f}")
check("P1 U(399, 399) over 25.3379 in [1.15, 1.30]",
1.15 <= U / R399 <= 1.30, f"U {U:.4f}, ratio {U / R399:.4f}, "
f"live cells {live}")
check("KILL K2 not fired: D0's U at or above 25.3379", U >= R399,
f"{U:.4f}")
for Pa, Pb, P, w in INSIDE:
t0 = time.time()
_, _, r, at, n, share = dominance(cells_at(Pa), Pa, Pb, P, w)
msg = f"largest ratio {r:.4f} at tau {at:.4f} over {n} cells"
ok = r <= 1
if P % 2:
PLANT["gamma"] = False
try:
rp = dominance(cells_at(Pa), Pa, Pb, P, w)[2]
finally:
PLANT["gamma"] = True
msg += f"; gamma plant {rp:.4f}"
ok = ok and rp > 1
check(f"D0 inside the box [{Pa}, {Pb}] at P = {P}: dominates"
+ (", gamma plant past 1" if P % 2 else ""), ok,
msg + f" ({time.time() - t0:.0f} s)")
if Pb > 1600:
print(f" the truth's largest cell share from classes past W "
f"{share:.4g}", flush=True)
return U
def section_cut_plant(n_cells):
"""(D3) PLANT: |o| read as 3|o| in the per-class cut; a ratio past 1 at
some control P (399 and INSIDE's 869)."""
OF[0] = 3.0
try:
rs = [dominance(n_cells)[2], dominance(n_cells, *INSIDE[1])[2]]
finally:
OF[0] = 1.0
check("D3 PLANT: |o| read as 3|o|, a ratio past 1 at P = 399 or 869",
max(rs) > 1, f"largest ratios {rs[0]:.4f}, {rs[1]:.4f}")
def section_coverage(boxes):
print("D9 the u-cells past W hold every class past W", flush=True)
bx = [(Pa, Pb) for _, Pa, Pb, _ in boxes if W / Pb < 0.25]
res = [coverage(Pa, Pb) for Pa, Pb in bx]
worst, whole = min(r[0] for r in res), all(r[1] for r in res)
PLANT["edge"] = False
try:
pwhole = all(coverage(Pa, Pb)[1] for Pa, Pb in bx)
finally:
PLANT["edge"] = True
check("D9 the u-cells cover [W/Pb, 1/4], each count within (u2 - u1)Pb "
"+ 1; PLANT (the last u-cell dropped) leaves a gap",
whole and worst >= 0 and not pwhole,
f"{len(bx)} boxes, 3 P each: least slack {worst:.3f}; plant "
f"{'covers' if pwhole else 'leaves a gap'}")
def section_boxes(bx):
"""Per box: c, F_box against F_edge at Pa, the middle and Pb, U, the R
bound U/Pa and the margin budget(Pb) - F_box - X - U/Pa; per region
the extremes."""
print(f"D9 {len(bx)} boxes from {bx[0][1]} to {bx[-1][2]:.5g}",
flush=True)
rows = {}
for reg, Pa, Pb, n_cells in bx:
t0 = time.time()
c = c_box(Pa, Pb)
Fb = F_box(c, Pa, Pb)
fe = max(F_edge(c, P)[0] for P in (Pa, (Pa + Pb) // 2, Pb))
U = box_bound(Pa, Pb, n_cells)[0]
m = budget(Pb) - Fb - X_MAX - U / Pa
rows.setdefault(reg, []).append(
(Pa, Pb, c, Fb / fe, U, U / Pa, m) + PW[(Pa, Pb)]
+ (time.time() - t0,))
for reg, rs in rows.items():
top = max(rs, key=lambda r: r[4])
low = min(rs, key=lambda r: r[6])
print(f" {reg} [{rs[0][0]}, {rs[-1][1]:.5g}): {len(rs)} boxes, c "
f"{min(r[2] for r in rs)}..{max(r[2] for r in rs)}; U max "
f"{top[4]:.6g} at [{top[0]}, {top[1]}], U/Pa max "
f"{max(r[5] for r in rs):.4g}; least margin {low[6]:.6f} at "
f"[{low[0]}, {low[1]}]; F_box/F_edge least "
f"{min(r[3] for r in rs):.15f}; past-W add max "
f"{max(r[7] for r in rs):.4g}, over Pa "
f"{max(r[7] / r[0] for r in rs):.4g} ({sum(r[8] for r in rs)} "
f"uncut) ({sum(r[9] for r in rs):.0f} s)", flush=True)
return rows
def section_predict(rows):
"""The predictions P2..P4, each a check."""
closed = [r for reg, rs in rows.items() if reg != "open" for r in rs]
top = max(closed, key=lambda r: r[4])
first = rows["range 0"][0]
check("P2 the largest U below 1.339e10 at range 0's first box, in [38, "
"55]", top is first and 38 <= top[4] <= 55,
f"{top[4]:.4f} at [{top[0]}, {top[1]}]")
low = min(closed, key=lambda r: r[6])
check("P2 the least margin below 1.339e10 at range 0's first box, in "
"[0.11, 0.16]", low is first and 0.11 <= low[6] <= 0.16,
f"{low[6]:.6f} at [{low[0]}, {low[1]}]")
add = max(r[7] for r in closed)
check("P3 past-W add at most 5 a box below 1.339e10", add <= 5,
f"largest {add:.4g}")
op = rows.get("open", [(1,) + (0,) * 7])
mo = min(r[6] for r in op)
check("P4 the open boxes' margins at least 0.237", mo >= 0.237,
f"least {mo:.6f}")
check("P4 the open boxes' least margin in [0.2398, 0.2400]",
0.2398 <= mo <= 0.2400, f"{mo:.6f}")
add = max(r[7] / r[0] for r in op)
check("P4 the open boxes' past-W add over Pa at most 1e-6", add <= 1e-6,
f"largest {add:.4g}")
def main():
t0 = time.time()
n_cells = cells_at(399)
section_control(n_cells)
if not (CHECKS[0] and CHECKS[1] and CHECKS[3]) and not REHEARSE:
print("control failed; stopping")
sys.exit(1)
section_cut_plant(n_cells)
bx = layout()
if REHEARSE:
bx = [b for i, b in enumerate(bx) if i == 0 or b[0] != bx[i - 1][0]]
section_coverage(bx)
rows = section_boxes(bx)
allr = [r for rs in rows.values() for r in rs]
check("D10 F_box at or above F_edge at Pa, the middle and Pb of every box,"
" to 1e-13 relative", min(r[3] for r in allr) >= 1 - 1e-13,
f"least ratio {min(r[3] for r in allr):.15f}")
section_predict(rows)
worst = min(r[6] for r in allr)
check("KILL not fired: every box's margin above 0", worst > 0,
f"least margin {worst:.6f}")
wall = time.time() - t0
check("P5 wall at most 600 s", wall <= 600, f"{wall:.0f} s")
print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed, {wall:.0f} s")
sys.exit(as_recorded(MISSED))
def as_recorded(missed):
"""0 iff the checks that failed are exactly the record's misses, by name."""
got = {n for n, ok in zip(NAMES, CHECKS) if not ok}
print("misses %s; the record's %s: %s" % (
"; ".join(sorted(got)) or "none", "; ".join(sorted(missed)) or "none",
"as recorded" if got == set(missed) else "NOT AS RECORDED"))
return 0 if got == set(missed) else 1
if __name__ == "__main__":
main()