compbound.py

"""compbound.py -- at the takeover theorem's depth J*, does a closed-form
lower bound put every non-champion pure product above F_e's circle mean?

QUESTION. At J* = (4e - 1) l - 1, l = r - 2, height(F_e) < H_ch (the takeover
theorem) and the cell (J* + r, J*) fails once every other
admissible pure product is above F_e. takeph.py read that by enumeration
to r = 25. A bound that sees a shape only through a few numbers would
carry it to every rank; at which orders s does this one, at e = 3, 4?

BY HAND (re-derived before the code). A shape: k >= 0, parts g_i >= 2,
k + sum(g_i - 1) <= r - 1; s = #even parts, o = #odd parts, Pi = prod
over even parts of g/2. After x -> -x, [g](-x) = (1 + x^g)/(1 + x) at odd
g and (1 - x)[g/2](x^2) at even g, so up to sign
  P = (1 + x)^(J + k - o) (1 - x)^s prod_even [g/2](x^2) prod_odd (1 + x^g).
On x = e^(2iu), P = e^(iuN) (-i)^s 2^(J+k+s) Pi cos^(J+k-o) u sin^s u Phi,
Phi = prod_even D_(g/2)(2u) prod_odd cos(gu), D_j(y) = sin(jy)/(j sin y).
The test w = x^m (1 + x)^p (1 - x)^q, 2m + p + q = N, q = s mod 2,
p = J + k mod 2 (the centring itself: 2m + p + q = N) gives
  |<P, w>| = (2^(J+k+s+p+q)/pi) Pi B(al, be) |E_nu Phi|,
  al = (s + q + 1)/2, be = (J + k - o + p + 1)/2,
nu ~ cos^(2be - 1) sin^(2al - 1) on [0, pi/2]. Each factor of Phi has
|f| <= 1, so 1 - Phi <= sum (1 - f); 1 - D_j(y) <= (j^2 - 1) y^2/6, 1 -
cos(gu) <= g^2 u^2/2, u <= tan u, E_nu tan^2 u = al/(be - 1). With T =
sum_even (g^2 - 4)/6 + sum_odd g^2/2,
  height(P) >= LB = (2^(J+k+s+p+q)/pi) Pi B(al, be) (1 - T al/(be - 1))_+
                    / |w|_1,
|w|_1 = 2^p (q = 0), 2 C(p, floor(p/2)) (q = 1). The bound must die at one
odd part g with g^2/2 >~ 4e - 1: a part 5 at e = 3, a part 7 at e = 4.

PREDICTIONS (frozen before the code). P1 at e = 3 (odd r = 21..33) every
order s >= r - 4 clears, least LB/mean(F_e) >= 1.79 (order r - 2). P2 at
e = 3 some order s <= r - 5 fails, the first failure (largest s) a shape
with a part 5. P3 at e = 4 (odd r = 35..51) the largest failing order is
r - 7 or r - 6, its shape with a part 7 or two parts 5. P4 the direct
scan (exact ph, r = 27 timed, then odd r to 35 if r = 27 runs under 30 s)
reads height(F_e) < ph at each rank run.
KILL. K1 at e = 3 a least LB/mean <= 1 at an order s >= r - 4. K2 the
orders failing, printed per rank of the e = 4 arm: non-empty kills the
bound as the whole proof. K3 a rank of P4's arm with height(F_e) >= ph.
CONTROLS, read first. C1 at r = 11, 13, 21, 25, every shape: LB <= exact
height; at T = 0 the closed form equals |<P, w>|/|w|_1 in integers at
p = p0, p0 + 2, p0 + 10 (relative 1e-9). PLANTS: Pi doubled must give LB
> height somewhere; m moved off centre by one must break the equality at
every T = 0 shape. (Frozen was p of the wrong parity; replaced before the
run: 2m + p + q = N makes p = J + k mod 2 the centring itself.) The T = 0
shapes have Pi = 1, o = 0, Phi = 1: Pi, o and the remainder are tested
only through LB <= height. C2 mean(F_e)/H_ch against takeph's READ
0.7370, 0.7016, 0.6709 (r = 21, 23, 25); mean >= exact height(F_e). C3 at
r = 21, 25 the frontier sweep (a key off the frontier is evaluated only
when the frontier key dominating it fails) against every key evaluated:
the same failing keys and per-order least.
Floats: lgamma in log space; a ratio within 1e-6 above 1 counts as
failing.

FINDINGS (copied from the printed output).
  (1) C1 holds at r = 11, 13, 21, 25: LB <= height at every shape, the
      T = 0 closed form equal to the integer pairing at all 195, 270, 690,
      972 (shape, p) pairs; the m plant caught at every T = 0 shape, the
      Pi plant at 391, 775, 8817, 21939 shapes. C2 reproduces 0.7370,
      0.7016, 0.6709, mean >= height(F_e) (0.5292, 0.5134, 0.4966).
  (2) K1 misses: at e = 3 (odd r = 21..33) orders r - 2, r - 3, r - 4 clear,
      least LB/mean 1.7226 (r = 21) .. 2.2225 (r = 33), always at order r - 2
      (k=0 2^(r-2), p = 10). P1's 1.79 misses at r = 21 (1.7226;
      1.8154 at r = 23): the centred test with q = s mod 2 sits below takeph's best,
      which also shifts off centre and tries the other q.
  (3) K2 fires, as P2 and P3 said: at e = 3 the first failure is order
      r - 5 at the part 5 (k=0 2^(r-5) 5, LB/mean 0.465 .. 0.166 over r = 21..
      33); order r - 6 clears to r = 29 (1.013) and fails from r = 31; at
      e = 4 (odd r = 35..51) orders r - 2 .. r - 6 clear (least 1.9656 ..
      2.4075, at r - 2) and the first failure is r - 7 at the part 7 (LB/mean
      0.014 .. 0.000). At each failing order printed (r - 5 .. r - 10) the
      least shape is the order-s champion (1 + x)^J (1 - x)^s with one odd
      part g = 5, 7 or 9 traded in.
  (4) K3 misses (P4 holds): the direct scan at odd r = 27..35 reads ph =
      H_ch, the champion, and ph/height(F_e) = 2.0854, 2.1613, 2.2402,
      2.3216, 1.7969, so J_hi(r) >= J* there, computed.
READS. The failing part sits where the hand's remainder put it (g^2/2 >~
4e - 1) and near the half wave of the champion's oscillation: u
concentrates at u0 = arcsin(1/(2 sqrt(e))): g u0 = 1.46 at g = 5, e = 3 and
1.77 at g = 7, e = 4, near pi/2 (cos 0.106, -0.197), so E cos(gu) itself
is small there, not only its bound (a TRANSPLANT from the Gaussian
concentration, unread; sin^2 u0 = 1/(4e) is U's mean at p = 0 only
nearly (0.0756 against 0.0833 at 2^16 5, r = 21), and a failing shape's
best p moves it).

  (5) The failing keys (a key failing when LB at its least Pi does)
      number 14, 26, 45, 87, 154, 260, 430 at odd r = 21..33 and 141, 221,
      371, 584, 895, 1376, 2035, 2994, 4332 at odd r = 35..51; C3 holds at
      r = 21, 25 (14, 45). The least-Pi shape carries an odd part >= 5 at
      all but 0, 0, 0, 1, 2, 3, 6 of them (r = 21..33) and 0, 0, 0, 0, 1,
      4, 7, 15, 32 (r = 35..51): parts 3 only (k=0 2^18 3 3 3 3 at r = 27)
      or a large even part (k=0 2^21 12 and 2^19 14 at r = 33, 2^27 16 at
      r = 43, 2^31 14 3 at r = 47).

RUN RECORD. 74.2 s (the scan 36.5 s of it), peak commit 71.0 MB (one k
of the key table at a time).
"""
import math
import os
import sys
import time

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

LN2 = math.log(2)
LNPI = math.log(math.pi)


def stats(k, parts):
    """(s, o, Pi, 6T) of a shape; 6T is an integer."""
    s = o = 0
    pi = 1
    t6 = 0
    for g in parts:
        if g % 2:
            o += 1
            t6 += 3 * g * g
        else:
            s += 1
            pi *= g // 2
            t6 += g * g - 4
    return s, o, pi, t6


def keys(r, k):
    """Every non-champion shape with k extra (x - 1)'s grouped by (s, o,
    6T), keeping the least Pi and one shape attaining it: the bound is Pi
    times a function of the key (s, k, o, 6T). One k at a time, so the
    table never holds every k at once."""
    top = r - 1
    out = {}
    stack = [((), top - k, top - k, 0, 0, 1, 0)]
    while stack:
        parts, left, cap, s, o, pi, t6 = stack.pop()
        if not (k == 0 and s == top):
            key = (s, o, t6)
            if key not in out or pi < out[key][0]:
                out[key] = (pi, parts)
        for x in range(min(left, cap), 0, -1):
            g = x + 1
            if g % 2:
                stack.append((parts + (g,), left - x, x, s, o + 1, pi,
                              t6 + 3 * g * g))
            else:
                stack.append((parts + (g,), left - x, x, s + 1, o,
                              pi * (g // 2), t6 + g * g - 4))
    return out


def log_w1(p, q):
    if q == 0:
        return p * LN2
    return LN2 + math.lgamma(p + 1) - math.lgamma(p // 2 + 1) \
        - math.lgamma(p - p // 2 + 1)


def log_lb(J, s, k, o, pi, t6, p, plant=False):
    """log LB at test exponent p (-inf where 1 - T tau <= 0)."""
    q = s % 2
    al = (s + q + 1) / 2
    be = (J + k - o + p + 1) / 2
    f = 1 - t6 / 6 * al / (be - 1)
    if f <= 0:
        return -math.inf
    lpi = math.log(2 * pi if plant else pi)
    return ((J + k + s + p + q) * LN2 - LNPI + lpi + math.lgamma(al)
            + math.lgamma(be) - math.lgamma(al + be) + math.log(f)
            - log_w1(p, q))


def best_lb(J, s, k, o, pi, t6, plant=False, patience=40):
    """(log LB, p) maximised over p = J + k mod 2, scanned upward until the
    bound has fallen for `patience` steps past its best."""
    q = s % 2
    al = (s + q + 1) / 2
    # 1 - T tau > 0 needs be > 1 + T al: start at the least such p
    pmin = math.floor(2 * (t6 / 6 * al + 1) - (J + k - o + 1)) + 1
    p = max((J + k) % 2, pmin)
    p += (p - J - k) % 2
    best = (-math.inf, p)
    since = 0
    while since < patience:
        v = log_lb(J, s, k, o, pi, t6, p, plant)
        if v > best[0]:
            best, since = (v, p), 0
        else:
            since += 1
        p += 2
    return best


def int_pairing(c, J, k, parts, p, shift=0):
    """(|<P, w>|, |w|_1) in integers, P the flipped coefficients of c, w =
    x^m (1 + x)^p (1 - x)^q centred (shift: the PLANT, m moved)."""
    s = takeph.order(parts)
    q = s % 2
    cf = [v if i % 2 == 0 else -v for i, v in enumerate(c)]
    n = len(cf) - 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)]
    m = (n - p - q) // 2 + shift
    tot = sum(cf[m + i] * a for i, a in enumerate(v) if 0 <= m + i <= n)
    return abs(tot), sum(abs(a) for a in v)


def c1(r):
    """LB <= height at every shape; T = 0 closed form against integers;
    the two plants."""
    J = takeph.theorem_J(r)
    worst = -math.inf
    planted = 0
    eq = 0
    eqbad = 0
    off = 0
    offcaught = 0
    for c, k, parts in takeph.walk_all(J, r):
        if (k, parts) == takeph.champion(r):
            continue
        s, o, pi, t6 = stats(k, parts)
        lh = math.log(flatten.height(c))
        lb = best_lb(J, s, k, o, pi, t6)[0]
        worst = max(worst, lb - lh)
        if best_lb(J, s, k, o, pi, t6, plant=True)[0] > lh + 1e-9:
            planted += 1
        if t6 == 0:
            p0 = (J + k) % 2
            for p in (p0, p0 + 2, p0 + 10):
                num, w1 = int_pairing(c, J, k, parts, p)
                want = math.exp(log_lb(J, s, k, o, pi, t6, p) - lh)
                got = num / w1 / flatten.height(c)
                eq += 1
                if abs(got - want) > 1e-9 * want:
                    eqbad += 1
            num, w1 = int_pairing(c, J, k, parts, p0, shift=1)
            want = math.exp(log_lb(J, s, k, o, pi, t6, p0) - lh)
            off += 1
            if abs(num / w1 / flatten.height(c) - want) > 1e-9 * want:
                offcaught += 1
    return worst, planted, eq, eqbad, off, offcaught


def bound_rank(r, lbar=None):
    """Per order s the least LB/mean(F_e) over non-champion shapes; a key
    fails at or below the log bar lbar (default the mean)."""
    J, e = takeph.theorem_J(r), takeph.theorem_e(r)
    lmean = float(circlemean.log_mean(r, J, e))
    lhch = math.log(circlemean.exact_Hch(r, J))
    bar = (lmean if lbar is None else lbar) + math.log(1 + 1e-6)
    # at fixed (s, k, o) LB rises with Pi and falls with T: the frontier,
    # each key with a Pi below every key of larger T, holds each order's
    # least; a key off it is above the frontier key that dominates it, so
    # it is evaluated only when that key fails
    nkeys = nfront = 0
    by_s = {}
    failing = []
    for k in range(r):
        groups = {}
        for (s, o, t6), (pi, parts) in keys(r, k).items():
            groups.setdefault((s, o), []).append((t6, pi, parts))
            nkeys += 1
        for (s, o), lst in groups.items():
            low, low_v = math.inf, None
            for t6, pi, parts in sorted(lst, reverse=True):
                if pi < low:
                    v, p = best_lb(J, s, k, o, pi, t6)
                    low, low_v = pi, v
                    nfront += 1
                    if s not in by_s or v < by_s[s][0]:
                        by_s[s] = (v, p, k, parts)
                elif low_v > bar:
                    continue
                else:
                    v, p = best_lb(J, s, k, o, pi, t6)
                if v <= bar:
                    failing.append((s, k, o, t6, pi, parts))
    rows = {}
    for s in sorted(by_s, reverse=True):
        v, p, k, parts = by_s[s]
        rows[s] = (math.exp(v - lmean) if v > -math.inf else 0.0, p,
                   takeph.name(k, parts))
    return J, e, math.exp(lmean - lhch), (nkeys, nfront), rows, failing


def main():
    t0 = time.time()
    ok = True
    for r in (11, 13, 21, 25):
        worst, planted, eq, eqbad, off, offc = c1(r)
        print("C1 r = %d: max log(LB/height) %.4f; T = 0 equality %d of %d"
              " off; PLANT m moved by one caught %d of %d; PLANT Pi doubled"
              " LB > height at %d shapes  [%.1f s]"
              % (r, worst, eqbad, eq, offc, off, planted,
                 time.time() - t0), flush=True)
        ok &= worst <= 1e-9 and eqbad == 0 and offc == off and planted > 0
    print("C1 %s" % ("holds" if ok else "FAILS"), flush=True)
    for r, want in ((21, 0.7370), (23, 0.7016), (25, 0.6709)):
        J, e = takeph.theorem_J(r), takeph.theorem_e(r)
        m = float(mp.exp(circlemean.log_mean(r, J, e))) \
            / circlemean.exact_Hch(r, J)
        hf = circlemean.exact_height_F(r, J, e) / circlemean.exact_Hch(r, J)
        print("C2 r = %d: mean(F_e)/H_ch %.4f (takeph %.4f), height(F_e)/H_ch"
              " %.4f, mean >= height %s" % (r, m, want, hf, m >= hf),
              flush=True)
        ok &= abs(m - want) < 6e-5 and m >= hf
    for r in (21, 25):
        J, e, mr, nk, rows, failing = bound_rank(r)
        lmean = float(circlemean.log_mean(r, J, e))
        bar = lmean + math.log(1 + 1e-6)
        every, least = 0, {}
        for k in range(r):
            for (s, o, t6), (pi, parts) in keys(r, k).items():
                v = best_lb(J, s, k, o, pi, t6)[0]
                every += v <= bar
                least[s] = min(least.get(s, math.inf), v)
        same = every == len(failing) and all(
            abs(math.exp(least[s] - lmean) - rows[s][0]) < 1e-12 for s in rows)
        print("C3 r = %d: every key evaluated, %d failing (frontier sweep %d),"
              " per-order least the same %s" % (r, every, len(failing), same),
              flush=True)
        ok &= same
    k1 = []
    for r in list(range(21, 34, 2)) + list(range(35, 52, 2)):
        J, e, mr, nk, rows, failing = bound_rank(r)
        fail = [s for s, (q, p, nm) in rows.items() if q <= 1 + 1e-6]
        top = [s for s in rows if s >= r - 4]
        lt = min(top, key=lambda s: rows[s][0])
        print("r = %2d e = %d J* = %3d: %s keys, frontier; mean/H_ch %.4f;"
              " least over"
              " s >= r-4 %.4f at s = r-%d (%s, p = %d); orders failing: %s"
              "  [%.1f s]"
              % (r, e, J, nk, mr, rows[lt][0], r - lt, rows[lt][2],
                 rows[lt][1], ["r-%d" % (r - s) for s in fail],
                 time.time() - t0), flush=True)
        print("        by order: %s" % ", ".join(
            "r-%d %.3f (%s)" % (r - s, q, nm)
            for s, (q, p, nm) in sorted(rows.items(), reverse=True)[:9]),
            flush=True)
        if fail:
            s = max(fail)
            print("        first failure r-%d: %s" % (r - s, rows[s][2]),
                  flush=True)
        rest = [f for f in failing if not any(g % 2 and g >= 5 for g in f[5])]
        print("        failing keys %d, %d of them without an odd part >= 5"
              " in the least-Pi shape%s" % (
                  len(failing), len(rest), (": " + ", ".join(
                      takeph.name(f[1], f[5]) for f in rest))
                  if 0 < len(rest) <= 8 else ""), flush=True)
        if e == 3 and rows[lt][0] <= 1 + 1e-6:
            k1.append(r)
    print("K1 %s" % ("fires at %s" % k1 if k1 else "misses"), flush=True)
    k3 = []
    for r in range(27, 36, 2):
        t1 = time.time()
        J, e = takeph.theorem_J(r), takeph.theorem_e(r)
        best, hch, by_s = takeph.read_cell(J, r)
        hf = circlemean.exact_height_F(r, J, e)
        print("SCAN r = %d e = %d J* = %d: ph/H_ch %.4f at %s; ph/height(F_e)"
              " %.4f  [%.1f s]"
              % (r, e, J, best[0] / hch, takeph.name(best[1], best[2]),
                 best[0] / hf, time.time() - t1), flush=True)
        if hf >= best[0]:
            k3.append(r)
        if r == 27 and time.time() - t1 > 30:
            print("SCAN stops: r = 27 ran over 30 s")
            break
    print("K3 %s" % ("fires at %s" % k3 if k3 else "misses"))
    print("%.1f s" % (time.time() - t0))
    return not ok or bool(k1) or bool(k3)


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