"""Numerical companion to Section 9 ("The circle method") of the paper "Analytic q-inequalities and a proof of Bergeron's conjecture".
Python 3 + numpy.  Floating point throughout: these are checks and experiments, not proofs.

C_n(q) = prod_{j=2}^n (1-q^{n+j})/(1-q^j) = sum_k c_k q^k,  mu = n(n-1)/2,  rho = e^{-s},
C_n(e^{-w}) = e^{-mu w} Chat(w),  L = log Chat,  h(y) = log(sinh(y)/y),
L^{(r)}(w) = sum_{j=2}^n [((n+j)/2)^r h^{(r)}((n+j)w/2) - (j/2)^r h^{(r)}(jw/2)],
kappa_r = (-1)^r L^{(r)}(s)  (cumulants of the tilted law),  mean = mu - L'(s).

  (A) the two constants used in Step 3 of the proof of Theorem C;
  (B) Lemma 9.5 (Bernoulli pairs) against exact values for [2n choose n]_q and [2n choose n-1]_q, n = 40;
  (C) Lemma 9.3 (derivative bounds) against the actual sup of |L^{(r)}(s - it)|, n = 60;
  (D) Lemma 9.2 (difference formula) against exact coefficients, n = 12;
  (E) Remark 9.8: at n = NE (default 1001), explicit Taylor bounds on [0,T] plus the exact modulus
      |psi| on [T,pi], relative to the Gaussian main term M, at 41 saddle points spanning n+2 <= k <= mu.

Usage:  python3 circle_check.py [NE]
"""
import sys, time
from fractions import Fraction
from math import pi, sqrt, exp, log, sin, cosh, factorial, comb
import numpy as np

# ---------- h and its derivatives ----------
def _bernoulli(m):
    B = [Fraction(1)]
    for k in range(1, m + 1):
        B.append(-sum(comb(k + 1, j) * B[j] for j in range(k)) / Fraction(k + 1))
    return B
_B = _bernoulli(80)
_NS = 40
_a = [float(_B[2 * k] * 2 ** (2 * k) / (2 * k * factorial(2 * k))) for k in range(1, _NS + 1)]  # h = sum a_k y^{2k}

def Li(r, x):
    """Li_{1-r}(x) = sum_l l^{r-1} x^l, r = 2..6 (x real or complex, |x| < 1)."""
    if r == 2: return x / (1 - x) ** 2
    if r == 3: return x * (1 + x) / (1 - x) ** 3
    if r == 4: return x * (1 + 4 * x + x * x) / (1 - x) ** 4
    if r == 5: return x * (1 + 11 * x + 11 * x * x + x ** 3) / (1 - x) ** 5
    if r == 6: return x * (1 + 26 * x + 66 * x * x + 26 * x ** 3 + x ** 4) / (1 - x) ** 6
    raise ValueError

def hder(r, y):
    """h^{(r)}(y) for r = 1..6, y real >= 0 or complex with Re y >= 0 (array)."""
    y = np.asarray(y, dtype=complex)
    out = np.zeros_like(y)
    small = np.abs(y) < 1.0
    ys = y[small]
    acc = np.zeros_like(ys)
    for k in range(1, _NS + 1):
        p = 2 * k
        if p >= r:
            acc += _a[k - 1] * factorial(p) / factorial(p - r) * ys ** (p - r)
    out[small] = acc
    yb = y[~small]
    if r == 1:
        out[~small] = 1 / np.tanh(yb) - 1 / yb
    else:
        out[~small] = (-1) ** r * (factorial(r - 1) / yb ** r - 2 ** r * Li(r, np.exp(-2 * yb)))
    return out

def Lder(n, w, r):
    j = np.arange(2, n + 1, dtype=float)
    m1 = n + j
    return complex(np.sum((m1 / 2) ** r * hder(r, m1 * w / 2) - (j / 2) ** r * hder(r, j * w / 2)))

def cumulants(n, s):
    return dict(k=n * (n - 1) / 2 - Lder(n, s, 1).real, k2=Lder(n, s, 2).real,
                k3=-Lder(n, s, 3).real, k4=Lder(n, s, 4).real)

def saddle(n, k):
    mu = n * (n - 1) / 2
    lo, hi = 0.0, 1.0
    while cumulants(n, hi)['k'] > k: hi *= 2
    for _ in range(100):
        mid = (lo + hi) / 2
        if cumulants(n, mid)['k'] > k: lo = mid
        else: hi = mid
    return (lo + hi) / 2

# ---------- Lemma 9.3: pair bounds ----------
def HB(r, u, v, L=400):
    """sum_{l != 0} (u^2 + (pi|l| - v)^2)^{-r/2}  >=  sup_{|v'| <= v} sum_{l != 0} |u - i v' - i pi l|^{-r}   (0 <= v < pi)."""
    u = np.asarray(u, dtype=float)[..., None]
    v = np.asarray(v, dtype=float)[..., None]
    l = np.arange(1, L + 1, dtype=float)
    return 2 * np.sum((u ** 2 + (pi * l - v) ** 2) ** (-r / 2), axis=-1) + 2 * (pi * L - v[..., 0]) ** (1 - r) / (pi * (r - 1))

def pair_bound(n, sp, T, r, u_zero=False):
    """Upper bound for sup_{|t|<=T} |L^{(r)}(sp - it)|: per pair (j, n+j) the smaller of the h-bound (needs
    (n+j)T/2 < pi) and the l-bound (needs sp > 0).  u_zero: h-bound with Re y = 0, valid for all 0 <= s' <= sp."""
    j = np.arange(2, n + 1, dtype=float)
    m1, m2 = n + j, j
    tot = np.full(j.shape, np.inf)
    ok = m1 * T / 2 < pi
    if np.any(ok):
        u1, u2 = (0 * m1, 0 * m2) if u_zero else (m1 * sp / 2, m2 * sp / 2)
        v1, v2 = np.minimum(m1 * T / 2, 0.999999 * pi), np.minimum(m2 * T / 2, 0.999999 * pi)
        hb = factorial(r - 1) * ((m1 / 2) ** r * HB(r, u1, v1) + (m2 / 2) ** r * HB(r, u2, v2))
        tot = np.where(ok, hb, np.inf)
    if sp > 0 and not u_zero:
        lb = m1 ** r * Li(r, np.exp(-m1 * sp)) + m2 ** r * Li(r, np.exp(-m2 * sp))
        tot = np.minimum(tot, lb)
    return float(np.sum(tot))

def logabs_psi(n, s, th):
    """log |C_n(rho e^{i th})| / C_n(rho), from the product formula (th > 0 array)."""
    th = np.asarray(th, dtype=float)
    j = np.arange(2, n + 1, dtype=float)
    out = np.zeros_like(th)
    for i0 in range(0, len(th), 256):
        t = th[i0:i0 + 256]
        if s > 0:
            for m, sg in ((n + j, 0.5), (j, -0.5)):
                X = np.sin(np.outer(t, m) / 2) ** 2 / np.sinh(m * s / 2)[None, :] ** 2
                out[i0:i0 + 256] += sg * np.sum(np.log1p(X), axis=1)
        else:
            a = np.abs(np.sin(np.outer(t, n + j) / 2)) / (n + j)[None, :]
            b = np.abs(np.sin(np.outer(t, j) / 2)) / j[None, :]
            out[i0:i0 + 256] += np.sum(np.log(a) - np.log(b), axis=1)
    return out

# ---------- exact polynomials (small n) ----------
def catalan_coeffs(n):
    C = [1]
    for m in range(2, n + 1):
        R = C + [0] * m
        for i, v in enumerate(C): R[i + m] += v
        R2 = R + [0] * (2 * m - 1)
        for i, v in enumerate(R): R2[i + 2 * m - 1] -= v
        Q = list(R2)
        for i in range(m + 1, len(Q)): Q[i] += Q[i - m - 1]
        C = Q[:len(Q) - m - 1]
    return C

def gauss_coeffs(M, r):
    P = [1]
    for i in range(1, r + 1):
        m = M - r + i
        R = P + [0] * m
        for a, v in enumerate(P): R[a + m] -= v
        for a in range(i, len(R)): R[a] += R[a - i]
        P = R[:len(R) - i]
    return P

def simpson(y, x):
    h = (x[-1] - x[0]) / (len(x) - 1)
    return h / 3 * (y[0] + y[-1] + 4 * y[1:-1:2].sum() + 2 * y[2:-1:2].sum())

# ---------- (A) ----------
def check_A():
    g = lambda x: (x / 2) ** 2 / np.sinh(x / 2) ** 2
    x = np.linspace(1e-4, 2, 200001)
    gp = np.gradient(g(x), x)
    c1 = np.min(-gp / x)
    x1 = np.linspace(1e-9, 1, 100001)
    c2 = np.min(g(x1) - g(x1 + 1))
    print("(A) g(x) = (x/2)^2/sinh^2(x/2):  min_{0<x<=2} -g'(x)/x = %.4f (>= 1/10);  min_{0<=x<=1} g(x)-g(x+1) = %.4f (>= 0.079)" % (c1, c2))
    assert c1 >= 0.1 and c2 >= 0.079

# ---------- (B) ----------
def check_B(n=40):
    M = 2 * n
    worst = 0.0
    for r in (n, n - 1):
        G = np.array(gauss_coeffs(M, r), dtype=float)
        kk = np.arange(len(G))
        for s in (0.0, 0.01, 0.05, 0.2):
            rho = exp(-s)
            i = np.arange(M, dtype=float)
            z0 = (M - 1) / 2
            p = 1 / (1 + np.exp(s * (i - z0)))          # z = rho^{-(M-1)/2}
            dist = np.array([1.0])
            for pp in p: dist = np.convolve(dist, [1 - pp, pp])
            PS = dist[r]
            for th in (0.03, 0.1, 0.4, 1.0, 2.0, 3.1):
                exact = abs(np.sum(G * rho ** kk * np.exp(1j * th * kk))) / np.sum(G * rho ** kk)
                best = np.inf
                for d in range(1, M):
                    first = np.array([a for a in range(M) if a % (2 * d) < d and a + d <= M - 1], dtype=int)
                    if len(first) == 0: continue
                    pi_i = p[first] * (1 - p[first + d]) + (1 - p[first]) * p[first + d]
                    R = sqrt(max(0.0, 1 - sin(th * d / 2) ** 2 / cosh(d * s / 2) ** 2))
                    best = min(best, float(np.prod(1 - pi_i * (1 - R))) / PS)
                worst = max(worst, exact / best)
    print("(B) Lemma 9.5, n=%d, r in {n, n-1}, 4 values of s, 6 of theta: max exact/bound = %.3f (<= 1)" % (n, worst))
    assert worst <= 1 + 1e-9

# ---------- (C) ----------
def check_C(n=60):
    worst = 0.0
    for s in (0.0, 0.002, 0.01, 0.05, 0.2):
        N = n if s <= 1 / n else 1 / s
        T = pi / (4 * N)
        for r in (3, 4, 5, 6):
            act = max(abs(Lder(n, s - 1j * t, r)) for t in np.linspace(-T, T, 201))
            worst = max(worst, act / pair_bound(n, s, T, r))
    # parity (Lemma 9.4): |Re L5(s-it)| <= s sup_{s'<=s} |L6(s'-it)|
    s = 0.004; T = pi / (4 * n); par = 0.0
    for t in np.linspace(-T, T, 101):
        m6 = max(abs(Lder(n, sp - 1j * t, 6)) for sp in np.linspace(0, s, 21))
        par = max(par, abs(Lder(n, s - 1j * t, 5).real) / (s * m6))
    print("(C) Lemma 9.3, n=%d: max over tests of sup|L^(r)| / bound = %.3f (<= 1);  Lemma 9.4: max ratio = %.3f (<= 1)" % (n, worst, par))
    assert worst <= 1 + 1e-9 and par <= 1 + 1e-6

# ---------- (D) ----------
def check_D(n=12):
    C = np.array(catalan_coeffs(n), dtype=float)
    kk = np.arange(len(C)); err = 0.0
    for k in (5, 20, 40, 66):
        for s in (0.0, 0.05):
            rho = exp(-s); th = np.linspace(0, pi, 20001)
            Cq = np.array([np.sum(C * (rho * np.exp(1j * t)) ** kk) for t in th])
            val = simpson(((Cq * np.exp(-1j * k * th) / np.sum(C * rho ** kk)) * (1 - rho * np.exp(1j * th))).real, th) / pi
            err = max(err, abs(val * np.sum(C * rho ** kk) * rho ** (-k) - (C[k] - C[k - 1])))
    print("(D) Lemma 9.2, n=%d: max |formula - (c_k - c_{k-1})| = %.1e" % (n, err))

# ---------- (E) ----------
def local_bounds(n, s, cu, lam, K5, K6, M6):
    k2, k3, k4 = cu['k2'], cu['k3'], cu['k4']
    rho = exp(-s)
    u = np.linspace(0, lam, 4001); th = u / sqrt(k2)
    r4 = abs(k4) / k2 ** 2 * u ** 4 / 24 + K5 / k2 ** 2.5 * u ** 5 / 120          # 4th order, 5th-derivative remainder
    r6lo = k4 / k2 ** 2 * u ** 4 / 24 - K6 / k2 ** 3 * u ** 6 / 720               # 6th order
    r6hi = k4 / k2 ** 2 * u ** 4 / 24 + K6 / k2 ** 3 * u ** 6 / 720
    lo = -u ** 2 / 2 + np.maximum(-r4, r6lo)                                        # Re f >= lo
    hi = -u ** 2 / 2 + np.minimum(r4, r6hi)                                         # Re f <= hi
    ph = abs(k3) * u ** 3 / (6 * k2 ** 1.5) + min(K5, s * M6) * u ** 5 / (120 * k2 ** 2.5)   # |Im f| <= ph
    w = 1 - rho * np.cos(th)
    pos = np.exp(lo) * w
    neg = np.exp(hi) * (w * np.minimum(ph ** 2 / 2, 2.0) + rho * np.abs(np.sin(th)) * np.minimum(ph, 1.0))
    c = 1 / (pi * sqrt(k2))
    return c * simpson(pos, u), c * simpson(neg, u)

def experiment_E(n):
    t0 = time.time()
    s0 = saddle(n, n + 2)
    sigs = [0.0] + list(np.geomspace(1e-3, s0 * n, 40))
    worst = None
    for sig in sigs:
        s = sig / n; rho = exp(-s)
        cu = cumulants(n, s); k2 = cu['k2']
        Mg = ((1 - rho) + rho / (2 * k2)) / sqrt(2 * pi * k2)
        best = None
        for lam in (3.0, 4.0, 5.0, 6.0, 8.0, 10.0):
            T = lam / sqrt(k2)
            K5 = pair_bound(n, s, T, 5); K6 = pair_bound(n, s, T, 6)
            M6 = pair_bound(n, s, T, 6, u_zero=True) if n * T < pi else np.inf
            pos, neg = local_bounds(n, s, cu, lam, K5, K6, M6)
            th = np.concatenate([np.geomspace(T, min(pi, 200 * T), 4000), np.linspace(min(pi, 200 * T), pi, 4000)[1:]])
            vals = np.exp(logabs_psi(n, s, th)) * np.sqrt((1 - rho) ** 2 + 4 * rho * np.sin(th / 2) ** 2)
            G = np.sum(np.maximum(vals[:-1], vals[1:]) * np.diff(th)) / pi
            m = (pos - neg - G) / Mg
            if best is None or m > best[0]: best = (m, lam)
        if worst is None or best[0] < worst[0]: worst = (best[0], cu['k'], best[1])
    print("(E) n=%d: min over 41 saddle points of (local lower bound - |rest|)/M = %.4f at k = %.0f (lambda = %g); %.0fs"
          % (n, worst[0], worst[1], worst[2], time.time() - t0))

if __name__ == "__main__":
    NE = int(sys.argv[1]) if len(sys.argv) > 1 else 1001
    check_A(); check_B(); check_C(); check_D(); experiment_E(NE)