erring_cells.py
Checks: SIZE: The cell tent bound theorem; property.
Cited on: SIZE.
To run it: save erring_cells.py, with erring_label.py, erring_lower.py, erring_widths.py, erring_zones.py beside it; needs numpy.
"""
erring_cells.py -- the cell tent bound: each cell pair errs at least the
larger of its two one-boundary tent integrals, and over nondecreasing
widths that bound puts the error at least (1 - o(1))/(2 log m), and the
sorted bound that would carry it to every partition, if true, is sharp
on a grid; the per-shift bound LB3 holds it where the tents fail.
QUESTION. erring_lower.py's LB charges a cell one label boundary at a
time and reaches liminf kappa_m log m >= 1/4. Charging the cell pair
{(a, b), (b, a)} the larger of its two directions, does the continuum
error E(P) reach (1 - o(1))/(2 log m), the upper side's constant? Over
which partitions, and does sorting the widths ever lower E? Does E
itself stay above the increasing bound of its sorted widths, and how
close does it come?
PROOF (P, E and LB as in erring_lower.py; widths w_a, boundaries c_a).
(1) T_ab is the integral over y in I_b of h_a(mid_a + y), h_a(t) the
largest (r - |t - beta|)_+ over boundaries beta, r = w_a/2: the
cell I_a x I_b errs at least T_ab (erring_lower.py's step (2)
with the boundary nearest the shifted midpoint at each y). The
cells (a, b) and (b, a) hold the same sums, so each errs at least
max(T_ab, T_ba), and E >= LB2 = sum over ordered (a, b) of
max(T_ab, T_ba). In closed form: a gap of length g holds r^2 of
H_a, the integral of h_a, if g >= 2r, else rg - g^2/4.
From here the widths are nondecreasing, w_1 <= ... <= w_m,
c_j = w_1 + ... + w_j, S_j = w_1^2 + ... + w_j^2.
(2) LB2 >= G2 = sum_a T_aa + (1/2) sum_j S'_j, S'_j the sum of w_b^2
over b <= j with c_(b-1) >= c_j/2. LB2 >= sum T_aa +
2 sum_b sum_{a<b} T_ba, and sum_{a<b} T_ba integrates I_b's tent
over t in [mid_b, 2c_(b-1) + r), r = w_b/2. Each boundary c_j,
j >= b, c_j <= 2c_(b-1), carries a whole tent there (its ends
c_j - r >= c_(b-1) + r and c_j + r <= 2c_(b-1) + r, a tie
included), the tents disjoint as later gaps are at least w_b:
w_b^2/4 each, and the pairs b <= j, c_(b-1) >= c_j/2 sum to
sum_j S'_j.
(3) J = the j with w_j <= c_(j-1); S'_j >= w_j^2 there. A j outside J
(a doubling step) is wider than all earlier widths together, so
the doubling widths past any point sum below twice the last.
(4) For j in J, rho_j = (S_j - S'_j)/S_j <= 3/4: S_j - S'_j = S_i,
i the least index with c_i >= c_j/2, S_i <= w_i c_i, S'_j >= w_i s with
s = c_j - c_i >= max(w_i, c_j/2 - w_i) >= c_j/4.
(5) LEMMA. If c_j > delta for every j > j0 >= 1, then
Lambda = sum_{j in J, j0 < j < m} w_j^2/S'_j
<= log(S_(m-1)/S_j0) + 8 W, W = ceil(log2(1/delta)). By (4),
w^2/S' <= w^2/S + 4 w^2 X/S^2, X = S - S'; the first telescopes;
in a window c_j in (U/2, U], U = 2^-k, with f its first index,
X_j <= S_f, the term j = f is at most 1 and the rest sum to at
most S_f sum (1/S_(j-1) - 1/S_j) <= 1.
(6) j0 the last j with c_j <= delta, A = c_j0. A < delta/2 gives
E >= LB >= S_m/8 > delta(2 - delta)/32; w_m >= eta gives
E >= S_m/8 >= eta^2/8; else S_j0 >= delta^2/(4m), the J widths
strictly between j0 and m sum past 1 - delta - 3 eta (w_m < eta,
and by (3) the doubling widths there sum below 2 w_(m-1) < 2 eta),
and Cauchy-Schwarz with (5) gives E >= G2 >= F3's third term:
F3 = min(delta(2 - delta)/32, eta^2/8,
max(0, 1 - delta - 3 eta)^2
/ (2 log(4m/delta^2) + 16 ceil(log2(1/delta)))).
(7) delta = 8/log m, eta = 2/sqrt(log m): F3 = (1 - o(1))/(2 log m).
So over nondecreasing widths E(P) >= (1 - o(1))/(2 log m). The bound
does not pass to kappa_m by sorting: sorting the widths increasing
can raise E (the run's rounding cut sets).
THE SORTED BOUND. E(P) >= G2(sort P), sort P the widths sorted
increasing, would carry (7) to every P. The tents alone cannot give it:
LB2(P) < G2(sort P) at a seven-interval P (section_sorted). It is sharp
if true, on the half-shifted grid H_n: widths w/2, w, ..., w, w/2 (n - 1
full ones, m = n + 1, w = 1/n), I_0 = [0, w/2), I_a = [(a - 1/2)w,
(a + 1/2)w), I_n = [(n - 1/2)w, 1), boundaries (a + 1/2)w and 1 = nw.
(H1) E(H_n) = (n^2 + 3n - 1)/(8n^2), n >= 2. A full cell (a, b) has
its sums on a tent over [(a + b - 1)w, (a + b + 1)w): for
a + b <= n - 1 the boundaries (a + b -+ 1/2)w leave 3/4 in the
middle label, error w^2/4, (n - 1)(n - 2)/2 cells; for a + b = n
the tent's peak sits at 1, half of it overflows, error w^2/2,
n - 1 cells; past n every sum overflows. I_0 with a full I_b (both
orders, 2(n - 1) cells) is a trapezoid of mass w^2/2 over
[(b - 1/2)w, (b + 1)w) with 3w^2/8 below the boundary
(b + 1/2)w: error w^2/8. I_n with a full interval overflows
whole. The half-half cells (0, 0), (0, n), (n, 0) split their
tents of mass w^2/4 at w/2 or at 1: w^2/8 each; (n, n) overflows.
Summed: (w^2/8)((n - 1)(n - 2) + 4(n - 1) + 2(n - 1) + 3).
(H2) G2(sort H_n) = (n^2 + n - 1/2)/(8n^2). Sorted, c_1 = w/2 and
c_j = (j - 1)w for j >= 2: S'_1 = 0, S'_2 = w^2/4 (the tie
c_1 = c_2/2) and for j >= 3 the b counted are b >= (j + 3)/2, so
S'_j = w^2 floor((j - 1)/2), summing with S'_2 to
w^2(floor(n^2/4) + 1/4), with a tie at b = (j + 3)/2 for every
odd j: the closed form is the inclusive G2's (ties out it is
smaller). T_11 = w^2/16 (the boundary w/2 at the
middle of mid_1 + I_1), T_22 = 0 (mid_2 + I_2 = [5w/4, 7w/4)
keeps w/4 from both w and 2w), and for a >= 3 mid_a + I_a centres
on (2a - 3)w, a whole tent w^2/4 while 2a - 3 <= n, nothing
past: floor((n - 1)/2) of them. At n = 2k and 2k + 1 the sum is
8G2 = n^2 + n - 1/2 times w^2.
So E/G2(sort H_n) = 1 + (2n - 1/2)/(n^2 + n - 1/2) -> 1: a constant
above 1 in front of G2 fails on H_n.
LB3, the per-shift bound. e_a(y) = w_a less the longest piece of
I_a + y cut at the boundaries; U_ab its integral over y in I_b;
LB3 = sum over ordered (a, b) of max(U_ab, U_ba). The cell's one
label holds at most the longest piece at each y, so E_ab >= U_ab,
U_ba; e_a(y) >= h_a(mid_a + y), equal with one boundary inside. So
LB2 <= LB3 <= E. Between breakpoints (the c_b and beta - c_(a-1),
beta - c_a) the longest piece is max(p - y, q + y, r), linear
between its crossings, so the trapezoid rule cut there is exact.
At the tie (2, 4, 2, 5)/13 by hand (lengths in 1/13, areas in
1/169): U rows (1, 1, 1,
1), (2, 6, 1/2, 7/2), (1, 1/2, 1/2, 0), (2, 17/4, 0, 0), LB3 = 27/169;
sorted (2, 2, 4, 5), T^s diagonal (1, 0, 7/2, 0), S' = (0, 4, 16,
25) with both ties in, G2 = 27/169: equal.
DESIGN. tents computes T exactly (partials of H_a through
Q(z) = rz - z^2/2, z <= r, r^2/2 past), tents_brute by a midpoint grid
of the largest tent. G2's membership c_(b-1) >= c_j/2 is read at
c_j/2 (1 - 1e-12), a tie included as step (2) allows. The general
population: erring_lower.py's partitions, its exponential cuts
t -> t + eps e^(K(t - 1)) at K = 3..8, eps in 1/K, 1/K^2 (m <= 3100),
and the eps = 1/K widths shuffled (seed 1616); E where m <= 2500. The
increasing population: the members with nondecreasing widths, the
exponential cuts (a clipped last width narrower than the one before it
merged into it), constant widths (m = 2, 10, 100, 1000) and geometric
widths w_j ~ q^j (q in 1.01, 1.1, 1.5, 1.9, 2, 2.1, 3; m in 10, 40,
200). Steps (5) and (6) are read at delta in {1/2, 1/4, 0.1, 8/log m
when m > e^8} and eta in {1/4, 1/8, 2/sqrt(log m)}. section_sorted
reads fixed partitions: the tent witness (a search's least point, its
widths to 10 places), H_n, and the top: H_n scaled onto [0, 1/2)
beside one interval [1/2, 1) (m = n + 2), with G2 ties in and out.
section_lb3: cell_errs, erring_zones.py's _rows per cell; shifts, the
exact U, the interior's longest gap from a sparse table; shifts_brute
a midpoint grid; read at the tent witness, the tie (2, 4, 2, 5)/13,
H_50, H_200, the top at n = 200 and the rearrangement at four places.
PREDICTIONS AND KILLS (written before the run).
C1 [0, 1]: T_11 = 1/4.
C2 closed-form T against the brute grid at 200,000 points, three
small partitions: |diff| < 1e-6.
C3 break: T + T^t summed in place of their max prints a ratio over E
above 1.
K0 LB <= sum T <= LB2 <= E (1e-9) wherever E is computed.
P0 E/LB2 on the exponential cuts about 1.03 to 1.10 (E/LB 1.6 to 1.9).
C4 break: G2 with the S' coefficient 1 prints G2/LB2 above 1.
K1 max G2/LB2 <= 1 + 1e-9 over the increasing population; kills (2).
C5 break: rho's bound 1/2 for 3/4 prints max rho above 1/2.
K2 max rho over J <= 3/4 + 1e-9; kills (4).
C6 break: the lemma's bound without 8W prints Lambda/bound above 1.
K3 max Lambda/(log(S_(m-1)/S_j0) + 8W) <= 1 + 1e-9 over (partition,
delta) pairs with A >= delta/2; kills (5).
C7 break: the doubling widths past j0 against their last, not twice
it, print a ratio above 1.
K4 max over those pairs of the doubling widths past j0 over twice the
last of them below 1; kills (3).
C8 break: the J widths strictly between j0 and m against 1 - delta
(no allowance for w_m or the doubling widths) print a shortfall at
some triple of (6)'s main case (first written "at every triple",
which the code never tests: restated after the run).
K5 min over (6)'s main case (A >= delta/2, w_m < eta) of the J widths
strictly between j0 and m less (1 - delta - 3 eta) above 0; kills
(6)'s count.
P1 F3 log m on the grid tiny at every m run; at the limit's choice in
closed form about 0.37 at m = 10^1024.
P2 E against the same widths sorted: E/E_up below 1 somewhere
(rounding), E_down/E_up at least 1; E(P)/G2(sort P) and
LB2(P)/G2(sort P) at least 1 over the general population.
The fold of three probes run before (their predictions met there):
C9 (H1) and (H2) against cont_error and g2 at n = 2..60, 100, 200,
400, 800: |diff| < 1e-12.
C10 break: the S' coefficient 1 prints E/G2(sort H_400) below 1.
P3 the tent witness: LB2(P)/G2(sort P) 0.960165 ties in and out and
by the brute grid, E/G2(sort P) 1.151538 there.
P4 m(E/G2(sort P) - 1) on H_n near 2 (1.990 at n = 50); on the top
0.414 to 0.402 at n = 100..800; m(LB2/G2(sort P) - 1) 0.490 to
0.499 on H_n at n = 50..400, about 0.10 on the top.
K6 min E(P)/G2(sort P) over section_sorted's partitions at least
1 - 1e-9.
The fold of a fourth probe (its predictions met there):
C11 exact U against the midpoint grid, four small partitions, one
with several boundaries inside a shifted interval: < 1e-8.
C12 break: U + U^t summed in place of their max prints a ratio over
E above 1.
K7 T <= U and max(U_ab, U_ba) <= E_ab (1e-12) at every cell read.
P5 LB3/G2(sort P) 1.015148 at the tent witness; 169 LB3 = 169 G2 =
27 at the tie, ties in (1.588235 out); m(LB3/G2(sort P) - 1) 0.980,
0.995 on H_50, H_200, 0.201 on the top at n = 200.
P6 LB3(P)/LB3(sort P) below 1 at the rearrangement (0.945197 at its
unrounded widths), so the bound does not pass through sorting.
K8 min LB3(P)/G2(sort P) over these at least 1 - 1e-9.
The fold of a fifth probe, the lattice (its predictions met there over
a wider range). The least of LB3/G2(sort P) sits ON ties, where G2
jumps up and a search lands only by luck; integer widths put every
tie of small denominator on the lattice exactly: every composition
(every order) of N <= 12, and those of at most five parts to N = 20.
A reading below 1 - 1e-9 is a witness only if its strict reading one
width up (a width 1e-7 of the total put first, ties out) is below 1.
C13 the tie (2, 4, 2, 5) reads LB3 - G2(sort P) within 1e-12 of 0.
C14 break: the kill run on LB2 in place of LB3 fires.
K9 a strict witness on the lattice kills LB3(P) >= G2(sort P) there
(after an audit it also fails on a ties-in reading below 1 - 1e-9,
which the rule at integer widths states).
FINDINGS (from the printed run; C1 0.250000000000, C2 max |diff|
5.31e-13 and C3 max 1.777 passed first; C4 max 2.000, C5 max 0.7497,
C6 max 2.506, C7 max 1.764 and C8 max shortfall 0.9000 before K1).
K0 silent over 1202 partitions (m up to 2984): largest ratio 1.000000.
E/LB2 is 1.976 to 2.000 on tiny then wide and 1.996 to 1.999 on
alternating, where comparable widths carry the mass; 1.344 to 1.474
on the shuffles.
P0 held: E/LB2 1.029 to 1.105 over the nine exponential cut sets with
E (m = 21 to 2418), where E/LB is 1.633 to 1.905. On this one
family LB2 holds nearly all of E.
K1 silent over 511 increasing partitions (m up to 2984): max G2/LB2
1.000000. K2 silent: max rho 0.7497 (3/4 is approached where c_j/2
falls just past a boundary among nearly equal widths). At q = 2
the doubling steps from j = 40 on (2^j >= 10^12) pass the 1e-12
guard and count in J, so that group reads G2/LB2 1.0000 and tests
the doubling regime only at its first steps. K3 silent
over 724 (partition, delta) pairs: max 0.2988. K4 silent: max
0.8818. K5 silent over 936 main-case triples: min margin 0.1783.
F3/max(G, G2) at most 0.0372.
P1 held: the best F3 log m on the grid is at most 0.0156 at the m run;
at the limit's choice it reads 0.0000 at m = 10^8 and 10^16
(1 - delta - 3 eta <= 0 there), 0.0772, 0.2458 and 0.3679 at
m = 10^64, 10^256 and 10^1024, against its limit 1/2: the 1/2 is
asymptotic.
P2 held: E/E_up reaches 0.5706 (rounding) and 0.9889 (Dirichlet), so
sorting increasing can raise E; E_down/E_up at least 1.0002 over
1200 partitions. min E(P)/G2(sort P) 1.0066 (sorted), 1.0173 at
rounding; min LB2(P)/G2(sort P) 1.0000 (log-uniform, tiny then
wide), 1.0025 at rounding.
The fold (C9 max |diff| 4.4e-14 and C10 0.5031 passed first).
P3 held: LB2/G2(sort P) 0.960165 ties in, out and by the brute grid;
E/G2(sort P) 1.151538. The tents lose 4% at a narrow interval just
after a wide one, where E keeps 15%.
P4 held: m(E/G2(sort P) - 1) 1.990, 1.995, 1.998, 1.999, 1.999 on H_n
at n = 50..800, as (H1)/(H2) give, and 0.414, 0.407, 0.404, 0.402 on
the top at n = 100..800; m(LB2/G2(sort P) - 1) 0.490 to 0.499 on H_n,
0.101 to 0.100 on the top. On both grids the tents hold G2 too, so
their only loss read is the witness's shape.
K6 silent: min E(P)/G2(sort P) 1.000501 (the top, n = 800).
The fourth fold (C11 max |diff| 5.31e-13, C12 max 1.909 and K7 max
excess 1.1e-16 passed first).
P5 held: LB3/G2(sort P) 1.015148 at the tent witness, ties in and
out; at the tie 169 LB3 = 169 G2 = 27.000000000 ties in (1.588235
out), the hand values; m(LB3/G2(sort P) - 1) 0.980, 0.995 on H_50,
H_200 and 0.201 on the top at n = 200, twice LB2's.
P6 held: LB3(P)/LB3(sort P) 0.945197 at the four-place widths
(LB3(P)/G2(sort P) 1.540309 there): sorting can raise LB3.
K8 silent: min LB3(P)/G2(sort P) 1.000000000, at the tie.
The fifth fold (C13 2.8e-17 and C14 20 strict LB2 witnesses, least
0.962617 at (1, 1, 4, 2, 5), passed first).
K9 silent over 24209 lattice points: least LB3/G2(sort P)
1.000000000. The primitive unsorted points within 1e-9 of 1 are
(2, 4, 2, 5) and (1, 1, 4, 2, 5), the tie and its first 2 split.
Run: 46.2 s, peak commit 162.4 MB (T and U are m x m), 24/24; K9 gating
the ties-in least too (a scratch break scaling it by 1 - 1e-6 fails K9
alone), 46.2 s, 162.8 MB, 24/24.
"""
import os
os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import itertools # noqa: E402
import math # noqa: E402
import numpy as np # noqa: E402
from erring_lower import (bounds, cuts_of, exp_cuts, gain, # noqa: E402
partitions, sums)
from erring_zones import below, cont_error # noqa: E402
CHECKS = []
GUARD = 1e-12
def check(name, ok, detail=""):
CHECKS.append(bool(ok))
tail = f" ({detail})" if detail else ""
print(f" {'PASS' if ok else 'FAIL'} {name}{tail}")
def Q(z, r):
z = np.clip(z, 0.0, None)
return np.where(z <= r, r * z - z * z / 2, r * r / 2)
def tents(cuts):
"""T[a, b]: the one-boundary tent bound of cell (a, b), exact."""
C = np.asarray(cuts, dtype=float)
m = len(C) - 1
B = C[1:]
g = np.diff(B)
T = np.empty((m, m))
for a in range(m):
r = (C[a + 1] - C[a]) / 2
J = 2 * Q(g / 2, r)
HB = r * r / 2 + np.concatenate([[0.0], np.cumsum(J)])
t = (C[a] + r) + C
i = np.searchsorted(B, t, side="right") - 1
H = np.empty_like(t)
lead = i < 0
H[lead] = r * r / 2 - Q(B[0] - t[lead], r)
tail = i == len(B) - 1
H[tail] = HB[-1] + Q(t[tail] - B[-1], r)
mid = ~lead & ~tail
k = i[mid]
s, gk = t[mid] - B[k], g[k]
H[mid] = (HB[k] + Q(np.minimum(s, gk / 2), r)
+ np.where(s > gk / 2, Q(gk / 2, r) - Q(gk - s, r), 0.0))
T[a] = np.diff(H)
return T
def tents_brute(cuts, n=200_000):
C = np.asarray(cuts, dtype=float)
m = len(C) - 1
B = C[1:]
T = np.empty((m, m))
for a in range(m):
r = (C[a + 1] - C[a]) / 2
for b in range(m):
y = C[b] + (np.arange(n) + 0.5) * (C[b + 1] - C[b]) / n
t = C[a] + r + y
h = np.clip(r - np.abs(t[:, None] - B[None, :]), 0, None).max(1)
T[a, b] = h.sum() * (C[b + 1] - C[b]) / n
return T
def lb2(T):
return np.maximum(T, T.T).sum()
def cell_errs(cuts):
"""E[a, b]: the continuum error of cell (a, b), overflow [1, 2)."""
C = np.asarray(cuts, dtype=float)
lab = np.append(C, 2.0)
A, W = C[:-1], np.diff(C)
m = len(W)
base = (A[:, None] + A[None, :]).ravel()
a = np.broadcast_to(W[:, None], (m, m)).ravel()
b = np.broadcast_to(W[None, :], (m, m)).ravel()
lo = np.searchsorted(lab, base, side="right") - 1
hi = np.searchsorted(lab, base + a + b, side="left") - 1
best = np.zeros_like(base)
K = 8
for k in range(K + 1):
L = np.minimum(lo + k, len(lab) - 2)
land = (below(a, b, lab[L + 1] - base) - below(a, b, lab[L] - base))
best = np.maximum(best, np.where(lo + k <= hi, land, 0.0))
for i in np.flatnonzero(hi - lo > K):
L = np.arange(lo[i], hi[i] + 1)
land = (below(a[i], b[i], lab[L + 1] - base[i])
- below(a[i], b[i], lab[L] - base[i]))
best[i] = land.max()
return (a * b - best).reshape(m, m)
def sparse(g):
st = [np.asarray(g, dtype=float)]
while 2 ** len(st) <= len(g):
h = 2 ** (len(st) - 1)
st.append(np.maximum(st[-1][:-h], st[-1][h:]))
return st
def range_max(st, lo, hi):
"""max g[lo:hi] where hi > lo, else -inf."""
out = np.full(len(lo), -np.inf)
ok = hi > lo
if ok.any():
n = hi[ok] - lo[ok]
k = np.floor(np.log2(n)).astype(int)
for kk in np.unique(k):
s = k == kk
idx = np.flatnonzero(ok)[s]
out[idx] = np.maximum(st[kk][lo[ok][s]],
st[kk][hi[ok][s] - 2 ** kk])
return out
def shifts(cuts):
"""U[a, b]: the integral over y in I_b of w_a less the longest piece
of I_a + y cut at the boundaries c_1..c_m strictly inside, exact."""
C = np.asarray(cuts, dtype=float)
m = len(C) - 1
B = C[1:]
st = sparse(np.diff(B)) if m > 1 else None
U = np.zeros((m, m))
for a in range(m):
L0, R0 = C[a], C[a + 1]
w = R0 - L0
y = np.unique(np.clip(np.concatenate([C, B - L0, B - R0]), 0, 1))
y0, y1 = y[:-1], y[1:]
keep = y1 - y0 > 0
y0, y1 = y0[keep], y1[keep]
ym = (y0 + y1) / 2
ilo = np.searchsorted(B, L0 + ym, side="right")
ihi = np.searchsorted(B, R0 + ym, side="left") - 1
has = ilo <= ihi
e = np.zeros(len(ym))
if has.any():
lo, hi = ilo[has], ihi[has]
p = B[lo] - L0
q = R0 - B[hi]
r = (range_max(st, lo, hi) if st is not None
else np.full(len(lo), -np.inf))
s0, s1 = y0[has], y1[has]
pts = np.stack([s0, (p - q) / 2, p - r, r - q, s1], axis=1)
pts = np.sort(np.clip(np.nan_to_num(pts, neginf=-1, posinf=2),
s0[:, None], s1[:, None]), axis=1)
f = np.maximum(np.maximum(p[:, None] - pts, q[:, None] + pts),
r[:, None])
longest = ((f[:, 1:] + f[:, :-1]) / 2 * np.diff(pts, axis=1)
).sum(1)
e[has] = w * (s1 - s0) - longest
bb = np.searchsorted(C, ym, side="right") - 1
U[a] = np.bincount(bb, weights=e, minlength=m)[:m]
return U
def shifts_brute(cuts, n=200_000):
C = np.asarray(cuts, dtype=float)
m = len(C) - 1
B = C[1:]
U = np.empty((m, m))
for a in range(m):
w = C[a + 1] - C[a]
for b in range(m):
y = C[b] + (np.arange(n) + 0.5) * (C[b + 1] - C[b]) / n
L, R = C[a] + y, C[a + 1] + y
ins = (B[None, :] > L[:, None]) & (B[None, :] < R[:, None])
pts = np.where(ins, B[None, :], np.nan)
pts = np.concatenate([L[:, None], pts, R[:, None]], axis=1)
pts = np.sort(pts, axis=1)
longest = np.nanmax(np.diff(pts, axis=1), axis=1)
U[a, b] = (w - longest).sum() * (C[b + 1] - C[b]) / n
return U
def lb3(U):
return np.maximum(U, U.T).sum()
def rising(cu):
return bool(np.all(np.diff(np.diff(np.asarray(cu, dtype=float)))
>= -1e-15))
def window(cu, incl=True):
"""S (S_0 = 0), S'_j and the J mask for j = 1..m; ties c_(b-1) =
c_j/2 in, or out."""
C = np.asarray(cu, dtype=float)
w = np.diff(C)
S = np.concatenate([[0.0], np.cumsum(w * w)])
j = np.arange(1, len(C))
i = np.searchsorted(C, C[1:] / 2 * (1 - GUARD if incl else 1 + GUARD),
side="left")
Sp = np.where(i <= j - 1, S[j] - S[np.minimum(i, j)], 0.0)
return S, Sp, i <= j - 1
def g2(cu, T, coef=0.5, incl=True):
return np.trace(T) + coef * window(cu, incl)[1].sum()
def split(cu, d):
"""j0 (the last j with c_j <= delta) and its J and doubling indices
strictly between j0 and m, or None when A < delta/2."""
C = np.asarray(cu, dtype=float)
m = len(C) - 1
j0 = int(np.searchsorted(C, d, side="right")) - 1
if j0 < 1 or C[j0] < d / 2:
return None
inJ = window(cu)[2]
idx = np.arange(j0 + 1, m)
return j0, idx[inJ[idx - 1]], idx[~inJ[idx - 1]]
def lam(cu, d, window_term=True):
"""Lambda and the lemma's bound at delta, or None if A < delta/2."""
got = split(cu, d)
if got is None:
return None
j0, idx, _ = got
S, Sp, _ = window(cu)
w = np.diff(np.asarray(cu, dtype=float))
L = float((w[idx - 1] ** 2 / Sp[idx - 1]).sum())
bound = math.log(S[len(w) - 1] / S[j0])
if window_term:
bound += 8 * math.ceil(math.log2(1 / d))
return L, bound
def f3(m, d, e):
W = math.ceil(math.log2(1 / d))
return min(d * (2 - d) / 32, e * e / 8,
max(0.0, 1 - d - 3 * e) ** 2
/ (2 * math.log(4 * m / d ** 2) + 16 * W))
def grid(m):
L = math.log(m)
ds = [0.5, 0.25, 0.1] + ([8 / L] if m > math.e ** 8 else [])
es = [0.25, 0.125, 2 / math.sqrt(L)] if m > 1 else [0.25, 0.125]
return ds, es
def general():
"""(group, cuts): erring_lower's partitions, the exponential cuts and
the eps = 1/K widths shuffled."""
rng = np.random.default_rng(1567)
for name, cu in partitions(rng):
key = name.split(",")[0].split(" m")[0]
yield key, np.asarray(cu, dtype=float)
shuf = np.random.default_rng(1616)
for K in (3, 4, 5, 6, 7, 8):
for eps in (1 / K, 1 / K ** 2):
cu = np.asarray(exp_cuts(K, eps))
if len(cu) - 1 > 3100:
continue
yield f"exp K={K} eps=1/{round(1 / eps)}", cu
if eps == 1 / K:
w = np.diff(cu)
shuf.shuffle(w)
yield f"shuf K={K} eps=1/{K}", cuts_of(w)
def rising_population():
rng = np.random.default_rng(1567)
for name, cu in partitions(rng):
if rising(cu):
yield name.split(",")[0].split(" m")[0], np.asarray(cu)
for K in (3, 4, 5, 6, 7, 8):
for eps in (1 / K, 1 / K ** 2):
cu = np.asarray(exp_cuts(K, eps))
if not rising(cu):
cu = np.delete(cu, -2)
if len(cu) - 1 <= 3100:
yield "exponential", cu
for m in (2, 10, 100, 1000):
yield "constant", cuts_of(np.ones(m))
for q in (1.01, 1.1, 1.5, 1.9, 2.0, 2.1, 3.0):
for m in (10, 40, 200):
yield f"geometric q={q}", cuts_of(q ** np.arange(1.0, m + 1))
def section_controls():
print("controls")
t11 = tents([0.0, 1.0])[0, 0]
check("C1: T_11 at [0, 1] is 1/4", abs(t11 - 0.25) < 1e-12,
f"{t11:.12f}")
worst = 0.0
for cu in ([0, 0.3, 1], [0, 0.1, 0.25, 0.6, 1], [0, 0.5, 0.55, 1]):
worst = max(worst, np.abs(tents(cu) - tents_brute(cu)).max())
check("C2: closed-form T against the brute grid", worst < 1e-6,
f"max |diff| {worst:.2e}")
def section_chain(pop):
"""K0 and the exponential cuts' E/LB2; returns E by index."""
print("K0: LB <= sum T <= LB2 <= E")
rows, k0, c3, Es, exp = {}, 0.0, 0.0, {}, []
for n, (key, cu) in enumerate(pop):
m = len(cu) - 1
T = tents(cu)
lb = bounds(cu)[0]
lbt, two = T.sum(), lb2(T)
del T
k0 = max(k0, lb / lbt, lbt / two)
group = key.split(" K=")[0]
r = rows.setdefault(group, [0, 0, math.inf, 0.0])
r[0] += 1
r[1] = max(r[1], m)
if m <= 2500:
E = Es[n] = cont_error(cu)
k0 = max(k0, two / E)
c3 = max(c3, 2 * lbt / E)
r[2] = min(r[2], E / two)
r[3] = max(r[3], E / two)
if group == "exp":
exp.append((key, m, E / two, E / lb))
check("C3: T + T^t in place of max prints a ratio over E above 1",
c3 > 1, f"max {c3:.3f}")
for key, (k, mm, lo, hi) in rows.items():
print(f" {key:16s} n = {k:3d}, m <= {mm:4d}: E/LB2 "
f"{lo:.3f}..{hi:.3f}")
check("K0: LB <= sum T <= LB2 <= E", k0 <= 1 + 1e-9,
f"over {len(pop)} partitions, largest ratio {k0:.6f}")
print("P0: the exponential cuts with E")
for key, m, a, b in exp:
print(f" {key:20s} m = {m:4d}: E/LB2 {a:.3f}, E/LB {b:.3f}")
return Es
def section_rise():
print("K1 to K5 over the increasing partitions")
pop = list(rising_population())
assert all(rising(cu) for _, cu in pop)
c4 = c5 = c6 = c7 = k1 = k2 = k3 = k4 = f3g = 0.0
c8, k5 = math.inf, math.inf
pairs = triples = main = 0
rows = {}
for key, cu in pop:
m = len(cu) - 1
w = np.diff(cu)
T = tents(cu)
two = lb2(T)
G2 = g2(cu, T)
c4 = max(c4, g2(cu, T, coef=1.0) / two)
del T
k1 = max(k1, G2 / two)
S, Sp, inJ = window(cu)
if inJ.any():
rho = float((1 - Sp[inJ] / S[1:][inJ]).max())
k2, c5 = max(k2, rho), max(c5, rho)
Sm, pair, _ = sums(cu)
top = max(gain(Sm, pair), G2)
ds, es = grid(m)
best = 0.0
for d in ds:
got = split(cu, d)
if got is not None:
pairs += 1
L, b = lam(cu, d)
k3 = max(k3, L / b if b > 0 else (math.inf if L else 0))
L0, b0 = lam(cu, d, window_term=False)
if b0 > 0:
c6 = max(c6, L0 / b0)
_, idx, dbl = got
if len(dbl):
D, last = float(w[dbl - 1].sum()), float(w[dbl[-1] - 1])
k4 = max(k4, D / (2 * last))
c7 = max(c7, D / last)
for e in es:
if w[-1] < e:
main += 1
sJ = float(w[idx - 1].sum())
k5 = min(k5, sJ - (1 - d - 3 * e))
c8 = min(c8, sJ - (1 - d))
for e in es:
triples += 1
v = f3(m, d, e)
best = max(best, v)
f3g = max(f3g, v / top)
r = rows.setdefault(key, [0, 0, 0.0, 0.0])
r[0] += 1
r[1] = max(r[1], m)
r[2] = max(r[2], G2 / two)
r[3] = max(r[3], best * math.log(m) if m > 1 else 0.0)
n0 = len(CHECKS)
check("C4: G2 with coefficient 1 prints G2/LB2 above 1", c4 > 1,
f"max {c4:.3f}")
check("C5: max rho over J above 1/2", c5 > 0.5, f"max {c5:.4f}")
check("C6: the lemma without 8W prints Lambda/bound above 1", c6 > 1,
f"max {c6:.3f}")
check("C7: the doubling widths against their last print above 1",
c7 > 1, f"max {c7:.3f}")
check("C8: the J widths against 1 - delta fall short", c8 < 0,
f"max shortfall {-c8:.4f}")
if not all(CHECKS[n0:]):
print(" CONTROL FAILED: K1 to K5 are not read.")
return
print(f" {len(pop)} partitions (m up to "
f"{max(len(cu) for _, cu in pop) - 1}), {pairs} (partition, "
f"delta) pairs with A >= delta/2, {main} triples in (6)'s main "
f"case, {triples} triples")
for key, (n, mm, a, b) in rows.items():
print(f" {key:20s} n = {n:3d}, m <= {mm:4d}: max G2/LB2 {a:.4f};"
f" best F3 log m {b:.4f}")
check("K1: G2 <= LB2", k1 <= 1 + 1e-9, f"max G2/LB2 {k1:.6f}")
check("K2: rho <= 3/4 over J", k2 <= 0.75 + 1e-9, f"max {k2:.4f}")
check("K3: Lambda <= log(S_(m-1)/S_j0) + 8W", k3 <= 1 + 1e-9,
f"max ratio {k3:.4f}")
check("K4: the doubling widths past j0 sum below twice the last",
k4 < 1, f"max D/(2 w_last) {k4:.4f}")
check("K5: the J widths between j0 and m sum past 1 - delta - 3 eta",
k5 > 0, f"min margin {k5:.4f}")
print(f" F3 <= max(G, G2): max ratio {f3g:.4f}")
print("P1: the limit's choice delta = 8/log m, eta = 2/sqrt(log m), "
"in closed form:")
for ex in (8, 16, 64, 256, 1024):
L = ex * math.log(10)
d, e = 8 / L, 2 / math.sqrt(L)
W = math.ceil(math.log2(1 / d))
v = min(d * (2 - d) / 32, e * e / 8, max(0.0, 1 - d - 3 * e) ** 2
/ (2 * (math.log(4) + L - 2 * math.log(d)) + 16 * W))
print(f" m = 10^{ex:<5d} F3 log m {v * L:.4f}")
def section_rearrange(pop, Es):
print("P2: E against the same widths sorted increasing")
rows = {}
for n, (key, cu) in enumerate(pop):
if n not in Es:
continue
w = np.diff(cu)
E = Es[n]
cs = cuts_of(np.sort(w))
up = cont_error(cs)
dn = cont_error(cuts_of(np.sort(w)[::-1]))
T = tents(cs)
G2s = g2(cs, T)
del T
two = lb2(tents(cu))
group = key.split(" K=")[0]
r = rows.setdefault(group, [0, math.inf, 0.0, math.inf, math.inf,
math.inf])
r[0] += 1
r[1] = min(r[1], E / up)
r[2] = max(r[2], E / up)
r[3] = min(r[3], dn / up)
r[4] = min(r[4], E / G2s)
r[5] = min(r[5], two / G2s)
for key, (n, lo, hi, dn, eg, lg) in rows.items():
print(f" {key:16s} n = {n:3d}: E/E_up {lo:.4f}..{hi:.4f}; min "
f"E_down/E_up {dn:.4f}, E/G2_up {eg:.4f}, LB2/G2_up {lg:.4f}")
print(f" over {sum(r[0] for r in rows.values())} partitions: min "
f"E/E_up {min(r[1] for r in rows.values()):.6f}, min "
f"E_down/E_up {min(r[3] for r in rows.values()):.6f}, min E/G2_up "
f"{min(r[4] for r in rows.values()):.4f}, min LB2/G2_up "
f"{min(r[5] for r in rows.values()):.4f}")
WITNESS = [0.0720326820, 0.0042024314, 0.0752904574, 0.3027208640,
0.0753392577, 0.0760274920, 0.3943868156]
TIE = [2 / 13, 4 / 13, 2 / 13, 5 / 13]
REARRANGED = [0.0821, 0.3561, 0.2057, 0.3561]
def h_grid(n):
return np.concatenate([[0.5], np.ones(n - 1), [0.5]])
def top(n):
return np.concatenate([h_grid(n) / (2 * n), [0.5]])
def against_sorted(w, brute=False):
"""E(P), LB2(P), G2(sort P) ties in and ties out."""
cu, cs = cuts_of(w), cuts_of(np.sort(w))
tf = tents_brute if brute else tents
Ts = tf(cs)
return (cont_error(cu), lb2(tf(cu)), g2(cs, Ts), g2(cs, Ts, incl=False))
def section_sorted():
print("E(P) and LB2(P) against G2(sort P), fixed partitions")
worst = 0.0
for n in list(range(2, 61)) + [100, 200, 400, 800]:
E, _, G, _ = against_sorted(h_grid(n))
worst = max(worst, abs(E - (n * n + 3 * n - 1) / (8 * n * n)),
abs(G - (n * n + n - 0.5) / (8 * n * n)))
check("C9: (H1) and (H2) at n = 2..60, 100, 200, 400, 800",
worst < 1e-12, f"max |diff| {worst:.1e}")
cs = cuts_of(np.sort(h_grid(400)))
one = cont_error(cuts_of(h_grid(400))) / g2(cs, tents(cs), coef=1.0)
check("C10: the S' coefficient 1 prints E/G2(sort H_400) below 1",
one < 1, f"{one:.4f}")
if not all(CHECKS[-2:]):
print(" CONTROL FAILED: the readings below are not read.")
return
E, two, G, Gx = against_sorted(WITNESS)
_, twob, _, Gxb = against_sorted(WITNESS, brute=True)
print(f" P3 the tent witness, m = 7: LB2/G2(sort P) {two / G:.6f} ties "
f"in, {two / Gx:.6f} out, {twob / Gxb:.6f} by the brute grid; "
f"E/G2(sort P) {E / G:.6f}")
low = E / G
print(" P4 m(ratio - 1), E and LB2 against G2(sort P), ties in")
for name, f, ns in (("H_n", h_grid, (50, 100, 200, 400, 800)),
("top", top, (100, 200, 400, 800))):
row = []
for n in ns:
w = f(n)
E, two, G, _ = against_sorted(w)
low = min(low, E / G)
row.append(f"n = {n}: {len(w) * (E / G - 1):.3f}, "
f"{len(w) * (two / G - 1):.3f}")
print(f" {name:3s} " + "; ".join(row))
check("K6: E(P) >= G2(sort P) at these partitions", low >= 1 - 1e-9,
f"min {low:.6f}")
section_lb3()
def section_lb3():
print("LB3: each cell charged w_a less the longest piece of I_a + y")
c11 = max(np.abs(shifts(cu) - shifts_brute(cu)).max() for cu in
([0, 0.3, 1], [0, 0.1, 0.25, 0.6, 1], [0, 0.5, 0.55, 1],
[0, 0.5, 0.52, 0.55, 0.6, 0.62, 1]))
check("C11: exact U against the midpoint grid", c11 < 1e-8,
f"max |diff| {c11:.2e}")
k7 = c12 = 0.0
for w in (WITNESS, h_grid(30), top(20), TIE, REARRANGED):
cu = cuts_of(w)
E, T, U = cell_errs(cu), tents(cu), shifts(cu)
k7 = max(k7, (T - U).max(), (np.maximum(U, U.T) - E).max())
c12 = max(c12, (U + U.T).sum() / E.sum())
check("C12 break: U + U^t in place of the max prints above E",
c12 > 1, f"max {c12:.3f}")
check("K7: T <= U and max(U_ab, U_ba) <= E_ab at every cell",
k7 <= 1e-12, f"max excess {k7:.1e}")
if not all(CHECKS[-3:]):
print(" CONTROL FAILED: the readings below are not read.")
return
ratio = []
for name, w in (("the tent witness", WITNESS), ("the tie", TIE)):
cu, cs = cuts_of(w), cuts_of(np.sort(w))
U, Ts = shifts(cu), tents(cs)
three, G, Gx = lb3(U), g2(cs, Ts), g2(cs, Ts, incl=False)
ratio.append(three / G)
print(f" P5 {name}: LB3/G2(sort P) {three / G:.6f} ties in, "
f"{three / Gx:.6f} out; 169 LB3 {169 * three:.9f}, "
f"169 G2 {169 * G:.9f}")
row = []
for name, f, n in (("H_n", h_grid, 50), ("H_n", h_grid, 200),
("top", top, 200)):
w = f(n)
cs = cuts_of(np.sort(w))
r = lb3(shifts(cuts_of(w))) / g2(cs, tents(cs))
ratio.append(r)
row.append(f"{name} n = {n}: {len(w) * (r - 1):.3f}")
print(" P5 m(LB3/G2(sort P) - 1): " + "; ".join(row))
w = np.asarray(REARRANGED)
r = lb3(shifts(cuts_of(w))) / lb3(shifts(cuts_of(np.sort(w))))
cs = cuts_of(np.sort(w))
print(f" P6 the rearrangement: LB3(P)/LB3(sort P) {r:.6f}, "
f"LB3(P)/G2(sort P) "
f"{lb3(shifts(cuts_of(w))) / g2(cs, tents(cs)):.6f}")
check("K8: LB3(P) >= G2(sort P) at these partitions",
min(ratio) >= 1 - 1e-9, f"min {min(ratio):.9f}")
def lattice():
"""Integer widths: every composition of N <= 12, and those of at most
five parts to N = 20."""
for N in range(1, 21):
for m in range(1, min(N, N if N <= 12 else 5) + 1):
for cut in itertools.combinations(range(1, N), m - 1):
yield np.diff((0,) + cut + (N,))
def sorted_g2(w, incl=True):
cs = cuts_of(np.sort(w))
return g2(cs, tents(cs), incl=incl)
def one_up(w, f):
"""f/G2(sort P) ties out, a width 1e-7 of the total put first."""
w = np.asarray(w, dtype=float)
v = np.concatenate([[1e-7 * w.sum()], w])
return f(cuts_of(v)) / sorted_g2(v, incl=False)
def section_lattice():
print("LB3/G2(sort P) over the integer lattice, ties in")
f2, f3 = (lambda c: lb2(tents(c))), (lambda c: lb3(shifts(c)))
d = f3(cuts_of(TIE)) - sorted_g2(np.asarray(TIE))
check("C13: the tie (2, 4, 2, 5) reads LB3 - G2(sort P) within 1e-12",
abs(d) < 1e-12, f"{d:.1e}")
n, low, eq, dead2, dead3 = 0, math.inf, [], [], []
for w in lattice():
G = sorted_g2(w)
if G <= 0:
continue
n += 1
cu = cuts_of(w)
r2, r3 = f2(cu) / G, f3(cu) / G
low = min(low, r3)
if r2 < 1 - 1e-9 and one_up(w, f2) < 1:
dead2.append((r2, tuple(w.tolist())))
if r3 < 1 - 1e-9 and one_up(w, f3) < 1:
dead3.append((r3, tuple(w.tolist())))
if abs(r3 - 1) <= 1e-9 and np.any(np.diff(w) < 0) \
and math.gcd(*w.tolist()) == 1:
eq.append(tuple(w.tolist()))
r, w = min(dead2, default=(math.inf, None))
check("C14 break: the kill run on LB2 in place of LB3 fires",
bool(dead2), f"{len(dead2)} strict witnesses, least {r:.6f} at {w}")
if not all(CHECKS[-2:]):
print(" CONTROL FAILED: the readings below are not read.")
return
print(f" {n} points; least LB3/G2(sort P) {low:.9f}; primitive "
f"unsorted points within 1e-9 of 1: {eq}")
check("K9: LB3(P) >= G2(sort P) over the lattice",
not dead3 and low >= 1 - 1e-9,
f"{len(dead3)} strict witnesses, least ties in {low:.9f}")
def main():
section_controls()
if not all(CHECKS):
print("\nCONTROL FAILED: nothing below is read.")
raise SystemExit(1)
pop = list(general())
Es = section_chain(pop)
section_rise()
section_rearrange(pop, Es)
section_sorted()
section_lattice()
print(f"{sum(CHECKS)}/{len(CHECKS)} checks passed")
raise SystemExit(0 if all(CHECKS) else 1)
if __name__ == "__main__":
main()