#!/usr/bin/env python3 """ Certificate for Theorems A and A' of the paper "Analytic q-inequalities and a proof of Bergeron's conjecture": Theorem A'. For 1 <= a <= 20 and a <= b <= c, write bc = aD + r with 0 <= r < a. Then a*[b+c choose b]_q >= (a-r)*[a+D choose a]_q + r*[a+D+1 choose a]_q in degrees 0..floor(bc/2). It follows from the fractional Pascal step FStep (Conjecture 8.10 of the paper), which this script proves for 3 <= a <= 20 and all a < b <= c. Theorem A. (a consequence) For 1 <= a <= 20, a <= b <= c and a*d <= b*c: [b+c choose b]_q >= [a+d choose a]_q (coefficientwise). In particular Bergeron's conjecture holds whenever the smallest side is <= 20. Notation (as in the paper): G(x,y) = [x+y choose x]_q with coefficients g_{x,y}(k); for a < b <= c D = floor(bc/a), D' = floor((b-1)c/a), D'' = floor(b(c-1)/a), k_max = floor(aD/2); T(b,c) is the statement G(b,c) >= G(a,D); FT(b,c) is the statement of Theorem A' for (b,c). FT(b,c) implies T(b,c) (degrees <= bc/2, then palindromy and nesting), and FT(a,c) is trivial. Proof structure (Sections 8.4 and 8.7 of the paper) * FT(b-1,c) and FT(b,c-1) [or FT(b-1,b) if b=c] give lower bounds for G(b-1,c), G(b,c-1) in every degree (palindromy and unimodality), and FT(b,c) follows by Pascal's rule if, for k <= bc/2, L_FT(k) := g_{a,D}(k) - g_{a,D'}(k) + (r/a) g_{a-1,D+1}(k-D-1) <= R(k) := g_{a,D''}(k-b), or if k lies in the range of the fractional matching identity (k <= min(2D'+2-b, (b-1)c/2)). * Proposition 8.4: for a 'cell' kl*bc/a <= k <= kh*bc/a, L_FT(k) <= R(k) follows from P_cell(b,c) >= 0, where P_cell is an explicit polynomial with rational coefficients (built below from Irwin-Hall densities; the bounds used are Lemmas 5.4-5.7 of the paper: smoothing identity, unimodality, and m! e_m <= p_1^m <= m! h_m). What this script does, for every a in [3, amax] (Proposition 8.5): 1. proves P_cell >= 0 for every cell on the regions {b fixed, c >= C_a(b)} (a < b < B1) and {b >= B1, c >= b}, in exact rational arithmetic: after the shift c = C_a(b) + y, resp. b = B1 + x, c = B1 + x + y, either all coefficients are >= 0, or (Sturm sequences, via sympy) the polynomial in y has no root in [0, oo), resp. all its coefficient polynomials in x (or all those in y) are >= 0 on [0, oo); The cells cover kappa in [min(2 - 2/b0 - a/c0, (a/2)(1 - 1/b0)), a/2], i.e. every k that the matching identity does not handle; 2. verifies FT(b,c) exactly on the finite complement F_a = {(b,c): a < b < B1, b <= c < C_a(b)}: a*G(b,c) >= (a-r)*G(a,D) + r*G(a,D+1) in degrees 0..floor(bc/2), by big-integer (Kronecker) arithmetic; 3. verifies FStep exactly on F_a (so, with 1, FStep holds for all a < b <= c): a*[q^k]H(bc) <= Lam_{b-1,c}(k) + Lam_{b,c-1}(k-b) for 0 <= k <= bc/2, with a*[q^k]H(n) = (a-r) g_{a,D}(k) + r g_{a,D+1}(k) for n = aD + r, and Lam_{m}(k) the maximum of a*[q^j]H(n') over j <= min(k, n'-k) (n' the product of the neighbour; Lam = 0 outside [0, n']). a = 1, 2 are proved by hand in the paper (a = 2 is spot-checked here; for a = 2 the extra term of FT lives above degree bc/2, so FT is T). Options: --amax A (default 20) --brute N independent exact check of FT for all a < b <= c, bc <= N --selftest compare every bound with exact coefficients on random small cases Requires: Python 3.8+, sympy. """ import argparse, random, sys, time from fractions import Fraction as Fr from functools import lru_cache from math import comb, factorial, floor import sympy as sp # (J, Jbi, B1, [C_a(a+1), C_a(a+2), ..., C_a(B1-1)]) -- thresholds found by bisection, proved below THRESHOLDS = { 3: (8, 4, 8, [21, 13, 10, 8]), 4: (8, 4, 11, [37, 25, 19, 16, 13, 12]), 5: (8, 4, 15, [59, 41, 32, 26, 22, 19, 17, 16, 15]), 6: (8, 4, 19, [85, 62, 49, 40, 34, 30, 27, 24, 22, 21, 20, 18]), 7: (8, 4, 23, [116, 88, 70, 58, 50, 44, 39, 35, 32, 29, 28, 26, 25, 24, 23]), 8: (8, 4, 27, [152, 118, 96, 81, 70, 61, 54, 49, 45, 41, 38, 36, 34, 32, 30, 29, 28, 27]), 9: (8, 4, 31, [192, 153, 127, 107, 93, 82, 73, 66, 60, 55, 51, 48, 45, 42, 40, 38, 37, 35, 34, 33, 32]), 10: (8, 4, 36, [236, 193, 161, 138, 120, 107, 96, 86, 79, 73, 67, 62, 58, 55, 52, 50, 47, 46, 44, 42, 41, 40, 38, 37, 36]), 11: (8, 4, 41, [286, 237, 201, 173, 152, 135, 122, 110, 101, 93, 86, 80, 75, 70, 66, 63, 60, 58, 55, 53, 51, 50, 48, 47, 45, 44, 43, 42, 41]), 12: (8, 4, 46, [339, 285, 244, 213, 188, 168, 151, 138, 126, 116, 108, 100, 94, 88, 83, 79, 75, 72, 69, 66, 64, 62, 60, 58, 56, 54, 53, 51, 50, 49, 48, 47, 46]), 13: (8, 4, 51, [398, 338, 293, 257, 228, 204, 185, 169, 155, 143, 133, 124, 116, 109, 103, 97, 92, 88, 84, 81, 78, 75, 73, 70, 68, 66, 64, 63, 61, 59, 58, 57, 55, 54, 53, 52, 51]), 14: (8, 4, 56, [461, 396, 346, 305, 272, 245, 223, 204, 188, 174, 162, 151, 141, 133, 125, 119, 113, 107, 102, 98, 94, 91, 88, 85, 82, 80, 78, 75, 73, 71, 70, 68, 66, 65, 64, 62, 61, 60, 59, 58, 57]), 15: (8, 8, 61, [528, 458, 403, 358, 321, 290, 265, 243, 224, 208, 193, 181, 170, 160, 151, 143, 135, 129, 123, 117, 113, 109, 105, 101, 98, 95, 92, 90, 87, 85, 83, 81, 79, 77, 76, 74, 72, 71, 70, 68, 67, 66, 65, 63, 62]), 16: (8, 8, 67, [600, 525, 465, 415, 374, 340, 311, 286, 264, 245, 229, 214, 201, 190, 179, 170, 161, 153, 146, 140, 134, 129, 124, 120, 116, 113, 109, 106, 103, 100, 98, 95, 93, 91, 89, 87, 85, 83, 82, 80, 79, 77, 76, 74, 73, 72, 71, 70, 69, 68]), 17: (8, 8, 73, [676, 596, 531, 477, 432, 394, 361, 332, 308, 287, 268, 251, 236, 223, 211, 200, 190, 181, 172, 165, 158, 151, 145, 141, 136, 132, 128, 124, 121, 118, 115, 112, 109, 106, 104, 102, 99, 97, 95, 93, 92, 90, 88, 87, 85, 84, 82, 81, 80, 78, 77, 76, 75, 74, 73]), 18: (8, 8, 79, [757, 672, 602, 543, 494, 451, 415, 383, 356, 332, 311, 292, 275, 259, 246, 233, 222, 211, 201, 193, 185, 177, 170, 164, 158, 153, 149, 144, 140, 137, 133, 130, 126, 123, 121, 118, 115, 113, 110, 108, 106, 104, 102, 100, 98, 97, 95, 94, 92, 91, 89, 88, 87, 85, 84, 83, 82, 81, 80, 79]), 19: (8, 8, 85, [843, 752, 677, 614, 560, 514, 474, 439, 408, 381, 357, 336, 317, 299, 284, 269, 256, 244, 233, 223, 214, 205, 197, 190, 183, 177, 172, 167, 162, 158, 153, 149, 146, 142, 139, 136, 133, 130, 127, 125, 122, 120, 117, 115, 113, 111, 109, 107, 106, 104, 102, 101, 99, 98, 96, 95, 94, 92, 91, 90, 89, 88, 87, 86, 85]), 20: (8, 8, 91, [933, 837, 757, 689, 631, 580, 536, 498, 464, 434, 408, 384, 362, 343, 325, 309, 294, 281, 268, 257, 246, 236, 227, 219, 211, 203, 197, 191, 186, 181, 176, 171, 167, 163, 159, 155, 152, 149, 146, 143, 140, 137, 134, 132, 129, 127, 125, 123, 121, 119, 117, 115, 113, 112, 110, 108, 107, 105, 104, 103, 101, 100, 99, 97, 96, 95, 94, 93, 92, 91]), } # ============================================================================ exact verification def _offset(h, B): return (1 << (B - 1)) * (((1 << ((h + 1) * B)) - 1) // ((1 << B) - 1)) def geq_upto(U, V, h, B): """U, V encode polynomials (base 2^B digits, coefficients < 2^(B-1)); True iff U_k >= V_k, k <= h.""" mask = (1 << ((h + 1) * B)) - 1 W = (U & mask) - (V & mask) + _offset(h, B) # digit k = U_k - V_k + 2^(B-1) in [0, 2^B) if W < 0: return False nbytes = (h + 1) * B // 8 bs = W.to_bytes(nbytes + 1, 'little') step = B // 8 tops = bs[step - 1:nbytes:step] # top byte of every digit return len(tops) == h + 1 and min(tops) >= 128 def verify_pairs(a, F): """Exactly verify FT(b,c): a*G(b,c) >= (a-r)*G(a,D) + r*G(a,D+1) in degrees 0..floor(bc/2) (bc = aD + r) for (b,c) in F. Returns the list of failures.""" F = sorted(set(F)) if not F: return [] bmax = max(b for b, c in F); cmax = max(c for b, c in F) Dmax = max((b * c) // a for b, c in F) + 1 bits = max(max((a * comb(b + c, b)).bit_length() for b, c in F), (a * comb(a + Dmax, a)).bit_length()) B = ((bits + 2 + 7) // 8) * 8 col = [1] * (a + 1); T = [1] for D in range(1, Dmax + 1): # G(i,D) = G(i,D-1) + q^D G(i-1,D) new = [1] for i in range(1, a + 1): new.append(col[i] + (new[i - 1] << (D * B))) col = new; T.append(col[a]) need = {} for b, c in F: need.setdefault(c, []).append(b) bneed = [0] * (cmax + 2) for b, c in F: bneed[c] = max(bneed[c], b) for cc in range(cmax - 1, -1, -1): bneed[cc] = max(bneed[cc], bneed[cc + 1]) L = max((b * c) // 2 + 1 for b, c in F) gmask = (1 << (L * B)) - 1 prev = [1] * (bmax + 1) fails = [] for c in range(1, cmax + 1): cur = [1] for b in range(1, bneed[c] + 1): # G(b,c) = G(b,c-1) + q^c G(b-1,c) cur.append((prev[b] + (cur[b - 1] << (c * B))) & gmask) for b in need.get(c, []): D, r = divmod(b * c, a) if not geq_upto(a * cur[b], (a - r) * T[D] + r * T[D + 1], (b * c) // 2, B): fails.append((b, c)) prev = cur return fails def _digits(X, upto, B): """The first upto+1 base-2^B digits of the nonnegative integer X.""" nb = B // 8 bs = X.to_bytes((X.bit_length() + 7) // 8 + nb, 'little') return [int.from_bytes(bs[k * nb:(k + 1) * nb], 'little') for k in range(upto + 1)] def verify_fstep(a, F): """Exactly verify FStep(b,c) for (b,c) in F. Returns the list of failures (b, c, first failing k).""" F = sorted(set(F)) if not F: return [] Dmax = max((b * c) // a for b, c in F) + 1 B = (((a * comb(a + Dmax, a)).bit_length() + 8 + 7) // 8) * 8 col = [1] * (a + 1); T = [1] for D in range(1, Dmax + 1): new = [1] for i in range(1, a + 1): new.append(col[i] + (new[i - 1] << (D * B))) col = new; T.append(col[a]) def aH(n, upto): # a * coefficients of H_a(n/a), degrees 0..upto D, r = divmod(n, a) return [(a - r) * x + r * y for x, y in zip(_digits(T[D], upto, B), _digits(T[D + 1], upto, B))] def lam(bb, cc): # prefix maxima over the lower half of (bb, cc) n = bb * cc; pref = []; m = 0 for x in aH(n, n // 2): m = max(m, x); pref.append(m) return pref, n def lv(L, k): pref, n = L return 0 if k < 0 or k > n else pref[min(k, n - k)] fails = [] for b, c in F: n = b * c H = aH(n, n // 2) L1 = lam(b - 1, c); L2 = lam(b, c - 1) if c - 1 >= b else lam(b - 1, b) for k in range(n // 2 + 1): if lv(L1, k) + lv(L2, k - b) < H[k]: fails.append((b, c, k)); break return fails # ============================================================================ Irwin-Hall density @lru_cache(maxsize=None) def f_IH(m, y): """Density of U_1+...+U_m (U_i uniform on [0,1)) at the rational point y.""" y = Fr(y) if y <= 0 or y >= m: return Fr(0) return sum((-1) ** i * comb(m, i) * (y - i) ** (m - 1) for i in range(floor(y) + 1)) / factorial(m - 1) def phi(m): return f_IH(m, Fr(m, 2)) # ============================================================================ exact bivariate polynomials def padd(P, Q, s=1): R = dict(P) for k, v in Q.items(): R[k] = R.get(k, 0) + s * v return {k: v for k, v in R.items() if v != 0} def pscale(P, s): return {k: s * v for k, v in P.items()} if s != 0 else {} def pmul(P, Q): R = {} for (i1, j1), v1 in P.items(): for (i2, j2), v2 in Q.items(): R[(i1 + i2, j1 + j2)] = R.get((i1 + i2, j1 + j2), 0) + v1 * v2 return {k: v for k, v in R.items() if v != 0} def ppow(P, n): R = {(0, 0): Fr(1)} for _ in range(n): R = pmul(R, P) return R def const(v): return {(0, 0): Fr(v)} if v != 0 else {} class Blocks: """Polynomials (in the shifted variables x, y) entering P_cell: 'uni': b = b0, c = c0 + y; 'bi': b = b0 + x, c = b0 + x + y. Dhi = bc/a, D1lo = ((b-1)c-(a-1))/a, D2lo = (b(c-1)-(a-1))/a, tau_j = D1lo + (j/J)(Dhi - D1lo), R = (D2lo+1)^(a-1), dH_j = H(tau_{j+1}) - H(tau_j) with H(x) = (x+a)^(a-1)/(a-1), Q = ((Dhi+a)^(a-2) - (D1lo+a)^(a-2))/2.""" def __init__(self, a, J, mode, b0, c0): self.a, self.mode = a, mode if mode == 'uni': b = const(b0); c = padd(const(c0), {(0, 1): Fr(1)}) else: b = padd(const(b0), {(1, 0): Fr(1)}); c = padd(const(b0), {(1, 0): Fr(1), (0, 1): Fr(1)}) bc = pmul(b, c) Dhi = pscale(bc, Fr(1, a)) D1lo = pscale(padd(padd(bc, c, -1), const(a - 1), -1), Fr(1, a)) D2lo = pscale(padd(padd(bc, b, -1), const(a - 1), -1), Fr(1, a)) diff = padd(Dhi, D1lo, -1) taus = [padd(D1lo, pscale(diff, Fr(j, J))) for j in range(J + 1)] Hs = [pscale(ppow(padd(t, const(a)), a - 1), Fr(1, a - 1)) for t in taus] self.R = ppow(padd(D2lo, const(1)), a - 1) self.Q = pscale(padd(ppow(padd(Dhi, const(a)), a - 2), ppow(padd(D1lo, const(a)), a - 2), -1), Fr(1, 2)) self.dH = [padd(Hs[j + 1], Hs[j], -1) for j in range(J)] def cell_poly(self, Phis, fb, fZ): """P_cell = fZ * R - a * sum_j Phi_j dH_j - a * fb * Q.""" P = pscale(self.R, fZ) for ph, dh in zip(Phis, self.dH): if ph: P = padd(P, dh, -self.a * ph) if fb: P = padd(P, self.Q, -self.a) return P def cell_consts(a, J, b0, c0, kl, kh): """Constants of Proposition 8.5 for the cell [kl, kh] and the region with corner (b0, c0).""" m = a - 1; C2 = comb(a - 1, 2) Phis = []; fb = False for j in range(J): beta = Fr(b0 - 1) + Fr(j, J) E1 = Fr((a - 1) * (J - j), J) + a * (C2 + a - 2) Ybar = kh * b0 / beta - 1 + E1 / (beta * c0) if Ybar <= Fr(m, 2) * (1 - Fr(a) / (beta * c0)): Phis.append(f_IH(m, max(Ybar, Fr(0)))) # smoothing upper bound, Lemma 5.5(b) else: Phis.append(phi(m)); fb = True # maximum bound, Lemma 5.5(c) for j in range(J - 2, -1, -1): # nonincreasing in j (automatic; kept as a guard) Phis[j] = max(Phis[j], Phis[j + 1]) return Phis, fb, f_IH(a, max(Fr(0), kl - Fr(a, c0))) def _univ_nonneg(coeffs): """coeffs: {degree: Fraction}. Sufficient test for 'polynomial >= 0 on [0, oo)'.""" coeffs = {d: v for d, v in coeffs.items() if v != 0} if not coeffs or all(v >= 0 for v in coeffs.values()): return True d0 = min(coeffs) coeffs = {d - d0: v for d, v in coeffs.items()} deg = max(coeffs) if coeffs[0] < 0 or coeffs[deg] < 0: return False t = sp.Symbol('t') p = sp.Poly([sp.Rational(coeffs.get(d, Fr(0)).numerator, coeffs.get(d, Fr(0)).denominator) for d in range(deg, -1, -1)], t, domain='QQ') return p.count_roots(0, None) == 0 # p(0) > 0 and no real root in [0, oo) def poly_nonneg(P): """Sufficient test for P >= 0 on the closed quadrant x, y >= 0.""" if P.get((0, 0), 0) <= 0: return False if all(v >= 0 for v in P.values()): return True if all(i == 0 for (i, j) in P): # univariate in y return _univ_nonneg({j: v for (i, j), v in P.items()}) byy, byx = {}, {} for (i, j), v in P.items(): byy.setdefault(j, {})[i] = v byx.setdefault(i, {})[j] = v # P = sum_j p_j(x) y^j = sum_i r_i(y) x^i ; nonnegative coefficient polynomials suffice return all(_univ_nonneg(d) for d in byy.values()) or all(_univ_nonneg(d) for d in byx.values()) def prove_region(a, J, mode, b0, c0, n0=16, maxdepth=12): """Prove P_cell >= 0 on the region for cells covering [kmin, a/2] (bisecting as needed), kmin = min(2 - 2/b0 - a/c0, (a/2)(1 - 1/b0)). Returns the number of cells used, or raises.""" blk = Blocks(a, J, mode, b0, c0) kmin = min(2 - Fr(2, b0) - Fr(a, c0), Fr(a, 2) * (1 - Fr(1, b0))); kmax = Fr(a, 2) stack = [(kmin + (kmax - kmin) * i / n0, kmin + (kmax - kmin) * (i + 1) / n0, 0) for i in range(n0)] ncells = 0 while stack: kl, kh, d = stack.pop() if poly_nonneg(blk.cell_poly(*cell_consts(a, J, b0, c0, kl, kh))): ncells += 1 continue if d >= maxdepth: raise RuntimeError("cell inequality not proved: a=%d mode=%s b0=%d c0=%d cell=[%s,%s]" % (a, mode, b0, c0, kl, kh)) mid = (kl + kh) / 2 stack += [(kl, mid, d + 1), (mid, kh, d + 1)] return ncells # ============================================================================ driver def certify(a): t0 = time.time() J, Jbi, B1, Cs = THRESHOLDS[a] assert len(Cs) == B1 - a - 1 and all(C >= b for b, C in zip(range(a + 1, B1), Cs)) cells = 0 for b, C in zip(range(a + 1, B1), Cs): cells += prove_region(a, J, 'uni', b, C) cells += prove_region(a, Jbi, 'bi', B1, B1) F = [(b, c) for b, C in zip(range(a + 1, B1), Cs) for c in range(b, C)] fails = verify_pairs(a, F) fails2 = verify_fstep(a, F) print("a=%2d: step proved on {c >= C_a(b)} for %d <= b < %d and on {b >= %d, c >= b} (%d cells); " "C_a(%d)=%d; FT on |F_a|=%d pairs (max area %d) verified exactly: %s; FStep on F_a: %s [%.1fs]" % (a, a + 1, B1, B1, cells, a + 1, Cs[0], len(F), max(b * c for b, c in F), "OK" if not fails else fails[:5], "OK" if not fails2 else fails2[:5], time.time() - t0), flush=True) return not fails and not fails2 def check_a2(): for b in range(3, 400): for c in range(b, 400): assert (b * c) // 2 <= 2 * (((b - 1) * c) // 2) + 2 - b, (b, c) assert (b * c) // 2 + 1 >= b * c - (b * c) // 2 # the extra term of FT starts above bc/2 print("a= 2: window [2D'+3-b, k_max] empty for all 3 <= b <= c (parity proof in the paper; " "spot-checked for c < 400); FT = T for a = 2") def brute(N): t0 = time.time(); tot = 0; bad = [] for a in range(1, N): F = [(b, c) for b in range(a + 1, N) for c in range(b, N // b + 1)] if not F: break tot += len(F) bad += [(a,) + f for f in verify_pairs(a, F)] print("brute force (FT, hence T): %d triples (a,b,c) with a 2500: continue D, r = divmod(b * c, a); D1 = ((b - 1) * c) // a; D2 = (b * (c - 1)) // a; khalf = (b * c) // 2 kmin = min(2 - Fr(2, b0) - Fr(a, c0), Fr(a, 2) * (1 - Fr(1, b0))); kmx = Fr(a, 2); Dhi = Fr(b * c, a) lo = min(2 * D1 + 3 - b, ((b - 1) * c) // 2 + 1) if khalf < lo: continue tests += 1 blk = Blocks(a, J, mode, b0, c0) xv, yv = (0, c - c0) if mode == 'uni' else (b - b0, c - b) ev = lambda P: sum(v * Fr(xv) ** i * Fr(yv) ** j for (i, j), v in P.items()) Rv, Qv, dHv = ev(blk.R), ev(blk.Q), [ev(p) for p in blk.dH] assert kmin * Dhi <= lo and khalf <= kmx * Dhi w = (kmx - kmin) / ncell for k in range(lo, khalf + 1): Lk = sum(g(a - 1, t, k - t) for t in range(D1 + 1, D + 1)) + Fr(r, a) * g(a - 1, D + 1, k - D - 1) Rk = g(a, D2, k - b) i = min(int((Fr(k) / Dhi - kmin) / w), ncell - 1) Phis, fb, fZ = cell_consts(a, J, b0, c0, kmin + i * w, kmin + (i + 1) * w) Lb = (sum(p * h for p, h in zip(Phis, dHv)) + (Qv if fb else 0)) / factorial(a - 1) Rb = fZ * Rv / factorial(a) assert Lk <= Lb and Rk >= Rb, ("bound violated", a, J, mode, b0, c0, b, c, k) nk += 1 print("selftest: %d random regions, %d degrees k, all bounds hold (with the fractional term) [%.1fs]" % (tests, nk, time.time() - t0)) if __name__ == "__main__": ap = argparse.ArgumentParser() ap.add_argument("--amax", type=int, default=20) ap.add_argument("--amin", type=int, default=3) ap.add_argument("--brute", type=int, default=0) ap.add_argument("--selftest", action="store_true") args = ap.parse_args() if args.selftest: selftest() ok = True if args.amin <= 2: check_a2() for a in range(max(3, args.amin), min(args.amax, max(THRESHOLDS)) + 1): ok = certify(a) and ok if args.brute: ok = brute(args.brute) and ok print("ALL CHECKS PASSED" if ok else "SOME CHECK FAILED")