# Fibre counts for the convolution map of X = SL2^3 \ Sp4 x Sp4 (Sakellaridis Ex. 8.0.1)
# Orbit label nu = (d1, d2, e); Hecke label lam = ((a1,b1),(a2,b2)), a_i >= b_i >= 0.
# N(lam, nu)(q) = #{(C,Q1,Q2) : pos_1(C+Q1) = lam1, pos_2(phiC+Q2) = lam2} = (T_lam * phi_0)(x_nu)
# Polynomials in q are dicts {deg: coeff}.
import itertools, sys

def padd(a, b):
    r = dict(a)
    for k, v in b.items():
        r[k] = r.get(k, 0) + v
        if r[k] == 0: del r[k]
    return r

def pmul(a, b):
    r = {}
    for i, x in a.items():
        for j, y in b.items():
            r[i+j] = r.get(i+j, 0) + x*y
    return {k: v for k, v in r.items() if v != 0}

ONE = {0: 1}
Q = {1: 1}
def qpow(n): return {n: 1}
def qm1(): return {1: 1, 0: -1}
def qp1(): return {1: 1, 0: 1}

def h(e, D1, D2):
    """# vertices C with dist(C,P1)=D1, dist(C,P2)=D2, dist(P1,P2)=e, as poly in q."""
    if (D1 + D2 - e) % 2: return {}
    hh = (D1 + D2 - e)//2
    x = (D1 - D2 + e)//2
    if hh < 0 or x < 0 or x > e: return {}
    if hh == 0: return dict(ONE)
    if e == 0: return pmul(qp1(), qpow(hh-1))
    if x == 0 or x == e: return qpow(hh)
    return pmul(qm1(), qpow(hh-1))

def pos(d, DC, DQ, sigma):
    u2 = DC - d   # = 2u
    w2 = DQ - d   # = 2w
    lam1 = (2*d + max(u2, w2))//2
    tot = (2*d + max(u2, 0) + max(w2, 0) + sigma)//2
    return (lam1, tot - lam1)

def g(d, DC, lam, Rmax):
    """# Q (dist(Q,R) = d mod 2) with pos(C+Q) = lam, given dist(C,P)=DC."""
    rho = min(DC, d)
    total = {}
    for DQ in range(d % 2, Rmax+1, 2):
        rhop = min(DQ, d)
        if DQ > d:
            ext = pmul(qp1() if d == 0 else Q, qpow(DQ - d - 1))
        else:
            ext = dict(ONE)
        for sigma in range(0, 2*d+1):
            cnt = h(rho, rhop, sigma)
            if not cnt: continue
            if pos(d, DC, DQ, sigma) != lam: continue
            total = padd(total, pmul(cnt, ext))
    return total

def N(lam, nu, Rmax=None):
    (a1, b1), (a2, b2) = lam
    d1, d2, e = nu
    if Rmax is None:
        Rmax = 2*(a1 + a2 + b1 + b2) + 2*(d1 + d2 + e) + 4
    total = {}
    for D1 in range(d1 % 2, Rmax+1, 2):
        g1 = g(d1, D1, (a1, b1), Rmax)
        if not g1: continue
        for D2 in range(d2 % 2, Rmax+1, 2):
            cnt = h(e, D1, D2)
            if not cnt: continue
            g2 = g(d2, D2, (a2, b2), Rmax)
            if not g2: continue
            total = padd(total, pmul(cnt, pmul(g1, g2)))
    return total

def deg(p): return max(p) if p else None
def lead(p): return p[max(p)] if p else 0
def two_rho(lam): (a1,b1),(a2,b2) = lam; return 4*a1+2*b1+4*a2+2*b2
def delta(nu): d1,d2,e = nu; return 3*d1+3*d2+e
def phi(lam): (a1,b1),(a2,b2) = lam; return (a1+b1, a2+b2, a1-b1+a2-b2)
def pstr(p):
    if not p: return "0"
    return " + ".join(f"{p[k]}q^{k}" for k in sorted(p, reverse=True))