"""
myarr.py -- verifier B's own library for central hyperplane arrangements.

Written from the definitions, BEFORE the programs of the claim were opened.

Conventions
-----------
* An arrangement is a list of n non-zero integer vectors a_H in Z^k (the linear forms alpha_H = sum a_j x_j),
  pairwise non-proportional.  The field is Q (p = 0) or F_p (p prime; the vectors are read mod p).
* A subset of the arrangement is a bitmask (bit h <-> hyperplane h).
* rk(mask) = dim span{alpha_h : h in mask} = codim of the intersection of these hyperplanes.
* A flat is a closed mask: F = closure(F) = {h : alpha_h in span(alpha_g : g in F)}.  The flat F corresponds to
  the subspace X_F = intersection of the hyperplanes in F (X of codimension rk(F)), and F = A_X.
  L(A) <-> set of flats; reverse inclusion of subspaces <-> inclusion of flats.
* Modularity is tested from the definition "X + Y in L(A) for all Y in L(A)":
      X_F + X_G in L(A)  <=>  X_F + X_G = intersection of the hyperplanes containing X_F + X_G
                         <=>  dim(X_F + X_G) = k - rk(F & G)                [F & G is the set of these hyperplanes]
  and dim(X_F + X_G) = dim X_F + dim X_G - dim(X_F cap X_G) = (k - rk F) + (k - rk G) - (k - rk(F | G)).
  `is_modular` uses these rank numbers; `is_modular_literal` builds the subspaces explicitly (bases of
  null spaces, rank of the stacked bases) and is used as a cross-check.
* Derivations: theta = sum_j P_j d/dx_j, P_j homogeneous of the same degree d ("degree d"; Euler has degree 1).
  theta in D(A)  <=>  theta(alpha_H) = sum_j a_j P_j is divisible by alpha_H for every H.

Freeness test (characteristic 0), see `free_with_exponents`:
  Let pt be a point with Q(pt) != 0 and e_d = dim_Q { theta(pt) : theta in D(A)_d } (a subspace of Q^k).
  CLAIM.  For a multiset d_1 <= ... <= d_k with sum d_i = n:
      A is free with exponents (d_i)   <=>   e_d >= #{i : d_i <= d} for every value d occurring among the d_i.
  PROOF.  (=>) If theta_1..theta_k is a homogeneous basis with deg theta_i = d_i, then by Saito's criterion
  det[theta_i(x_j)] = c*Q with c != 0, so the vectors theta_i(pt) are linearly independent; multiplying theta_i by
  a suitable power of a linear form not vanishing at pt shows theta_i(pt) in E_d for every d >= d_i, so e_d >=
  #{i : d_i <= d}.  (<=) E_d is increasing in d (multiply by a linear form l with l(pt) != 0).  Choose, degree by
  degree, elements theta_i in D(A)_{d_i} whose values at pt are linearly independent (possible because
  e_d >= #{i : d_i <= d}).  Then det[theta_i(x_j)] does not vanish at pt, so theta_1..theta_k are linearly
  independent over S; they are homogeneous elements of D(A) with sum of degrees n = |A|, hence a basis of D(A)
  ([OT, Theorem 4.23]; equivalently: det is a multiple of Q of degree n, so det = c*Q, c != 0, and Saito applies).
  Exponents of a free arrangement are unique, so the exponents are exactly (d_i).
  (If sum d_i != n the answer is "no", because the exponents of a free arrangement sum to n.)

  All dimensions over Q are computed exactly: the null space of the integer constraint matrix is computed mod
  P31 = 2^31-1, every basis vector is lifted by rational reconstruction and VERIFIED exactly over Z; if all
  nullity_p vectors verify, then nullity_Q >= nullity_p, and nullity_Q <= nullity_p holds always (the rank of an
  integer matrix cannot go up mod p), so the lifted vectors are an exact Q-basis.  If a lift fails, an exact
  Fraction elimination is used instead.  The early exit "nullity_p < required dimension => not free with these
  exponents" is rigorous for the same reason.
"""
import itertools
import math
import random
from fractions import Fraction
from functools import lru_cache

import numpy as np

P31 = 2147483647  # 2^31 - 1, prime


# --------------------------------------------------------------------------------------------------
# small exact linear algebra on integer vectors (field Q if p == 0, else F_p)
# --------------------------------------------------------------------------------------------------

def _norm_row(v, p):
    if p:
        return [x % p for x in v]
    g = 0
    for x in v:
        g = math.gcd(g, x)
    if g > 1:
        return [x // g for x in v]
    return list(v)


def reduce_vec(v, basis, p=0):
    """Reduce the integer vector v by the echelon basis [(pivot_col, row), ...] (built by `add_to_basis`).
    The result is 0 iff v lies in the span."""
    v = list(v)
    for c, row in basis:
        b = v[c]
        if (b % p if p else b):
            a = row[c]
            v = [a * x - b * y for x, y in zip(v, row)]
            if p:
                v = [x % p for x in v]
    return _norm_row(v, p)


def add_to_basis(basis, v, p=0):
    """Try to enlarge the echelon basis by v.  Returns True iff v was independent."""
    w = reduce_vec(v, basis, p)
    for c, x in enumerate(w):
        if x:
            basis.append((c, w))
            return True
    return False


def rank_of(vectors, p=0):
    basis = []
    r = 0
    for v in vectors:
        if add_to_basis(basis, v, p):
            r += 1
    return r


def rref_frac(rows, ncols, p=0):
    """Reduced row echelon form over Q (Fractions) or F_p.  Returns (R, pivots)."""
    if p:
        M = [[x % p for x in r] for r in rows]
    else:
        M = [[Fraction(x) for x in r] for r in rows]
    piv = []
    r = 0
    for c in range(ncols):
        pr = None
        for i in range(r, len(M)):
            if M[i][c] != 0:
                pr = i
                break
        if pr is None:
            continue
        M[r], M[pr] = M[pr], M[r]
        pv = M[r][c]
        if p:
            iv = pow(pv, -1, p)
            M[r] = [(x * iv) % p for x in M[r]]
        else:
            M[r] = [x / pv for x in M[r]]
        for i in range(len(M)):
            if i != r and M[i][c] != 0:
                f = M[i][c]
                if p:
                    M[i] = [(a - f * b) % p for a, b in zip(M[i], M[r])]
                else:
                    M[i] = [a - f * b for a, b in zip(M[i], M[r])]
        piv.append(c)
        r += 1
        if r == len(M):
            break
    return M[:r], piv


def nullspace_frac(rows, ncols, p=0):
    """Basis of {x in K^ncols : rows * x = 0} (Fractions over Q, ints mod p over F_p)."""
    if not rows:
        R, piv = [], []
    else:
        R, piv = rref_frac(rows, ncols, p)
    pivset = set(piv)
    basis = []
    for f in range(ncols):
        if f in pivset:
            continue
        v = [0] * ncols
        v[f] = 1
        for i, c in enumerate(piv):
            v[c] = (-R[i][f]) % p if p else -R[i][f]
        basis.append(v)
    return basis


def clear_denominators(v):
    l = 1
    for x in v:
        x = Fraction(x)
        l = l * x.denominator // math.gcd(l, x.denominator)
    w = [int(Fraction(x) * l) for x in v]
    g = 0
    for x in w:
        g = math.gcd(g, x)
    if g > 1:
        w = [x // g for x in w]
    return w


def normalize_vec(v, p=0):
    """Primitive integer representative (Q) or representative with first non-zero entry 1 (F_p) of the form."""
    if p:
        v = [x % p for x in v]
        for x in v:
            if x:
                iv = pow(x, -1, p)
                return tuple((y * iv) % p for y in v)
        raise ValueError("zero vector")
    v = clear_denominators(v)
    for x in v:
        if x:
            if x < 0:
                v = [-y for y in v]
            return tuple(v)
    raise ValueError("zero vector")


def popcount(m):
    return bin(m).count("1")


def bits(m):
    out = []
    i = 0
    while m:
        if m & 1:
            out.append(i)
        m >>= 1
        i += 1
    return out


# --------------------------------------------------------------------------------------------------
# the arrangement class: lattice, L_2, formality, irreducibility, modular elements, supersolvability
# --------------------------------------------------------------------------------------------------

class Arr:
    def __init__(self, vecs, p=0, check=True):
        self.p = p
        self.vecs = [tuple(int(x) for x in v) for v in vecs]
        self.n = len(self.vecs)
        self.k = len(self.vecs[0]) if self.vecs else 0
        self.full = (1 << self.n) - 1
        self._rk = {0: 0}
        self._flats = None
        self._flatrank = None
        if check:
            for v in self.vecs:
                assert any((x % p if p else x) for x in v), "zero form"
            for i in range(self.n):
                for j in range(i):
                    assert rank_of([self.vecs[i], self.vecs[j]], p) == 2, "proportional forms %d %d" % (j, i)

    # ---- rank and closure --------------------------------------------------------------------------
    def rk(self, mask):
        r = self._rk.get(mask)
        if r is None:
            r = rank_of([self.vecs[h] for h in bits(mask)], self.p)
            self._rk[mask] = r
        return r

    def rank(self):
        return self.rk(self.full)

    def is_essential(self):
        return self.rank() == self.k

    def closure(self, mask):
        basis = []
        for h in bits(mask):
            add_to_basis(basis, self.vecs[h], self.p)
        out = mask
        for h in range(self.n):
            if not (mask >> h) & 1:
                w = reduce_vec(self.vecs[h], basis, self.p)
                if not any(w):
                    out |= 1 << h
        return out

    # ---- the lattice of flats ------------------------------------------------------------------------
    def flats(self):
        """dict rank -> sorted list of flats (closed masks); also fills self._flatrank (flat -> rank)."""
        if self._flats is not None:
            return self._flats
        r_top = self.rank()
        levels = {0: [0]}
        frank = {0: 0}
        cur = [0]
        for r in range(1, r_top + 1):
            nxt = set()
            for F in cur:
                covered = F
                for h in range(self.n):
                    if not (covered >> h) & 1:
                        G = self.closure(F | (1 << h))
                        covered |= G
                        nxt.add(G)
            cur = sorted(nxt)
            levels[r] = cur
            for G in cur:
                frank[G] = r
                self._rk[G] = r
        self._flats = levels
        self._flatrank = frank
        return levels

    def all_flats(self):
        self.flats()
        return self._flatrank

    def L2(self):
        """list of the rank-2 flats (masks)."""
        out = set()
        for i in range(self.n):
            for j in range(i):
                m = (1 << i) | (1 << j)
                if not any((m & Z) == m for Z in out):
                    out.add(self.closure(m))
        return sorted(out)

    def multiple_flats(self):
        return [Z for Z in self.L2() if popcount(Z) >= 3]

    def degrees(self):
        """deg(H) = number of rank-2 flats of multiplicity >= 3 containing H."""
        deg = [0] * self.n
        for Z in self.multiple_flats():
            for h in bits(Z):
                deg[h] += 1
        return deg

    def uv(self):
        mult = [popcount(Z) for Z in self.multiple_flats()]
        return mult.count(3), mult.count(4), max(mult) if mult else 2

    def sum_binom(self):
        """sum over rank-2 flats Z of C(m_Z - 1, 2)."""
        return sum((popcount(Z) - 1) * (popcount(Z) - 2) // 2 for Z in self.L2())

    # ---- relation space and formality ------------------------------------------------------------------
    def relation_space_dim(self):
        return self.n - self.rank()

    def local_relations(self, Z):
        """basis of F_Z: the relations among the forms of the rank-2 flat Z (m_Z - 2 vectors in K^n)."""
        hs = bits(Z)
        if len(hs) < 3:
            return []
        rows = [list(col) for col in zip(*[self.vecs[h] for h in hs])]  # k rows, m columns
        ns = nullspace_frac(rows, len(hs), self.p)
        out = []
        for v in ns:
            if not self.p:
                v = clear_denominators(v)
            r = [0] * self.n
            for h, x in zip(hs, v):
                r[h] = x
            out.append(r)
        assert len(out) == len(hs) - 2
        return out

    def formal_rank(self):
        """dim of sum_Z F_Z."""
        rel = []
        for Z in self.multiple_flats():
            rel.extend(self.local_relations(Z))
        return rank_of(rel, self.p)

    def is_formal(self):
        return self.formal_rank() == self.relation_space_dim()

    def formal_rank_by_triples(self):
        """dim of the span of all relations among three forms (the definition of formality, literally)."""
        rel = []
        for a, b, c in itertools.combinations(range(self.n), 3):
            m = (1 << a) | (1 << b) | (1 << c)
            if self.rk(m) == 2:
                rows = [list(col) for col in zip(self.vecs[a], self.vecs[b], self.vecs[c])]
                ns = nullspace_frac(rows, 3, self.p)
                assert len(ns) == 1
                v = ns[0] if self.p else clear_denominators(ns[0])
                r = [0] * self.n
                r[a], r[b], r[c] = v
                rel.append(r)
        return rank_of(rel, self.p)

    # ---- irreducibility ---------------------------------------------------------------------------------
    def components(self):
        """connected components of the matroid of the forms, via the fundamental graph of a basis
        (Oxley, Matroid Theory, Prop. 4.3.2).  Returns a list of masks."""
        cols = self.vecs
        rows = [list(r) for r in zip(*cols)]  # k x n matrix with the forms as columns
        R, piv = rref_frac(rows, self.n, self.p)
        parent = list(range(self.n))

        def find(x):
            while parent[x] != x:
                parent[x] = parent[parent[x]]
                x = parent[x]
            return x

        for i, c in enumerate(piv):
            for e in range(self.n):
                if e != c and R[i][e] != 0:
                    parent[find(e)] = find(c)
        comp = {}
        for h in range(self.n):
            comp.setdefault(find(h), 0)
            comp[find(h)] |= 1 << h
        return sorted(comp.values())

    def is_irreducible(self):
        return len(self.components()) == 1

    def is_irreducible_bruteforce(self):
        """from the definition: no partition A = A' u A'' (both non-empty) with rk A' + rk A'' = rk A."""
        r = self.rank()
        for m in range(1, 1 << (self.n - 1)):
            if self.rk(m) + self.rk(self.full ^ m) == r:
                return False
        return True

    # ---- modular elements and supersolvability -----------------------------------------------------------
    def _flat_basis_bits(self, F):
        b = []
        basis = []
        for h in bits(F):
            if add_to_basis(basis, self.vecs[h], self.p):
                b.append(h)
        return b

    def is_modular(self, F, flatlist=None):
        """F (a flat) is modular iff for every flat G:  rk F + rk G - rk(F | G) = rk(F & G)
        (this is 'X_F + X_G in L(A)', see the module docstring)."""
        frank = self.all_flats()
        rF = frank[F]
        bF = self._bb(F)
        if flatlist is None:
            flatlist = self._flat_order()
        for G in flatlist:
            if (F & G) == G or (F & G) == F:
                continue  # comparable: the sum is the larger subspace
            rG = frank[G]
            meet = F & G                      # hyperplanes containing X_F + X_G ; a flat
            rmeet = frank.get(meet)
            if rmeet is None:
                rmeet = self.rk(meet)
            rjoin = rank_of([self.vecs[h] for h in bF] + [self.vecs[h] for h in self._bb(G)], self.p)
            if rF + rG - rjoin != rmeet:
                return False
        return True

    def _bb(self, F):
        d = self.__dict__.setdefault("_bbcache", {})
        b = d.get(F)
        if b is None:
            b = self._flat_basis_bits(F)
            d[F] = b
        return b

    def _flat_order(self):
        o = self.__dict__.get("_forder")
        if o is None:
            lev = self.flats()
            o = []
            for r in sorted(lev):
                o.extend(lev[r])
            self._forder = o
        return o

    def subspace_basis(self, F):
        """basis of X_F = intersection of the hyperplanes in F (list of vectors in K^k)."""
        rows = [list(self.vecs[h]) for h in bits(F)]
        return nullspace_frac(rows, self.k, self.p)

    def is_modular_literal(self, F):
        """the definition with explicit subspaces: for all flats G, X_F + X_G is an element of L(A)."""
        XF = self.subspace_basis(F)
        for G in self._flat_order():
            XG = self.subspace_basis(G)
            if XF or XG:
                if self.p:
                    dsum = rank_of([[int(x) for x in v] for v in XF + XG], self.p)
                else:
                    dsum = rank_of([clear_denominators(v) for v in XF + XG], 0)
            else:
                dsum = 0
            # the hyperplanes of A containing X_F + X_G are exactly those in F & G
            W = self.subspace_basis(F & G)
            if dsum != len(W):
                return False
        return True

    def modular_flats(self):
        return [F for F in self._flat_order() if self.is_modular(F)]

    def modular_chain(self):
        """a maximal chain of modular flats 0 = F_0 < F_1 < ... < F_r = full (r = rank), or None."""
        lev = self.flats()
        r = self.rank()
        memo = {}

        def mod(F):
            v = memo.get(F)
            if v is None:
                v = self.is_modular(F)
                memo[F] = v
            return v

        def down(F, rF):
            # F is modular of rank rF; find modular chain below it
            if rF == 0:
                return [0]
            for G in lev[rF - 1]:
                if (G & F) == G and mod(G):
                    c = down(G, rF - 1)
                    if c is not None:
                        return c + [F]
            return None

        return down(self.full, r)

    def is_supersolvable(self):
        return self.modular_chain() is not None

    # ---- characteristic polynomial ---------------------------------------------------------------------
    def charpoly(self):
        """coefficients c[i] of chi(t) = sum_F mu(F) t^(k - rk F), as list indexed by rank (c[r] multiplies
        t^(k-r))."""
        lev = self.flats()
        mu = {0: 1}
        coeffs = [0] * (self.k + 1)
        coeffs[0] = 1
        order = self._flat_order()
        for F in order:
            if F == 0:
                continue
            s = 0
            for G in order:
                if G != F and (G & F) == G and G in mu:
                    s += mu[G]
            mu[F] = -s
            coeffs[self._flatrank[F]] += mu[F]
        return coeffs

    def charpoly_whitney(self):
        """Whitney: chi(t) = sum over subsets S of (-1)^|S| t^(k - rk S)."""
        coeffs = [0] * (self.k + 1)
        for m in range(1 << self.n):
            coeffs[self.rk(m)] += -1 if popcount(m) & 1 else 1
        return coeffs

    # ---- sub-arrangements, localization, restriction ----------------------------------------------------
    def sub(self, mask):
        return Arr([self.vecs[h] for h in bits(mask)], self.p, check=False)

    def quotient_by_center(self, mask):
        """the essential arrangement A_mask / X, X = centre of the sub-arrangement: coordinates of the forms
        in a basis of their span."""
        hs = bits(mask)
        vs = [self.vecs[h] for h in hs]
        basis = []
        bidx = []
        for i, v in enumerate(vs):
            if add_to_basis(basis, v, self.p):
                bidx.append(i)
        # express every form in the basis (vs[i], i in bidx)
        r = len(bidx)
        rows = [list(row) for row in zip(*[vs[i] for i in bidx])]  # k x r
        out = []
        for v in vs:
            aug = [row + [x] for row, x in zip(rows, v)]
            R, piv = rref_frac(aug, r + 1, self.p)
            assert r not in piv, "form not in span"
            coord = [0] * r
            for i, c in enumerate(piv):
                coord[c] = R[i][r]
            out.append(normalize_vec(coord, self.p))
        return Arr(out, self.p, check=False)

    def restriction(self, h):
        """A^H: the arrangement induced on the hyperplane H = ker(alpha_h) (forms restricted to a basis of H)."""
        B = self.subspace_basis(1 << h)            # k-1 vectors spanning H
        if not self.p:
            B = [clear_denominators(v) for v in B]
        seen = {}
        for g in range(self.n):
            if g == h:
                continue
            w = [sum(a * b for a, b in zip(self.vecs[g], v)) for v in B]
            if any((x % self.p if self.p else x) for x in w):
                seen.setdefault(normalize_vec(w, self.p), g)
        return Arr(list(seen.keys()), self.p, check=False)


# --------------------------------------------------------------------------------------------------
# class membership and the lemmas of the claim (stated here from RESULT.md, tested on examples)
# --------------------------------------------------------------------------------------------------

def numerics_C(A):
    """(C3)-(C5): n = 2k, all m_Z <= 4, sum C(m_Z-1,2) = k+1   (A assumed essential of rank k)."""
    u, v, mx = A.uv()
    return A.n == 2 * A.k and mx <= 4 and A.sum_binom() == A.k + 1


def numerics_Cprime(A):
    """(C3')-(C5'): n = 2k-1, all m_Z <= 3, number of triple points = k-1."""
    u, v, mx = A.uv()
    return A.n == 2 * A.k - 1 and mx <= 3 and u == A.k - 1


def good_triples(A):
    """all (Z0, i, j): Z0 a triple point, H_i, H_j in Z0 with deg = 1."""
    deg = A.degrees()
    out = []
    for Z in A.multiple_flats():
        if popcount(Z) == 3:
            ones = [h for h in bits(Z) if deg[h] == 1]
            for i, j in itertools.combinations(ones, 2):
                out.append((Z, i, j))
    return out


def hypergraph_connected(A):
    """connectedness of the hypergraph (vertices: hyperplanes; hyperedges: multiple rank-2 flats)."""
    parent = list(range(A.n))

    def find(x):
        while parent[x] != x:
            parent[x] = parent[parent[x]]
            x = parent[x]
        return x

    for Z in A.multiple_flats():
        hs = bits(Z)
        for h in hs[1:]:
            parent[find(h)] = find(hs[0])
    return len({find(h) for h in range(A.n)}) == 1


def check_lemma3(A, Z0, i, j, deep=True):
    """Test Lemma 3 (a)-(d) of the claim for the triple point Z0 and its hyperplanes i, j of degree 1.
    Returns a dict of booleans (a_line, a_AX, b_modular, c_modular_lift, d_flats, d_formal, d_irreducible)."""
    res = {}
    k = A.rank()
    AX = A.full & ~((1 << i) | (1 << j))
    res["a_line"] = (A.rk(AX) == k - 1)
    res["a_AX"] = (A.closure(AX) == AX)
    frank = A.all_flats()
    res["a_inL"] = AX in frank
    if not res["a_inL"]:
        return res
    res["b_modular"] = A.is_modular(AX)
    # (d) multiple flats of A_X = multiple flats of A other than Z0, same multiplicities
    multA = sorted(Z for Z in A.multiple_flats() if Z != Z0)
    B = A.sub(AX)
    hs = bits(AX)

    def lift(m):
        out = 0
        for t in bits(m):
            out |= 1 << hs[t]
        return out

    multB = sorted(lift(Z) for Z in B.multiple_flats())
    res["d_flats"] = (multA == multB)
    res["d_formal"] = B.is_formal()
    Bq = A.quotient_by_center(AX)
    res["d_irreducible"] = (Bq.is_irreducible() if A.is_irreducible() else None)
    res["quotient"] = Bq
    if deep:
        # (c) every flat of A_X that is modular in L(A_X) is modular in L(A)
        ok = True
        for F in B._flat_order():
            if B.is_modular(F):
                if not A.is_modular(lift(F)):
                    ok = False
                    break
        res["c_modular_lift"] = ok
    return res


# --------------------------------------------------------------------------------------------------
# logarithmic derivations, degree by degree (characteristic 0)
# --------------------------------------------------------------------------------------------------

@lru_cache(maxsize=None)
def monomials(k, d):
    """exponent tuples of the monomials of degree d in k variables (fixed order) and the index dict."""
    out = []

    def rec(pos, left, cur):
        if pos == k - 1:
            out.append(tuple(cur + [left]))
            return
        for e in range(left, -1, -1):
            rec(pos + 1, left - e, cur + [e])

    if k == 0:
        return ((),) if d == 0 else (), {(): 0} if d == 0 else {}
    rec(0, d, [])
    return tuple(out), {m: i for i, m in enumerate(out)}


def restriction_matrix(a, d):
    """Matrix (rows: monomials of degree d in the k-1 variables y; columns: monomials of degree d in x) of the
    ring homomorphism phi: Q[x] -> Q[y], x_j -> a_piv*y_j (j != piv), x_piv -> -sum_{j != piv} a_j y_j,
    where piv is the first index with a_piv != 0.  phi is surjective with kernel (alpha), alpha = sum a_j x_j,
    so f is divisible by alpha iff phi(f) = 0."""
    k = len(a)
    piv = next(j for j in range(k) if a[j] != 0)
    others = [j for j in range(k) if j != piv]
    monsX, _ = monomials(k, d)
    monsY, idxY = monomials(k - 1, d)
    # powers of L = -sum a_j y_j as dicts {exponent tuple of length k-1 : coefficient}
    L = {}
    for t, j in enumerate(others):
        if a[j] != 0:
            e = [0] * (k - 1)
            e[t] = 1
            L[tuple(e)] = -a[j]
    Lpow = [{tuple([0] * (k - 1)): 1}]
    for _ in range(d):
        nxt = {}
        for e1, c1 in Lpow[-1].items():
            for e2, c2 in L.items():
                e = tuple(x + y for x, y in zip(e1, e2))
                nxt[e] = nxt.get(e, 0) + c1 * c2
        Lpow.append(nxt)
    R = [[0] * len(monsX) for _ in range(len(monsY))]
    ap = a[piv]
    for c, m in enumerate(monsX):
        e = m[piv]
        rest = tuple(m[j] for j in others)
        scale = ap ** (d - e)
        for ey, coef in Lpow[e].items():
            tgt = tuple(x + y for x, y in zip(ey, rest))
            R[idxY[tgt]][c] += scale * coef
    return R


def constraint_block(a, d):
    """integer matrix of the conditions 'theta(alpha) divisible by alpha' on the coefficient vector of
    theta = sum_j P_j d/dx_j, deg P_j = d.  Unknown (j, m) has index j*N + idx(m), N = #monomials."""
    R = np.array(restriction_matrix(a, d), dtype=np.int64)
    return np.hstack([int(aj) * R for aj in a])


def rref_mod(M, p=P31, want_rref=True):
    """row reduction mod p of an int64 matrix (entries reduced mod p).  Returns (R, pivots);
    R is the reduced row echelon form restricted to its non-zero rows if want_rref."""
    M = np.array(M, dtype=np.int64) % p
    nr, nc = M.shape
    piv = []
    r = 0
    for c in range(nc):
        if r == nr:
            break
        col = M[r:, c]
        nz = np.flatnonzero(col)
        if nz.size == 0:
            continue
        pr = r + int(nz[0])
        if pr != r:
            M[[r, pr]] = M[[pr, r]]
        iv = pow(int(M[r, c]), -1, p)
        M[r] = (M[r] * iv) % p
        colv = M[:, c].copy()
        colv[r] = 0
        idx = np.flatnonzero(colv)
        if idx.size:
            M[idx] = (M[idx] - (colv[idx, None] * M[r][None, :]) % p) % p
        piv.append(c)
        r += 1
    return M[:r], piv


def rank_mod(M, p=P31):
    return len(rref_mod(M, p)[1])


_RR_BOUND = int(math.isqrt((P31 - 1) // 2))


def rat_recon(a, p=P31, cache={}):
    """rational reconstruction: n/d = a mod p with |n|, d <= sqrt((p-1)/2); None if impossible."""
    a %= p
    if a in cache:
        return cache[a]
    if a == 0:
        res = (0, 1)
    else:
        r0, r1 = p, a
        s0, s1 = 0, 1
        while r1 > _RR_BOUND:
            q = r0 // r1
            r0, r1 = r1, r0 - q * r1
            s0, s1 = s1, s0 - q * s1
        if s1 == 0 or abs(s1) > _RR_BOUND or math.gcd(r1, abs(s1)) != 1:
            res = None
        else:
            res = (r1, s1) if s1 > 0 else (-r1, -s1)
    if len(cache) < 200000:
        cache[a] = res
    return res


STATS = {"lift_fail": 0, "crt3": 0, "exact_fallback": 0, "nullspaces": 0}
P_B = 2147483629   # 2^31 - 19, prime
P_C = 2147483587   # 2^31 - 61, prime


def _rat_recon_big(a, M, bound):
    a %= M
    if a == 0:
        return (0, 1)
    r0, r1 = M, a
    s0, s1 = 0, 1
    while r1 > bound:
        q = r0 // r1
        r0, r1 = r1, r0 - q * r1
        s0, s1 = s1, s0 - q * s1
    if s1 == 0 or abs(s1) > bound or math.gcd(r1, abs(s1)) != 1:
        return None
    return (r1, s1) if s1 > 0 else (-r1, -s1)


def _lift_try(M, primes):
    """null space of the integer matrix M from its reduced echelon forms modulo the given primes (CRT +
    rational reconstruction); every vector is verified exactly over Z.  Returns the list of integer vectors
    (one per free column of the echelon form modulo primes[0]) or None if a lift or a verification fails."""
    nr, nc = M.shape
    forms = [rref_mod(M, p) for p in primes]
    piv = forms[0][1]
    for R, pv in forms[1:]:
        if pv != piv:
            return None
    pivset = set(piv)
    free = [c for c in range(nc) if c not in pivset]
    maxM = int(np.abs(M).max()) if M.size else 0
    if len(primes) == 1:
        p = primes[0]
        recon = lambda res: rat_recon(res[0], p)
    else:
        mod = 1
        for p in primes:
            mod *= p
        bound = math.isqrt(mod // 2)
        coef = []
        for p in primes:
            mi = mod // p
            coef.append(mi * pow(mi, -1, p))

        def recon(res):
            x = sum(c * r for c, r in zip(coef, res)) % mod
            return _rat_recon_big(x, mod, bound)
    Mobj = None
    basis = []
    for f in free:
        v = [Fraction(0)] * nc
        v[f] = Fraction(1)
        for i, c in enumerate(piv):
            res = [int((-R[i, f]) % p) for (R, _), p in zip(forms, primes)]
            if any(res):
                rr = recon(res)
                if rr is None:
                    return None
                v[c] = Fraction(rr[0], rr[1])
        w = clear_denominators(v)
        mx = max(abs(x) for x in w)
        if mx * maxM * nc < 2 ** 62:
            prod = M @ np.array(w, dtype=np.int64)
            zero = not prod.any()
        else:
            if Mobj is None:
                Mobj = M.astype(object)
            prod = Mobj @ np.array(w, dtype=object)
            zero = all(int(x) == 0 for x in prod)
        if not zero:
            return None
        basis.append(w)
    return basis


def exact_nullspace_int(M):
    """Exact Q-basis (integer vectors) of the null space of the integer matrix M.
    If all nullity_p vectors of the echelon form modulo p = 2^31-1 lift to exact null vectors, they are a basis
    over Q (nullity_Q <= nullity_p always, and the lifted vectors are independent).  Otherwise an exact
    elimination with Fractions is used."""
    M = np.array(M, dtype=np.int64)
    nr, nc = M.shape
    STATS["nullspaces"] += 1
    res = _lift_try(M, [P31])
    if res is not None:
        return res
    STATS["lift_fail"] += 1
    res = _lift_try(M, [P31, P_B, P_C])
    if res is not None:
        STATS["crt3"] += 1
        return res
    res = _lift_try(M, [P31, P_B, P_C, 2147483579, 2147483563, 2147483549])
    if res is not None:
        STATS["crt6"] = STATS.get("crt6", 0) + 1
        return res
    STATS["exact_fallback"] += 1
    ns = nullspace_frac([[int(x) for x in row] for row in M.tolist()], nc, 0)
    return [clear_denominators(v) for v in ns]


class DerSystem:
    """constraint matrices of D(A)_d for an arrangement given by integer vectors (characteristic 0)."""

    def __init__(self, vecs, block_cache=None):
        self.vecs = [tuple(int(x) for x in v) for v in vecs]
        self.n = len(self.vecs)
        self.k = len(self.vecs[0])
        self.cache = block_cache if block_cache is not None else {}

    def block(self, a, d):
        key = (a, d)
        b = self.cache.get(key)
        if b is None:
            b = constraint_block(a, d)
            self.cache[key] = b
        return b

    def matrix(self, d):
        return np.vstack([self.block(a, d) for a in self.vecs])

    def nullity_mod(self, d, p=P31):
        M = self.matrix(d)
        return M.shape[1] - rank_mod(M, p)

    def basis(self, d):
        """exact Q-basis of D(A)_d: integer coefficient vectors, index j*N + idx(m)."""
        if d == 0:
            ns = nullspace_frac([list(a) for a in self.vecs], self.k, 0)
            return [clear_denominators(v) for v in ns]
        return exact_nullspace_int(self.matrix(d))

    def evaluate(self, basis, d, pt):
        """values theta(pt) in Z^k of the derivations given by the coefficient vectors in `basis`."""
        mons, _ = monomials(self.k, d)
        N = len(mons)
        mv = []
        for m in mons:
            x = 1
            for c, e in zip(pt, m):
                x *= c ** e
            mv.append(x)
        vals = []
        for w in basis:
            vals.append([sum(w[j * N + t] * mv[t] for t in range(N) if w[j * N + t]) for j in range(self.k)])
        return vals


def point_off(vecs, rng=None, lo=-40, hi=40):
    rng = rng or random.Random(12345)
    k = len(vecs[0])
    while True:
        pt = [rng.randint(lo, hi) for _ in range(k)]
        if all(sum(a * b for a, b in zip(v, pt)) != 0 for v in vecs):
            return pt


def required_dim(exps, k, d):
    """dim of the degree-d part of the free module with generators in degrees exps."""
    return sum(math.comb(d - e + k - 1, k - 1) for e in exps if d >= e)


def free_with_exponents(vecs, exps, block_cache=None, pt=None, info=None):
    """Exact decision: is the arrangement (integer vectors, char 0) free with exponents exactly `exps`?
    See the module docstring for the proof that the test decides in both directions."""
    vecs = [tuple(int(x) for x in v) for v in vecs]
    n, k = len(vecs), len(vecs[0])
    exps = sorted(exps)
    assert len(exps) == k
    if sum(exps) != n:
        return False
    ds = DerSystem(vecs, block_cache)
    degs = sorted(set(exps))
    # cheap rigorous rejections first: nullity mod p is an upper bound for dim_Q D(A)_d
    for d in degs:
        if d == 0:
            continue
        if ds.nullity_mod(d) < required_dim(exps, k, d):
            if info is not None:
                info["reject"] = ("nullity_mod", d)
            return False
    if pt is None:
        pt = point_off(vecs)
    for d in degs:
        B = ds.basis(d)
        if len(B) < required_dim(exps, k, d):
            if info is not None:
                info["reject"] = ("dimQ", d)
            return False
        vals = ds.evaluate(B, d, pt)
        e_d = rank_of(vals, 0)
        if e_d < sum(1 for e in exps if e <= d):
            if info is not None:
                info["reject"] = ("eval_rank", d, e_d)
            return False
    return True


def free_exponents(vecs, dmax=None, block_cache=None, pt=None):
    """General exact freeness test: returns the exponents if A is free, None if A is not free, and
    'undecided' if the degree bound dmax was reached before the evaluation rank reached k.
    A is free iff sum_d d*(e_d - e_{d-1}) = n, e_d = dim of the values at pt of D(A)_d (see docstring)."""
    vecs = [tuple(int(x) for x in v) for v in vecs]
    n, k = len(vecs), len(vecs[0])
    ds = DerSystem(vecs, block_cache)
    if pt is None:
        pt = point_off(vecs)
    if dmax is None:
        dmax = n
    exps = []
    prev = 0
    total = 0
    for d in range(0, dmax + 1):
        B = ds.basis(d)
        vals = ds.evaluate(B, d, pt)
        e_d = rank_of(vals, 0)
        exps.extend([d] * (e_d - prev))
        total += d * (e_d - prev)
        prev = e_d
        if total > n:
            return None
        if e_d == k:
            return exps if total == n else None
        # lower bound on the final total: the remaining k - e_d generators have degree >= d+1
        if total + (k - e_d) * (d + 1) > n:
            return None
    return "undecided"


def saito_certificate(vecs, exps, block_cache=None, pt=None):
    """explicit derivations theta_1..theta_k (coefficient vectors with their degrees) with independent values
    at pt, or None.  Used for the symbolic re-check with sympy."""
    vecs = [tuple(int(x) for x in v) for v in vecs]
    n, k = len(vecs), len(vecs[0])
    exps = sorted(exps)
    if sum(exps) != n:
        return None
    ds = DerSystem(vecs, block_cache)
    if pt is None:
        pt = point_off(vecs)
    chosen = []
    basis = []
    for d in sorted(set(exps)):
        B = ds.basis(d)
        vals = ds.evaluate(B, d, pt)
        need = sum(1 for e in exps if e == d)
        got = 0
        for w, val in zip(B, vals):
            if got == need:
                break
            if add_to_basis(basis, val, 0):
                chosen.append((d, w))
                got += 1
        if got < need:
            return None
    return chosen


# --------------------------------------------------------------------------------------------------
# ground sets and symmetry
# --------------------------------------------------------------------------------------------------

def ground_pm1(k):
    """all non-zero vectors in {0,+-1}^k up to sign (first non-zero entry +1)."""
    out = []
    for v in itertools.product((0, 1, -1), repeat=k):
        nz = [x for x in v if x]
        if nz and nz[0] == 1:
            out.append(tuple(v))
    return sorted(out)


def ground_B(k):
    out = []
    for i in range(k):
        v = [0] * k
        v[i] = 1
        out.append(tuple(v))
    for i in range(k):
        for j in range(i + 1, k):
            for s in (1, -1):
                v = [0] * k
                v[i] = 1
                v[j] = s
                out.append(tuple(v))
    return sorted(out)


def ground_D(k):
    return [v for v in ground_B(k) if sum(1 for x in v if x) == 2]


def signed_perm_group(ground):
    """the permutations of the ground set induced by the signed coordinate permutations
    (distinct permutations only), as an int array of shape (|G|, N)."""
    k = len(ground[0])
    idx = {v: i for i, v in enumerate(ground)}
    perms = set()
    for sigma in itertools.permutations(range(k)):
        for eps in itertools.product((1, -1), repeat=k):
            img = []
            for v in ground:
                w = [eps[i] * v[sigma[i]] for i in range(k)]
                for x in w:
                    if x:
                        if x < 0:
                            w = [-y for y in w]
                        break
                img.append(idx[tuple(w)])
            perms.add(tuple(img))
    return np.array(sorted(perms), dtype=np.int64)


def burnside_orbit_count(perms, size):
    """number of orbits of `size`-subsets under the permutation group `perms` (Burnside's lemma)."""
    total = 0
    N = perms.shape[1]
    for g in perms:
        seen = [False] * N
        poly = [1] + [0] * size
        for s in range(N):
            if not seen[s]:
                L = 0
                x = s
                while not seen[x]:
                    seen[x] = True
                    x = int(g[x])
                    L += 1
                if L <= size:
                    new = poly[:]
                    for t in range(L, size + 1):
                        new[t] += poly[t - L]
                    poly = new
        total += poly[size]
    assert total % len(perms) == 0
    return total // len(perms)


def _rev_weights(N):
    return (np.int64(1) << (N - 1 - np.arange(N, dtype=np.int64)))


def is_canonical_batch(idx, perms, N):
    """idx: int array (B, size) of sorted element lists.  A set S is CANONICAL if its sorted element tuple is
    lexicographically minimal in its orbit; with the weights w[x] = 2^(N-1-x) this means
    sum w[g(S)] <= sum w[S] for all g.  Returns a boolean array."""
    w = _rev_weights(N)
    B = idx.shape[0]
    own = w[idx].sum(axis=1)
    alive = np.arange(B)
    cur = idx
    curown = own
    for g in perms:
        img = w[g[cur]].sum(axis=1)
        bad = img > curown
        if bad.any():
            keep = ~bad
            alive = alive[keep]
            cur = cur[keep]
            curown = curown[keep]
            if len(alive) == 0:
                break
    out = np.zeros(B, dtype=bool)
    out[alive] = True
    return out


def orbit_representatives(perms, size, verbose=False, seed=1, return_levels=False):
    """Orderly generation (Read/Faradzev): the canonical `size`-subsets of {0..N-1} (lexicographically minimal
    sorted tuples in their orbits) are exactly the sets R + {x} with R a canonical (size-1)-subset, x > max R,
    that are canonical themselves (prefixes of lexicographically minimal tuples are lexicographically minimal).
    Returns an int array (number of orbits, size) of sorted element lists."""
    N = perms.shape[1]
    rng = np.random.default_rng(seed)
    order = rng.permutation(len(perms))
    P = perms[order]
    reps = np.zeros((1, 0), dtype=np.int64)
    levels = []
    for level in range(1, size + 1):
        parts = []
        for x in range(N):
            if level == 1:
                sel = reps
            else:
                sel = reps[reps[:, -1] < x]
            if len(sel):
                parts.append(np.hstack([sel, np.full((len(sel), 1), x, dtype=np.int64)]))
        cand = np.vstack(parts)
        ok = np.zeros(len(cand), dtype=bool)
        chunk = 400000
        for s in range(0, len(cand), chunk):
            ok[s:s + chunk] = is_canonical_batch(cand[s:s + chunk], P, N)
        reps = cand[ok]
        levels.append(len(reps))
        if verbose:
            print("level", level, "candidates", len(cand), "orbits", len(reps), flush=True)
    if return_levels:
        return reps, levels
    return reps


def canonical_form(idx_list, perms, N):
    """the canonical (lexicographically minimal) image of one set, as a sorted tuple."""
    w = _rev_weights(N)
    idx = np.array(sorted(idx_list), dtype=np.int64)
    imgs = perms[:, idx]
    vals = w[imgs].sum(axis=1)
    g = int(np.argmax(vals))
    return tuple(sorted(int(x) for x in imgs[g]))


def orbit_size(idx_list, perms):
    N = perms.shape[1]
    idx = np.array(list(idx_list), dtype=np.int64)
    pow2 = (np.int64(1) << np.arange(N, dtype=np.int64))
    imgs = pow2[perms[:, idx]].sum(axis=1)
    return len(np.unique(imgs))


# --------------------------------------------------------------------------------------------------
# additions (still before the claim's programs were opened): other ground sets, freeness over F_p
# --------------------------------------------------------------------------------------------------

def ground_bin(k):
    """all non-zero 0/1 vectors: the points of PG(k-1, 2) (to be used with p = 2)."""
    return sorted(v for v in itertools.product((0, 1), repeat=k) if any(v))


def ground_box(k, c):
    """all non-zero integer vectors with entries in [-c, c], primitive, up to sign."""
    out = set()
    for v in itertools.product(range(-c, c + 1), repeat=k):
        if any(v):
            out.add(normalize_vec(v))
    return sorted(out)


def perm_group(ground):
    """permutations of the ground set induced by coordinate permutations only."""
    k = len(ground[0])
    idx = {v: i for i, v in enumerate(ground)}
    perms = set()
    for sigma in itertools.permutations(range(k)):
        perms.add(tuple(idx[tuple(v[sigma[i]] for i in range(k))] for v in ground))
    return np.array(sorted(perms), dtype=np.int64)


def constraint_block_modp(a, d, p):
    a = tuple(int(x) % p for x in a)
    R = restriction_matrix(a, d)
    R = np.array([[x % p for x in row] for row in R], dtype=np.int64)
    return np.hstack([(aj * R) % p for aj in a])


def nullspace_mod(M, ncols, p):
    """basis (list of int lists) of the null space over F_p of the int matrix M."""
    if M is None or len(M) == 0:
        return [[1 if i == j else 0 for i in range(ncols)] for j in range(ncols)]
    R, piv = rref_mod(M, p)
    pivset = set(piv)
    basis = []
    for f in range(ncols):
        if f in pivset:
            continue
        v = [0] * ncols
        v[f] = 1
        for i, c in enumerate(piv):
            v[c] = int((-R[i, f]) % p)
        basis.append(v)
    return basis


def free_with_exponents_modp(vecs, exps, p, cache=None, rng=None):
    """the same criterion as free_with_exponents, over the field F_p (all linear algebra exact in F_p).
    Returns True/False, or None if no point of F_p^k off the arrangement was found (field too small)."""
    vecs = [tuple(int(x) % p for x in v) for v in vecs]
    n, k = len(vecs), len(vecs[0])
    exps = sorted(exps)
    if sum(exps) != n:
        return False
    cache = cache if cache is not None else {}
    rng = rng or random.Random(4321)
    pt = None
    for _ in range(400):
        c = [rng.randrange(p) for _ in range(k)]
        if all(sum(a * b for a, b in zip(v, c)) % p for v in vecs):
            pt = c
            break
    if pt is None:
        return None
    for d in sorted(set(exps)):
        if d == 0:
            B = nullspace_mod(np.array(vecs, dtype=np.int64), k, p)
        else:
            blocks = []
            for a in vecs:
                key = (a, d, p)
                b = cache.get(key)
                if b is None:
                    b = constraint_block_modp(a, d, p)
                    cache[key] = b
                blocks.append(b)
            M = np.vstack(blocks)
            B = nullspace_mod(M, M.shape[1], p)
        if len(B) < required_dim(exps, k, d):
            return False
        mons, _ = monomials(k, d)
        N = len(mons)
        mv = []
        for m in mons:
            x = 1
            for c, e in zip(pt, m):
                x = (x * pow(c, e, p)) % p
            mv.append(x)
        vals = [[sum(w[j * N + t] * mv[t] for t in range(N)) % p for j in range(k)] for w in B]
        if rank_of(vals, p) < sum(1 for e in exps if e <= d):
            return False
    return True


def orbit_representatives_all(perms, size, seed=1):
    """orderly generation as in orbit_representatives, returning the representatives of every level 1..size."""
    N = perms.shape[1]
    rng = np.random.default_rng(seed)
    P = perms[rng.permutation(len(perms))]
    reps = np.zeros((1, 0), dtype=np.int64)
    out = []
    for level in range(1, size + 1):
        parts = []
        for x in range(N):
            sel = reps if level == 1 else reps[reps[:, -1] < x]
            if len(sel):
                parts.append(np.hstack([sel, np.full((len(sel), 1), x, dtype=np.int64)]))
        cand = np.vstack(parts)
        ok = np.zeros(len(cand), dtype=bool)
        chunk = 400000
        for s in range(0, len(cand), chunk):
            ok[s:s + chunk] = is_canonical_batch(cand[s:s + chunk], P, N)
        reps = cand[ok]
        out.append(reps)
    return out
