#!/usr/bin/env python3
"""Check of Proposition 8.19 (the band) in "Analytic q-inequalities and a proof of Bergeron's conjecture", Section 8.8.  Python 3 with python-flint.

Notation.  G(m, n) = [m+n choose m]_q.  For x >= 3 and A = xD + r = (x-1)E + s (0 <= r < x, 0 <= s < x-1),
    D_A = (x-1)[(x-r) G(x,D) + r G(x,D+1)] - x[(x-1-s) G(x-1,E) + s G(x-1,E+1)],
and R(x, A) says that D_A has nonnegative coefficients in degrees k <= A/2.  The base series is
    Gamma_{x,N} = prod_{m=2}^{x} (1 - q^{N+m}) / (1 - q^m)          (a power series in q),
which is the rational q-Catalan polynomial C_{x,N+1} when gcd(x, N+1) = 1.

Proposition 8.19 writes D_A = (1-q) Gamma_{x,N} Y at every breakpoint A of the band x^2 <= A <= 2x(x-1):
    rows   A = x(x+j),       0 <= j <= x-2:  N = x+j,  Y = x q^x W_j / (1 - q^{x+j+2}),
           W_j = (x-1)[j+1][x+j+2] - (j+1) q^{j+2} [x][x-1];
    duals  A = (x-1)(x+1+i), 1 <= i <= x-2:  N = x+i,  Y = (x-1) q^x V_i / (1 - q^{2x+i}),
           V_i = (x-1-i)[i+1][2x+i] + (i+1)[i][x+i+1] - (i+1) q^{i+1} [x][x-1].
With gamma_m the coefficients of Gamma_{x,N} (gamma_0 = 1, gamma_1 = 0, gamma_2 = 1), Abel summation gives
    [q^k] D_A = (Y_k - Y_{k-1} + Y_{k-2}) + sum_{m >= 3} (gamma_m - gamma_{m-1}) Y_{k-m},
and R follows on the band from
    (II)  Y >= 0,   (III)  Y_t - Y_{t-1} + Y_{t-2} >= 0,   both proved in the paper, and
    (M)   gamma_{m-1} <= gamma_m for 3 <= m <= floor((x-1)N/2) - 1,  for x <= N <= 2x-2
          (the case of Conjecture 8.21 that the band needs).
Each breakpoint needs these up to degree K = floor(A'/2), A' the next breakpoint (K = A/2 at the top A = 2x(x-1)).

For 3 <= x <= X the script checks (M) exactly, and (II), (III) up to degree K.  For x <= XI it also checks the
identity D_A = (1-q) Gamma Y modulo q^(K+1) and, directly, that D_A >= 0 in degrees <= K.

Usage:  python3 band_check.py [X] [XI] [XS]   defaults X = 60, XI = 30, XS = 3 (under a minute);
        XS resumes a long run at x = XS.
"""
import sys, time
from math import gcd
import flint
from flint import fmpz_poly, fmpz_series


def inv_prod(ms, K):
    """prod_{m in ms} 1/(1 - q^m) modulo q^(K+1)"""
    flint.ctx.cap = K + 1
    P = fmpz_poly([1])
    for m in ms:
        P = (P - P.left_shift(m)).truncate(K + 1)
    return fmpz_poly((1 / fmpz_series(P.coeffs(), prec=K + 1)).coeffs())


def num_prod(es, K):
    """prod_{e in es} (1 - q^e) modulo q^(K+1)"""
    P = fmpz_poly([1])
    for e in es:
        if e <= K:
            P = (P - P.left_shift(e)).truncate(K + 1)
    return P


def coeffs(P, K):
    c = [int(v) for v in P.coeffs()][:K + 1]
    return c + [0] * (K + 1 - len(c))


def qint(n):
    return fmpz_poly([1] * n) if n > 0 else fmpz_poly([])


def breakpoints(x):
    """(kind, t, A, N, Y-constant, weight polynomial, period L, K) for the band of x"""
    out = []
    for j in range(0, x - 1):
        A = x * (x + j)
        K = A // 2 if j == x - 2 else ((x - 1) * (x + j + 2)) // 2
        W = (x - 1) * qint(j + 1) * qint(x + j + 2) - (j + 1) * (qint(x) * qint(x - 1)).left_shift(j + 2)
        out.append(('row', j, A, x + j, x, W, x + j + 2, K))
    for i in range(1, x - 1):
        A = (x - 1) * (x + 1 + i)
        K = (x * (x + i)) // 2
        V = (x - 1 - i) * qint(i + 1) * qint(2 * x + i) + (i + 1) * qint(i) * qint(x + i + 1) \
            - (i + 1) * (qint(x) * qint(x - 1)).left_shift(i + 1)
        out.append(('dual', i, A, x + i, x - 1, V, 2 * x + i, K))
    return out


def Y_series(x, c, Wp, L, K):
    """Y = c q^x W / (1 - q^L) modulo q^(K+1), as a list"""
    W = [int(v) for v in Wp.coeffs()]
    Y = [0] * (K + 1)
    for l in range(0, (K - x) // L + 1):
        base = x + l * L
        for d, w in enumerate(W):
            if base + d > K:
                break
            Y[base + d] += c * w
    return Y


def check_x(x, identities):
    bps = breakpoints(x)
    Kmax = max(b[-1] for b in bps)
    # (M): gamma nondecreasing on [3, floor(mu)-1] for x <= N <= 2x-2
    Zx = inv_prod(range(2, x + 1), Kmax)
    gam, mfail = {}, []
    for N in range(x, 2 * x - 1):
        top = ((x - 1) * N) // 2 - 1
        Kn = max(top + 1, max(b[-1] for b in bps if b[3] == N))
        Gm = num_prod([N + m for m in range(2, x + 1)], Kn).mul_low(Zx, Kn + 1)
        dG = (Gm - Gm.left_shift(1)).truncate(top + 1).coeffs()       # gamma(m) - gamma(m-1), m <= top
        dG += [0] * (top + 1 - len(dG))
        bad = next((m for m in range(3, top + 1) if dG[m] < 0), None)
        if bad is not None:
            mfail.append((N, gcd(x, N + 1), bad))
        if identities:
            gam[N] = coeffs(Gm, Kn)
    # (II), (III) and, for small x, the identity and D_A >= 0
    yfail, sfail, ifail, dfail = [], [], [], []
    if identities:
        Z1 = inv_prod(range(1, x + 1), Kmax)
        Zm = inv_prod(range(1, x), Kmax)
        cache = {}
        def G(m, n):
            if (m, n) not in cache:
                cache[(m, n)] = coeffs(num_prod([n + t for t in range(1, m + 1)], Kmax).mul_low(Z1 if m == x else Zm, Kmax + 1), Kmax)
            return cache[(m, n)]
    for kind, t, A, N, c, Wp, L, K in bps:
        Y = Y_series(x, c, Wp, L, K)
        if any(v < 0 for v in Y):
            yfail.append((kind, t))
        if any(Y[u] - (Y[u - 1] if u >= 1 else 0) + (Y[u - 2] if u >= 2 else 0) < 0 for u in range(K + 1)):
            sfail.append((kind, t))
        if identities:
            D, r = divmod(A, x)
            E, s = divmod(A, x - 1)
            G0, G1, H0, H1 = G(x, D), G(x, D + 1), G(x - 1, E), G(x - 1, E + 1)
            DA = [(x - 1) * ((x - r) * G0[k] + r * G1[k]) - x * ((x - 1 - s) * H0[k] + s * H1[k]) for k in range(K + 1)]
            g = gam[N]
            dg = fmpz_poly([g[0]] + [g[m] - g[m - 1] for m in range(1, K + 1)])
            rec = coeffs(dg.mul_low(fmpz_poly(Y), K + 1), K)     # (1-q) Gamma Y modulo q^(K+1)
            if rec != DA:
                ifail.append((kind, t))
            if any(v < 0 for v in DA):
                dfail.append((kind, t))
    return len(bps), mfail, yfail, sfail, ifail, dfail


if __name__ == "__main__":
    X = int(sys.argv[1]) if len(sys.argv) > 1 else 60
    XI = int(sys.argv[2]) if len(sys.argv) > 2 else 30
    XS = int(sys.argv[3]) if len(sys.argv) > 3 else 3          # optional first x (to resume a long run)
    t0 = time.time()
    all_ok = True
    for x in range(XS, X + 1):
        nb, mf, yf, sf, idf, df = check_x(x, x <= XI)
        ok = not (mf or yf or sf or idf or df)
        extra = " identity and D_A >= 0 checked;" if x <= XI else ""
        print("x=%3d: %3d breakpoints;%s Y >= 0: %s; small-m: %s; (M) for x <= N <= 2x-2: %s" % (
            x, nb, extra, "ok" if not yf else "FAILS %s" % yf[:3], "ok" if not sf else "FAILS %s" % sf[:3],
            "ok" if not mf else "fails at (N, gcd(x,N+1), degree) %s" % mf[:4]) +
            ("" if not idf else "; IDENTITY FAILS %s" % idf[:3]) + ("" if not df else "; D_A < 0 at %s" % df[:3]), flush=True)
        if yf or sf or idf:
            all_ok = False
        if mf and x >= 13:
            all_ok = False
    print("done (%.0fs): %s" % (time.time() - t0, "identities, (II) and (III) hold; (M) holds for every %d <= x <= %d"
                                  % (max(XS, 13), X) if all_ok else "SOMETHING FAILED"))