"""exacttest.py -- at the takeover theorem's depth J*, does the exact
centred test put every shape compbound's bound misses above F_e's circle
mean?

QUESTION. compbound.py bounds a pure product's centred pairing with
x^m (1 + x)^p (1 - x)^q by Pi B(al, be) (1 - T al/(be - 1)) and fails at
its failing keys. Every shape whose bound fails (its Pi below its key's
threshold, not only the least-Pi one): does the exact centred test put it
above F_e's mean, at odd r = 21..51 (e = 3, 4)?

BY HAND (re-derived before the code). D_j(2u) = sin(2ju)/(j sin 2u) =
U_(2j-1)(c)/(2jc), c = cos u, since sin(2ju)/sin u = U_(2j-1)(cos u)
(j = 1: 1; j = 2: U_3/(4c) = 2c^2 - 1 = cos 2u). cos(gu) = T_g(c). So
Pi Phi = 2^-s Psi(c), Psi = prod_even U_(g-1)(c)/c prod_odd T_g(c), an
integer polynomial with Psi(-c) = (-1)^o Psi(c). With A = J + k - o + p,
S = s + q, I(a, b) = int_0^(pi/2) cos^a u sin^b u du,
  |<P, w>| = (2^(J+k+p+q+1)/pi) |sum_i psi_i I(A + i, S)|.
p = J + k mod 2 makes A = o mod 2 and every i with psi_i != 0 too, so A + i
and S are even, I(2a, 2b) = (pi/2)(2a - 1)!!(2b - 1)!!/(2a + 2b)!!, and the
pairing is the integer 2^(J+k+p+q) |sum_i psi_i (A+i-1)!!(S-1)!!/(A+i+S)!!|
((1 + x)^2, p = q = 0: 4 (1)/(2) = 2 = C(2, 1)). The sum runs by integer
Horner on I(A+i0+2t, S)/I(A+i0, S) = prod_(v<t) (A+i0+1+2v)/(A+i0+S+2+2v):
psi's coefficients reach ~2^r and cancel to O(1), so floats die there;
lgamma carries only the prefactor. X = |<P, w>|/|w|_1, maximised over
p = J + k mod 2. A shape fails the bound when Pi F(key) <= mean, F the
bound at Pi = 1, so its key's least Pi fails and Pi <= Pi_min exp(bar -
log LB_min).

PREDICTIONS (frozen before the code). P1 at e = 3 (odd r = 21..33) every
failing shape with an odd part >= 5 or parts 3 only clears, least X/mean
>= 1.5. P2 where a failing shape has a large even part (r = 33, 43..51)
the rank's least X/mean is at one, in 1.0 .. 1.5. P3 at e = 4 (odd r =
35..51) the least >= 1.2.
KILL. Printed per rank: the least X/mean over every failing shape, with
its shape, e = 3 and e = 4 arms. Any <= 1 + 1e-6 fires: the centred test
is not the whole proof there, and those shapes owe the off-centre test.
CONTROLS, read first. C1 at r = 11, 13, 21, 25, every shape with a part
above 2, p = p0, p0 + 2, p0 + 10: the closed form equals compbound's
int_pairing as integers, exactly. PLANTS: the largest odd part's T_g as
T_(g-2) must break the equality at every shape with an odd part >= 5; the
largest even part's U_(g-1) as U_(g-3) at every shape with an even part
>= 4. C2 X >= LB at the bound's own p at every failing shape. C3 at
r = 21 every T = 0 shape: the best X equals compbound's best LB (Phi = 1,
the scan and the prefactor read against a known answer). C4 (added after
the run: C3's shapes never enter the Horner loop) at r = 11, 13, every
shape with a part above 2 and the three p: log X, the verdict's own
function, against log(exact pairing/|w|_1), zero pairings matched as
zero. Frozen C2 also had X <= height at every shape; C1's equality
carries it, so it does not run.

FINDINGS (copied from the printed output).
  (1) C1 holds at r = 11, 13, 21, 25: the closed form equals the integer
      pairing exactly at all 1071, 2397, 32247, 95310 (shape, p) pairs;
      the odd plant is caught at 108, 288, 5922, 19811 of as many shapes,
      the even plant at 174, 431, 7509, 23993; the enumeration's shapes
      are takeph.walk_all's (422, 889, 10979, 32094). Read at p0 alone,
      as frozen (the first run's print), the plants missed 6, 1, 0, 4
      shapes (k=0 2^3 7, 2^4 6,
      k=1 2^1 17, ...): there w has one or two coefficients, and the
      planted twin shares them. Every miss is caught at p0 + 2 and
      p0 + 10, so a plant now reads over the three p (changed after the
      run). C3 holds at r = 21 (230 T = 0 shapes). C4 holds: 3468 (shape,
      p), 7 zero, max |difference| 3.6e-14; with the Horner step's
      numerator factor off by one, C4 reads inf and C3 still holds.
      C2 holds: X >= LB at
      the bound's p at every failing shape.
  (2) The KILL misses: every failing shape clears at every odd r =
      21..51. At e = 3 (r = 21..33) the least X/mean is 2.6738, 2.6630,
      2.6617, 2.6673, 2.6780, 2.6925, 2.7100, always at the part-5 shape
      k=0 2^(r-5) 5 with p = 0, its central coefficient. At e = 4 (r =
      35..51) it is 7.0171 (k=0 2^26 7 3, p = 44), 7.6582 (k=0 2^28 7 3),
      then 9.8599 .. 14.0320 at k=0 2^(r-7) 7 with p = 0. The failing
      shapes number 14 .. 442 at e = 3 and 142 .. 5046 at e = 4. With
      compbound's bound clearing every other shape, every non-champion
      pure product is above F_e's circle mean at J*, so J_hi(r) >= J* at
      every odd r = 21..51, computed.
  (3) P1 holds: least 2.66 >= 1.5; the unrecorded quadrature's 2.67 at
      r = 21 reads 2.6738. P2 fails: the large even parts stand far
      above, 287.75 (k=0 2^21 12, r = 33) and 1144 .. 31076 at r = 43..51.
      Their bound failed on the remainder (g^2 - 4)/6, not on the test.
      Parts 3 only: least 20.36 (r = 29). P3 holds: 7.0171 >= 1.2.
READS. The e = 3 least stays flat in r (2.66 .. 2.71), at the p = 0
central coefficient of one shape. compbound's READ (E cos(gu) small near
g u0 = pi/2) does not bind: the exact test clears there by 2.66 of the
mean.

RUN RECORD. 46.2 s (C1 4.6 s of it), peak commit 71.3 MB.
"""
import math
import os
import sys
import time
from fractions import Fraction

os.environ.setdefault("OPENBLAS_NUM_THREADS", "1")
import circlemean  # noqa: E402
import takeph  # noqa: E402
import compbound  # noqa: E402

LN2 = math.log(2)
LNPI = math.log(math.pi)
_T = {0: [1], 1: [0, 1]}
_U = {0: [1], 1: [0, 2]}


def _cheb(tab, n):
    while n not in tab:
        m = max(tab)
        a, b = tab[m], tab[m - 1]
        nxt = [0] + [2 * v for v in a]
        for i, v in enumerate(b):
            nxt[i] -= v
        tab[m + 1] = nxt
    return tab[n]


def mul(a, b):
    out = [0] * (len(a) + len(b) - 1)
    for i, x in enumerate(a):
        if x:
            for j, y in enumerate(b):
                out[i + j] += x * y
    return out


def psi(parts, plant=None):
    """Psi = prod_even U_(g-1)(c)/c prod_odd T_g(c), low to high. plant
    'odd': the largest odd part's T_g as T_(g-2); 'even': the largest even
    part's U_(g-1) as U_(g-3)."""
    odd = [g for g in parts if g % 2]
    even = [g for g in parts if g % 2 == 0]
    out = [1]
    for i, g in enumerate(odd):
        n = g - 2 if plant == "odd" and i == 0 and g >= 5 else g
        out = mul(out, _cheb(_T, n))
    for i, g in enumerate(even):
        n = g - 3 if plant == "even" and i == 0 and g >= 4 else g - 1
        out = mul(out, _cheb(_U, n)[1:])
    return out


_DF = {-1: 1, 0: 1}


def dfact(n):
    while n not in _DF:
        m = max(_DF)
        _DF[m + 1] = (m + 1) * _DF[m - 1]
    return _DF[n]


def exact_pairing(J, k, parts, p, plant=None):
    """|<P, w>| by the closed form, as an exact integer (C1)."""
    s = takeph.order(parts)
    o = len(parts) - s
    q = s % 2
    A, S = J + k - o + p, s + q
    ps = psi(parts, plant)
    M = A + len(ps) - 1 + S
    top = dfact(M)
    num = sum(v * dfact(A + i - 1) * (top // dfact(A + i + S))
              for i, v in enumerate(ps) if v)
    return abs(Fraction(num * dfact(S - 1) * 2 ** (J + k + p + q), top))


def log_x(J, k, s, o, ps, p):
    """log X at test exponent p (-inf where the pairing vanishes)."""
    q = s % 2
    A, S = J + k - o + p, s + q
    i0 = o % 2
    co = ps[i0::2]
    H = len(co) - 1
    N, D = co[H], 1
    for t in range(H - 1, -1, -1):
        y = A + i0 + S + 2 + 2 * t
        N = co[t] * D * y + N * (A + i0 + 1 + 2 * t)
        D *= y
    if N == 0:
        return -math.inf
    a = A + i0
    li = (math.lgamma((a + 1) / 2) + math.lgamma((S + 1) / 2)
          - math.lgamma((a + S) / 2 + 1) - LN2)
    return ((J + k + p + q + 1) * LN2 - LNPI + li + math.log(abs(N))
            - math.log(D) - compbound.log_w1(p, q))


def best_x(J, k, parts, patience=40):
    """(log X, p) maximised over p = J + k mod 2, scanned upward until X
    has fallen for `patience` steps past its best."""
    s = takeph.order(parts)
    o = len(parts) - s
    ps = psi(parts)
    p = (J + k) % 2
    best = (-math.inf, p)
    since = 0
    while since < patience:
        v = log_x(J, k, s, o, ps, p)
        if v > best[0]:
            best, since = (v, p), 0
        else:
            since += 1
        p += 2
    return best


def shapes(r, k):
    """Every non-champion shape with k extra (x - 1)'s: (parts, s, o, Pi,
    6T), parts in falling order."""
    top = r - 1
    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):
            yield parts, s, o, pi, t6
        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))


def failing_shapes(r, lbar=None):
    """Every shape whose bound fails at the log bar lbar (default the
    mean): (k, parts, s, o, pi, t6, p_LB)."""
    J, e, mr, nk, rows, failing = compbound.bound_rank(r, lbar)
    lmean = float(circlemean.log_mean(r, J, e))
    bar = (lmean if lbar is None else lbar) + math.log(1 + 1e-6)
    thr = {}
    for s, k, o, t6, pi, parts in failing:
        v, p = compbound.best_lb(J, s, k, o, pi, t6)
        thr[(s, k, o, t6)] = (math.log(pi) + bar - v, p)
    out = []
    for k in range(r):
        for parts, s, o, pi, t6 in shapes(r, k):
            t = thr.get((s, k, o, t6))
            if t and math.log(pi) <= t[0] + 1e-12:
                out.append((k, parts, s, o, pi, t6, t[1]))
    return J, e, lmean, len(failing), out


def kind(parts):
    if any(g % 2 and g >= 5 for g in parts):
        return "odd>=5"
    if any(g % 2 == 0 and g >= 4 for g in parts):
        return "even>=4"
    return "3s"


def c1(r):
    J = takeph.theorem_J(r)
    n = bad = 0
    po = pon = pe = pen = 0
    walked = 0
    for c, k, parts in takeph.walk_all(J, r):
        if (k, parts) == takeph.champion(r):
            continue
        walked += 1
        if max(parts, default=2) <= 2:
            continue
        p0 = (J + k) % 2
        want = {}
        for p in (p0, p0 + 2, p0 + 10):
            got = exact_pairing(J, k, parts, p)
            want[p] = compbound.int_pairing(c, J, k, parts, p)[0]
            n += 1
            bad += got != want[p]
        # a plant is caught when the pairing moves at one of the three p:
        # at p0, w has one or two coefficients, and a few shapes share
        # them with their planted twin
        if any(g % 2 and g >= 5 for g in parts):
            po += 1
            pon += any(exact_pairing(J, k, parts, p, "odd") != want[p]
                       for p in want)
        if any(g % 2 == 0 and g >= 4 for g in parts):
            pe += 1
            pen += any(exact_pairing(J, k, parts, p, "even") != want[p]
                       for p in want)
    enum = sum(1 for k in range(r) for _ in shapes(r, k))
    return n, bad, (pon, po), (pen, pe), (walked, enum)


def c4(rs):
    """log X, the verdict's own function, against log(exact pairing/|w|_1)
    at every shape with a part above 2 and p = p0, p0 + 2, p0 + 10, zero
    pairings matched as zero: (pairs, zeros, max |difference|)."""
    n = zero = 0
    worst = 0.0
    for r in rs:
        J = takeph.theorem_J(r)
        for k in range(r):
            for parts, s, o, pi, t6 in shapes(r, k):
                if max(parts, default=2) <= 2:
                    continue
                ps = psi(parts)
                p0 = (J + k) % 2
                for p in (p0, p0 + 2, p0 + 10):
                    got = log_x(J, k, s, o, ps, p)
                    ex = exact_pairing(J, k, parts, p)
                    n += 1
                    if ex == 0 or got == -math.inf:
                        zero += 1
                        worst = max(worst, 0.0 if ex == 0 and
                                    got == -math.inf else math.inf)
                        continue
                    want = math.log(ex) - compbound.log_w1(p, s % 2)
                    worst = max(worst, abs(got - want))
    return n, zero, worst


def main():
    t0 = time.time()
    ok = True
    for r in (11, 13, 21, 25):
        n, bad, (pon, po), (pen, pe), (walked, enum) = c1(r)
        print("C1 r = %d: closed form against the integer pairing, %d of %d"
              " (shape, p) differ; PLANT odd part caught %d of %d; PLANT"
              " even part caught %d of %d; shapes enumerated %d, walked %d"
              "  [%.1f s]"
              % (r, bad, n, pon, po, pen, pe, enum, walked, time.time() - t0),
              flush=True)
        ok &= (bad == 0 and pon == po and pen == pe and po and pe
               and enum == walked)
    r = 21
    J = takeph.theorem_J(r)
    n = bad = 0
    for k in range(r):
        for parts, s, o, pi, t6 in shapes(r, k):
            if t6:
                continue
            v = best_x(J, k, parts)[0]
            w = compbound.best_lb(J, s, k, o, pi, t6)[0]
            n += 1
            bad += abs(v - w) > 1e-9
    print("C3 r = 21: best X against best LB at the %d T = 0 shapes, %d off"
          % (n, bad), flush=True)
    ok &= bad == 0 and n > 0
    n, zero, worst = c4((11, 13))
    print("C4 r = 11, 13: log X by Horner against the exact pairing at %d"
          " (shape, p), %d of them zero, max |difference| %.1e"
          % (n, zero, worst), flush=True)
    ok &= worst < 1e-9 and n > 0
    print("C1, C3, C4 %s" % ("hold" if ok else "FAIL"), flush=True)
    k1 = []
    c2bad = 0
    for r in list(range(21, 34, 2)) + list(range(35, 52, 2)):
        t1 = time.time()
        J, e, lmean, nkey, fs = failing_shapes(r)
        least = {}
        below = []
        for k, parts, s, o, pi, t6, plb in fs:
            ps = psi(parts)
            if log_x(J, k, s, o, ps, plb) < \
                    compbound.log_lb(J, s, k, o, pi, t6, plb) - 1e-9:
                c2bad += 1
            v, p = best_x(J, k, parts)
            ratio = math.exp(v - lmean)
            kd = kind(parts)
            if kd not in least or ratio < least[kd][0]:
                least[kd] = (ratio, p, takeph.name(k, parts))
            if ratio <= 1 + 1e-6:
                below.append((ratio, takeph.name(k, parts), p))
        lo = min(least.values())
        print("r = %2d e = %d: %d failing keys, %d failing shapes; least"
              " X/mean %.4f at %s (p = %d); by kind %s; at or below 1: %d"
              "  [%.1f s]"
              % (r, e, nkey, len(fs), lo[0], lo[2], lo[1],
                 ", ".join("%s %.4f (%s)" % (kd, least[kd][0], least[kd][2])
                           for kd in sorted(least)), len(below),
                 time.time() - t1), flush=True)
        if below:
            below.sort()
            print("        below: %s" % ", ".join(
                "%s %.4f p = %d" % (nm, q, p) for q, nm, p in below[:8]),
                flush=True)
            k1.append(r)
    print("C2 X >= LB at the bound's p: %d failing shapes off" % c2bad)
    print("KILL %s" % ("fires at %s" % k1 if k1 else "misses"))
    print("%.1f s" % (time.time() - t0))
    return not ok or c2bad > 0 or bool(k1)


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