#!/usr/bin/env python3
"""verify_lattice.py -- stand-alone checker for the lattice data in this directory.

Python 3 standard library only (exact integer / rational arithmetic throughout); if PARI/GP (`gp`) is on the
PATH it is used in addition to recompute automorphism-group orders. Run in this directory:

    python3 verify_lattice.py            # full run (all 43 charge-bound certificates; about three minutes on one core)
    python3 verify_lattice.py --quick    # the same checks with only the five Table-3 certificates (about one minute)
    python3 verify_lattice.py --no-gp    # do not call gp even if present

Exit status 0 if and only if every check passes. What is checked:

 rank19_lattices.json   for W512, W512', W1024, Lambda_19: symmetry, evenness, determinant, Smith invariants of the
                        discriminant group, minimum 4 and root-freeness and the kissing number by exhaustive exact
                        enumeration of all vectors of norm <= 4, |Aut| (gp only); Lambda_19: T G0 T^T = G, det T = +-1;
                        for each gluing: G even with det -1 and signature (3,19), X G X^T = -W, K G K^T = C, X G K^T = 0,
                        |det[X;K]| = glue index, det W det C = index^2 (so the overlattice is unimodular).
 genera54.json          54 genera, 98 complement classes, 1764 certificates; per genus: the ternary witness has the stated
                        determinant and discriminant group; per class: det, discriminant group, root count s_K, number of
                        dual shells, |Aut| (gp only); sum over classes of 1/|Aut| equals the stated mass of the complement
                        genus; every Farkas certificate evaluates to a strictly negative rational (exact).
 floor_certificates/    every certificate: B equals the named Gram matrix, beta^2 < 12, negative weights sit only on unit
                        vectors (and never on the released column), R is positive definite (exact LDL^T), the value V
                        recomputed from the certificate equals value_frac exactly, and V clears its threshold (48 for
                        row 1, rows 4-5 and all released-column certificates; 46 for rows 2-3); the 156 local cases: the
                        congruences q(M) = 48, det M = d to the stated p-adic precision and the spanning condition;
                        the released-pair coverage counts 55/171 and 40/171.
 rootfree_det160_224.json  for the seven lattices: det, evenness, minimum 4, root count 0, kissing number, Smith
                        invariants, |Aut| (gp only); every isometry U^T B U = A with det U = +-1.
"""
import json, gzip, os, sys, math, time, shutil, subprocess
from fractions import Fraction as Fr

HERE = os.path.dirname(os.path.abspath(__file__))
ARGS = set(sys.argv[1:])
QUICK = "--quick" in ARGS
USE_GP = ("--no-gp" not in ARGS) and shutil.which("gp") is not None
FAILS = []; NPASS = 0
def check(name, ok, detail=""):
    global NPASS
    if ok: NPASS += 1
    else: FAILS.append(name)
    print(("PASS  " if ok else "FAIL  ") + name + (("  [" + str(detail) + "]") if detail != "" else ""), flush=True)

# ---------------------------------------------------------------- exact linear algebra
def T(M): return [list(r) for r in zip(*M)]
def mul(A, B):
    Bt = T(B); return [[sum(a * b for a, b in zip(r, c)) for c in Bt] for r in A]
def is_sym(M): return all(M[i][j] == M[j][i] for i in range(len(M)) for j in range(i))
def det_int(M):
    """Bareiss fraction-free determinant of an integer matrix."""
    A = [list(map(int, r)) for r in M]; n = len(A); sign = 1; prev = 1
    for k in range(n - 1):
        if A[k][k] == 0:
            piv = next((r for r in range(k + 1, n) if A[r][k] != 0), None)
            if piv is None: return 0
            A[k], A[piv] = A[piv], A[k]; sign = -sign
        for i in range(k + 1, n):
            for j in range(k + 1, n):
                A[i][j] = (A[i][j] * A[k][k] - A[i][k] * A[k][j]) // prev
            A[i][k] = 0
        prev = A[k][k]
    return sign * A[n - 1][n - 1]
def ldl(Q):
    """exact LDL^T of a symmetric matrix of Fractions; returns pivots d and unit upper factor mu (Q = U^T D U). Stops at a non-positive pivot."""
    n = len(Q); A = [[Fr(x) for x in r] for r in Q]
    d = [None] * n; mu = [[Fr(0)] * n for _ in range(n)]
    for i in range(n):
        d[i] = A[i][i] - sum(d[k] * mu[k][i] * mu[k][i] for k in range(i))
        if d[i] <= 0: return d[:i + 1], None
        for j in range(i + 1, n):
            mu[i][j] = (A[i][j] - sum(d[k] * mu[k][i] * mu[k][j] for k in range(i))) / d[i]
    return d, mu
def inertia(Q):
    """(n+, n-, n0) of a rational symmetric matrix by symmetric Gaussian elimination with 1x1/2x2 pivots avoided via diagonal shifts:
    simple exact congruence diagonalization."""
    n = len(Q); A = [[Fr(x) for x in r] for r in Q]; pos = neg = zero = 0; active = list(range(n))
    while active:
        i = next((k for k in active if A[k][k] != 0), None)
        if i is None:
            # all remaining diagonals zero: find off-diagonal nonzero (i,j), replace row/col i by i+j
            pair = next(((a, b) for a in active for b in active if a < b and A[a][b] != 0), None)
            if pair is None:
                zero += len(active); break
            a, b = pair
            for k in range(n): A[a][k] += A[b][k]
            for k in range(n): A[k][a] += A[k][b]
            continue
        p = A[i][i]
        if p > 0: pos += 1
        else: neg += 1
        row = A[i][:]
        for r in active:
            if r == i or A[r][i] == 0: continue
            f = A[r][i] / p
            for c in active: A[r][c] -= f * row[c]
        for c in active: A[i][c] = Fr(0)
        for r in active: A[r][i] = Fr(0)
        active.remove(i)
    return (pos, neg, zero)
def smith_invariants(M):
    """invariant factors > 1 of an integer matrix (abelian group Z^n / row space)."""
    A = [list(map(int, r)) for r in M]; n = len(A); m = len(A[0]); inv = []
    r0 = 0
    for c0 in range(min(n, m)):
        # bring smallest nonzero to (r0,r0) repeatedly until it divides its row and column
        while True:
            piv = None
            for i in range(r0, n):
                for j in range(r0, m):
                    if A[i][j] != 0 and (piv is None or abs(A[i][j]) < abs(A[piv[0]][piv[1]])): piv = (i, j)
            if piv is None:
                return sorted([x for x in inv if x > 1])
            i, j = piv
            A[r0], A[i] = A[i], A[r0]
            for row in A: row[r0], row[j] = row[j], row[r0]
            p = A[r0][r0]; done = True
            for i in range(r0 + 1, n):
                q = A[i][r0] // p
                if q:
                    for j in range(r0, m): A[i][j] -= q * A[r0][j]
                if A[i][r0] != 0: done = False
            for j in range(r0 + 1, m):
                q = A[r0][j] // p
                if q:
                    for i in range(r0, n): A[i][j] -= q * A[i][r0]
                if A[r0][j] != 0: done = False
            if done:
                # ensure p divides all remaining entries (else add a row and repeat)
                bad = next(((i, j) for i in range(r0 + 1, n) for j in range(r0 + 1, m) if A[i][j] % p != 0), None)
                if bad is None:
                    inv.append(abs(p)); r0 += 1; break
                for j in range(r0, m): A[r0][j] += A[bad[0]][j]
    return sorted([x for x in inv if x > 1])
def inv_frac(M):
    n = len(M); A = [[Fr(x) for x in r] + [Fr(int(i == j)) for j in range(n)] for i, r in enumerate(M)]
    for c in range(n):
        p = next(r for r in range(c, n) if A[r][c] != 0); A[c], A[p] = A[p], A[c]
        piv = A[c][c]; A[c] = [x / piv for x in A[c]]
        for r in range(n):
            if r != c and A[r][c] != 0:
                f = A[r][c]; A[r] = [a - f * b for a, b in zip(A[r], A[c])]
    return [r[n:] for r in A]
def isqrt_fr(x):
    if x < 0: return -1
    v = math.isqrt(x.numerator // x.denominator)
    while (v + 1) * (v + 1) <= x: v += 1
    while v * v > x: v -= 1
    return v
def short_vectors(Q, bound, strict=False):
    """all nonzero integer x with x^T Q x <= bound (< bound if strict), both signs; Q positive definite (int or Fraction). Exact."""
    n = len(Q); d, mu = ldl(Q); assert mu is not None, "not positive definite"
    bound = Fr(bound); out = []; x = [0] * n
    def rec(i, rem):
        c = sum(mu[i][j] * x[j] for j in range(i + 1, n))
        r = isqrt_fr(rem / d[i]) + 1
        for xi in range(math.floor(-c - r), math.ceil(-c + r) + 1):
            t = d[i] * (xi + c) ** 2
            if t > rem: continue
            x[i] = xi
            if i == 0:
                if any(x):
                    nrm = bound - (rem - t)
                    if (not strict) or nrm < bound: out.append((list(x), nrm))
            else:
                rec(i - 1, rem - t)
        x[i] = 0
    rec(n - 1, bound)
    return out
def norm(Q, v): return sum(Q[i][j] * v[i] * v[j] for i in range(len(v)) for j in range(len(v)))
GP_OK = USE_GP
def gp_aut_order(G):
    global GP_OK
    if not GP_OK: return None
    s = "[" + ";".join(",".join(str(x) for x in r) for r in G) + "]"
    try:
        r = subprocess.run(["gp", "-q", "--default", "parisizemax=1G"], input=f"print(qfauto({s})[1]);\nquit\n", capture_output=True, text=True, timeout=600)
        return int(r.stdout.strip().splitlines()[-1])
    except Exception as e:
        print("note: gp call failed (" + str(e)[:60] + "); |Aut| checks skipped from here on"); GP_OK = False; return None

def load(rel):
    p = os.path.join(HERE, rel)
    if rel.endswith(".gz"):
        with gzip.open(p, "rt") as f: return json.load(f)
    return json.load(open(p))

t0 = time.time()
print("verify_lattice.py  mode=%s  gp=%s" % ("quick" if QUICK else "full", "yes" if USE_GP else "no"))
# ================================================================ 1. rank-19 lattices and gluings
D = load("rank19_lattices.json")
GRAM = {}
for name, d in D["lattices"].items():
    G = d["gram"]; GRAM[name] = G; n = len(G)
    check(f"{name}: {n}x{n} symmetric integer Gram, rank {d['rank']}", n == d["rank"] and all(len(r) == n for r in G) and is_sym(G))
    check(f"{name}: even", all(G[i][i] % 2 == 0 for i in range(n)))
    check(f"{name}: det = {d['det']}", det_int(G) == d["det"])
    check(f"{name}: discriminant group invariants {d['discriminant_group_invariants']}", smith_invariants(G) == sorted(d["discriminant_group_invariants"]))
    sv = short_vectors(G, 4)
    n2 = sum(1 for v, nn in sv if nn == 2); n4 = sum(1 for v, nn in sv if nn == 4); nother = len(sv) - n2 - n4
    check(f"{name}: minimum 4, root-free (vectors of norm 2: {n2}; of norm 1 or 3: {nother})", n2 == 0 and nother == 0 and d["minimum"] == 4 and d["n_roots"] == 0)
    check(f"{name}: kissing number {d['kissing_number']} (vectors of norm 4, both signs)", n4 == d["kissing_number"], n4)
    a = gp_aut_order(G)
    if a is not None: check(f"{name}: |Aut| = {d['aut_order']} (gp qfauto)", a == d["aut_order"], a)
d = D["lattices"]["Lambda19"]
check("Lambda19: T * G_construction * T^T = G_LLL and det T = %d" % d["det_T"], mul(mul(d["T"], d["gram_construction_basis"]), T(d["T"])) == d["gram"] and det_int(d["T"]) == d["det_T"] and abs(d["det_T"]) == 1)
for name, g in D["gluings"].items():
    G = g["overlattice_gram"]; X = g["embedding_rows"]; K = g["complement_basis_rows"]; C = g["complement_gram"]; W = GRAM[name]
    check(f"gluing {name}: overlattice Gram 22x22 symmetric, even, det {g['overlattice_det']}", len(G) == 22 and is_sym(G) and all(G[i][i] % 2 == 0 for i in range(22)) and det_int(G) == g["overlattice_det"] == -1)
    check(f"gluing {name}: signature (3,19)", inertia(G) == (3, 19, 0) and tuple(g["overlattice_inertia"]) == (3, 19, 0))
    check(f"gluing {name}: X G X^T = -{name}", mul(mul(X, G), T(X)) == [[-x for x in r] for r in W])
    check(f"gluing {name}: K G K^T = diag{tuple(C[i][i] for i in range(3))} and X G K^T = 0", mul(mul(K, G), T(K)) == C and all(x == 0 for r in mul(mul(X, G), T(K)) for x in r))
    idx = abs(det_int(X + K))
    check(f"gluing {name}: glue index |det[X;K]| = {g['glue_index']}, det W * det C = index^2 (unimodular overlattice)", idx == g["glue_index"] and det_int(W) * det_int(C) == idx * idx)

# ================================================================ 2. genera
GE = load("genera54.json"); H = GE["coxeter_numbers_H"]
IJ = [(i, j) for i in range(5) for j in range(i, 5)]
check("genera: 54 genera, 98 complement classes, 1764 certificates, 18 Coxeter numbers",
      len(GE["genera"]) == GE["n_genera"] == 54 and sum(len(g["complement_genus"]["classes"]) for g in GE["genera"]) == GE["n_complement_classes"] == 98
      and sum(len(c["farkas_certificates"]) for g in GE["genera"] for c in g["complement_genus"]["classes"]) == GE["n_certificates"] == 1764 and len(H) == 18)
from collections import Counter
check("genera: per-determinant counts consistent", {str(k): v for k, v in Counter(g["det"] for g in GE["genera"]).items()} == GE["per_det_count"])
ncert_ok = 0; ncert = 0; weakest = None; mass_ok = 0; cls_ok = 0; aut_gp_ok = 0; aut_gp_n = 0; tern_ok = 0
for g in GE["genera"]:
    a, b, c, d_, e, f = map(int, g["ternary_witness"].split())
    Tm = [[2 * a, f, e], [f, 2 * b, d_], [e, d_, 2 * c]]
    grp = sorted(int(x) for x in g["discriminant_group"].split("x") if int(x) > 1)
    if det_int(Tm) == g["det"] and smith_invariants(Tm) == grp and ldl(Tm)[1] is not None: tern_ok += 1
    else: print("  ternary witness mismatch in genus", g["n"])
    s = Fr(0)
    for ci, cl in enumerate(g["complement_genus"]["classes"], 1):
        GK = cl["gram"]; s += Fr(1, cl["aut_order"])
        ok = is_sym(GK) and all(GK[i][i] % 2 == 0 for i in range(5)) and det_int(GK) == g["det"] and smith_invariants(GK) == grp
        roots = [v for v, nn in short_vectors(GK, 2) if nn == 2]
        sK = len(roots); ok = ok and sK == cl["n_roots"]
        S = [[0] * 5 for _ in range(5)]
        for v in roots:
            u = [sum(GK[i][j] * v[j] for j in range(5)) for i in range(5)]
            for i in range(5):
                for j in range(5): S[i][j] += u[i] * u[j]
        Ginv = inv_frac(GK)
        shells = []
        for u, t in short_vectors(Ginv, 2, strict=True):
            if next(x for x in u if x != 0) < 0: continue
            cap = 2 * ((Fr(2) / t).numerator // (Fr(2) / t).denominator)
            shells.append((u, cap))
        ok = ok and len(shells) == cl["n_dual_shells"]
        if ok: cls_ok += 1
        else: print("  class data mismatch genus", g["n"], "class", ci)
        ag = gp_aut_order(GK)
        if ag is not None:
            aut_gp_n += 1; aut_gp_ok += (ag == cl["aut_order"])
        for h in H:
            y = [Fr(x) for x in cl["farkas_certificates"][str(h)]]; ncert += 1
            bvec = [2 * h * GK[i][j] - S[i][j] for (i, j) in IJ] + [24 * h - sK]
            val = sum(yi * bi for yi, bi in zip(y, bvec))
            for (u, cap) in shells:
                col = sum(y[r] * u[i] * u[j] for r, (i, j) in enumerate(IJ)) + y[15]
                if col < 0: val += cap * (-col)
            if val < 0:
                ncert_ok += 1; weakest = val if (weakest is None or val > weakest) else weakest
            else: print("  FARKAS CERTIFICATE NOT NEGATIVE: genus", g["n"], "class", ci, "h", h, "value", val)
    if s == Fr(g["complement_genus"]["mass"]): mass_ok += 1
    else: print("  mass mismatch genus", g["n"], s, g["complement_genus"]["mass"])
check("genera: every ternary witness has the stated determinant and discriminant group and is positive definite", tern_ok == 54, tern_ok)
check("genera: every complement class has the stated det, discriminant group, root count and dual-shell count", cls_ok == 98, cls_ok)
check("genera: sum of 1/|Aut| over the classes equals the mass of the complement genus, all 54 genera", mass_ok == 54, mass_ok)
if aut_gp_n: check(f"genera: |Aut(K)| recomputed by gp for {aut_gp_n} classes", aut_gp_ok == aut_gp_n, aut_gp_ok)
check(f"genera: all {ncert} Farkas certificates strictly negative ({ncert_ok} negative; weakest negative value {weakest}; offenders, if any, printed above)", ncert_ok == ncert == 1764)

# ================================================================ 3. charge-bound certificates (Theorem 3.5, Table 3)
IDX = load("floor_certificates/index.json")
def introot(N, k):
    if N <= 0: return 0
    x = 1 << ((N.bit_length() + k - 1) // k)
    while True:
        y = ((k - 1) * x + N // x ** (k - 1)) // k
        if y >= x: break
        x = y
    while x ** k > N: x -= 1
    while (x + 1) ** k <= N: x += 1
    return x
def trunc(V, k=11):
    q = (V.numerator * 10 ** k) // V.denominator; s = str(q)
    return s[:-k] + "." + s[-k:] + "..."
def verify_cert(c):
    n = 19; B = lat_gram = GRAM[c["multiplier"]]
    beta = Fr(c["beta"]); D_ = c["D"]
    if not (beta * beta < 12 and beta > 0): return False, "beta"
    R = [[Fr(B[i][j]) for j in range(n)] for i in range(n)]
    sumw = Fr(0); neg_cols = set(); unit_only = True
    for y, w in zip(c["rank1_y"], c["rank1_w"]):
        w = Fr(w); sumw += w
        nz = [i for i in range(n) if y[i]]
        if w < 0:
            if not (len(nz) == 1 and abs(y[nz[0]]) == 1): unit_only = False
            neg_cols.add(nz[0])
        for i in nz:
            for j in nz: R[i][j] -= w * y[i] * y[j]
    pairsum = Fr(0)
    P = c["pair_terms"]
    if P:
        for yi, yj, l1, l2, m in zip(P["yi"], P["yj"], P["l1"], P["l2"], P["m"]):
            l1 = Fr(l1); l2 = Fr(l2); m = Fr(m)
            if l1 < 0 or l2 < 0: return False, "l<0"
            if not any(yi[i] * yj[j] - yi[j] * yj[i] for i in range(n) for j in range(i + 1, n)): return False, "dependent pair"
            a_ = l1 * l1; b_ = m * m + l2 * l2; c_ = l1 * m; pairsum += 2 * l1 * l2
            nzi = [i for i in range(n) if yi[i]]; nzj = [i for i in range(n) if yj[i]]
            for i in nzi:
                for j in nzi: R[i][j] -= a_ * yi[i] * yi[j]
            for i in nzj:
                for j in nzj: R[i][j] -= b_ * yj[i] * yj[j]
            for i in set(nzi) | set(nzj):
                for j in set(nzi) | set(nzj):
                    cr = yi[i] * yj[j] + yj[i] * yi[j]
                    if cr: R[i][j] -= c_ * cr
    if not is_sym(R): return False, "R not symmetric"
    piv, mu = ldl(R)
    if mu is None: return False, "R not positive definite"
    detR = Fr(1)
    for p in piv: detR *= p
    X = D_ * detR
    t = introot((X.numerator * 10 ** 57) // X.denominator, 19)
    if not (Fr(t, 1000) ** 19 <= X): return False, "t"
    V = 4 * sumw + pairsum * beta + 19 * Fr(t, 1000)
    if f"{V.numerator}/{V.denominator}" != c["value_frac"]: return False, f"V={V} != value_frac"
    cond = c["conditional_on_minimal_columns"]
    if not cond and neg_cols: return False, "negative weight in an unconditional certificate"
    if cond and not unit_only: return False, "negative weight on a non-unit vector"
    J = c.get("released_column")
    if J is not None and J in neg_cols: return False, "released column carries a negative weight"
    thr = 46 if c.get("table3_row") in (2, 3) else 48
    if not V > thr: return False, f"V <= {thr}"
    return True, f"V = {trunc(V)} > {thr}; rank-1 terms {len(c['rank1_w'])}, pair terms {len(P['l1']) if P else 0}; " + ("conditional on G_jj = 4" if cond else "unconditional")
files = [f for f in IDX["files"] if (not QUICK) or f["table3_row"] is not None]
for f in files:
    c = load("floor_certificates/" + f["file"]); tt = time.time()
    ok, msg = verify_cert(c)
    check(f"certificate {f['label']} ({f['file']}): {msg}  ({time.time()-tt:.0f}s)", ok)
check("certificates: 5 Table-3 rows + 38 released-column certificates indexed", IDX["n_certificates"] == 43 and sum(1 for f in IDX["files"] if f["table3_row"]) == 5)
if not QUICK:
    mins = {}
    for f in IDX["files"]:
        if f["table3_row"] is None: mins[f["multiplier"]] = min(mins.get(f["multiplier"], Fr(10 ** 9)), Fr(f["value_frac"]))
    check("released-column certificates: minima %s (W512) and %s (W512') both > 48 (decimals truncated)" % (trunc(mins["W512"]), trunc(mins["W512p"])), all(v > 48 for v in mins.values()))
cov = IDX["two_released_columns_surviving_position_pairs"]
check("two released columns: 55 of 171 (W512) and 40 of 171 (W512') position pairs listed", len(cov["W512"]) == 55 and len(cov["W512p"]) == 40)
# local solvability
LOC = load("floor_certificates/local_solvability_156.json")["cases"]
nloc = 0
for key, v in LOC.items():
    pname, ps, ds = key.split("|"); p = int(ps[1:]); dd = int(ds[1:]); Kk = v["K"]
    An, Bn = pname.split("x"); A = GRAM[An]; B = GRAM[Bn]; M = v["final_M"]
    G = mul(mul(M, B), T(M)); q = sum(G[i][j] * A[j][i] for i in range(19) for j in range(19)); Dm = det_int(M)
    mq, md = ((1 << (Kk + 2), 1 << (Kk + 1)) if p == 2 else (3 ** (Kk + 1), 3 ** (Kk + 1)))
    ok = (q - 48) % mq == 0 and (Dm - dd) % md == 0
    Pm = mul(mul(A, M), B); Minv = inv_frac(M)
    Cad = T([[int(Minv[i][j] * Dm) for j in range(19)] for i in range(19)])  # cofactor matrix cof(M) = adj(M)^T = (det M * M^-1)^T
    if p == 2:
        cls = {(Pm[i][j] % 2, Cad[i][j] % 2) for i in range(19) for j in range(19)}; cls.discard((0, 0)); ok = ok and len(cls) >= 2
    else:
        pts = {((2 * Pm[i][j]) % 3, Cad[i][j] % 3) for i in range(19) for j in range(19)}; pts.discard((0, 0))
        ok = ok and any(a1 * b2 - a2 * b1 for (a1, b1) in pts for (a2, b2) in pts)
    nloc += ok
check(f"local solvability: {nloc}/156 cases verified (q = 48 and det = d solvable over Z_p, p in {{2,3}}, d in +-1..13, three unordered det-512 pairs)", nloc == 156 == len(LOC))

# ================================================================ 4. the seven root-free lattices at det 160..224 and their identifications
C7 = load("rootfree_det160_224.json")
G7 = {}
for name, d in C7["lattices"].items():
    G = d["gram"]; G7[name] = G
    ok = is_sym(G) and all(G[i][i] % 2 == 0 for i in range(19)) and det_int(G) == d["det"] and smith_invariants(G) == sorted(d["smith_invariants"])
    sv = short_vectors(G, 4); n2 = sum(1 for v, nn in sv if nn == 2); n4 = sum(1 for v, nn in sv if nn == 4)
    ok = ok and n2 == 0 == d["n_roots"] and len(sv) == n4 == d["kissing_number"] and d["minimum"] == 4
    a = gp_aut_order(G)
    check(f"{name}: det {d['det']}, even, root-free, minimum 4, kissing {d['kissing_number']}, Smith {d['smith_invariants']}" + (f", |Aut| = {d['aut_order']} (gp)" if a is not None else ""), ok and (a is None or a == d["aut_order"]))
REF = {k: v["gram"] for k, v in C7["references"].items()}
niso = 0
for iso in C7["isometries"]:
    A = G7[iso["A"]] if iso["A_is_project"] else REF[iso["A_ref"]]; B = REF[iso["B_ref"]]; U = iso["U"]
    ok = mul(mul(T(U), B), U) == A and abs(det_int(U)) == 1 and det_int(U) == iso["det_U"]
    niso += ok
    if not ok: print("  ISOMETRY FAILS: A =", iso["A"], " B_ref =", iso["B_ref"])
check(f"isometries: {niso}/{len(C7['isometries'])} satisfy U^T B U = A exactly with det U = +-1 (failing pairs, if any, printed above)", niso == len(C7["isometries"]) == 18)

print("\n%d checks passed, %d failed, %.0f s" % (NPASS, len(FAILS), time.time() - t0))
if FAILS:
    print("FAILED:"); [print("  " + f) for f in FAILS]
    sys.exit(1)
print("ALL CHECKS PASSED"); sys.exit(0)