"""Direct lattice count for GL_{2m}/Sp_{2m}, independent of the Lagrangian reformulation.

N(lam, c) = #{Lambda'' containing O^n of type lam : x_c unimodular on Lambda''}, n = 2m,
x_c = (+)_i t^{c_i} J_2.  Enumerate Lambda'' as duals L^vee = H^{-1} O^n of sublattices
L = (row space of H) of O^n, H in Hermite normal form over O = F_p[[t]]:
upper triangular, H_jj = t^{d_j}, H_ij (i<j) a polynomial of degree < d_j.
Lambda''/O^n is dual to O^n/L, so has the same type.  x_c is unimodular on Lambda''
iff the Gram matrix H^{-T} x_c H^{-1} is integral (val Pf x_c = |c| = sum d_j).

Usage: python3 -I symp_hnf.py p m "lam1;lam2;..." "c1;c2;..."
"""
import sys, itertools


def ladd(a, b, p):
    r = dict(a)
    for e, x in b.items():
        r[e] = (r.get(e, 0) + x) % p
        if r[e] == 0:
            del r[e]
    return r


def lmul(a, b, p):
    r = {}
    for e1, x1 in a.items():
        for e2, x2 in b.items():
            e = e1 + e2
            r[e] = (r.get(e, 0) + x1 * x2) % p
            if r[e] == 0:
                del r[e]
    return r


def lneg(a, p):
    return {e: (-x) % p for e, x in a.items()}


def val(a):
    return min(a) if a else 10 ** 9


def det(M, p):
    n = len(M)
    if n == 1:
        return M[0][0]
    total = {}
    for j in range(n):
        if not M[0][j]:
            continue
        minor = [row[:j] + row[j + 1:] for row in M[1:]]
        term = lmul(M[0][j], det(minor, p), p)
        if j % 2:
            term = lneg(term, p)
        total = ladd(total, term, p)
    return total


def smith_exponents(M, p):
    """Sorted Smith exponents e_1<=...<=e_n of a square matrix of Laurent polys (via minors)."""
    n = len(M)
    dk = [0]
    for k in range(1, n + 1):
        best = 10 ** 9
        for rows in itertools.combinations(range(n), k):
            for cols in itertools.combinations(range(n), k):
                sub = [[M[r][c] for c in cols] for r in rows]
                best = min(best, val(det(sub, p)))
        dk.append(best)
    return [dk[k] - dk[k - 1] for k in range(1, n + 1)]


def hnf_iter(n, D, p):
    for d in itertools.product(range(D + 1), repeat=n):
        if sum(d) != D:
            continue
        slots = [(i, j) for j in range(n) for i in range(j)]
        choices = []
        for (i, j) in slots:
            choices.append(list(itertools.product(range(p), repeat=d[j])))
        for pick in itertools.product(*choices):
            H = [[{} for _ in range(n)] for _ in range(n)]
            for j in range(n):
                H[j][j] = {d[j]: 1}
            for (i, j), coeffs in zip(slots, pick):
                H[i][j] = {e: c for e, c in enumerate(coeffs) if c}
            yield d, H


def inv_upper(H, p):
    n = len(H)
    X = [[{} for _ in range(n)] for _ in range(n)]
    for j in range(n):
        dj = next(iter(H[j][j]))
        X[j][j] = {-dj: 1}
        for i in range(j - 1, -1, -1):
            s = {}
            for k in range(i + 1, j + 1):
                if H[i][k] and X[k][j]:
                    s = ladd(s, lmul(H[i][k], X[k][j], p), p)
            di = next(iter(H[i][i]))
            X[i][j] = lneg(lmul({-di: 1}, s, p), p)
    return X


def xform(c, p):
    n = 2 * len(c)
    x = [[{} for _ in range(n)] for _ in range(n)]
    for i, ci in enumerate(c):
        x[2 * i][2 * i + 1] = {ci: 1}
        x[2 * i + 1][2 * i] = {ci: (-1) % p}
    return x


def gram_integral(Hinv, x, p):
    n = len(x)
    # G = Hinv^T x Hinv
    xH = [[{} for _ in range(n)] for _ in range(n)]
    for i in range(n):
        for j in range(n):
            s = {}
            for k in range(n):
                if x[i][k] and Hinv[k][j]:
                    s = ladd(s, lmul(x[i][k], Hinv[k][j], p), p)
            xH[i][j] = s
    for i in range(n):
        for j in range(i + 1, n):
            s = {}
            for k in range(n):
                if Hinv[k][i] and xH[k][j]:
                    s = ladd(s, lmul(Hinv[k][i], xH[k][j], p), p)
            if s and val(s) < 0:
                return False
    return True


if __name__ == '__main__':
    p = int(sys.argv[1]); m = int(sys.argv[2])
    lams = [tuple(int(a) for a in s.split(',')) for s in sys.argv[3].split(';')]
    cs = [tuple(int(a) for a in s.split(',')) for s in sys.argv[4].split(';')]
    n = 2 * m
    for lam in lams:
        D = sum(lam)
        target = sorted(lam)
        target = [0] * (n - len(target)) + target
        target = sorted(target)
        row = []
        cnt_type = 0
        xs = [xform(c, p) for c in cs]
        counts = [0] * len(cs)
        for d, H in hnf_iter(n, D, p):
            e = smith_exponents(H, p)
            if sorted(e) != target:
                continue
            cnt_type += 1
            Hinv = inv_upper(H, p)
            for k, (c, x) in enumerate(zip(cs, xs)):
                if sum(c) != D:
                    continue
                if gram_integral(Hinv, x, p):
                    counts[k] += 1
        print('lam', lam, 'lattices of this type', cnt_type, 'N(lam, c) for c in', cs, '=', counts)