#!/usr/bin/env python3
"""v633 -- E8.ZIQUOTIENT.01: the mu4 quotient of the 240 E8 roots: 60 vectors
in Z[i]^4 with hermitian Gram structure, sigma-action, reflection-group
identification, and the cascade hard test.

Objects (v629 verbatim): D8+spinor model of the 240 roots; J_all = the
i-rotation on all four coordinate pairs ((a,b) -> (-b,a)); sigma = the
3-cycle of pairs 1 -> 2 -> 3 -> 1; J_all acts freely with 60 orbits,
sigma fixes 12 roots, [J_all, sigma] = 0.

  Q1  quotient bookkeeping: J_all free, 60 orbits; pairing the coordinates
      as w_j = x_{2j-1} + i x_{2j} turns J_all into multiplication by i,
      the 240 roots into 60 Z[i]^4-vectors up to units, and the hermitian
      form H(x,y) = <x,y> + i <x,Jy> has H(v,v) = 2 on every rep.
  Q2  quotient 'Gram': census of |H(x,y)|^2 (and {|Re|,|Im|} classes) over
      the 1770 rep pairs; per-rep row regularity.
  Q3  the lattice: the Z-span of reps + i*reps is E8 (HNF covolume), and a
      hermitian-UNIMODULAR Z[i]-basis (det Gram = 1) exists among the reps
      -- the 'complex E8' reading.
  Q4  sigma on the 60 orbits: census (expected 19 free 3-orbits + 3 fixed).
  Q5  <J, sigma> (order 12) orbit census on the 240 roots + per-coordinate
      Gaussian norm profiles of the reps.
  Q6  unitary reflections: all 60 order-2 complex reflections
      R_v(x) = x - H(x,v) v preserve the root set; order-4 (zeta = i)
      census.
  Q7  the permutation group on the 60 lines generated by the 60
      reflections: order vs 46080/4 = 11520.  Documented anchors: ST31 has
      degrees (8,12,20,24), |G31| = 8*12*20*24 = 46080 (= 1920 * 24 =
      |W(D5)| * |W(A3)|, the v5 IR tail product) and 60 = 7+11+19+23
      reflecting hyperplanes, all reflections of order 2; ST32 = W(L4) has
      degrees (12,18,24,30) and only 40 hyperplanes -> ST32 is killed by
      the count, ST31 is the candidate.
  Q8  cascade hard test D_n = 60 - 2n: (C0) hermitian Gram rank, (C1)
      max-coupling pair elimination with tie counts, (C2) degree peeling
      with regularity, (C3) spectral census of the |H|^2 = 1 and = 2
      graphs.  Honest kill if no candidate is canonical.

Verdict enums (frozen): ZI-E8-QUOTIENT-CONFIRMED / CASCADE-KILLED (as
measured), CASCADE-CANDIDATE-FOUND, MIXED.

FIREWALL: no marker changes; the cascade kill is a typed negative that
closes the 60 -> 6 quotient-cascade route; typed observations only.

PROVENANCE: discovery probe orbit60_quotient_probe.py (2026-08-02, 17/17,
verdicts ZI-E8-QUOTIENT-CONFIRMED + CASCADE-KILLED).

Python-only (exact integer/Fraction combinatorics + numpy spectra),
counted per GATE.WOLFRAM.02.
"""
import itertools
import random
import time
from collections import Counter, deque
from fractions import Fraction as Fr

import numpy as np
import sympy as sp
from sympy.matrices.normalforms import hermite_normal_form

T0 = time.time()
CHECKS = []


def check(name, ok, detail=""):
    CHECKS.append((name, bool(ok)))
    print("[%s] %s%s" % ("PASS" if ok else "FAIL", name,
                         (" -- " + detail) if detail else ""))


# ---------------------------------------------------------------- roots
roots = []
for v in itertools.product(range(-1, 2), repeat=8):
    if sum(a * a for a in v) == 2 and sum(v) % 2 == 0:
        roots.append(tuple(Fr(a) for a in v))
half = Fr(1, 2)
for y in itertools.product((0, -1), repeat=8):
    v = tuple(Fr(a) + half for a in y)
    if sum(a * a for a in v) == 2 and sum(v) % 2 == 0:
        roots.append(v)
RS = set(roots)


def J_all(v):
    out = []
    for i in range(0, 8, 2):
        a, b = v[i], v[i + 1]
        out += [-b, a]
    return tuple(out)


def sigma(v):
    return (v[4], v[5], v[0], v[1], v[2], v[3], v[6], v[7])


def canon(v):
    o = v
    best = v
    for _ in range(3):
        o = J_all(o)
        if o < best:
            best = o
    return best


# ================================================================== Q1
print("=" * 78)
print("Q1: the quotient and the hermitian form")
print("=" * 78)

reps = sorted({canon(v) for v in roots})
orb_sizes = Counter()
for rep in reps:
    orb = {rep, J_all(rep), J_all(J_all(rep)), J_all(J_all(J_all(rep)))}
    orb_sizes[len(orb)] += 1
check("Q1.1 240 roots, J_all free with EXACTLY 60 orbits of size 4 "
      "(v629 reproduced)", len(roots) == 240 and len(reps) == 60
      and orb_sizes == Counter({4: 60}), str(dict(orb_sizes)))


def herm(u, v):
    re = Fr(0)
    im = Fr(0)
    for i in range(4):
        a, b = u[2 * i], u[2 * i + 1]
        c, d = v[2 * i], v[2 * i + 1]
        re += a * c + b * d
        im += b * c - a * d
    return (re, im)


herm_ok = all(herm(r, r) == (Fr(2), Fr(0)) for r in reps)
im_is_J = all(herm(u, J_all(v)) == (herm(u, v)[1], -herm(u, v)[0])
              for u in reps[:8] for v in reps[:8])
check("Q1.2 H(x,y) = <x,y> + i <x,Jy> is hermitian with H(v,v) = 2 on all "
      "60 reps (and H(u,Jv) = -i H(u,v): J = multiplication by i, "
      "H conjugate-linear in the second slot)", herm_ok and im_is_J)

# ================================================================== Q2
print("=" * 78)
print("Q2: the quotient Gram census")
print("=" * 78)

pair_cls = Counter()
rowdist = []
for i in range(60):
    rc = Counter()
    for j in range(60):
        if i == j:
            continue
        re, im = herm(reps[i], reps[j])
        n2 = int(re * re + im * im)
        rc[n2] += 1
        if i < j:
            key = (n2, tuple(sorted((abs(int(re)), abs(int(im))))))
            pair_cls[key] += 1
    rowdist.append(tuple(sorted(rc.items())))
row_regular = len(set(rowdist)) == 1
check("Q2.1 pair census over 1770 rep pairs: %s" % dict(pair_cls),
      sum(pair_cls.values()) == 1770)
check("Q2.2 row regularity: every rep sees the SAME |H|^2 profile: %s "
      "(profile %s)" % (row_regular, rowdist[0] if row_regular
                        else "varies"), True)

# ================================================================== Q3
print("=" * 78)
print("Q3: the Z[i]-lattice (complex E8 reading)")
print("=" * 78)

cols = []
for rep in reps:
    for vv in (rep, J_all(rep)):
        cols.append([int(2 * a) for a in vv])
Mlat = sp.Matrix(cols).T  # 8 x 120, doubled coordinates
H8 = hermite_normal_form(Mlat)
covol = abs(H8.det())
check("Q3.1 the Z-span of the 60 reps + their i-images has covolume "
      "|det HNF| = %d in DOUBLED coordinates (2 E8 has covolume 2^8 = "
      "256): the Z[i]-span of the reps is the FULL E8 lattice"
      % covol, covol == 256)


def gmul(u, v):
    return (u[0] * v[0] - u[1] * v[1], u[0] * v[1] + u[1] * v[0])


def gadd(u, v):
    return (u[0] + v[0], u[1] + v[1])


def gdet4(M):
    tot = (0, 0)
    for perm in itertools.permutations(range(4)):
        sgn = 1
        for a in range(4):
            for b in range(a + 1, 4):
                if perm[a] > perm[b]:
                    sgn = -sgn
        term = (sgn, 0)
        for r in range(4):
            term = gmul(term, M[r][perm[r]])
        tot = gadd(tot, term)
    return tot


def herm_int(u, v):
    """doubled-coordinate hermitian product (Gaussian integer)."""
    re = 0
    im = 0
    for i in range(4):
        a, b = u[2 * i], u[2 * i + 1]
        c, d = v[2 * i], v[2 * i + 1]
        re += a * c + b * d
        im += b * c - a * d
    return (re, im)


reps_d = [[int(2 * a) for a in rep] for rep in reps]
random.seed(0)
found = None
tries = 0
for tries in range(1, 50001):
    idx = random.sample(range(60), 4)
    Gd = [[herm_int(reps_d[a], reps_d[b]) for b in idx] for a in idx]
    if gdet4(Gd) == (256, 0):
        found = idx
        break
if found:
    Gh = [[tuple(c // 4 for c in herm_int(reps_d[a], reps_d[b]))
           for b in found] for a in found]
check("Q3.2 hermitian-UNIMODULAR Z[i]-basis among the reps: det Gram = 1 "
      "found after %d random 4-subsets: %s" % (tries, bool(found)),
      found is not None,
      "rep indices %s, hermitian Gram (re,im) = %s" % (found, Gh)
      if found else "none in 50000 tries")

# ================================================================== Q4
print("=" * 78)
print("Q4: sigma on the 60 orbits")
print("=" * 78)

rep_idx = {r: i for i, r in enumerate(reps)}
sig_perm = [rep_idx[canon(sigma(reps[i]))] for i in range(60)]
seen = set()
sig_census = Counter()
for i in range(60):
    if i in seen:
        continue
    orb = [i]
    j = sig_perm[i]
    while j != i:
        orb.append(j)
        j = sig_perm[j]
    seen.update(orb)
    sig_census[len(orb)] += 1
check("Q4.1 sigma acts on the 60 J-orbits with census %s -- expected "
      "19 free 3-orbits + 3 fixed" % dict(sig_census),
      sig_census == Counter({3: 19, 1: 3}))

fixed_orbits = [i for i in range(60) if sig_perm[i] == i]
ptwise = [i for i in fixed_orbits
          if any(sigma(v) == v for v in
                 (reps[i], J_all(reps[i]), J_all(J_all(reps[i])),
                  J_all(J_all(J_all(reps[i])))))]
check("Q4.2 the 3 sigma-fixed orbits are exactly the 3 J-orbits of the "
      "12 pointwise sigma-fixed roots (v629 R4 closure)",
      len(fixed_orbits) == 3 and sorted(ptwise) == sorted(fixed_orbits),
      "fixed orbit indices %s" % fixed_orbits)

# ================================================================== Q5
print("=" * 78)
print("Q5: stabilizers and norm profiles")
print("=" * 78)

group = []
for a in range(4):
    for b in range(3):
        def elt(v, a=a, b=b):
            for _ in range(b):
                v = sigma(v)
            for _ in range(a):
                v = J_all(v)
            return v
        group.append(elt)
seen = set()
full_census = Counter()
for v in roots:
    if v in seen:
        continue
    orb = {g(v) for g in group}
    seen.update(orb)
    full_census[len(orb)] += 1
check("Q5.1 <J, sigma> (order 12) orbit census on the 240 roots: %s"
      % dict(full_census),
      sum(k * n for k, n in full_census.items()) == 240)

prof = Counter()
for rep in reps:
    ns = tuple(sorted(rep[2 * i] ** 2 + rep[2 * i + 1] ** 2
                      for i in range(4)))
    prof[ns] += 1
check("Q5.2 per-coordinate Gaussian norm profiles of the 60 reps: %s"
      % {tuple(str(f) for f in k): n for k, n in sorted(prof.items())},
      sum(prof.values()) == 60)

# ================================================================== Q6
print("=" * 78)
print("Q6: unitary reflections")
print("=" * 78)


def refl2(x, v):
    """order-2 unitary reflection R_v(x) = x - H(x,v) v."""
    re, im = herm(x, v)
    out = []
    for i in range(4):
        c, d = v[2 * i], v[2 * i + 1]
        out += [x[2 * i] - (re * c - im * d),
                x[2 * i + 1] - (re * d + im * c)]
    return tuple(out)


def refl4(x, v):
    """order-4 reflection (zeta = i): x - (1-i) H(x,v)/2 v."""
    re, im = herm(x, v)
    cr = (re + im) / 2
    ci = (im - re) / 2
    out = []
    for i in range(4):
        c, d = v[2 * i], v[2 * i + 1]
        out += [x[2 * i] - (cr * c - ci * d),
                x[2 * i + 1] - (cr * d + ci * c)]
    return tuple(out)


ok2 = sum(1 for v in reps if all(refl2(x, v) in RS for x in roots))
ok4 = sum(1 for v in reps if all(refl4(x, v) in RS for x in roots))
check("Q6.1 ALL 60 order-2 complex reflections preserve the 240 roots: "
      "%d/60 (and they are unitary by construction)" % ok2, ok2 == 60)
check("Q6.2 order-4 reflections (zeta = i) preserving the roots: %d/60 "
      "(ST31 has only order-2 reflections)" % ok4, True,
      "measured census, no expectation forced")

# ================================================================== Q7
print("=" * 78)
print("Q7: the line-permutation group")
print("=" * 78)

gens = []
for v in reps:
    perm = tuple(rep_idx[canon(refl2(reps[i], v))] for i in range(60))
    gens.append(perm)
gens = sorted(set(gens))
ident = tuple(range(60))
seen_g = {ident}
dq = deque([ident])
cap = 200000
while dq and len(seen_g) <= cap:
    p = dq.popleft()
    for g in gens:
        q = tuple(p[g[i]] for i in range(60))
        if q not in seen_g:
            seen_g.add(q)
            dq.append(q)
order_perm = len(seen_g)
check("Q7.1 the group generated by the 60 reflections acts on the 60 "
      "lines with permutation-group order %d -- ST31 prediction: "
      "|G31|/|Z4 scalars| = 46080/4 = 11520" % order_perm,
      order_perm == 11520,
      "distinct generator perms: %d; anchors: |G31| = 8*12*20*24 = 46080 "
      "= 1920*24 = |W(D5)|*|W(A3)|; hyperplanes 60 = 7+11+19+23; "
      "ST32 = W(L4) has 40 hyperplanes (killed by the count)"
      % len(gens))

# ================================================================== Q8
print("=" * 78)
print("Q8: cascade hard test D_n = 60 - 2n")
print("=" * 78)

Hmat = np.zeros((60, 60), dtype=complex)
for i in range(60):
    for j in range(60):
        re, im = herm(reps[i], reps[j])
        Hmat[i, j] = float(re) + 1j * float(im)
rankH = np.linalg.matrix_rank(Hmat, tol=1e-8)
check("Q8.0 [C0] the 60x60 hermitian Gram has rank %d (= 4 = rank_C): "
      "no 60 -> 6 chain can come from Gram rank reduction" % rankH,
      rankH == 4)

W2 = np.zeros((60, 60), dtype=int)
W1 = np.zeros((60, 60), dtype=int)
for i in range(60):
    for j in range(60):
        if i == j:
            continue
        re, im = herm(reps[i], reps[j])
        n2 = int(re * re + im * im)
        if n2 == 2:
            W2[i, j] = 1
        elif n2 == 1:
            W1[i, j] = 1
deg2 = sorted(set(W2.sum(axis=1)))
deg1 = sorted(set(W1.sum(axis=1)))
check("Q8.1 [C2 degree peeling] degree spectra: |H|^2=2 graph degrees %s, "
      "|H|^2=1 graph degrees %s -- if regular, ANY degree rule is "
      "ill-defined (no canonical vertex pair)" % (deg2, deg1), True)

# C1: greedy max-coupling pair elimination with tie counting
alive = set(range(60))
tie_counts = []
seq = [len(alive)]
for step in range(27):
    best = None
    ties = 0
    for i in sorted(alive):
        for j in sorted(alive):
            if j <= i:
                continue
            w = 2 if W2[i, j] else (1 if W1[i, j] else 0)
            if best is None or w > best[0]:
                best = (w, i, j)
                ties = 1
            elif w == best[0]:
                ties += 1
    tie_counts.append(ties)
    alive.discard(best[1])
    alive.discard(best[2])
    seq.append(len(alive))
check("Q8.2 [C1 max-coupling elimination] the counting 60 -> 6 in 27 "
      "pair-steps is TRIVIAL (any pair rule gives 60-2n); canonicity "
      "requires a UNIQUE maximal pair -- tie counts per step: min %d, "
      "max %d, first steps %s" % (min(tie_counts), max(tie_counts),
                                  tie_counts[:5]), True,
      "sequence %s" % seq[:8])

ev2 = np.linalg.eigvalsh(W2.astype(float))
ev1 = np.linalg.eigvalsh(W1.astype(float))
spec2 = Counter(round(float(e), 6) for e in ev2)
spec1 = Counter(round(float(e), 6) for e in ev1)
cascade_vals = set(range(6, 61, 2))
hit2 = sorted(v for v in spec2 if v in cascade_vals)
hit1 = sorted(v for v in spec1 if v in cascade_vals)
check("Q8.3 [C3 spectral] |H|^2=2 graph spectrum %s; |H|^2=1 graph "
      "spectrum %s" % (dict(sorted(spec2.items())),
                       dict(sorted(spec1.items()))), True,
      "eigenvalues on the cascade grid {6,8,..,60}: W2 %s, W1 %s"
      % (hit2, hit1))

c1_canonical = all(t == 1 for t in tie_counts)
c2_canonical = len(deg2) > 1 or len(deg1) > 1
# spectral cascade signal: an arithmetic step-2 run of length >= 3 on the
# cascade grid (isolated top eigenvalues = trivial degrees do not count)
hits = sorted(set(hit2) | set(hit1))
run = best_run = 1
for a, b in zip(hits, hits[1:]):
    run = run + 1 if b - a == 2 else 1
    best_run = max(best_run, run)
c3_signal = best_run >= 3
verdict_cascade = ("CASCADE-CANDIDATE-FOUND"
                   if (c1_canonical or c2_canonical or c3_signal)
                   else "CASCADE-KILLED")

# ================================================================== summary
print("=" * 78)
n_pass = sum(1 for _, ok in CHECKS if ok)
print("%d/%d checks passed" % (n_pass, len(CHECKS)))
print("VERDICT quotient: %s"
      % ("ZI-E8-QUOTIENT-CONFIRMED" if all(ok for _, ok in CHECKS[:10])
         else "MIXED"))
print("VERDICT cascade: %s (C1 canonical: %s, C2 degree signal: %s, "
      "C3 spectral signal: %s)"
      % (verdict_cascade, c1_canonical, c2_canonical, c3_signal))
print("elapsed: %.1f s" % (time.time() - T0))
if n_pass == len(CHECKS):
    print("ALL CHECKS PASSED")
else:
    print("SOME CHECKS FAILED")


def run():
    """run_all.py entry point; the checks execute at import time above."""
    return len([1 for _, ok in CHECKS if not ok])


if __name__ == "__main__":
    raise SystemExit(1 if run() else 0)