"""prunerr.py -- does the enumeration's squared-radius slack carry its float error at every cell the band scans read, heights to rank 14's, and on the chart? QUESTION. flatten.enumerate_min prunes a node when its float partial squared norm exceeds R2 = (s - 1)^2 M (1 + 1e-9) + 1e-9; it misses a vector of height at most s - 1 only if the float sum passes R2 at a node whose true sum is at most R = (s - 1)^2 M. gap.py reads it at heights to 6e44 (rank 10). The band scans also read ranks 11..14 to J_c - 1 (bandgrowth.py), heights near 1e86. Does the slack hold there, and at every other cell? BY HAND. n = r; the reduced basis's B_i = d_(i+1)/d_i and mu_ji = lam_ji/d_(j+1) are correctly rounded int/int quotients (relative u = 2^-53 however large the ints, the quotient below 1.8e308). At level i, c_i = -sum_(j>i) mu_ji x_j, y_i = x_i - c_i, p = sum y_i^2 B_i. (a) Inside the true ball every ancestor has y_j^2 B_j <= R, so X_(n-1) = sqrt(R/B_(n-1)), X_i = S_i + sqrt(R/B_i), S_i = sum_(j>i) |mu_ji| X_j bound |x_i|. (b) |c~_i - c_i| <= g_(n+1) S_i (g_k = k u/(1 - k u)), so dy_i = g_(n+1) S_i + u (sqrt(R/B_i) + g_(n+1) S_i) bounds |y~_i - y_i|. (c) y~^2 B rounds three times and the partial sums n - 1 times, all of nonnegative terms, so |p~ - p| <= E R, E = (1 + g_3) E_0 + (g_3 + g_n)(1 + E_0), E_0 = sum_i (2 sqrt(B_i/R) dy_i + (B_i/R) dy_i^2). Every X, S, dy is sqrt(R) times its value at R = 1: E depends on neither R nor the lattice's scale (B scales, mu does not). The heights enter only through overflow, R2 below 1.8e308. (d) The zigzag orders a level's points by distance from c~, not c; a point v after the one w that breaks the level has |y_v| >= |y_w| - 2 dc_i, so w's true sum passes v's by at most 4 sqrt(B_i R) dc_i + 4 B_i dc_i^2, twice level i's term of E_0 but for 2 (B_i/R) dc_i^2 of R, at most E^2/2 (z_i = sqrt(B_i/R) dc_i <= sqrt(B_i/R) dy_i, sum z_i <= E_0/2, so 2 z_i^2 <= E_0^2/2). w lies outside the true ball, |y_w| up to sqrt(R/B_i) + 2 dc_i, so its float sum errs by at most E R plus about 4 (B_i/R) dc_i dy_i <= E_0^2 and terms of order u E <= E^2/4 (E >= 4u): a missed v costs E + 2E + 2E^2. SAFE: 3E + 2E^2 + 4u < 1e-9 (4u for R2's own roundings). (b) also needs every x_j exact in a float: |x_j| <= sqrt(R2) X_j(R = 1) < 2^53. THE CELLS. Every cell band.py, gap.py and bandgrowth.py read with flatten.least_height: ranks 1..10 at depths 2 .. max(140, J_c - 1), ranks 11..14 at 2 .. J_c - 1, ranks 15..22 at 2 .. 60 (band.py stops there at its first failure or 60), J_c = 227, 241, 305, 319 (bandgrowth.py FINDINGS (2)) and 47, 57, 97, 108, 157, 170 at ranks 5..10 (gap.py). Each cell's basis is flatten.lll(flatten.basis(M, J)), least_height's. PREDICTIONS (frozen before the engine). PS1 E < 1e-11 at every cell. PS2 the largest H_ch read is rank 14's at J = 318, in [1e86, 4e86] (I_14(318) = 3.9e86 bounds it above). PS3 the largest float R2 below 1e180. KILL (printed): a cell with 2E + 4u >= 1e-9, or an R2 at or past 1e300; amended after review, before the rerun, to 3E + 4u (step (d)'s second term) and an x bound at or past 2^53, over ranks 1..22; after a later review, to 3E + 2E^2 + 4u ((d)'s second-order terms), every print but the slack line's label unchanged. CONTROL, read first. Ranks 5..8, J = 10..30: every node with true p <= R, R = h^2 M (h the least height), walked exactly in Fractions beside flatten's float order of operations; the largest |p~ - p|/R must be at most E at every cell. PLANT: u = 0 in E, caught if any read error is nonzero; replaced after review (caught by construction, it read nothing) by E/30, chosen after the first run read 0.0487 E, so it is a regression check, not a positive control: it fails only if the walk or E changes. Read, not gated: E with (b)'s c~ term dropped, and with its own rounding of y dropped. FINDINGS (copied from the printed output). (1) Control: 84 cells, 6554 nodes walked exactly, no read error above E, the largest read 0.0487 E; the plant E/30 caught at 4 cells (so a bound understated about 21 times, derived, would show). Neither term-drop plant is caught: each term of E alone covers every error the walk reads, and (b), (c), (d) rest on their proof. (2) PS2 and PS3 hold; PS1 holds at ranks 1..21 and fails at rank 22 (1.158e-11; ranks 15..22 were added after the prediction, which the print's 'PS1 (ranks 1..14)' reads). The KILL misses. Over 2991 cells E is at most 1.158e-11, at rank 22, depth 3; per rank its largest is 6.661e-16 at rank 1, 1.613e-13 at rank 10, 8.586e-13 at rank 14 (at depth 2 there, printed by the first run: E is set by the basis's shape, not the deep cells' heights), so at ranks 1..14 3E + 2E^2 + 4u is below 1e-9 by about 388 (derived from the printed E), and at rank 22 by 28.8 (printed), the least over the cells. The largest H_ch read is 3.197e86 at rank 14, J = 318, its float R2 3.394e175, far below overflow; the largest bound on a coefficient |x_j| is 76.04 (rank 22, J = 7), far below 2^53. (3) So the enumeration's minimality at every cell band.py, gap.py and bandgrowth.py read with flatten.least_height rests on this bound, not on the slack's size alone; E is at most 1.613e-13 over rank 10's cells, gap.py's. THE CHART ARM. The chart and its J = 1 column, M = 4..40, J = 1..min(30, M - 1), 732 cells, read by flatten.stage_chart. A cell where the bound fails is still exact if its least height is 1: no nonzero integer vector is lower, so nothing can be missed. NOT BLIND. A one-off read came first: the bound failed at 7 cells, ranks 33..38 at depths 2..4, h = 1 at each. The arm re-reads that with the witness checked in integers; its prints are a record, not a test of a prediction. KILL (printed): a cell where 3E + 2E^2 + 4u >= 1e-9 or the coefficient bound reaches 2^53, and the returned witness is not of height 1 with J exact divisions by x - 1. CONTROL: at the escaped cells the witness is divided by x - 1 J times in integers, the remainder read zero each time, the division itself checked on (x - 1)^3 first. FINDINGS (copied from the printed output). (4) The division control is ok; the KILL misses. (5) Of the 732 cells the bound covers 725; the largest E is 1.246e-9 at (40, 3). The 7 others, (M, J) = (36, 3), (37, 3), (38, 3), (39, 3), (40, 2), (40, 3), (40, 4) at ranks 33..38, have h = 1, a height-1 witness dividing by (x - 1)^J. So the chart's least heights are exact at every cell. RUN RECORD. Timed first: one cell 0.04 s (rank 14, J = 318). 15.2 s in all before the chart arm, 16.5 s with it, peak commit 14.2 MB under a 512 MB guard. Two first runs died in the control (a Fraction compared with an mpf, then a float zero with no numerator). A first full run read ranks 1..14 against 2E + 4u with the plant u = 0 (caught at 66 cells by construction); review found (d)'s second E to be 2E, ranks 15..22 unread and the plant empty: the KILL, the cells and the plant were amended before this run, the ranks 1..14 prints unchanged. Exit red on a control or a KILL. """ import os import sys import time from fractions import Fraction from math import isqrt os.environ.setdefault("OPENBLAS_NUM_THREADS", "1") import mpmath as mp # noqa: E402 from flatten import lll, basis, least_height # noqa: E402 from bandend import champ_height # noqa: E402 mp.mp.dps = 40 U = mp.mpf(2) ** -53 JC = {5: 47, 6: 57, 7: 97, 8: 108, 9: 157, 10: 170, 11: 227, 12: 241, 13: 305, 14: 319} def g(k, u): return k * u / (1 - k * u) def err_bound(d, lam, n, u=U, cterm=True, yterm=True): """E of (c) and the largest X_i, at R = 1, from the exact d and lam; cterm=False drops (b)'s c~ term, yterm=False its own rounding of y (the PLANTS).""" B = [mp.mpf(d[i + 1]) / d[i] for i in range(n)] mu = [[mp.mpf(lam[k][j]) / d[j + 1] for j in range(k)] for k in range(n)] X = [mp.mpf(0)] * n E0 = mp.mpf(0) gn1 = g(n + 1, u) if cterm else 0 for i in range(n - 1, -1, -1): S = sum((abs(mu[j][i]) * X[j] for j in range(i + 1, n)), mp.mpf(0)) r = 1 / mp.sqrt(B[i]) X[i] = S + r dy = gn1 * S + (u * (r + gn1 * S) if yterm else 0) E0 += 2 * mp.sqrt(B[i]) * dy + B[i] * dy ** 2 g3, gn = g(3, u), g(n, u) return (1 + g3) * E0 + (g3 + gn) * (1 + E0), max(X) def walk(b, d, lam, R): """Every node with true p <= R, walked exactly beside flatten's float order: the largest |p~ - p|/R read, and the node count.""" n = len(b) Bf = [d[i + 1] / d[i] for i in range(n)] muf = [[lam[k][j] / d[j + 1] for j in range(k)] for k in range(n)] Bq = [Fraction(d[i + 1], d[i]) for i in range(n)] muq = [[Fraction(lam[k][j], d[j + 1]) for j in range(k)] for k in range(n)] x = [0] * n worst = [Fraction(0), 0] def rec(i, pf, pq): cf = 0.0 cq = Fraction(0) for j in range(i + 1, n): if x[j]: cf -= muf[j][i] * x[j] cq -= muq[j][i] * x[j] # integer xi with (xi - cq)^2 Bq <= R - pq room = (R - pq) / Bq[i] span = isqrt(int(room) + 1) + 2 for xi in range(int(cq) - span, int(cq) + span + 1): yq = xi - cq q = pq + yq * yq * Bq[i] if q > R: continue yf = xi - cf f = pf + yf * yf * Bf[i] worst[1] += 1 worst[0] = max(worst[0], abs(Fraction(f) - q) / R) x[i] = xi if i > 0: rec(i - 1, f, q) x[i] = 0 rec(n - 1, 0.0, Fraction(0)) return worst def control(): off = offp = offq = offs = cells = nodes = 0 worst_ratio = 0.0 for r in range(5, 9): for J in range(10, 31): M = r + J b, d, lam = lll(basis(M, J)) h = least_height(M, J)[0] R = h * h * M e, nn = walk(b, d, lam, R) e = mp.mpf(e.numerator) / e.denominator E = err_bound(d, lam, r)[0] Ep = err_bound(d, lam, r, cterm=False)[0] Eq = err_bound(d, lam, r, yterm=False)[0] cells += 1 nodes += nn off += e > E offp += e > Ep offq += e > Eq offs += e > E / 30 worst_ratio = max(worst_ratio, float(e / E)) print("CONTROL ranks 5..8, J = 10..30: %d cells, %d nodes walked exactly," " read error above E at %d; largest read/E %.3g; PLANT the c~ term" " dropped: above at %d; PLANT the y~ term dropped: above at %d;" " PLANT E/30: above at %d" % (cells, nodes, off, worst_ratio, offp, offq, offs), flush=True) return off == 0 and offs > 0 def divides(v, J): """J exact divisions of v by x - 1 (synthetic division at 1).""" p = list(v) for _ in range(J): q, acc = [], 0 for c in reversed(p): acc += c q.append(acc) if q[-1] != 0: return False p = list(reversed(q[:-1])) return True def chart(): """The chart arm: the bound per cell, else a height-1 witness.""" t0 = time.time() ctl = divides([-1, 3, -3, 1], 3) and not divides([-1, 3, -3, 1], 4) print("CONTROL divides: (x - 1)^3 to depth 3 and not 4: %s" % ("ok" if ctl else "OFF")) cells = covered = 0 escaped, bad = [], [] worst = (0, None) for M in range(4, 41): for J in range(1, min(30, M - 1) + 1): r = M - J b, d, lam = lll(basis(M, J)) s0 = min(max(abs(c) for c in v) for v in b) R2 = (s0 - 1) ** 2 * M * (1 + 1e-9) + 1e-9 E, X = err_bound(d, lam, r) cells += 1 if E > worst[0]: worst = (E, (M, J)) if (3 * E + 2 * E * E + 4 * U < mp.mpf(10) ** -9 and mp.sqrt(R2) * X < 2 ** 53): covered += 1 continue h, v = least_height(M, J)[:2] ok = h == 1 and max(abs(c) for c in v) == 1 and divides(v, J) escaped.append((M, J, r, mp.nstr(E, 4), h)) if not ok: bad.append((M, J)) print("CHART %d cells: the bound covers %d; largest E %s at %s" % (cells, covered, mp.nstr(worst[0], 4), worst[1])) print("escaped by h = 1 (M, J, r, E, h): %s" % escaped) print("CHART KILL %s %s (%.1f s)" % (bad, "fires" if bad else "misses", time.time() - t0)) return ctl and not bad def main(): t0 = time.time() ctl = control() print("controls: %s (%.1f s)" % ("ok" if ctl else "OFF", time.time() - t0), flush=True) worstE, worstR2, worstH = (0, None), (0.0, None), (0, None) E14, xmax = 0, (0, None) bad, cells = [], 0 for r in range(1, 23): if r >= 15: top = 60 elif r >= 11: top = JC[r] - 1 else: top = max(140, JC.get(r, 0) - 1) rE = 0 for J in range(2, top + 1): M = r + J b, d, lam = lll(basis(M, J)) s0 = min(max(abs(c) for c in v) for v in b) R2 = (s0 - 1) ** 2 * M * (1 + 1e-9) + 1e-9 E, X = err_bound(d, lam, r) xb = mp.sqrt(R2) * X cells += 1 rE = max(rE, E) if r <= 14: E14 = max(E14, E) if xb > xmax[0]: xmax = (xb, (r, J)) if E > worstE[0]: worstE = (E, (r, J)) if R2 > worstR2[0]: worstR2 = (R2, (r, J)) if J >= r - 1: H = champ_height(r, J) if H > worstH[0]: worstH = (H, (r, J)) if (3 * E + 2 * E * E + 4 * U >= mp.mpf(10) ** -9 or R2 >= 1e300 or xb >= 2 ** 53): bad.append((r, J)) print(" r = %2d depths 2..%d largest E %s" % (r, top, mp.nstr(rE, 4)), flush=True) print("CELLS %d; largest E %s at %s; largest float R2 %.3e at %s;" " largest H_ch read %.3e at %s (%.1f s)" % (cells, mp.nstr(worstE[0], 4), worstE[1], worstR2[0], worstR2[1], float(worstH[0]), worstH[1], time.time() - t0)) wE = worstE[0] print("largest x bound %s at %s; slack over 3E + 2E^2 + 4u at the" " largest E: %.1f" % (mp.nstr(xmax[0], 4), xmax[1], float(mp.mpf(10) ** -9 / (3 * wE + 2 * wE * wE + 4 * U)))) ps1 = E14 < mp.mpf(10) ** -11 ps2 = worstH[1] == (14, 318) and 1e86 <= float(worstH[0]) <= 4e86 ps3 = worstR2[0] < 1e180 print("PS1 (ranks 1..14) %s; PS2 %s; PS3 %s" % tuple("holds" if v else "fails" for v in (ps1, ps2, ps3))) print("KILL %s %s" % (bad[:5], "fires" if bad else "misses")) ok = chart() return not ctl or bool(bad) or not ok if __name__ == "__main__": sys.exit(1 if main() else 0)