#!/usr/bin/env python3
"""verify_charge22.py -- Python checker for charge22.json (the charge-22 K3 x K3 flux; main text Section 3.4, SM Section S6).

Python 3 standard library only; exact integer/rational arithmetic in Q[x]/(Phi_66) and on integer matrices. The signs of the
real embeddings are evaluated in floating point (mpmath at 50 digits if installed, else double precision) against the
exact margin 8.9e-6; if sympy is installed, irreducibility of m(x) is also checked directly. Run in this directory:

    python3 verify_charge22.py        (well under a minute)

Exit status 0 iff every check passes. verify_charge22.gp performs the same checks in PARI/GP.
"""
import json, os, sys, math, cmath, time
from fractions import Fraction as Fr
HERE = os.path.dirname(os.path.abspath(__file__))
D = json.load(open(os.path.join(HERE, "charge22.json")))
FAILS = []; NP = 0
def check(name, ok, detail=""):
    global NP
    NP += ok
    if not ok: FAILS.append(name)
    print(("PASS  " if ok else "FAIL  ") + name + (f"  [{detail}]" if detail != "" else ""), flush=True)
try:
    import mpmath; mpmath.mp.dps = 50; HAVE_MP = True
except Exception: HAVE_MP = False
try:
    import sympy; HAVE_SYMPY = True
except Exception: HAVE_SYMPY = False

# ---------------------------------------------------------------- polynomials over Q as coefficient lists c_0..c_n
def ptrim(p):
    p = list(p)
    while len(p) > 1 and p[-1] == 0: p.pop()
    return p
def padd(a, b): n = max(len(a), len(b)); return ptrim([(a[i] if i < len(a) else 0) + (b[i] if i < len(b) else 0) for i in range(n)])
def pscale(a, c): return ptrim([c * x for x in a])
def pmul(a, b):
    r = [0] * (len(a) + len(b) - 1)
    for i, x in enumerate(a):
        if x:
            for j, y in enumerate(b): r[i + j] += x * y
    return ptrim(r)
def pdivmod(a, b):
    a = [Fr(x) for x in a]; b = ptrim([Fr(x) for x in b]); q = [Fr(0)] * max(1, len(a) - len(b) + 1)
    while len(ptrim(a)) >= len(b) and any(a):
        a = ptrim(a); k = len(a) - len(b); c = a[-1] / b[-1]; q[k] = c
        for i, y in enumerate(b): a[i + k] -= c * y
        a = ptrim(a)
        if len(a) < len(b): break
    return ptrim(q), ptrim(a)
def pmod(a, m): return pdivmod(a, m)[1]
def pgcd(a, b):
    a = ptrim([Fr(x) for x in a]); b = ptrim([Fr(x) for x in b])
    while any(b): a, b = b, pmod(a, b)
    return pscale(a, 1 / a[-1])
def pderiv(a): return ptrim([i * a[i] for i in range(1, len(a))]) if len(a) > 1 else [0]
def peval(a, x):
    r = 0
    for c in reversed(a): r = r * x + c
    return r
def cyclotomic(n):
    # Phi_n = (x^n - 1) / prod_{d | n, d < n} Phi_d
    num = [-1] + [0] * (n - 1) + [1]
    for d_ in range(1, n):
        if n % d_ == 0:
            q, r = pdivmod(num, cyclotomic(d_)); assert not any(r); num = [int(c) for c in q]
    return num
PHI = D["field"]["Phi66"]
check("Phi66 as printed equals the 66th cyclotomic polynomial", cyclotomic(66) == PHI)
# arithmetic in K = Q[x]/(Phi66): elements are coefficient lists of length <= 20
def kmul(a, b): return [Fr(c) for c in pmod(pmul(a, b), PHI)]
def kpow(a, e):
    r = [Fr(1)]; base = a
    while e:
        if e & 1: r = kmul(r, base)
        base = kmul(base, base); e >>= 1
    return r
def kinv(a):
    # extended Euclid in Q[x]
    r0, r1 = ptrim([Fr(c) for c in PHI]), ptrim([Fr(c) for c in a]); s0, s1 = [Fr(0)], [Fr(1)]
    while any(r1):
        q, r = pdivmod(r0, r1); r0, r1 = r1, r; s0, s1 = s1, padd(s0, pscale(pmul(q, s1), -1))
    return kmul(pscale(s0, 1 / r0[0]), [1]) if len(r0) == 1 else None
def kint(a): a = [Fr(c) for c in a] + [Fr(0)] * (20 - len(a)); assert all(c.denominator == 1 for c in a); return [int(c) for c in a[:20]]
Z = [0, 1]
def zp(e): return kpow(Z, e % 66)
def conjk(a):  # zeta -> zeta^-1
    r = [Fr(0)]
    for j, c in enumerate(a):
        if c: r = padd(r, pscale(zp(-j), Fr(c)))
    return kmul(r, [1])
# power sums / traces via Newton's identities on Phi66 (monic, degree 20)
n = 20; e = PHI  # e[0..20], e[20] = 1
# Newton: p_k + e_{n-1} p_{k-1} + ... + e_{n-k+1} p_1 + k e_{n-k} = 0 (k <= n); for k > n: p_k + sum_{i=1}^{n} e_{n-i} p_{k-i} = 0
psum = [Fr(20)]
for k in range(1, 200):
    s = Fr(0)
    for i in range(1, min(k, n) + 1):
        if i == k: s -= k * e[n - k]
        else: s -= e[n - i] * psum[k - i]
    psum.append(s)
def ktrace(a): return sum(Fr(c) * psum[j] for j, c in enumerate(a))
check("trace check: Tr(1) = 20, Tr(zeta) = mu(66) = -1, Tr(zeta^33) = -20", psum[0] == 20 and psum[1] == -1 and psum[33] == -20 and psum[66] == 20)

# ---------------------------------------------------------------- delta, eps_S, the trace form
delta = kmul(padd(zp(22), pscale(zp(-22), -1)), kpow(padd(zp(6), pscale(zp(-6), -1)), 9))
check("delta = (zeta^22 - zeta^-22)(zeta^6 - zeta^-6)^9 equals the printed power-basis element", kint(delta) == D["delta"]["power_basis"])
check("delta is real", kint(conjk(delta)) == D["delta"]["power_basis"])
def u(a):  # (zeta^a - zeta^-a)/(zeta - zeta^-1) = sum_{i=0}^{a-1} zeta^(2i-(a-1))
    r = [Fr(0)]
    for i in range(a): r = padd(r, zp(2 * i - (a - 1)))
    return kmul(r, [1])
eps = [Fr(D["eps_S"]["overall_sign"])]
for a, ex in D["eps_S"]["u_exponents"].items():
    if ex: eps = kmul(eps, kpow(u(int(a)), ex))
check("eps_S = -u7 u17 u19 u23 u29 equals the printed power-basis element", kint(eps) == D["eps_S"]["power_basis"])
epsinv = kinv(eps)
check("eps_S is a real unit (eps_S^-1 integral)", kint(conjk(eps)) == D["eps_S"]["power_basis"] and epsinv is not None and all(Fr(c).denominator == 1 for c in epsinv))
dinv = kinv(delta); w = kmul(eps, dinv)
tau = [ktrace(kmul(w, zp(k))) for k in range(20)]
check("tau(k) = Tr(eps_S delta^-1 zeta^k), k = 0..19, integral and equal to the printed first row", all(t.denominator == 1 for t in tau) and [int(t) for t in tau] == D["L"]["tau_first_row"])
tau33 = [ktrace(kmul(w, zp(k))) for k in range(34)]
check("tau(33-k) = -tau(k)", all(tau33[33 - k] == -tau33[k] for k in range(34)))
B = D["L"]["gram_B"]
check("B is the symmetric Toeplitz matrix of tau (B_ij = Tr(eps delta^-1 zeta^(i-j)))", all(B[i][j] == tau[abs(i - j)] for i in range(20) for j in range(20)))
# ---------------------------------------------------------------- integer matrices
def T(M): return [list(r) for r in zip(*M)]
def mul(A, Bm): Bt = T(Bm); return [[sum(a * b for a, b in zip(r, c)) for c in Bt] for r in A]
def eye(n): return [[int(i == j) for j in range(n)] for i in range(n)]
def madd(A, Bm, s=1): return [[a + s * b for a, b in zip(r, q)] for r, q in zip(A, Bm)]
def mpow(M, e):
    R = eye(len(M))
    while e:
        if e & 1: R = mul(R, M)
        M = mul(M, M); e >>= 1
    return R
def det_int(M):
    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 inertia(Q):
    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:
            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]; pos += p > 0; neg += p < 0; row = A[i][:]
        for r in active:
            if r != i and A[r][i] != 0:
                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 charpoly(M):
    """Faddeev-LeVerrier, exact; returns c_0..c_n of det(x I - M)."""
    n = len(M); c = [Fr(0)] * (n + 1); c[n] = Fr(1); Mk = [[Fr(0)] * n for _ in range(n)]; I = eye(n)
    for k in range(1, n + 1):
        Mk = mul(M, madd(Mk, [[c[n - k + 1] * I[i][j] for j in range(n)] for i in range(n)]))
        c[n - k] = -sum(Mk[i][i] for i in range(n)) / k
    return [int(x) for x in c]
def blockdiag(A, Bm):
    n, m = len(A), len(Bm)
    return [A[i] + [0] * m for i in range(n)] + [[0] * n + Bm[i] for i in range(m)]
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]

check("B: even (diagonal 12726), det B = 1, signature (2,18)", all(B[i][i] == 12726 for i in range(20)) and det_int(B) == 1 and inertia(B) == (2, 18, 0) and D["L"]["signature"] == [2, 18])
Mz = [[int(Fr(c)) for c in (zp(j + 1) + [Fr(0)] * 20)[:20]] for j in range(20)]; Mz = T(Mz)   # column j = zeta^(j+1)
check("M_zeta (multiplication by zeta) equals the printed matrix; isometry of B; order 66; charpoly Phi66",
      Mz == D["L"]["M_zeta"] and mul(mul(T(Mz), B), Mz) == B and mpow(Mz, 66) == eye(20) and all(mpow(Mz, 66 // p) != eye(20) for p in (2, 3, 11)) and charpoly(Mz) == PHI)
dU = [[0, 1], [1, 0]]; d = blockdiag(dU, B)
check("d = d_U (+) B equals the printed Gram of Lambda; even, det -1, signature (3,19), rank 22", d == D["Lambda"]["gram_d"] and det_int(d) == -1 and inertia(d) == (3, 19, 0) and len(d) == 22)
F = D["flux"]; UB = D["U_block"]
Mzinv = mpow(Mz, 65)
check("M_(1+zeta) = 1 + M_zeta, M_(1+zeta^-1) = 1 + M_zeta^-1 as printed", F["M_1plus_zeta"] == madd(eye(20), Mz) and F["M_1plus_zeta_inverse"] == madd(eye(20), Mzinv))
g = blockdiag(UB["g_U"], madd(eye(20), Mzinv)); gt = blockdiag(UB["gtilde_U"], madd(eye(20), Mz))
check("g = g_U (+) M_(1+zeta^-1) and gtilde = gtilde_U (+) M_(1+zeta) equal the printed 22x22 matrices", g == F["g"] and gt == F["gtilde"])
check("U block: gtilde_U^T d_U = d_U g_U; S_U = g_U gtilde_U = [[3,4],[2,3]], charpoly x^2-6x+1", mul(T(UB["gtilde_U"]), dU) == mul(dU, UB["g_U"]) and mul(UB["g_U"], UB["gtilde_U"]) == UB["S_U"] == [[3, 4], [2, 3]] and charpoly(UB["S_U"]) == [1, -6, 1] == UB["charpoly"])
check("g and gtilde integral and mutually d-adjoint: gtilde^T d = d g; det g = det gtilde = 1", mul(T(gt), d) == mul(d, g) and det_int(g) == 1 and det_int(gt) == 1)
S = mul(g, gt); St = mul(gt, g)
check("S = g gtilde, Stilde = gtilde g equal the printed matrices and are d-self-adjoint", S == F["S"] and St == F["Stilde"] and mul(T(S), d) == mul(d, S) and mul(T(St), d) == mul(d, St))
trS = sum(S[i][i] for i in range(22))
check("tr S = 44 = 6 (U) + 38 (L); det S = 1; N_flux = 22 <= 24; n_M2 = 2", trS == 44 == F["tr_S"] and sum(S[i][i] for i in range(2)) == 6 and det_int(S) == 1 and F["N_flux"] == trS // 2 == 22 and F["n_M2"] == 24 - 22 == 2)
check("38 = Tr(2 + zeta + zeta^-1) = 40 + 2 mu(66)", ktrace(padd([2], padd(zp(1), zp(-1)))) == 38)
dinvM = inv_frac(d); N = mul([[Fr(x) for x in r] for r in g], dinvM)
check("N = g d^-1 is integral and equals the printed N; gtilde = N^T d; (1/2) tr(N d N^T d) = 22; max |N_ij| and nonzero count as printed",
      all(x.denominator == 1 for r in N for x in r) and [[int(x) for x in r] for r in N] == F["N"] and mul(T(F["N"]), d) == gt
      and sum(mul(mul(mul(F["N"], d), T(F["N"])), d)[i][i] for i in range(22)) == 44 and max(abs(x) for r in F["N"] for x in r) == F["N_max_abs_entry"] and sum(1 for r in F["N"] for x in r if x) == F["N_nonzero_entries"])
CP = D["characteristic_polynomial"]; m = CP["m"]; q2 = [1, -6, 1]
chi = charpoly(S)
check("charpoly(S) equals the printed chi_S and (x^2 - 6x + 1) m(x)^2; charpoly(Stilde) the same", chi == CP["chi_S"] == pmul(q2, pmul(m, m)) and charpoly(St) == chi)
alpha_ab = padd([2], padd(zp(1), zp(-1)))
mval = [Fr(0)]
for c in reversed(m): mval = padd(kmul(mval, alpha_ab), [c])
check("m(2 + zeta + zeta^-1) = 0 in K and deg m = 10 = [Q(zeta+zeta^-1):Q], so m is its minimal polynomial (irreducible)", not any(kmul(mval, [1])) and len(m) == 11)
if HAVE_SYMPY:
    x = sympy.Symbol('x'); check("m irreducible over Q (sympy)", sympy.Poly(list(reversed(m)), x).is_irreducible)
def disc_from_roots_free(p):  # discriminant via resultant of p and p' (Sylvester determinant)
    dp = pderiv(p); n_, m_ = len(p) - 1, len(dp) - 1; N_ = n_ + m_
    Sy = []
    for i in range(m_): Sy.append([0] * i + list(reversed(p)) + [0] * (N_ - n_ - 1 - i))
    for i in range(n_): Sy.append([0] * i + list(reversed(dp)) + [0] * (N_ - m_ - 1 - i))
    res = det_int(Sy); return (-1) ** (n_ * (n_ - 1) // 2) * res // p[-1]
check("disc(m) = 3^5 11^9 = %s as printed" % CP["disc_m"], disc_from_roots_free(m) == 3 ** 5 * 11 ** 9 == int(CP["disc_m"]))
check("gcd(m, m') = 1 and gcd(x^2-6x+1, m) = 1: minimal polynomial (x^2-6x+1) m is squarefree, S diagonalizable; sqrt 2 not rational (32 not a square)",
      pgcd(m, pderiv(m)) == [1] and pgcd(q2, m) == [1] and math.isqrt(32) ** 2 != 32)
# Sturm: m has 10 real roots in (0, 4]
def sturm_count(p, a, b):
    seq = [ptrim([Fr(c) for c in p]), ptrim([Fr(c) for c in pderiv(p)])]
    while len(seq[-1]) > 1 or seq[-1][0] != 0:
        r = pmod(seq[-2], seq[-1])
        if not any(r): break
        seq.append(pscale(r, -1))
    def V(x):
        vals = [peval(s, Fr(x)) for s in seq]; vals = [v for v in vals if v != 0]
        return sum(1 for i in range(len(vals) - 1) if (vals[i] > 0) != (vals[i + 1] > 0))
    return V(a) - V(b)
check("all twelve eigenvalues real and positive: m has 10 roots in (0,4), x^2-6x+1 has roots 3 +- 2 sqrt 2 > 0", sturm_count(m, 0, 4) == 10 and peval(m, 0) != 0)
# ---------------------------------------------------------------- real places, signs, eigenplanes (floating point against exact margins)
places = D["field"]["real_places"]
check("real places: k in [1,32] coprime to 66", places == [k for k in range(1, 33) if math.gcd(k, 66) == 1])
def emb(a, k):
    if HAVE_MP:
        zz = mpmath.exp(2j * mpmath.pi * k / 66); return sum(mpmath.mpf(Fr(c).numerator) / Fr(c).denominator * zz ** j for j, c in enumerate(a))
    zz = cmath.exp(2j * math.pi * k / 66); return sum(float(Fr(c)) * zz ** j for j, c in enumerate(a))
re = (lambda v: float(mpmath.re(v))) if HAVE_MP else (lambda v: v.real)
sw = [re(emb(eps, k) / emb(delta, k)) for k in places]; sd = [re(1 / emb(delta, k)) for k in places]
check("signs of sigma_k(delta^-1) at the real places as printed", [1 if v > 0 else -1 for v in sd] == D["delta"]["signs_of_delta_inverse_at_real_places"])
check("eps_S delta^-1 positive at exactly one real place (k = 1), negative at the nine others; min |sigma_k| = %.6e matches the printed margin" % min(map(abs, sw)),
      [1 if v > 0 else -1 for v in sw] == D["eps_S"]["signs_of_eps_delta_inverse_at_real_places"] and [k for k, v in zip(places, sw) if v > 0] == [D["eps_S"]["positive_place"]] == [1]
      and abs(min(map(abs, sw)) - float(D["eps_S"]["min_abs_eps_delta_inverse_at_real_places"])) < 1e-9 * (1 if HAVE_MP else 100))
ok = True
for row in D["spectrum"]["L_eigenplanes"]:
    k = row["k"]; lam = 2 + 2 * math.cos(2 * math.pi * k / 66)
    ok = ok and abs(lam - float(row["lambda_k"])) < 1e-12 and row["sign_of_b_on_P_k"] == (1 if sw[places.index(k)] > 0 else -1)
check("eigenplane table: lambda_k = 2 + 2 cos(2 pi k/66) and sign of b on P_k = sign sigma_k(eps_S delta^-1), for all ten k", ok)
# structural identity: B = sum_k 2 sigma_k(w) Re(v_k v_k^dagger), v_k = (sigma_k(zeta)^j)_j  => b|P_k definite with sign sigma_k(w)
err = 0.0
for i in range(20):
    for j in range(20):
        s = sum(2 * sw[t] * math.cos(2 * math.pi * places[t] * (i - j) / 66) for t in range(10))
        err = max(err, abs(s - B[i][j]))
check("B = sum over real places of 2 sigma_k(eps delta^-1) Re(v_k v_k^dagger) (max deviation %.1e): every P_k is b-definite with sign sigma_k" % err, err < 1e-6)
# U-line: exact in Z[sqrt2]: represent a + b sqrt2 as (a,b)
def q2mul(p, q): return (p[0] * q[0] + 2 * p[1] * q[1], p[0] * q[1] + p[1] * q[0])
v = [(0, 1), (1, 0)]  # (sqrt2, 1)
Sv = [tuple(sum(x) for x in zip(*[q2mul((UB["S_U"][i][j], 0), v[j]) for j in range(2)])) for i in range(2)]
lv = [q2mul((3, 2), v[i]) for i in range(2)]
normv = q2mul(v[0], v[1]); normv = (2 * normv[0], 2 * normv[1])   # v^T d_U v = 2 v0 v1
check("U: S_U (sqrt2,1) = (3+2sqrt2)(sqrt2,1) exactly and its d_U-norm is 2 sqrt2 > 0; with P_1 this gives exactly three positive directions = n_+(Lambda)",
      Sv == lv and normv == (0, 2) and sum(1 for s in sw if s > 0) * 2 + 1 == 3 == D["spectrum"]["n_positive_directions"] == inertia(d)[0])
check("Lambda cap Sigma^perp = 0: Phi66 = charpoly(M_zeta) is irreducible (cyclotomic), so sigma_1(x) = 0 forces x = 0; x^2-6x+1 irreducible, so u perp (sqrt2,1) forces u = 0", charpoly(Mz) == cyclotomic(66) and math.isqrt(32) ** 2 != 32)
# ---------------------------------------------------------------- Kondo / K3Groups
KD = D["kondo"]; G = KD["G_BH"]; M = KD["M_BH"]; A = T(M); P = KD["P"]
h = blockdiag(eye(2), Mz)
check("G_BH: even, det -1, signature (3,19)", all(G[i][i] % 2 == 0 for i in range(22)) and det_int(G) == -1 and inertia(G) == (3, 19, 0))
check("M_BH G_BH M_BH^T = G_BH (row convention); det 1; order 66; charpoly (x-1)^2 Phi66", mul(mul(M, G), T(M)) == G and det_int(M) == 1 and mpow(M, 66) == eye(22) and all(mpow(M, 66 // p) != eye(22) for p in (2, 3, 11)) and charpoly(M) == pmul([1, -2, 1], PHI))
check("h = 1_U (+) M_zeta is an isometry of d of order 66 with charpoly (x-1)^2 Phi66", mul(mul(T(h), d), h) == d and mpow(h, 66) == eye(22) and charpoly(h) == pmul([1, -2, 1], PHI))
check("det P = 1, P^T G_BH P = d, max |P_ij| = 30", det_int(P) == 1 == KD["det_P"] and mul(mul(T(P), G), P) == d and max(abs(x) for r in P for x in r) == 30)
A5 = mpow(A, 5)
check("A^5 P = P h and A P = P h^53 (A = M_BH^T)", mul(A5, P) == mul(P, h) and mul(A, P) == mul(P, mpow(h, 53)))
check("g h = h g and gtilde h = h gtilde (the flux is invariant under the order-66 isometry)", mul(g, h) == mul(h, g) and mul(gt, h) == mul(h, gt))
Cc = T([[int(Fr(c)) for c in (zp(-j) + [Fr(0)] * 20)[:20]] for j in range(20)])
Pc = mul(P, blockdiag(eye(2), Cc))
check("complex conjugation C on Z[zeta]: C^T B C = B, C M_zeta C^-1 = M_zeta^-1; P (1 (+) C) is an isometry carrying h to A^61", mul(mul(T(Cc), B), Cc) == B and mul(mul(Cc, Mz), Cc) == Mzinv and mul(mul(T(Pc), G), Pc) == d and mul(mpow(A, 61), Pc) == mul(Pc, h) and KD["admissible_powers"] == [5, 61])
Pinv = inv_frac(P)
gB = mul([[Fr(x) for x in r] for r in mul(P, g)], Pinv); gtB = mul([[Fr(x) for x in r] for r in mul(P, gt)], Pinv)
check("transported flux P g P^-1, P gtilde P^-1: integral, commute with A; charpoly of the product is (x^2-6x+1) m^2",
      all(x.denominator == 1 for r in gB + gtB for x in r) and mul(gB, A) == mul(A, gB) and charpoly([[int(x) for x in r] for r in mul(gB, gtB)]) == chi)
print("\n%d checks passed, %d failed (%s, %s)" % (NP, len(FAILS), "mpmath" if HAVE_MP else "float", "sympy" if HAVE_SYMPY else "no sympy"))
if FAILS:
    print("FAILED:"); [print("  " + f) for f in FAILS]; sys.exit(1)
print("ALL CHECKS PASSED"); sys.exit(0)