#!/usr/bin/env python3
"""Executable four-dimensional spinor identities with explicit conventions.

Python + NumPy, offline and standalone. Run --help for deterministic CSV/JSON
output options. Natural units c=hbar=1; eta=(+---), epsilon^(0123)=+1;
slash(p)=gamma^0 E-gamma^i p_i; gamma5=i gamma0 gamma1 gamma2 gamma3;
sigma^(mu nu)=i[gamma^mu,gamma^nu]/2. The initial matrices use the Dirac basis.
u(p)e^(-ip.x) and v(p)e^(+ip.x) have future-directed LABEL p=(E,p).
Thus v(p) has canonical momentum -p. Normalizations are ubar u=2m,
vbar v=-2m and u^dagger u=v^dagger v=2E. Hilbert orthogonality is
u(p)^dagger v(-p)=0, not u(p)^dagger v(p)=0.

The discrete maps act on COMMUTING solution columns:
P=gamma0 (linear), C=B K with B=i gamma2, T=U_T K with
U_T=-gamma1 gamma3=diag(-i sigma2,-i sigma2). K conjugates components.
P also reflects space; T reverses time; their coordinate arguments are not
replaced by component matrices. Ordered classical C P T has matrix -gamma5
and is linear after its two conjugations cancel. This computation is NOT a
proof of the quantum-field CPT theorem or a determination of quantum CPT's
square. Fermionic bilinear C signs have an additional statistics sign and
are not tested by commuting-number arrays.

Under psi'=V psi, gamma'=V gamma V^dagger, but B'=V B V^T and
U_T'=V U_T V^T. Rebuilding i gamma'^2 or using similarity for an antilinear
matrix generally fails in a complex basis. Dirac, chiral and a genuinely
complex unitary basis are tested, including transformed adjoints and norms.

Tests cover Clifford/traces/chirality, u/v spin sums, energy projectors,
polarization and all 16 bilinears, Gordon/Ward identities, discrete-map
squares/covariance, Majorana reality and basis congruence. Exact Fraction-
complex fixtures independently verify selected identities without NumPy.
Massless nonzero-momentum fixtures skip division-by-m polarization/Gordon
formulas; the zero-energy apex is excluded. Float residuals include absolute
errors AND declared scales; cancellation at high momentum is not hidden by
a claim of relative accuracy in a tiny invariant. Fixed tolerances apply.

Owners: /relativistic-qm/gamma-matrices/, /relativistic-qm/free-dirac-spinors/,
/relativistic-qm/bilinear-covariants/, /relativistic-qm/gamma-matrix-conventions/,
/relativistic-qm/symmetry-conventions/.
References: Bjorken & Drell, Relativistic Quantum Mechanics (1964);
Dreiner, Haber & Martin, Physics Reports 494, 1-196 (2010).
Tested environment: Python 3.12.14, NumPy 2.3.5; no plotting dependency.
"""
from __future__ import annotations

import argparse
import csv
from fractions import Fraction
import hashlib
from itertools import product
import json
import math
from pathlib import Path
import platform
import sys
import tempfile

import numpy as np

METRIC = np.diag([1.0, -1, -1, -1])
I4 = np.eye(4, dtype=complex)
I2 = np.eye(2, dtype=complex)
ZERO2 = np.zeros((2, 2), dtype=complex)
PAULI = [np.array([[0, 1], [1, 0]], dtype=complex),
         np.array([[0, -1j], [1j, 0]], dtype=complex),
         np.diag([1, -1]).astype(complex)]
ALGEBRA_TOL = 5e-13
KINEMATIC_TOL = 3e-9


def dagger(matrix):
    return np.asarray(matrix).conj().T


def base_package():
    gamma = [np.block([[I2, ZERO2], [ZERO2, -I2]])]
    gamma += [np.block([[ZERO2, s], [-s, ZERO2]]) for s in PAULI]
    g5 = 1j * gamma[0] @ gamma[1] @ gamma[2] @ gamma[3]
    return {"gamma": gamma, "gamma5": g5, "P": gamma[0], "B": 1j * gamma[2],
            "C": 1j * gamma[2] @ gamma[0], "T": -gamma[1] @ gamma[3]}


def transform_package(original, change):
    return {"gamma": [change @ g @ dagger(change) for g in original["gamma"]],
            "gamma5": change @ original["gamma5"] @ dagger(change),
            "P": change @ original["P"] @ dagger(change),
            **{name: change @ original[name] @ change.T for name in ("B", "C", "T")}}


def sigma(gamma, mu, nu):
    return 0.5j * (gamma[mu] @ gamma[nu] - gamma[nu] @ gamma[mu])


def slash(gamma, four_vector):
    return sum((g * value for g, value in zip(gamma, METRIC @ four_vector)), np.zeros((4, 4), complex))


def minkowski(a, b):
    return float(np.asarray(a) @ METRIC @ np.asarray(b))


def epsilon(indices):
    if len(set(indices)) < 4:
        return 0
    return (-1) ** sum(indices[i] > indices[j] for i in range(4) for j in range(i + 1, 4))


def spinors(momentum, mass, spin_basis=I2):
    """Columns u_s,v_s; v_s is labelled by future p but has canonical -p."""
    p = np.asarray(momentum, float)
    if p.shape != (3,) or not np.all(np.isfinite(p)) or not math.isfinite(mass) or mass < 0:
        raise ValueError("finite three-momentum and nonnegative mass required")
    spin_basis = np.asarray(spin_basis, complex)
    if spin_basis.shape != (2, 2) or not np.all(np.isfinite(spin_basis)) or not np.allclose(
            dagger(spin_basis) @ spin_basis, I2, atol=1e-12, rtol=0):
        raise ValueError("spin_basis must contain two orthonormal complex columns")
    e = math.hypot(mass, *p)
    if e == 0:
        raise ValueError("the massless zero-energy apex is excluded")
    pauli_p = sum((component * s for component, s in zip(p, PAULI)), np.zeros((2, 2), complex))
    upper = math.sqrt(e + mass) * spin_basis
    lower = pauli_p @ spin_basis / math.sqrt(e + mass)
    return np.vstack([upper, lower]), np.vstack([lower, upper]), np.r_[e, p]


def bilinear_matrices(package):
    gamma, g5 = package["gamma"], package["gamma5"]
    # label, Gamma, parity sign, time sign, classical commuting C sign.
    result = [("scalar", I4, 1, 1, -1), ("pseudoscalar_i_gamma5", 1j * g5, -1, -1, -1)]
    for mu in range(4):
        result += [(f"vector_{mu}", gamma[mu], 1 if mu == 0 else -1, 1 if mu == 0 else -1, 1),
                   (f"axial_{mu}", gamma[mu] @ g5, -1 if mu == 0 else 1, 1 if mu == 0 else -1, -1)]
    for mu in range(4):
        for nu in range(mu + 1, 4):
            result.append((f"tensor_{mu}{nu}", sigma(gamma, mu, nu), -1 if mu == 0 else 1, 1 if mu == 0 else -1, 1))
    return result


class Checks:
    def __init__(self):
        self.rows = []

    def close(self, name, case, actual, expected, scale=None, tolerance=ALGEBRA_TOL):
        absolute = float(np.max(np.abs(np.asarray(actual) - np.asarray(expected))))
        if scale is None:
            scale = max(1.0, float(np.max(np.abs(actual))), float(np.max(np.abs(expected))))
        residual = absolute / max(float(scale), 1e-300)
        self.rows.append({"name": name, "case": case, "absolute_error": absolute, "scale": float(scale),
                          "scaled_error": residual, "tolerance": tolerance,
                          "passed": math.isfinite(residual) and residual <= tolerance})

    def condition(self, name, case, condition):
        self.rows.append({"name": name, "case": case, "absolute_error": 0 if condition else 1,
                          "scale": 1, "scaled_error": 0 if condition else 1, "tolerance": 0, "passed": bool(condition)})

    @property
    def passed(self):
        return all(row["passed"] for row in self.rows)


def algebra_checks(package, basis_name, checks):
    gamma, g5 = package["gamma"], package["gamma5"]
    beta, b, c, t = gamma[0], package["B"], package["C"], package["T"]
    for mu, nu in product(range(4), repeat=2):
        checks.close("Clifford", f"{basis_name}/{mu}{nu}", gamma[mu] @ gamma[nu] + gamma[nu] @ gamma[mu], 2 * METRIC[mu, nu] * I4)
        checks.close("two-gamma trace", f"{basis_name}/{mu}{nu}", np.trace(gamma[mu] @ gamma[nu]), 4 * METRIC[mu, nu])
    for mu in range(4):
        checks.close("gamma adjoint", f"{basis_name}/{mu}", dagger(gamma[mu]), beta @ gamma[mu] @ beta)
        checks.close("gamma5 anticommutes", f"{basis_name}/{mu}", g5 @ gamma[mu] + gamma[mu] @ g5, 0 * I4)
        checks.close("C transpose identity", f"{basis_name}/{mu}", c @ gamma[mu].T @ dagger(c), -gamma[mu])
        checks.close("B conjugate identity", f"{basis_name}/{mu}", b @ gamma[mu].conj() @ dagger(b), -gamma[mu])
    for indices in product(range(4), repeat=4):
        mu, nu, rho, lam = indices
        matrix = gamma[mu] @ gamma[nu] @ gamma[rho] @ gamma[lam]
        expected = 4 * (METRIC[mu, nu] * METRIC[rho, lam] - METRIC[mu, rho] * METRIC[nu, lam] + METRIC[mu, lam] * METRIC[nu, rho])
        checks.close("four-gamma trace", f"{basis_name}/{indices}", np.trace(matrix), expected)
        checks.close("gamma5 epsilon trace", f"{basis_name}/{indices}", np.trace(g5 @ matrix), -4j * epsilon(indices))
    checks.close("gamma5 square", basis_name, g5 @ g5, I4)
    checks.close("gamma5 Hermitian", basis_name, dagger(g5), g5)
    left, right = (I4 - g5) / 2, (I4 + g5) / 2
    checks.close("chiral projector", basis_name, left @ left, left)
    checks.close("chiral orthogonality", basis_name, left @ right, 0 * I4)
    checks.close("adjoint swaps chirality", basis_name, dagger(left) @ beta, beta @ right)
    checks.close("parity square", basis_name, package["P"] @ package["P"], I4)
    checks.close("classical C square", basis_name, b @ b.conj(), I4)
    checks.close("T antiunitary square", basis_name, t @ t.conj(), -I4)
    checks.close("C relates to B and adjoint", basis_name, c, b @ beta.conj())
    cpt = b @ package["P"].conj() @ t.conj()
    checks.close("ordered classical CPT=-gamma5", basis_name, cpt, -g5)
    checks.close("classical CPT square only", basis_name, cpt @ cpt, I4)
    for label, matrix, psign, tsign, csign in bilinear_matrices(package):
        observable = beta @ matrix
        checks.close("bilinear reality", f"{basis_name}/{label}", dagger(observable), observable)
        checks.close("parity bilinear sign", f"{basis_name}/{label}", dagger(package["P"]) @ observable @ package["P"], psign * observable)
        checks.close("time bilinear sign", f"{basis_name}/{label}", (dagger(t) @ observable @ t).T, tsign * observable)
        checks.close("commuting C bilinear sign", f"{basis_name}/{label}", (dagger(b) @ observable @ b).T, csign * observable)


def kinematic_checks(package, change, basis_name, momentum, mass, spin_basis, label, checks):
    gamma, beta = package["gamma"], package["gamma"][0]
    base_u, base_v, p = spinors(momentum, mass, spin_basis)
    u, v = change @ base_u, change @ base_v
    ubar, vbar = dagger(u) @ beta, dagger(v) @ beta
    pslash, e = slash(gamma, p), p[0]
    case = basis_name + "/" + label
    def verify(name, actual, expected, scale=None):
        checks.close(name, case, actual, expected, scale=scale, tolerance=KINEMATIC_TOL)
    verify("u on shell", (pslash - mass * I4) @ u, 0 * u, 2 * e * math.sqrt(2 * e))
    verify("v on shell", (pslash + mass * I4) @ v, 0 * v, 2 * e * math.sqrt(2 * e))
    verify("u Hilbert normalization", dagger(u) @ u, 2 * e * I2, 2 * e)
    verify("v Hilbert normalization", dagger(v) @ v, 2 * e * I2, 2 * e)
    verify("u covariant normalization", ubar @ u, 2 * mass * I2, 2 * e)
    verify("v covariant normalization", vbar @ v, -2 * mass * I2, 2 * e)
    verify("covariant uv orthogonality", ubar @ v, 0 * I2, 2 * e)
    verify("u spin sum", u @ ubar, pslash + mass * I4, 2 * e)
    verify("v spin sum", v @ vbar, pslash - mass * I4, 2 * e)
    v_opposite = change @ spinors(-np.asarray(momentum), mass, spin_basis)[1]
    verify("fixed-momentum Hilbert uv orthogonality", dagger(u) @ v_opposite, 0 * I2, 2 * e)
    h = sum((p[i + 1] * beta @ gamma[i + 1] for i in range(3)), mass * beta)
    verify("u Hilbert energy projector", u @ dagger(u) / (2 * e), (I4 + h / e) / 2)
    verify("v(-p) Hilbert energy projector", v_opposite @ dagger(v_opposite) / (2 * e), (I4 - h / e) / 2)
    vector = u[:, 0]
    for name, operator, anti, target_sign in (("P", package["P"], False, 1), ("T", package["T"], True, 1), ("C", package["B"], True, -1)):
        transformed = operator @ (vector.conj() if anti else vector)
        reflected_h = sum((-p[i + 1] * beta @ gamma[i + 1] for i in range(3)), mass * beta)
        verify(name + " frequency/momentum covariance", reflected_h @ transformed, target_sign * e * transformed, 2 * e * math.sqrt(2 * e))
        verify(name + " Hilbert norm", np.vdot(transformed, transformed), np.vdot(vector, vector), 2 * e)
    if mass > 0:
        for spin_index in range(2):
            chi = spin_basis[:, spin_index]
            direction = np.array([np.vdot(chi, s @ chi).real for s in PAULI])
            projection = np.dot(momentum, direction)
            polarization = np.r_[projection / mass, direction + np.asarray(momentum) * projection / (mass * (e + mass))]
            spin_projector = (I4 + package["gamma5"] @ slash(gamma, polarization)) / 2
            column = u[:, spin_index]
            adjoint = column.conj() @ beta
            verify("polarization p dot s", minkowski(p, polarization), 0, 2 * e * max(1, float(np.max(np.abs(polarization)))))
            verify("polarization norm", minkowski(polarization, polarization), -1, max(1, float(np.max(np.abs(polarization))) ** 2))
            verify("polarized u outer product", np.outer(column, adjoint), (pslash + mass * I4) @ spin_projector, 2 * e)
            for mu in range(4):
                verify("diagonal vector", adjoint @ gamma[mu] @ column, 2 * p[mu], 2 * e)
                verify("diagonal axial", adjoint @ gamma[mu] @ package["gamma5"] @ column, 2 * mass * polarization[mu], 2 * e)
                for nu in range(mu + 1, 4):
                    p_lower, s_lower = METRIC @ p, METRIC @ polarization
                    expected = -2 * sum(epsilon((mu, nu, rho, lam)) * p_lower[rho] * s_lower[lam]
                                        for rho, lam in product(range(4), repeat=2))
                    verify("polarized tensor epsilon sign", adjoint @ sigma(gamma, mu, nu) @ column, expected, 2 * e)
            verify("diagonal pseudoscalar", adjoint @ package["gamma5"] @ column, 0, 2 * e)
    # Off-diagonal on-shell Ward identity and massive Gordon identity.
    outgoing_momentum = np.asarray(momentum) + max(mass, 1.0) * np.array([0.2, -0.4, 0.1])
    base_out, _, pp = spinors(outgoing_momentum, mass, spin_basis)
    outgoing = change @ base_out[:, 1]
    outbar, incoming = outgoing.conj() @ beta, u[:, 0]
    q = pp - p
    overlap_scale = 2 * math.sqrt(p[0] * pp[0])
    verify("on-shell Ward identity", outbar @ slash(gamma, q) @ incoming, 0, overlap_scale * max(1, float(np.max(np.abs(q)))))
    if mass > 0:
        for mu in range(4):
            left = 2 * mass * (outbar @ gamma[mu] @ incoming)
            operator = (pp[mu] + p[mu]) * I4 + 1j * sum((sigma(gamma, mu, nu) * (METRIC @ q)[nu] for nu in range(4)), 0 * I4)
            verify("massive Gordon identity", left, outbar @ operator @ incoming,
                   overlap_scale * max(2 * mass, abs(pp[mu] + p[mu]), float(np.max(np.abs(q)))))
    spinor_rows, bilinear_rows = [], []
    for frequency, columns in (("u_positive_frequency", u), ("v_negative_frequency_future_label", v)):
        for spin_index in range(2):
            column = columns[:, spin_index]
            for component, value in enumerate(column):
                spinor_rows.append({"basis": basis_name, "case": label, "mass": mass, "E_label": e,
                    "px_label": p[1], "py_label": p[2], "pz_label": p[3], "frequency": frequency,
                    "spin_column": spin_index, "component": component, "real": value.real, "imag": value.imag})
            for name, matrix, _, _, _ in bilinear_matrices(package):
                value = column.conj() @ beta @ matrix @ column
                bilinear_rows.append({"basis": basis_name, "case": label, "frequency": frequency,
                    "spin_column": spin_index, "bilinear": name, "real": value.real, "imag": value.imag})
    return spinor_rows, bilinear_rows


class RationalComplex:
    """Minimal exact complex arithmetic; a pair of Fraction values."""
    def __init__(self, real=0, imag=0):
        if isinstance(real, RationalComplex):
            self.real, self.imag = real.real, real.imag
        else:
            self.real, self.imag = Fraction(real), Fraction(imag)
    def __add__(self, other):
        z = RationalComplex(other)
        return RationalComplex(self.real + z.real, self.imag + z.imag)
    __radd__ = __add__
    def __neg__(self):
        return RationalComplex(-self.real, -self.imag)
    def __sub__(self, other):
        return self + -RationalComplex(other)
    def __mul__(self, other):
        z = RationalComplex(other)
        return RationalComplex(self.real * z.real - self.imag * z.imag, self.real * z.imag + self.imag * z.real)
    __rmul__ = __mul__
    def conjugate(self):
        return RationalComplex(self.real, -self.imag)
    def __eq__(self, other):
        z = RationalComplex(other)
        return self.real == z.real and self.imag == z.imag


def exact_checks(checks):
    """Small rational-complex matrix engine independent of NumPy arithmetic."""
    q, f = RationalComplex, Fraction
    def identity(n):
        return [[q(int(i == j)) for j in range(n)] for i in range(n)]
    def zeros(n, m):
        return [[q() for _ in range(m)] for _ in range(n)]
    def add(a, b):
        return [[x + y for x, y in zip(ar, br)] for ar, br in zip(a, b)]
    def scale(a, factor):
        return [[x * factor for x in row] for row in a]
    def mul(a, b):
        return [[sum((a[i][k] * b[k][j] for k in range(len(b))), q()) for j in range(len(b[0]))] for i in range(len(a))]
    def trans(a):
        return [list(row) for row in zip(*a)]
    def conj(a):
        return [[x.conjugate() for x in row] for row in a]
    def dag(a):
        return trans(conj(a))
    def block(a, b, c, d):
        return [ar + br for ar, br in zip(a, b)] + [cr + dr for cr, dr in zip(c, d)]
    def linear(matrices, values):
        result = zeros(4, 4)
        for matrix, value in zip(matrices, values):
            result = add(result, scale(matrix, value))
        return result
    def exact(name, actual, expected):
        checks.condition("exact " + name, "Fraction complex", actual == expected)
    i2, i4, z2 = identity(2), identity(4), zeros(2, 2)
    imaginary = q(0, 1)
    pauli = [[[q(0), q(1)], [q(1), q(0)]], [[q(0), -imaginary], [imaginary, q(0)]], [[q(1), q(0)], [q(0), q(-1)]]]
    gamma = [block(i2, z2, z2, scale(i2, -1))] + [block(z2, s, scale(s, -1), z2) for s in pauli]
    g5 = scale(mul(mul(mul(gamma[0], gamma[1]), gamma[2]), gamma[3]), imaginary)
    for mu, nu in product(range(4), repeat=2):
        exact(f"Clifford {mu}{nu}", add(mul(gamma[mu], gamma[nu]), mul(gamma[nu], gamma[mu])),
              scale(i4, 2 * (1 if mu == 0 else -1) if mu == nu else 0))
    exact("gamma5 square", mul(g5, g5), i4)
    exact("gamma5 Hermitian", dag(g5), g5)
    b, p = scale(gamma[2], imaginary), gamma[0]
    t = scale(mul(gamma[1], gamma[3]), -1)
    exact("classical C square", mul(b, conj(b)), i4)
    exact("T antiunitary square", mul(t, conj(t)), scale(i4, -1))
    exact("aligned Pauli T phase", t, block(scale(pauli[1], -imaginary), z2, z2, scale(pauli[1], -imaginary)))
    exact("ordered classical CPT", mul(mul(b, conj(p)), conj(t)), scale(g5, -1))
    # m=1,p=(3/4,0,0),E=5/4 makes sqrt(E+m)=3/2 exactly rational.
    chi = [[q(f(3, 5))], [q(f(12, 25), f(16, 25))]]
    u = scale(chi, f(3, 2)) + scale(mul(pauli[0], chi), f(1, 2))
    ubar = mul(dag(u), gamma[0])
    exact("polarized covariant norm", mul(ubar, u), [[q(2)]])
    exact("polarized Hilbert norm", mul(dag(u), u), [[q(f(5, 2))]])
    pslash = linear(gamma, [f(5, 4), f(-3, 4), 0, 0])
    spin = [f(54, 125), f(18, 25), f(96, 125), f(-7, 25)]
    sslash = linear(gamma, [spin[0], -spin[1], -spin[2], -spin[3]])
    projector = scale(add(i4, mul(g5, sslash)), f(1, 2))
    exact("polarized density matrix", mul(u, ubar), mul(add(pslash, i4), projector))
    exact("on-shell polarized spinor", mul(add(pslash, scale(i4, -1)), u), zeros(4, 1))
    for mu in range(4):
        exact(f"axial polarization {mu}", mul(mul(mul(ubar, gamma[mu]), g5), u), [[q(2 * spin[mu])]])
    # A complex global basis phase supplies an exact congruence fixture.
    phase = q(f(3, 5), f(4, 5))
    changed_b = scale(b, phase * phase)
    exact("complex-phase B congruence square", mul(changed_b, conj(changed_b)), i4)
    changed_u = scale(u, phase)
    exact("complex-phase C action", mul(changed_b, conj(changed_u)), scale(mul(b, conj(u)), phase))


def run(args):
    checks = Checks()
    exact_checks(checks)
    rng = np.random.default_rng(args.seed)
    original = base_package()
    chiral = np.block([[I2, -I2], [I2, I2]]) / math.sqrt(2)
    dft = np.array([[1j ** (row * column) for column in range(4)] for row in range(4)], complex) / 2
    complex_change = dft @ np.diag([1, (3 + 4j) / 5, 1j, (-5 + 12j) / 13])
    changes = [("Dirac", I4), ("chiral", chiral), ("complex_unitary", complex_change)]
    cases = [("rest", args.mass, np.zeros(3)), ("tiny", args.mass, args.mass * np.array([1e-10, -2e-10, 3e-10])),
             ("oblique", args.mass, args.mass * np.array([.3, -.4, .7])),
             ("high_ratio", args.mass, args.mass * args.max_ratio * np.array([2, -3, 6]) / 7),
             ("massless_nonzero", 0.0, np.array([.3, -.4, .5]))]
    for index in range(args.samples):
        direction = rng.normal(size=3)
        direction /= np.linalg.norm(direction)
        ratio = 10 ** rng.uniform(-6, math.log10(args.max_ratio))
        cases.append((f"seeded_{index}", args.mass, args.mass * ratio * direction))
    spin_bases = {}
    for label, _, _ in cases:
        matrix = rng.normal(size=(2, 2)) + 1j * rng.normal(size=(2, 2))
        spin_bases[label] = np.linalg.qr(matrix)[0]
    matrix_rows, spinor_rows, bilinear_rows = [], [], []
    reference_bilinears = {}
    for basis_name, change in changes:
        package = original if basis_name == "Dirac" else transform_package(original, change)
        checks.close("basis is unitary", basis_name, dagger(change) @ change, I4)
        algebra_checks(package, basis_name, checks)
        exported = {f"gamma{mu}": matrix for mu, matrix in enumerate(package["gamma"])}
        exported.update({key: package[key] for key in ("gamma5", "P", "B", "C", "T")})
        for name, matrix in exported.items():
            for row, column in product(range(4), repeat=2):
                matrix_rows.append({"basis": basis_name, "matrix": name, "row": row, "column": column,
                                    "real": matrix[row, column].real, "imag": matrix[row, column].imag})
        for label, mass, momentum in cases:
            srows, brows = kinematic_checks(package, change, basis_name, momentum, mass, spin_bases[label], label, checks)
            spinor_rows.extend(srows)
            bilinear_rows.extend(brows)
            for row in brows:
                key = (row["case"], row["frequency"], row["spin_column"], row["bilinear"])
                value = complex(row["real"], row["imag"])
                if basis_name == "Dirac":
                    reference_bilinears[key] = value
                else:
                    scale = 2 * math.hypot(mass, *momentum)
                    checks.close("bilinear basis invariance", basis_name + "/" + str(key), value, reference_bilinears[key], scale=scale, tolerance=KINEMATIC_TOL)
        # A generic commuting column tests Majorana reality and map covariance.
        psi = np.array([1 + 2j, -3 + .5j, .7 - 1j, 2 - .3j])
        changed = change @ psi
        other = change @ np.array([-.2 + .9j, 1 - 3j, -2 + .1j, .5 + .7j])
        for name in ("B", "T"):
            checks.close("antilinear basis covariance", basis_name + "/" + name, package[name] @ changed.conj(), change @ original[name] @ psi.conj())
            checks.close("antilinear inner product", basis_name + "/" + name,
                np.vdot(package[name] @ changed.conj(), package[name] @ other.conj()), np.vdot(other, changed))
        cpt = package["B"] @ package["P"].conj() @ package["T"].conj()
        checks.close("classical CPT norm", basis_name, np.vdot(cpt @ changed, cpt @ changed), np.vdot(changed, changed))
        checks.close("classical CPT linearity", basis_name, cpt @ (1j * changed), 1j * (cpt @ changed))
        majorana = changed + package["B"] @ changed.conj()
        checks.close("Majorana real projection", basis_name, package["B"] @ majorana.conj(), majorana)
        parity_majorana = 1j * package["P"] @ majorana
        checks.close("Majorana imaginary parity phase", basis_name, package["B"] @ parity_majorana.conj(), parity_majorana)
        time_majorana = package["T"] @ majorana.conj()
        checks.close("Majorana aligned time phase", basis_name, package["B"] @ time_majorana.conj(), time_majorana)
        checks.close("norm is basis invariant", basis_name, np.vdot(changed, changed), np.vdot(psi, psi))
    expected_chiral = [np.block([[ZERO2, I2], [I2, ZERO2]])] + [np.block([[ZERO2, s], [-s, ZERO2]]) for s in PAULI]
    for mu in range(4):
        checks.close("explicit chiral table", str(mu), chiral @ original["gamma"][mu] @ dagger(chiral), expected_chiral[mu])
    checks.close("explicit chiral gamma5", "table", chiral @ original["gamma5"] @ dagger(chiral), np.diag([-1, -1, 1, 1]))
    t_common = 1j * original["gamma"][1] @ original["gamma"][3]
    checks.close("Dirac C matrix square distinct from full map", "Dirac only", original["C"] @ original["C"], -I4)
    checks.close("aligned T ordinary matrix square", "Dirac only", original["T"] @ original["T"], -I4)
    checks.close("aligned T=i times common T", "Dirac", original["T"], 1j * t_common)
    checks.close("common T matrix square only", "Dirac", t_common @ t_common, I4)
    checks.close("common T antiunitary square", "Dirac", t_common @ t_common.conj(), -I4)
    checks.close("alternative classical CPT phase", "Dirac", original["B"] @ original["P"].conj() @ t_common.conj(), -1j * original["gamma5"])
    wrong_b = complex_change @ original["B"] @ dagger(complex_change)
    right_b = complex_change @ original["B"] @ complex_change.T
    wrong_gap = float(np.max(np.abs(wrong_b - right_b)))
    checks.condition("negative control: similarity is wrong for B", "complex unitary", wrong_gap > 0.1)
    rebuilt_gap = float(np.max(np.abs(right_b - 1j * transform_package(original, complex_change)["gamma"][2])))
    checks.condition("negative control: i gamma2 cannot be blindly rebuilt", "complex unitary", rebuilt_gap > 0.1)
    u, v, _ = spinors([.8, .2, -.3], 1)
    label_gap = float(np.max(np.abs(dagger(u) @ v)))
    checks.condition("negative control: wrong v momentum label is not Hilbert orthogonal", "equal future labels", label_gap > 0.1)
    try:
        spinors([0, 0, 0], 0)
    except ValueError:
        checks.condition("massless zero-energy apex is rejected", "E=0", True)
    else:
        checks.condition("massless zero-energy apex is rejected", "E=0", False)
    return checks, {"matrices": matrix_rows, "spinors": spinor_rows, "bilinears": bilinear_rows}, {
        "wrong_similarity_B_gap": wrong_gap, "blind_rebuild_B_gap": rebuilt_gap, "wrong_v_label_overlap": label_gap}


def parser():
    p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    p.add_argument("--output-dir", type=Path, help="new CSV/JSON directory; omitted: new temporary directory")
    p.add_argument("--json-output", type=Path, help="override summary JSON path")
    p.add_argument("--checks-csv", type=Path, help="override detailed checks CSV path")
    p.add_argument("--samples", type=int, default=12, help="seeded cases beyond five fixed fixtures, 0..1000 (default 12)")
    p.add_argument("--seed", type=int, default=20261002, help="nonnegative PCG64 seed")
    p.add_argument("--mass", type=float, default=1.0, help="positive massive-fixture mass, 1e-6..1e6 (default 1)")
    p.add_argument("--max-ratio", type=float, default=1000.0, help="largest massive |p|/m, 1..1e6 (default 1000)")
    return p


def main(argv=None):
    p = parser()
    args = p.parse_args(argv)
    if not 0 <= args.samples <= 1000 or args.seed < 0:
        p.error("samples must be 0..1000 and seed nonnegative")
    if not math.isfinite(args.mass) or not 1e-6 <= args.mass <= 1e6:
        p.error("mass must be finite and in [1e-6,1e6]; an additional massless nonzero fixture is always included")
    if not math.isfinite(args.max_ratio) or not 1 <= args.max_ratio <= 1e6:
        p.error("max-ratio must be finite and in [1,1e6]")
    output = args.output_dir or Path(tempfile.mkdtemp(prefix="spinor-algebra-"))
    paths = {name: output / (name + ".csv") for name in ("matrices", "spinors", "bilinears")}
    paths["checks"] = args.checks_csv or output / "checks.csv"
    paths["summary"] = args.json_output or output / "summary.json"
    if len({path.resolve() for path in paths.values()}) != len(paths):
        p.error("output paths must be distinct")
    for path in paths.values():
        if path.exists() or path.is_symlink():
            p.error(f"existing outputs are never overwritten: {path}")
        if ".git" in [part.lower() for part in path.resolve().parts]:
            p.error(f"output cannot be inside .git: {path}")
    checks, data, negative_controls = run(args)
    data["checks"] = checks.rows
    for name, rows in data.items():
        paths[name].parent.mkdir(parents=True, exist_ok=True)
        with paths[name].open("x", encoding="utf8", newline="") as stream:
            writer = csv.DictWriter(stream, fieldnames=list(rows[0]))
            writer.writeheader()
            writer.writerows(rows)
    summary_by_name = {}
    for row in checks.rows:
        item = summary_by_name.setdefault(row["name"], {"count": 0, "maximum_absolute_error": 0.0, "maximum_scaled_error": 0.0, "passed": True})
        item["count"] += 1
        item["maximum_absolute_error"] = max(item["maximum_absolute_error"], row["absolute_error"])
        item["maximum_scaled_error"] = max(item["maximum_scaled_error"], row["scaled_error"])
        item["passed"] = item["passed"] and row["passed"]
    report = {"program": "spinor-algebra.py", "passed": checks.passed, "check_count": len(checks.rows),
        "script_sha256": hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
        "environment": {"python": platform.python_version(), "numpy": np.__version__, "platform": platform.platform(), "rng": "PCG64"},
        "parameters": {name: getattr(args, name) for name in ("samples", "seed", "mass", "max_ratio")},
        "conventions": {"metric": "+---", "epsilon_0123_upper": 1, "units": "c=hbar=1",
            "u_v_norm": "ubar u=2m; vbar v=-2m; udag u=vdag v=2E; v future label has canonical -p",
            "T_phase": "U_T=-gamma1 gamma3 in Dirac basis; transform by congruence",
            "classical_CPT": "C P T=-gamma5 times spacetime inversion; linear commuting-column map only"},
        "tolerances": {"matrix_algebra_scaled": ALGEBRA_TOL, "kinematic_scaled": KINEMATIC_TOL},
        "summary": summary_by_name, "failed_checks": [row for row in checks.rows if not row["passed"]],
        "negative_controls": negative_controls,
        "limitations": ["Finite exact and floating fixtures do not prove identities for all representations or momenta.",
            "Absolute errors and scales are exported; high-momentum cancellation can lose digits in small invariants.",
            "Massless nonzero modes skip massive rest-polarization and Gordon formulas; E=0 is excluded.",
            "C/P/T are commuting solution maps with declared phases, not quantum Fock-space operators or a CPT theorem test.",
            "Only constant unitary component-basis changes are tested; Lorentz boosts are a different, generally nonunitary representation."],
        "output_files": {name: str(path) for name, path in paths.items()},
        "output_sha256": {name: hashlib.sha256(path.read_bytes()).hexdigest() for name, path in paths.items() if name != "summary"}}
    paths["summary"].parent.mkdir(parents=True, exist_ok=True)
    with paths["summary"].open("x", encoding="utf8") as stream:
        json.dump(report, stream, indent=2, allow_nan=False)
        stream.write("\n")
    print(json.dumps({"passed": checks.passed, "checks": len(checks.rows), "failed_checks": report["failed_checks"][:12],
                      "failed_count": len(report["failed_checks"]), "summary_json": str(paths["summary"])}, indent=2))
    return 0 if checks.passed else 1


if __name__ == "__main__":
    sys.exit(main())