takeph.py

"""takeph.py -- at the takeover theorem's depth, is the least pure product
the champion, and does F_e beat it?

QUESTION. At J* = (4e - 1) l - 1, l = r - 2, e = floor(sqrt(l/2)) at odd r
and floor(sqrt(l)/2) at even r, the takeover theorem gives height(F_e) <
H_ch. The cell (J* + r, J*) fails when height(F_e) < ph. Is ph = H_ch
there, which shape attains ph, and is F_e below it? r = 5..25.

BY HAND (re-derived before the code). An admissible pure product is
(x - 1)^(J+k) prod [g_i], [g] = 1 + x + ... + x^(g-1), parts g_i >= 2,
k + sum(g_i - 1) <= r - 1, k >= 0, at most J + k parts >= 2. Its order at
-1 is s = #even parts <= r - 1. s = r - 3 has the five cofactors 1, 1 + x,
(1 + x)^2, 1 + x^2, 1 - x + x^2 (after x -> -x); the last is F_1, so at
e = 1 the theorem's witness is itself pure and the cell need not fail.

PREDICTIONS (frozen before the code). P1 at the e = 1 ranks (r = 5, 7, 9
and even 6..16) ph < H_ch (proved: F_1 is pure), ph attained by the part-3
shape at most of them. P2 at the e >= 2 ranks (odd 11..25, even 18..24)
height(F_e) < ph, ph/height(F_e) >= 1.05 at each. P3 there ph = H_ch at
most ranks, else the part-3 shape, ph/H_ch in [0.9, 1]. P4 the least
order-(r - 2) shape over H_ch within 30% of sqrt(J*/(4(r - 1)))
(TRANSPLANT: Gaussian, M_(s-1)/M_s ~ s^(-1/2)).
KILL. K1 an e >= 2 rank where height(F_e) >= ph is printed.
CONTROLS, read first. C1 the restricted ph against flatten.ph_of(
pure_table(23)) at every cell M <= 24; its PLANT: k pinned to 0 (frozen),
and the parts capped at 3 (added after the run, the k plant read
nothing). C2 the champion in the scan against circlemean's
exact_Hch. C3 the lattice route at (J* + r, J*): the witness clears J, h
<= height(F_e), the enumeration complete by prunerr's bound (3E + 2E^2 +
2^-51 < 1e-9; its coordinate and overflow conditions are not printed).
A failure h < ph stands on its cleared witness and the scan's ph (C1);
completeness carries only h = ph.

FINDINGS (copied from the printed output).
  (1) C1 0 of 276 cells differ. The parts-at-3 plant is caught at 80 of
      276; the k plant at none: no cell of width <= 24 needs k > 0, so C1
      reads nothing for k (k > 0 attains no ph and no least non-champion
      below either). C2 holds. C3 holds at every rank, E <= 9.48e-14.
  (2) K1 misses: ph/height(F_e) >= 1.8895 (r = 21) over the e >= 2 arm
      (P2 holds). The lattice route reads h < ph at the cell (J* + r, J*)
      at r = 11, 13..25: J_hi(r) >= J* there, computed. At r = 5..10 and
      12 (e = 1) the cell does not fail, h = ph; at r = 14, 16 (e = 1) it
      fails, h/ph = 0.7407 and 0.6655.
  (3) P1: ph/H_ch = 0.1103 .. 0.4000 at the e = 1 ranks; ph is the part-3
      shape F_1 at r = 5, 7, 9 only; one part 3 with an excess unspent
      (order r - 4) at r = 6, and two or three parts 3 at r = 8, 10, 12,
      14, 16: P1's shape half fails.
  (4) P3 fails: ph = H_ch at r = 11, 13 (e = 2, the part 3 tying it, e -
      1 = 1) and 21, 23, 25 (e = 3);
      at r = 15 and 17..24 (e = 2) ph is the part-5 shape k = 0 2^(r-5) 5,
      of order r - 5, ph/H_ch = 0.9724, 0.8890, 0.8621, 0.8393, 0.8403,
      0.8236, 0.8109 at r = 15, 17, 18, 19, 20, 22, 24.
  (5) P4 holds: the order-(r - 2) shape over H_ch reads 1.2169 .. 1.3423
      at e = 2 and 1.6594 .. 1.6701 at e = 3, within 5.1% of the Gaussian
      1.2450 .. 1.6202; at e = 1 0.7857 .. 0.9667 against 0.7071 .. 0.8266.
  (6) At odd r the five order-(r - 3) shapes stand at exactly e - 1
      (the part 3), e (cofactor 1), 2e - 1 (the part 4), 2e - 1/2
      (k = 1) and 4e - 1 (k = 2) times H_ch, exact Fractions, at every
      odd rank read. By hand, Q = (1 + x)^J (1 - x)^(r-3): J* makes the
      Beta mean exactly 1/(4e), so the central coefficient of x Q is
      -e times the champion's, F_t's is (t - e) times it, and F_e's
      vanishes; the k = 1 shape's middle pair is 2e(1 - mu) = 2e - 1/2.
      That each height sits at its centre is read off the run.
READS (not verdicts). The middle coefficient bounds no class: it is 0 at some
shape at every rank read. The test family (p <= 16, q = 0, 1,
k within one of centring) puts every non-champion pure product above F_e's
circle mean at e = 3 (least test/H_ch 1.3174, 1.3189, 1.3201 at r = 21,
23, 25 against means 0.7370, 0.7016, 0.6709, the least at the order-(r - 2)
shape) and fails at every e = 2 rank read (0.1350 .. 0.4928 against 0.4581
.. 0.6905, the least at part-5 shapes).

RUN RECORD. 27.3 s, peak commit 23.5 MB (the reads 21 s of it).
"""
import math
import os
import sys
import time
from fractions import Fraction

os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import mpmath as mp  # noqa: E402
import flatten  # noqa: E402
import circlemean  # noqa: E402
import prunerr  # noqa: E402


def mul_block(p, g):
    """p * (1 + x + ... + x^(g-1)), by a running window."""
    n = len(p) + g - 1
    out = [0] * n
    acc = 0
    for i in range(n):
        if i < len(p):
            acc += p[i]
        if i - g >= 0:
            acc -= p[i - g]
        out[i] = acc
    return out


def pure_scan(J, r, kmax=None, xmax=None):
    """Every admissible pure product at (J + r, J): yields (height, k, parts)
    with parts the multiset of g_i >= 2, nonincreasing (kmax, xmax: the
    PLANTS, k and the part excess g - 1 capped)."""
    top = r - 1
    for k in range(0, top + 1 if kmax is None else min(top, kmax) + 1):
        base = flatten.pow_xm1(J + k)
        room = top - k
        stack = [(base, (), room, room if xmax is None else xmax)]
        while stack:
            p, parts, left, cap = stack.pop()
            yield flatten.height(p), k, parts
            if len(parts) >= J + k:
                continue
            for x in range(min(left, cap), 0, -1):
                stack.append((mul_block(p, x + 1), parts + (x + 1,),
                              left - x, x))


def order(parts):
    return sum(1 for g in parts if g % 2 == 0)


def champion(r):
    return (0, (2,) * (r - 1))


def read_cell(J, r, kmax=None, xmax=None):
    """ph with its shape, the champion's height, and per order s the least
    height with its shape."""
    best = None
    by_s = {}
    hch = None
    for h, k, parts in pure_scan(J, r, kmax, xmax):
        if best is None or h < best[0]:
            best = (h, k, parts)
        s = order(parts)
        if s not in by_s or h < by_s[s][0]:
            by_s[s] = (h, k, parts)
        if (k, parts) == champion(r):
            hch = h
    return best, hch, by_s


def theorem_e(r):
    l = r - 2
    return math.isqrt(l // 2) if r % 2 else math.isqrt(l) // 2


def theorem_J(r):
    return (4 * theorem_e(r) - 1) * (r - 2) - 1


def name(k, parts):
    """The shape in words: k extra (x - 1)'s and the parts above 2."""
    big = [g for g in parts if g != 2]
    twos = len(parts) - len(big)
    return "k=%d 2^%d%s" % (k, twos, "".join(" %d" % g for g in big))


def c1(kmax=None, xmax=None):
    """The restricted ph against flatten's full table, every M <= 24."""
    table = flatten.pure_table(23)
    diff = 0
    cells = 0
    for M in range(2, 25):
        for J in range(1, M):
            r = M - J
            best, _, _ = read_cell(J, r, kmax, xmax)
            want, _ = flatten.ph_of(table, M, J)
            cells += 1
            if best[0] != want:
                diff += 1
    return cells, diff


def walk_all(J, r):
    """Every admissible pure product's coefficients with its shape."""
    for k in range(r):
        room = r - 1 - k
        stack = [(flatten.pow_xm1(J + k), (), room, room)]
        while stack:
            p, parts, left, cap = stack.pop()
            yield p, k, parts
            if len(parts) >= J + k:
                continue
            for x in range(min(left, cap), 0, -1):
                stack.append((mul_block(p, x + 1), parts + (x + 1,),
                              left - x, x))


def test_family(c, pmax=16):
    """The best of |<c, x^k w>|/|w|_1 over w = (1 + x)^p (1 - x)^q, q = 0, 1,
    p <= pmax, k within one of centring: a lower bound on height(c)."""
    best = Fraction(0)
    n = len(c)
    for q in (0, 1):
        for p in range(pmax + 1):
            v = [math.comb(p, i) for i in range(p + 1)]
            if q:
                v = [a - b for a, b in zip(v + [0], [0] + v)]
            l1 = sum(abs(a) for a in v)
            k0 = (n - len(v)) // 2
            for k in (k0 - 1, k0, k0 + 1):
                if 0 <= k and k + len(v) <= n:
                    t = Fraction(abs(sum(c[k + i] * a
                                         for i, a in enumerate(v))), l1)
                    best = max(best, t)
    return best


def reads():
    """READS, not verdicts: the middle coefficient and the test family as
    lower bounds over every non-champion pure product, against F_e's
    circle mean, at J*."""
    for r in (11, 13, 15, 17, 18, 19, 21, 23, 24, 25):
        J, e = theorem_J(r), theorem_e(r)
        hch = circlemean.exact_Hch(r, J)
        mean = float(mp.exp(circlemean.log_mean(r, J, e))) / hch
        mid = tst = None
        o3 = []
        for c, k, parts in walk_all(J, r):
            if order(parts) == r - 3 and r % 2:
                o3.append((Fraction(flatten.height(c), hch), name(k, parts)))
            if (k, parts) == champion(r):
                continue
            n = len(c) - 1
            m = Fraction(max(abs(c[n // 2]), abs(c[(n + 1) // 2])), hch)
            t = test_family(c) / hch
            if mid is None or m < mid[0]:
                mid = (m, name(k, parts))
            if tst is None or t < tst[0]:
                tst = (t, name(k, parts), flatten.height(c) / hch)
        print("READ r = %2d e = %d: mean(F_e)/H_ch %.4f; least middle"
              " coefficient/H_ch %.4f (%s); least test/H_ch %.4f at %s"
              " (its height/H_ch %.4f)"
              % (r, e, mean, mid[0], mid[1], tst[0], tst[1], tst[2]),
              flush=True)
        if o3:
            print("        order r - 3 heights/H_ch: %s"
                  % ", ".join("%s %s" % (q, n) for q, n in sorted(o3)))


def main():
    t0 = time.time()
    cells, diff = c1()
    print("C1 restricted ph against pure_table(23): %d cells, %d differ"
          % (cells, diff), flush=True)
    pc, pdiff = c1(kmax=0)
    print("C1 PLANT (k pinned to 0): %d cells, %d differ  [%.1f s]"
          % (pc, pdiff, time.time() - t0), flush=True)
    qc, qdiff = c1(xmax=2)
    print("C1 PLANT (parts at most 3): %d cells, %d differ  [%.1f s]"
          % (qc, qdiff, time.time() - t0), flush=True)
    ok = diff == 0 and qdiff > 0
    c2 = c3 = True
    worst = None
    rows = []
    for r in range(5, 26):
        e = theorem_e(r)
        J = theorem_J(r)
        best, hch, by_s = read_cell(J, r)
        c2 &= hch == circlemean.exact_Hch(r, J)
        hf = circlemean.exact_height_F(r, J, e)
        fe = {t: circlemean.exact_height_F(r, J, t)
              for t in range(max(0, e - 2), e + 4)}
        eb = min(fe, key=fe.get)
        o2 = by_s.get(r - 2)
        nonch = min((v for v in by_s.values()
                     if (v[1], v[2]) != champion(r)), key=lambda v: v[0])
        print("r = %2d e = %d J* = %3d: ph/H_ch %.4f at %s (s = %d);"
              " F_e/H_ch %.4f, ph/F_e %.4f; best e %d (F/H_ch %.4f);"
              " least non-champion %s /H_ch %.4f; order r-2 %s /H_ch %.4f"
              " (Gauss %.4f)  [%.1f s]"
              % (r, e, J, best[0] / hch, name(best[1], best[2]),
                 order(best[2]), hf / hch, best[0] / hf, eb, fe[eb] / hch,
                 name(nonch[1], nonch[2]), nonch[0] / hch,
                 name(o2[1], o2[2]), o2[0] / hch,
                 math.sqrt(J / (4 * (r - 1))), time.time() - t0),
              flush=True)
        print("        least per order s/H_ch: %s"
              % ", ".join("s=%d %.3f" % (s, by_s[s][0] / hch)
                          for s in sorted(by_s, reverse=True)[:5]),
              flush=True)
        h, v, nodes = flatten.least_height(J + r, J)
        b, d, lam = flatten.lll(flatten.basis(J + r, J))
        E = prunerr.err_bound(d, lam, r)[0]
        whole = 3 * E + 2 * E ** 2 + mp.mpf(2) ** -51 < mp.mpf(10) ** -9
        c3 &= (flatten.clears(v, J) and flatten.height(v) == h
               and h <= hf and whole)
        print("        C3 lattice h/H_ch %.4f, h/ph %.4f, E %.2e (complete %s),"
              " %d nodes" % (h / hch, h / best[0], float(E), whole, nodes),
              flush=True)
        rows.append((r, e, J, best[0], hch, hf, h))
        if e >= 2:
            q = best[0] / hf
            if worst is None or q < worst[0]:
                worst = (q, r)
    print("C2 the champion's height equals circlemean.exact_Hch: %s" % c2)
    print("C3 the lattice witness clears J, its height is h <= height(F_e),"
          " the enumeration complete by its bound, every rank: %s" % c3)
    fails = [r for r, e, J, ph, hch, hf, h in rows if h < ph]
    print("cells (J* + r, J*) with h < ph: r = %s" % fails)
    k1 = [r for r, e, J, ph, hch, hf, h in rows if e >= 2 and hf >= ph]
    print("K1 least ph/height(F_e) over the e >= 2 arm: %.4f at r = %d;"
          " K1 %s" % (worst[0], worst[1],
                      "fires at %s" % k1 if k1 else "misses"))
    reads()
    print("%.1f s" % (time.time() - t0))
    return not (ok and c2 and c3) or bool(k1)


if __name__ == "__main__":
    sys.exit(1 if main() else 0)