"""Corner check for the paper "Analytic q-inequalities and a proof of Bergeron's conjecture", Section 8.8.  Pure Python 3, exact integers.

C_n(q) = [2n choose n]/[n+1] is MacMahon's q-Catalan polynomial, built from
    C_n = C_{n-1} (1 + q^n)(1 - q^{2n-1}) / (1 - q^{n+1}).
Write c_k = [q^k] C_n and y_k = c_k + c_{k-1}; the corner is (b, c) = (a+1, a+1) with n = a+1.
For n <= N the script checks:
  (A) c_k >= c_{k-2} and c_k >= c_{k-3} for all k <= (n^2-2n)/2.  With Proposition 8.12 this makes
      every term of the corner identity nonnegative, so FStep holds at the corner for a = n-1
      (Corollary 8.13; Stanley's theorem is not needed for these n).
  (B) Stanley's theorem: (1+q) C_n / (1+q^n) is a palindromic unimodal polynomial.
  (C) Where C_n fails to be nondecreasing: on degrees 2..n(n-1)/2, and on the corner range
      2..(n^2-2n)/2 (with the first failing degree).
  (F) (1+q) C_n is unimodal (parity unimodality, Corollary 9.14; for n >= 16 it follows from Theorem C').
For 2 <= a <= A it also checks, against independently computed Gaussian coefficients:
  (D) Proposition 8.12:  a c_k - [q^{k-2}] [2n choose n-2]
                           = sum_{i=2}^{a+1} (c_k - c_{k-i}) + sum_{j>=1} (y_{k-j(a+3)} - y_{k-j(a+3)-a}).
  (E) Lemma 8.11: with FStep's two sides computed from their definitions,
      a (RHS_k - LHS_k) = [q^{k-a-1}] (a C_n - q^2 [2n choose n-2])  for 0 <= k <= (a+1)^2/2.

Usage:  python3 corner_check.py [N] [A]      defaults N = 200, A = 40 (about a second).
        N = 1000 takes a few minutes.
"""
import sys, time


def mul_1pqm(P, m):                  # P * (1 + q^m)
    R = P + [0] * m
    for i, v in enumerate(P):
        R[i + m] += v
    return R


def mul_1mqm(P, m):                  # P * (1 - q^m)
    R = P + [0] * m
    for i, v in enumerate(P):
        R[i + m] -= v
    return R


def div_1mqm(P, m):                  # P / (1 - q^m), which must be a polynomial
    Q = list(P)
    for i in range(m, len(Q)):
        Q[i] += Q[i - m]
    assert all(v == 0 for v in Q[len(Q) - m:]), "not divisible"
    return Q[:len(Q) - m]


def div_1pqm(P, m):                  # P / (1 + q^m), which must be a polynomial
    Q = list(P)
    for i in range(m, len(Q)):
        Q[i] -= Q[i - m]
    assert all(v == 0 for v in Q[len(Q) - m:]), "not divisible"
    return Q[:len(Q) - m]


def gauss(N, k):                     # [N choose k]_q
    P = [1]
    for i in range(1, k + 1):
        P = div_1mqm(mul_1mqm(P, N - k + i), i)
    return P


def coef(P, k):
    return P[k] if 0 <= k < len(P) else 0


def palindromic_unimodal(P):
    if P != P[::-1]:
        return False
    i = 0
    while i + 1 < len(P) and P[i + 1] >= P[i]:
        i += 1
    while i + 1 < len(P) and P[i + 1] <= P[i]:
        i += 1
    return i == len(P) - 1


def check_corner(a, C):
    """(D) and (E) for n = a+1; C = C_n."""
    n = a + 1
    U = gauss(2 * n, n - 2)                                   # [2n choose n-2] = G(a-1, a+3)
    V = lambda m: a * coef(C, m) - coef(U, m - 2)             # coefficients of a C_n - q^2 U
    y = lambda m: coef(C, m) + coef(C, m - 1)
    for k in range((a * a - 1) // 2 + 1):                     # (D)
        rhs = sum(coef(C, k) - coef(C, k - i) for i in range(2, a + 2))
        j = 1
        while k - j * (a + 3) >= 0:
            rhs += y(k - j * (a + 3)) - y(k - j * (a + 3) - a)
            j += 1
        assert V(k) == rhs, ("Proposition 8.12 fails", a, k)
        assert V(k) >= 0, ("corner inequality fails", a, k)
    # (E)  FStep at (b, c) = (a+1, a+1):  [q^k] H_a(x) <= Lam(k) + Lam(k - b),  x = (a+1)^2/a.
    G1 = gauss(2 * a + 1, a)                                  # H_a((b-1)c/a) = H_a(b(c-1)/a) = G(a, a+1)
    G2 = gauss(2 * a + 2, a)                                  # G(a, a+2),  floor(x) = a+2,  {x} = 1/a
    top = a * (a + 1)                                         # degree of G(a, a+1)

    def lam(k):                                               # largest coefficient of G1 in degrees <= min(k, top-k)
        if k < 0 or k > top:
            return 0
        return max(G1[:min(k, top - k) + 1])

    for k in range((a + 1) ** 2 // 2 + 1):
        a_lhs = a * coef(G2, k) + coef(U, k - a - 3)          # a * [q^k] H_a(x)
        a_rhs = a * (lam(k) + lam(k - a - 1))
        assert a_rhs - a_lhs == V(k - a - 1), ("Lemma 8.11 fails", a, k)


def main():
    N = int(sys.argv[1]) if len(sys.argv) > 1 else 200
    A = int(sys.argv[2]) if len(sys.argv) > 2 else 40
    t0 = time.time()
    C = [1]                                                   # C_1
    step_fail, stanley_fail, mono_half, mono_corner, parity_fail = [], [], [], [], []
    for n in range(2, max(N, A + 1) + 1):
        C = div_1mqm(mul_1mqm(mul_1pqm(C, n), 2 * n - 1), n + 1)
        a = n - 1
        if 2 <= a <= A:
            check_corner(a, C)
        if n > N:
            continue
        top = (n * n - 2 * n) // 2
        bad = [k for k in range(2, top + 1) if C[k] < C[k - 2] or (k >= 3 and C[k] < C[k - 3])]
        if bad:
            step_fail.append((n, bad[:5]))
        if not palindromic_unimodal(div_1pqm(mul_1pqm(C, 1), n)):
            stanley_fail.append(n)
        if not palindromic_unimodal(mul_1pqm(C, 1)):
            parity_fail.append(n)
        dips = [k for k in range(2, n * (n - 1) // 2 + 1) if C[k] < C[k - 1]]
        if dips:
            mono_half.append(n)
            if dips[0] <= top:
                mono_corner.append((n - 1, dips[0]))
    print("(A) c_k >= c_{k-2} and c_k >= c_{k-3} for k <= (n^2-2n)/2:",
          "FAILS for %s" % step_fail if step_fail else "true for all 2 <= n <= %d" % N)
    if not step_fail:
        print("    => FStep holds at the corner (a+1, a+1) for every 2 <= a <= %d" % (N - 1))
    print("(B) Stanley's theorem:", "FAILS for %s" % stanley_fail if stanley_fail
          else "(1+q)C_n/(1+q^n) palindromic and unimodal for all 2 <= n <= %d" % N)
    print("(C) C_n not nondecreasing on [2, n(n-1)/2] exactly for n in", mono_half)
    print("    on the corner range [2, (n^2-2n)/2]: (a, first failing degree) =", mono_corner)
    print("(F) (1+q)C_n unimodal:", "FAILS for %s" % parity_fail if parity_fail
          else "true for all 2 <= n <= %d" % N)
    print("(D) Proposition 8.12 identity and corner inequality: verified for 2 <= a <= %d" % A)
    print("(E) Lemma 8.11 (FStep at the corner <=> q-Catalan form): verified for 2 <= a <= %d" % A)
    print("time %.1fs" % (time.time() - t0))


if __name__ == "__main__":
    main()