# G2 audit: (1) orbit sizes of alternating forms mod t^3 on F_q[t]^4 by elementary-divisor type,
# classified by the minimal entry valuation and v(Pf) (q = 2 brute force; closed formulas for all q);
# (2) dimension of the stabiliser Lie algebra {X : X w + w X^T = 0 mod t^N} of w = (+) t^{a_j} J_2
# by linear algebra over F_p, compared with (2m^2+m)N + sum_j (4j-3) a_j;
# (3) consequence: finite-type orbit dimensions and codimensions in EVERY component, compared with
# <rho, mu - nu>; the scalar shift sigma(w) = t w moves the Pfaffian valuation by m.
# python3 -I, standard library only.
import itertools
from collections import Counter

# ---------------- (1) orbit sizes mod t^3 over F_2 ----------------
def orbit_sizes_mod_t3(q=2, N=3):
    pairs = [(0,1),(0,2),(0,3),(1,2),(1,3),(2,3)]
    def mul(a, b):
        c = [0]*N
        for i in range(N):
            if a[i]:
                for j in range(N - i):
                    c[i+j] = (c[i+j] + a[i]*b[j]) % q
        return c
    def val(a):
        for i, x in enumerate(a):
            if x: return i
        return N
    cnt = Counter()
    polys = list(itertools.product(range(q), repeat=N))
    for w in itertools.product(polys, repeat=6):
        w01, w02, w03, w12, w13, w23 = w
        pf = [(x - y + z) % q for x, y, z in zip(mul(w01, w23), mul(w02, w13), mul(w03, w12))]
        vP = val(pf)
        b = min(val(x) for x in w)          # smallest elementary divisor exponent
        if vP < N:
            cnt[(vP - b, b)] += 1           # type (a, b) with a + b = v(Pf)
    return cnt

c = orbit_sizes_mod_t3()
q = 2
print("(1) q=2, N=3: #type(2,0) =", c[(2,0)], " formula q^10(q-1)(q^3-1)(q^2+1) =", q**10*(q-1)*(q**3-1)*(q**2+1))
print("               #type(1,1) =", c[(1,1)], " formula q^8(q-1)(q^3-1) =", q**8*(q-1)*(q**3-1))
# closed-form derivation used in the audit (valid for every q): type (1,1): w0 = 0, w1 nondegenerate, w2 free;
# type (2,0): w0 of rank 2 [(q^2+1)(q^3-1) choices], B(w0,w1)=0 [q^5], Pf(w1)+B(w0,w2) != 0 [q^6-q^5].
for qq in (2, 3, 5, 7):
    nondeg = qq**2*(qq-1)*(qq**3-1)
    t11 = qq**6 * nondeg
    t20 = (qq**2+1)*(qq**3-1) * qq**5 * (qq**6 - qq**5)
    assert t11 == qq**8*(qq-1)*(qq**3-1) and t20 == qq**10*(qq-1)*(qq**3-1)*(qq**2+1)
print("    closed-form counts agree with the stated polynomials for q = 2,3,5,7")

# ---------------- (2) stabiliser Lie algebra dimension ----------------
def rank_mod_p(rows, p):
    rows = [r[:] for r in rows]
    rank = 0; ncols = len(rows[0]) if rows else 0
    for col in range(ncols):
        piv = next((i for i in range(rank, len(rows)) if rows[i][col] % p), None)
        if piv is None: continue
        rows[rank], rows[piv] = rows[piv], rows[rank]
        inv = pow(rows[rank][col], p - 2, p)
        rows[rank] = [(x * inv) % p for x in rows[rank]]
        for i in range(len(rows)):
            if i != rank and rows[i][col] % p:
                f = rows[i][col]
                rows[i] = [(x - f * y) % p for x, y in zip(rows[i], rows[rank])]
        rank += 1
    return rank

def stab_dim(a, N, p):
    m = len(a); n = 2*m
    # w as an n x n matrix of polynomials mod t^N: w[i][j][e]
    w = [[[0]*N for _ in range(n)] for _ in range(n)]
    for j, e in enumerate(a):
        w[2*j][2*j+1][e] = 1
        w[2*j+1][2*j][e] = p - 1
    # unknowns X[r][s][e], index ((r*n)+s)*N+e ; equations: (X w + w X^T)[i][k][deg] = 0 for i<k, deg<N
    nv = n*n*N
    eqs = []
    for i in range(n):
        for k in range(i+1, n):
            for deg in range(N):
                row = [0]*nv
                # (X w)[i][k] = sum_l X[i][l] w[l][k]
                for l in range(n):
                    for e2 in range(N):
                        if w[l][k][e2]:
                            e1 = deg - e2
                            if 0 <= e1 < N:
                                row[((i*n)+l)*N+e1] = (row[((i*n)+l)*N+e1] + w[l][k][e2]) % p
                # (w X^T)[i][k] = sum_l w[i][l] X[k][l]
                for l in range(n):
                    for e2 in range(N):
                        if w[i][l][e2]:
                            e1 = deg - e2
                            if 0 <= e1 < N:
                                row[((k*n)+l)*N+e1] = (row[((k*n)+l)*N+e1] + w[i][l][e2]) % p
                eqs.append(row)
    return nv - rank_mod_p(eqs, p)

def formula(a, N):
    m = len(a)
    return (2*m*m + m)*N + sum((4*(j+1) - 3)*aj for j, aj in enumerate(a))

print("(2) stabiliser Lie algebra dimensions (p = 101 and p = 3):")
cases = [((2,0),3),((1,1),3),((2,1),4),((3,0),4),((1,0),3),((0,0),2),((2,2,0),3),((1,1,1),3),((2,1,1),3),((2,1,0),3),((1,0,0),2),((2,2,1),3)]
for a, N in cases:
    d1 = stab_dim(a, N, 101); d2 = stab_dim(a, N, 3)
    print(f"    a={a} N={N}: dim={d1} (p=101), {d2} (p=3), formula={formula(a,N)}", "OK" if d1 == d2 == formula(a, N) else "MISMATCH")

# ---------------- (3) codimensions in every component ----------------
def orbit_dim(a, N):
    m = len(a)
    return 4*m*m*N - formula(a, N)

def rho_pair(m, j):  # <rho, e_{2j-1} + e_{2j}> for GL_{2m}, j = 1..m
    return 2*m + 2 - 4*j

def rho_pairing(m, a):
    return sum(rho_pair(m, j+1)*a[j] for j in range(m))

def labels(m, s, lo, hi):
    return [a for a in itertools.product(range(hi, lo-1, -1), repeat=m) if sum(a) == s and list(a) == sorted(a, reverse=True)]

bad = 0; checked = 0
for m in (2, 3, 4):
    for s in range(0, 2*m + 2):
        labs = labels(m, s, 0, s)    # entries >= 0 (shift by sigma otherwise); finite-type model valid for N > a_1
        N = s + 2
        for mu in labs:
            for nu in labs:
                # dominance nu <= mu: partial sums
                if all(sum(mu[:i]) >= sum(nu[:i]) for i in range(1, m+1)):
                    codim = orbit_dim(mu, N) - orbit_dim(nu, N)
                    checked += 1
                    if codim != rho_pairing(m, mu) - rho_pairing(m, nu): bad += 1
print(f"(3) codim_(closure O_mu) O_nu = <rho, mu-nu> from the Lie-algebra count: {checked} pairs in components s=0..2m+1, m=2,3,4; failures: {bad}")
print("    sigma(w) = t w multiplies Pf by t^m: components s and s+m are identified; e.g. m=2: s odd never reaches s=0")