# S1 engine: enumerate lattices M in F^n (F = F_p((t))) in Hermite normal form,
# keep those on which a given alternating form x is unimodular, and record
# the elementary-divisor type inv(O^n, M).  Method: column-by-column HNF with
# the integrality conditions x(h_i, h_j) in O solved as affine-linear systems
# over F_p for each new column.  Independent of the second audit's scripts
# (which count Lagrangian submodules of an abstract finite symplectic module).
#
# Conventions.  A Laurent polynomial is a dict {exponent: coeff mod p}.
# x is given by block exponents a = (a_1,...,a_m): x = sum_p t^{a_p} J_2 on
# the coordinate pairs (2p, 2p+1) (0-indexed), J_2(e,f) = 1.
# HNF: M has a unique basis h_j = t^{d_j} e_j + sum_{i<j} H_ij e_i with
# H_ij a Laurent polynomial supported in [lo_i, d_i - 1], where lo_i is the
# a-priori lower bound M subset of (+)_i t^{lo_i} O e_i.

import itertools, sys
from collections import Counter

def lmul(f, g, p):
    r = {}
    for e1, c1 in f.items():
        for e2, c2 in g.items():
            e = e1 + e2
            r[e] = (r.get(e, 0) + c1 * c2) % p
    return {e: c for e, c in r.items() if c}

def ladd(f, g, p, s=1):
    r = dict(f)
    for e, c in g.items():
        r[e] = (r.get(e, 0) + s * c) % p
    return {e: c for e, c in r.items() if c}

def form_row(a, u, p):
    """w = x(u, .) as a list of Laurent polys: w[k] = x(u, e_k)."""
    n = 2 * len(a)
    w = [dict() for _ in range(n)]
    for q_, ap in enumerate(a):
        i, j = 2 * q_, 2 * q_ + 1
        # x(u, e_j) = t^ap * u_i ; x(u, e_i) = -t^ap * u_j
        w[j] = {e + ap: c % p for e, c in u[i].items() if c % p}
        w[i] = {e + ap: (-c) % p for e, c in u[j].items() if c % p}
    return w

def pair(w, v, p):
    r = {}
    for k in range(len(v)):
        if w[k] and v[k]:
            r = ladd(r, lmul(w[k], v[k], p), p)
    return r

def solve_affine(rows, rhs, nvar, p):
    """Solve rows * z = rhs over F_p.  Return (particular, nullbasis) or None."""
    A = [list(r) + [b % p] for r, b in zip(rows, rhs)]
    piv_cols = []
    r = 0
    for c in range(nvar):
        pr = None
        for i in range(r, len(A)):
            if A[i][c] % p:
                pr = i; break
        if pr is None:
            continue
        A[r], A[pr] = A[pr], A[r]
        inv = pow(A[r][c], p - 2, p)
        A[r] = [(v * inv) % p for v in A[r]]
        for i in range(len(A)):
            if i != r and A[i][c] % p:
                f = A[i][c]
                A[i] = [(vi - f * vr) % p for vi, vr in zip(A[i], A[r])]
        piv_cols.append(c)
        r += 1
        if r == len(A):
            break
    for i in range(r, len(A)):
        if A[i][nvar] % p:
            return None
    free = [c for c in range(nvar) if c not in piv_cols]
    part = [0] * nvar
    for i, c in enumerate(piv_cols):
        part[c] = A[i][nvar]
    basis = []
    for fcol in free:
        v = [0] * nvar
        v[fcol] = 1
        for i, c in enumerate(piv_cols):
            v[c] = (-A[i][fcol]) % p
        basis.append(v)
    return part, basis

def snf_type(cols, n, p, shift, prec):
    """Elementary-divisor exponents of the lattice spanned by cols (Laurent),
    sorted descending.  shift: multiply by t^shift to make integral."""
    # matrix entries as truncated power series lists of length prec
    Mx = [[[0] * prec for _ in range(n)] for _ in range(n)]  # Mx[row][col]
    for j, h in enumerate(cols):
        for i in range(n):
            for e, c in h[i].items():
                ee = e + shift
                assert ee >= 0
                if ee < prec:
                    Mx[i][j][ee] = c % p
    def val(f):
        for k, c in enumerate(f):
            if c:
                return k
        return prec
    def mul(f, g):
        r = [0] * prec
        for i, a in enumerate(f):
            if a:
                for j in range(prec - i):
                    if g[j]:
                        r[i + j] = (r[i + j] + a * g[j]) % p
        return r
    def inv_unit(f):
        # inverse of a unit power series
        g = [0] * prec
        g[0] = pow(f[0], p - 2, p)
        for k in range(1, prec):
            s = 0
            for i in range(1, k + 1):
                s += f[i] * g[k - i]
            g[k] = (-s * g[0]) % p
        return g
    exps = []
    rows = list(range(n)); colsI = list(range(n))
    A = Mx
    active_r = list(range(n)); active_c = list(range(n))
    while active_r:
        best = None
        for i in active_r:
            for j in active_c:
                v = val(A[i][j])
                if best is None or v < best[0]:
                    best = (v, i, j)
        v, pi, pj = best
        if v >= prec:
            raise ValueError('precision too low')
        exps.append(v - shift)
        # pivot f = t^v * u
        f = A[pi][pj]
        u = f[v:] + [0] * v
        ui = inv_unit(u)
        # eliminate other rows in column pj
        for i in active_r:
            if i == pi:
                continue
            g = A[i][pj]
            vg = val(g)
            if vg >= prec:
                continue
            # factor = (g / t^v) * u^{-1}
            q = mul(g[v:] + [0] * v, ui)
            for j in active_c:
                if any(A[pi][j]):
                    prod = mul(q, A[pi][j])
                    A[i][j] = [(a - b) % p for a, b in zip(A[i][j], prod)]
        active_r.remove(pi); active_c.remove(pj)
    return tuple(sorted(exps, reverse=True))

def enumerate_selfdual(a, p, lo, hi, total_d, callback, x_shift=None):
    """Enumerate lattices M with HNF diagonal exponents d_i in [lo_i, hi_i],
    sum d = total_d, M subset of (+) t^{lo_i} O e_i, on which x_a is integral.
    Since sum d = -val Pf(x), integral == unimodular.  callback(cols) for each."""
    n = 2 * len(a)
    cols = []
    ws = []
    def rec(j, dsum):
        if j == n:
            if dsum == total_d:
                callback(cols)
            return
        # remaining bounds for feasibility of total_d
        rem_lo = sum(lo[k] for k in range(j + 1, n))
        rem_hi = sum(hi[k] for k in range(j + 1, n))
        for dj in range(lo[j], hi[j] + 1):
            if not (rem_lo <= total_d - dsum - dj <= rem_hi):
                continue
            # variables: coefficients of H_ij for i<j, exponents lo_i .. d_i - 1
            var = []
            for i in range(j):
                di = DIAG[i]
                for e in range(lo[i], di):
                    var.append((i, e))
            nvar = len(var)
            base = [dict() for _ in range(n)]
            base[j] = {dj: 1}
            rows, rhs = [], []
            for ip in range(j):
                w = ws[ip]
                # constant part: x(h_ip, t^dj e_j)
                const = {e + dj: c for e, c in w[j].items()}
                # linear part: coefficient maps
                lin = {}
                for vi, (i, e) in enumerate(var):
                    for ee, c in w[i].items():
                        E = ee + e
                        if E < 0:
                            lin.setdefault(E, {})[vi] = (lin.get(E, {}).get(vi, 0) + c) % p
                Es = set(E for E in const if E < 0) | set(lin.keys())
                for E in Es:
                    row = [0] * nvar
                    for vi, c in lin.get(E, {}).items():
                        row[vi] = c
                    rows.append(row)
                    rhs.append((-const.get(E, 0)) % p)
            sol = solve_affine(rows, rhs, nvar, p) if rows else ([0] * nvar, [[1 if k == l else 0 for k in range(nvar)] for l in range(nvar)])
            if sol is None:
                continue
            part, basis = sol
            DIAG.append(dj)
            for coeffs in itertools.product(range(p), repeat=len(basis)):
                z = list(part)
                for cf, b in zip(coeffs, basis):
                    if cf:
                        z = [(zi + cf * bi) % p for zi, bi in zip(z, b)]
                h = [dict() for _ in range(n)]
                h[j] = {dj: 1}
                for vi, (i, e) in enumerate(var):
                    if z[vi]:
                        h[i][e] = z[vi]
                cols.append(h)
                ws.append(form_row(a, h, p))
                rec(j + 1, dsum + dj)
                cols.pop(); ws.pop()
            DIAG.pop()
    DIAG = []
    rec(0, 0)

def count_types(a, p, lo, hi, total_d, shift, prec, filt=None):
    n = 2 * len(a)
    cnt = Counter()
    def cb(cols):
        ty = snf_type(cols, n, p, shift, prec)
        if filt is None or filt(ty):
            cnt[ty] += 1
    enumerate_selfdual(a, p, lo, hi, total_d, cb)
    return cnt