prunerr.py

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