#!/usr/bin/env python3
"""Check of Conjecture 8.16 (the fractional first row) in "Analytic q-inequalities and a proof of Bergeron's conjecture", Section 8.8.  Pure Python 3.

G(m, n) = [m+n choose m]_q.  For real y >= 0 and integer m >= 1 the paper interpolates linearly in y:
    H_m(y) = G(m, floor y) + {y} q^{floor y + 1} G(m-1, floor y + 1),
so that m H_m(A/m) = (m - r) G(m, D) + r G(m, D+1) when A = mD + r, 0 <= r < m.

R(x, A), for integers x >= 2 and A >= x^2:   H_x(A/x) >=_q H_{x-1}(A/(x-1))   in degrees k <= A/2.

By Proposition 8.15, R(x, bc) for a < x <= b gives G(b,c) >=_q H_a(bc/a) in degrees <= bc/2 (Conjecture FT),
which implies Troyka's form and Bergeron's conjecture.  The script checks R(x, A) for all integers A in
[x^2, M x^2] and all degrees k <= A/2, in exact integer arithmetic, and reports the smallest relative margin.
In degrees k <= A/x both sides are partition counts, x p_{<=x}(k) and (x-1) p_{<=x-1}(k) after scaling,
so the comparison there is immediate; the script checks every degree anyway.

Usage:  python3 row_check.py [X] [M]        default X = 20, M = 3 (a few seconds); x runs over 2..X.
"""
import sys, time
from functools import lru_cache


@lru_cache(maxsize=None)
def G(m, n):
    """coefficient tuple of [m+n choose m]_q = prod_{j=1}^m (1 - q^{n+j})/(1 - q^j)"""
    if m < 0 or n < 0:
        return (0,)
    L = m * n + 1
    P = [0] * L
    P[0] = 1
    for j in range(1, m + 1):
        e = n + j
        for i in range(L - 1, e - 1, -1):      # times (1 - q^e)
            P[i] -= P[i - e]
        for i in range(j, L):                  # divided by (1 - q^j)
            P[i] += P[i - j]
    return tuple(P)


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


def check(x, Amax):
    """R(x, A) for x^2 <= A <= Amax.  Returns (failures, (smallest relative margin, A, k))."""
    fails, worst = [], None
    for A in range(x * x, Amax + 1):
        D, r = divmod(A, x)
        E, s = divmod(A, x - 1)
        L0, L1, R0, R1 = G(x, D), G(x - 1, D + 1), G(x - 1, E), G(x - 2, E + 1)
        for k in range(A // 2 + 1):
            lhs = x * coef(L0, k) + r * coef(L1, k - D - 1)          # x H_x(A/x)
            rhs = (x - 1) * coef(R0, k) + s * coef(R1, k - E - 1)    # (x-1) H_{x-1}(A/(x-1))
            d = (x - 1) * lhs - x * rhs                                # sign of H_x - H_{x-1}
            if d < 0:
                fails.append((A, k))
                break
            if rhs and k > D:
                rel = d / (x * rhs)
                if worst is None or rel < worst[0]:
                    worst = (rel, A, k)
    return fails, worst


if __name__ == "__main__":
    X = int(sys.argv[1]) if len(sys.argv) > 1 else 20
    M = float(sys.argv[2]) if len(sys.argv) > 2 else 3.0
    t0 = time.time()
    total = 0
    for x in range(2, X + 1):
        Amax = int(M * x * x)
        f, w = check(x, Amax)
        total += len(f)
        print("x=%3d, %d <= A <= %d: %s; smallest relative margin %s" % (
            x, x * x, Amax, "FAILURES %s" % f[:3] if f else "holds",
            "%.3g (A=%d, k=%d)" % w if w else "-"), flush=True)
        G.cache_clear()
    print("done: %d failures, %.0fs" % (total, time.time() - t0))