#!/usr/bin/env python3
"""Standalone cross-checker for OMEGA Paper IX Appendix A.3-A.4:
- 120-element local stabilizer Stab_{Aut(G)}(0) via exact Z[phi] inner products
- Orbit partitions of r-subsets of W = {0, 33, ..., 44} for r = 0..7
- Compressed-generator subset counts (5,117,618 for r=6; 45,111,870 and
  2,733,293,660 for r=7 across 30 row-subset representatives)
- Exact projected zonotope minor sums a_0 .. a_7 via Sylvester-Bareiss tree
- Exact permanent tail p_8 .. p_13 and capacity bounds L_7, U (Table 3, Eq. 34-35)

Usage:
  python3 verify_orbit_minor_sums.py            # Pure Python check (a_0..a_5 + counts)
  python3 verify_orbit_minor_sums.py --cpp      # Auto-compiles C++17 __int128 for a_0..a_7
  python3 verify_orbit_minor_sums.py --max-r 6  # Pure Python through a_6
"""

from hashlib import sha256
from itertools import combinations, permutations, product
from math import comb
import argparse
import json
import subprocess
import tempfile
from pathlib import Path

PIN = "40bc6a1d52efab3d75cf343e65ac06b314f9c9c375fc02bf6d45bf6bb21f9529"
DELTA = (2 ** 80) * (5 ** 27) * (11 ** 36) * (13 ** 25) * (41 ** 9) * (211 ** 4)

EXPECTED_A = {
    0: 1,
    1: 325,
    2: 45408,
    3: 3647440,
    4: 188983720,
    5: 6697090056,
    6: 167581390636,
    7: 3007006391716,
}

EXPECTED_ORBITS = {0: 1, 1: 2, 2: 4, 3: 8, 4: 15, 5: 22, 6: 30, 7: 30}

EXPECTED_PERM_TAIL = {
    8: 82326531807696,
    9: 879734614473680,
    10: 6510322192690880,
    11: 31533042290406720,
    12: 89385221732384160,
    13: 111664804452677280,
}


def require(ok: bool, message: str) -> None:
    if not ok:
        raise ArithmeticError(f"Check failed: {message}")


def build_600cell():
    """Reconstructs the 120 vertices, adjacency lists, and matrix M = 13I + A."""
    vertices = set()
    for i in range(4):
        for s in (-1, 1):
            v = [(0, 0)] * 4
            v[i] = (2 * s, 0)
            vertices.add(tuple(v))
    for signs in product((-1, 1), repeat=4):
        vertices.add(tuple((s, 0) for s in signs))
    even = [
        p for p in permutations(range(4))
        if sum(p[i] > p[j] for i in range(4) for j in range(i + 1, 4)) % 2 == 0
    ]
    for a, b, c in product((-1, 1), repeat=3):
        v = ((0, 0), (a, 0), (0, b), (-c, c))
        vertices.update(tuple(v[i] for i in p) for p in even)
    vertices = sorted(vertices)
    neighbors = [[] for _ in vertices]
    for i, j in combinations(range(len(vertices)), 2):
        x = y = 0
        for (a, b), (c, d) in zip(vertices[i], vertices[j]):
            u, v = a - c, b - d
            x += u * u + v * v
            y += 2 * u * v + v * v
        if (x, y) == (8, -4):
            neighbors[i].append(j)
            neighbors[j].append(i)

    fixture = {
        "coordinate_rule": "vertex coordinate is (a+b*phi)/2 for pair [a,b]",
        "neighbors": neighbors,
        "schema": "OMEGA_600CELL_GRAPH_FIXTURE/1",
        "vertex_order": "Python lexicographic order of tuples of four integer pairs",
        "vertices": vertices,
    }
    raw = (json.dumps(fixture, sort_keys=True, separators=(",", ":")) + "\n").encode()
    require(sha256(raw).hexdigest() == PIN, "graph SHA-256 fixture identity")
    require(len(vertices) == 120 and all(len(nb) == 12 for nb in neighbors), "12-regularity")

    M = [
        [13 if i == j else int(j in neighbors[i]) for j in range(120)]
        for i in range(120)
    ]
    W = [0] + list(range(33, 45))
    require(sorted([0] + neighbors[0]) == W, "W == {0, 33, ..., 44}")
    return vertices, neighbors, M, W


def zphi_dot4(u, v):
    """Exact 4 * <u, v> in Z[phi] as pair (x, y) representing x + y*phi."""
    x = y = 0
    for (a, b), (c, d) in zip(u, v):
        x += a * c + b * d
        y += a * d + b * c + b * d
    return (x, y)


def build_local_stabilizer(vertices, neighbors, W):
    """Constructs the 120-element stabilizer of vertex 0 on W and verifies that
    every element lifts to an exact automorphism of the full 120-vertex graph.
    """
    N0 = W[1:]
    adj_N0 = {u: set(neighbors[u]).intersection(N0) for u in N0}
    require(all(len(adj_N0[u]) == 5 for u in N0), "icosahedral neighborhood degree 5")

    u0 = N0[0]
    u1 = min(adj_N0[u0])
    u2 = min(adj_N0[u0].intersection(adj_N0[u1]))

    dot_table = [[zphi_dot4(vertices[v], vertices[w]) for w in W] for v in range(120)]
    sig_to_vertex = {tuple(dot_table[v]): v for v in range(120)}
    require(len(sig_to_vertex) == 120, "W spans R^4 so inner-product signatures are unique")

    adj_set = [set(nb) for nb in neighbors]
    w_index = {w: idx for idx, w in enumerate(W)}

    stab_on_W = []
    for v0 in N0:
        for v1 in adj_N0[v0]:
            for v2 in adj_N0[v0].intersection(adj_N0[v1]):
                basis = (0, u0, u1, u2)
                target_basis = (0, v0, v1, v2)
                basis_sigs = {
                    tuple(zphi_dot4(vertices[w], vertices[b]) for b in target_basis): w
                    for w in W
                }
                perm_W_vertices = []
                for w in W:
                    sig = tuple(zphi_dot4(vertices[w], vertices[b]) for b in basis)
                    if sig in basis_sigs:
                        perm_W_vertices.append(basis_sigs[sig])
                if len(set(perm_W_vertices)) != 13:
                    continue

                full_perm = []
                for v in range(120):
                    inv_sig = [None] * 13
                    for k, w_img in enumerate(perm_W_vertices):
                        inv_sig[w_index[w_img]] = dot_table[v][k]
                    v_img = sig_to_vertex.get(tuple(inv_sig))
                    if v_img is None:
                        break
                    full_perm.append(v_img)
                require(len(set(full_perm)) == 120, "bijective lift to 120 vertices")

                for i in range(120):
                    mapped_nb = {full_perm[j] for j in neighbors[i]}
                    require(mapped_nb == adj_set[full_perm[i]], "exact graph automorphism")

                stab_on_W.append(tuple(w_index[w_img] for w_img in perm_W_vertices))

    stab_on_W = sorted(set(stab_on_W))
    require(len(stab_on_W) == 120, "local stabilizer order == 120")
    return stab_on_W


def partition_row_orbits(stab_on_W, r: int):
    """Partitions all r-subsets of {0..12} into orbits under stab_on_W."""
    unvisited = set(combinations(range(13), r))
    orbits = []
    while unvisited:
        rep = min(unvisited)
        orb = {tuple(sorted(g[i] for i in rep)) for g in stab_on_W}
        require(orb.issubset(unvisited), "orbit closure")
        unvisited -= orb
        orbits.append((rep, len(orb)))
    require(len(orbits) == EXPECTED_ORBITS[r], f"orbit count for r={r}")
    require(sum(sz for _, sz in orbits) == comb(13, r), f"orbit size sum for r={r}")
    return orbits


def compress_columns(M, W, rep_local_indices):
    """Groups identical nonzero columns into (vector, multiplicity) pairs."""
    rows = [W[idx] for idx in rep_local_indices]
    r = len(rows)
    counts = {}
    for j in range(120):
        col = tuple(M[row][j] for row in rows)
        if any(col):
            counts[col] = counts.get(col, 0) + 1
    items = sorted(counts.items())
    vecs = [list(v) for v, _ in items]
    mults = [mu for _, mu in items]
    require(all(mu <= 98 for mu in mults), "generator multiplicity <= 98")
    return vecs, mults


def minor_sum_sylvester_bareiss(vecs, mults, r: int) -> int:
    """Weighted absolute minor sum using incremental fraction-free elimination."""
    m = len(vecs)
    if r == 0:
        return 1
    if m < r:
        return 0
    if r == 1:
        return sum(mults[i] * abs(vecs[i][0]) for i in range(m))
    if r == 2:
        total = 0
        for a in range(m - 1):
            ma, (xa, ya) = mults[a], vecs[a]
            for b in range(a + 1, m):
                det2 = xa * vecs[b][1] - ya * vecs[b][0]
                if det2:
                    total += ma * mults[b] * abs(det2)
        return total

    def walk(depth: int, candidates, prev_piv: int, weight_acc: int) -> int:
        k = r - depth
        need = k
        n_cand = len(candidates)
        if n_cand < need:
            return 0
        subtotal = 0
        for pos in range(n_cand - need + 1):
            mu_p, pvec = candidates[pos]
            piv_col = -1
            for c in range(k):
                if pvec[c] != 0:
                    piv_col = c
                    break
            if piv_col < 0:
                continue
            piv = pvec[piv_col]
            rem_cols = [c for c in range(k) if c != piv_col]

            if depth == r - 3:
                c0, c1 = rem_cols[0], rem_cols[1]
                p0, p1 = pvec[c0], pvec[c1]
                red2 = []
                for pos2 in range(pos + 1, n_cand):
                    mu2, v2 = candidates[pos2]
                    vc = v2[piv_col]
                    x2 = (v2[c0] * piv - vc * p0) // prev_piv
                    y2 = (v2[c1] * piv - vc * p1) // prev_piv
                    if x2 or y2:
                        red2.append((mu2, x2, y2))
                n_red = len(red2)
                if n_red >= 2:
                    pair_sum = 0
                    for a in range(n_red - 1):
                        mu_a, xa, ya = red2[a]
                        for b in range(a + 1, n_red):
                            mu_b, xb, yb = red2[b]
                            cross = xa * yb - ya * xb
                            if cross:
                                pair_sum += mu_a * mu_b * abs(cross)
                    q, rem = divmod(pair_sum, abs(piv))
                    require(rem == 0, "exact Sylvester 2x2 division")
                    subtotal += (weight_acc * mu_p) * q
            else:
                next_candidates = []
                for pos2 in range(pos + 1, n_cand):
                    mu2, v2 = candidates[pos2]
                    vc = v2[piv_col]
                    new_v = [(v2[c] * piv - vc * pvec[c]) // prev_piv for c in rem_cols]
                    if any(new_v):
                        next_candidates.append((mu2, new_v))
                if len(next_candidates) >= k - 1:
                    subtotal += walk(depth + 1, next_candidates, piv, weight_acc * mu_p)
        return subtotal

    initial_candidates = list(zip(mults, vecs))
    return walk(0, initial_candidates, 1, 1)


CPP_KERNEL = r"""
#include <iostream>
#include <vector>
#include <cmath>
#include <cstdint>

using i128 = __int128_t;
using u128 = __uint128_t;

static inline i128 iabs128(i128 x) { return x < 0 ? -x : x; }

struct Cand {
    i128 mu;
    i128 v[7];
};

int R_GLOB;

i128 walk(int depth, const std::vector<Cand>& cands, i128 prev_piv, i128 weight_acc) {
    int k = R_GLOB - depth;
    int n = (int)cands.size();
    if (n < k) return 0;
    i128 subtotal = 0;

    for (int pos = 0; pos <= n - k; ++pos) {
        const Cand& cp = cands[pos];
        int piv_col = -1;
        for (int c = 0; c < k; ++c) {
            if (cp.v[c] != 0) { piv_col = c; break; }
        }
        if (piv_col < 0) continue;
        i128 piv = cp.v[piv_col];
        int rem_cols[7], rc = 0;
        for (int c = 0; c < k; ++c) if (c != piv_col) rem_cols[rc++] = c;

        if (depth == R_GLOB - 3) {
            int c0 = rem_cols[0], c1 = rem_cols[1];
            i128 p0 = cp.v[c0], p1 = cp.v[c1];
            struct Vec2 { i128 mu, x, y; };
            std::vector<Vec2> red2;
            red2.reserve(n - pos - 1);
            for (int pos2 = pos + 1; pos2 < n; ++pos2) {
                const Cand& c2 = cands[pos2];
                i128 vc = c2.v[piv_col];
                i128 x2 = (c2.v[c0] * piv - vc * p0) / prev_piv;
                i128 y2 = (c2.v[c1] * piv - vc * p1) / prev_piv;
                if (x2 != 0 || y2 != 0) red2.push_back({c2.mu, x2, y2});
            }
            int nr = (int)red2.size();
            if (nr >= 2) {
                i128 pair_sum = 0;
                for (int a = 0; a < nr - 1; ++a) {
                    i128 mua = red2[a].mu, xa = red2[a].x, ya = red2[a].y;
                    for (int b = a + 1; b < nr; ++b) {
                        i128 cross = xa * red2[b].y - ya * red2[b].x;
                        if (cross != 0) pair_sum += mua * red2[b].mu * iabs128(cross);
                    }
                }
                subtotal += (weight_acc * cp.mu) * (pair_sum / iabs128(piv));
            }
        } else {
            std::vector<Cand> next_cands;
            next_cands.reserve(n - pos - 1);
            for (int pos2 = pos + 1; pos2 < n; ++pos2) {
                const Cand& c2 = cands[pos2];
                i128 vc = c2.v[piv_col];
                Cand nc;
                nc.mu = c2.mu;
                bool nz = false;
                for (int j = 0; j < k - 1; ++j) {
                    int col = rem_cols[j];
                    nc.v[j] = (c2.v[col] * piv - vc * cp.v[col]) / prev_piv;
                    if (nc.v[j] != 0) nz = true;
                }
                if (nz) next_cands.push_back(nc);
            }
            if ((int)next_cands.size() >= k - 1) {
                subtotal += walk(depth + 1, next_cands, piv, weight_acc * cp.mu);
            }
        }
    }
    return subtotal;
}

void print_i128(i128 x) {
    if (x == 0) { std::cout << "0\n"; return; }
    std::string s;
    while (x > 0) { s.push_back('0' + (int)(x % 10)); x /= 10; }
    for (int i = (int)s.size() - 1; i >= 0; --i) std::cout << s[i];
    std::cout << "\n";
}

int main() {
    int num_tasks;
    if (!(std::cin >> num_tasks)) return 0;
    for (int t = 0; t < num_tasks; ++t) {
        int r, m, mult_orb;
        std::cin >> r >> m >> mult_orb;
        R_GLOB = r;
        std::vector<Cand> cands(m);
        for (int i = 0; i < m; ++i) {
            long long mu;
            std::cin >> mu;
            cands[i].mu = mu;
            for (int j = 0; j < r; ++j) {
                long long val;
                std::cin >> val;
                cands[i].v[j] = val;
            }
        }
        i128 val = walk(0, cands, 1, 1) * (i128)mult_orb;
        print_i128(val);
    }
    return 0;
}
"""


def evaluate_via_cpp(tasks):
    """Compiles the embedded C++17 __int128 kernel and evaluates orbit tasks."""
    with tempfile.TemporaryDirectory() as tmp:
        src = Path(tmp) / "minor_enum.cpp"
        exe = Path(tmp) / "minor_enum"
        src.write_text(CPP_KERNEL)
        subprocess.run(
            ["g++", "-O3", "-std=c++17", str(src), "-o", str(exe)],
            check=True,
        )
        lines = [str(len(tasks))]
        for r, vecs, mults, orb_sz in tasks:
            lines.append(f"{r} {len(vecs)} {orb_sz}")
            for mu, v in zip(mults, vecs):
                lines.append(f"{mu} " + " ".join(str(x) for x in v))
        proc = subprocess.run(
            [str(exe)],
            input="\n".join(lines) + "\n",
            text=True,
            capture_output=True,
            check=True,
        )
        return [int(line.strip()) for line in proc.stdout.strip().splitlines()]


def verify_permanent_tail_and_bounds(M, W):
    """Verifies the permanent dynamic program and capacity arithmetic."""
    dp = [1] + [0] * ((1 << 13) - 1)
    for j in range(120):
        nz = [(1 << r, M[W[r]][j]) for r in range(13) if M[W[r]][j]]
        for mask in range((1 << 13) - 1, -1, -1):
            base = dp[mask]
            if base:
                for bit, weight in nz:
                    if not (mask & bit):
                        dp[mask | bit] += base * weight
    perm = [0] * 14
    for mask, val in enumerate(dp):
        perm[mask.bit_count()] += val

    for r in range(8, 14):
        require(perm[r] == EXPECTED_PERM_TAIL[r], f"permanent p_{r}")

    L7 = sum(EXPECTED_A[r] * (73 ** (13 - r)) for r in range(8))
    U = L7 + sum(perm[r] * (73 ** (13 - r)) for r in range(8, 14))
    require(L7 == 57869157529977646724169246, "exact L_7")
    require(U == 58067516542393048415542654, "exact U")
    require(U ** 120 < (2 ** 352 * DELTA) ** 13, "352-bit upper bound U^120 < (2^352 * Delta)^13")
    require(L7 ** 120 > (2 ** 351 * DELTA) ** 13, "351-bit lower barrier L7^120 > (2^351 * Delta)^13")
    return perm, L7, U


def main():
    parser = argparse.ArgumentParser()
    parser.add_argument(
        "--max-r", type=int, default=5,
        help="Maximum degree r to evaluate in pure Python (default: 5; use 7 for all)"
    )
    parser.add_argument(
        "--cpp", action="store_true",
        help="Compile embedded C++17 __int128 kernel to evaluate r=6 and r=7 rapidly"
    )
    args = parser.parse_args()

    vertices, neighbors, M, W = build_600cell()
    stab_on_W = build_local_stabilizer(vertices, neighbors, W)

    orbit_summary = {}
    computed_a = {}

    for r in range(8):
        orbits = partition_row_orbits(stab_on_W, r)
        rep_subsets = 0
        equiv_subsets = 0
        tasks = []
        for rep, orb_sz in orbits:
            vecs, mults = compress_columns(M, W, rep)
            sub_cnt = comb(len(vecs), r) if r > 0 else 1
            rep_subsets += sub_cnt
            equiv_subsets += orb_sz * sub_cnt
            tasks.append((r, vecs, mults, orb_sz))

        if r == 6:
            require(rep_subsets == 4800238, "r=6 representative subsets == 5,117,618")
        if r == 7:
            require(len(orbits) == 30, "r=7 orbit count == 30")
            require(rep_subsets == 40449768, "r=7 representative subsets == 45,111,870")
            require(equiv_subsets == 2474330048, "r=7 equivalent subsets == 2,733,293,660")

        orbit_summary[r] = {
            "orbits": len(orbits),
            "representative_subsets": rep_subsets,
            "equivalent_subsets": equiv_subsets,
        }

        if r <= args.max_r:
            total_r = sum(
                orb_sz * minor_sum_sylvester_bareiss(vecs, mults, r)
                for _, vecs, mults, orb_sz in tasks
            )
            require(total_r == EXPECTED_A[r], f"exact a_{r} match")
            computed_a[r] = total_r
        elif args.cpp and r >= 6:
            total_r = sum(evaluate_via_cpp(tasks))
            require(total_r == EXPECTED_A[r], f"exact a_{r} C++ match")
            computed_a[r] = total_r

    perm, L7, U = verify_permanent_tail_and_bounds(M, W)

    report = {
        "schema": "OMEGA_PAPER_IX_ORBIT_MINOR_CERTIFICATE/1",
        "status": "PASS",
        "graph_sha256": PIN,
        "local_stabilizer_order": len(stab_on_W),
        "orbit_and_subset_counts": orbit_summary,
        "verified_a_r": computed_a,
        "verified_p_8_to_13": {r: perm[r] for r in range(8, 14)},
        "L7": str(L7),
        "U": str(U),
    }
    print(json.dumps(report, indent=2))


if __name__ == "__main__":
    main()
