#!/usr/bin/env python3
"""
Computation 66 -- Lorentzian Dirac sector: Cl(1,3) ~= M_2(H) and internal SU(2)
================================================================================
The peer-review of v20.0 flagged that A_F = C oplus H oplus M_3(C) was
attributed to Furey 2014, but Furey 2014 only derives SU(3) x U(1) on one
generation under SU(3)_c.  The SU(3)-commutant inside Cl(0,6) is
C oplus C oplus M_2(C) by Schur, NOT A_F.

This computation does NOT close the full A_F derivation gap.  What it DOES
verify -- in clean form, as the necessary starting point of any
substrate-based SU(2)_L derivation -- is the existence of an internal SU(2)
structure on the spinor space of the emergent 4-d Lorentzian spacetime
M = R x S^3:

  P1-P3 + Mosco convergence (Comp 4)  =>  M = R x S^3 (4-d Lorentzian)
  Cl(1,3) ~= M_2(H)                    (Cartan classification, structurally
                                        forced by the Lorentzian (1,3) signature
                                        of the emergent spacetime)
  Right-multiplication by H on the
  M_2(H) module H^2                    =>  Internal SU(2), commutes with the
                                          external Lorentz-Clifford action
                                          (Dixon 2010, p. 2)

The internal SU(2) generated by right-H-multiplication is an honest
structural feature of any Lorentzian spinor framework.  Whether this SU(2)
can be IDENTIFIED with the SU(2)_L factor of A_F is a separate question
requiring either (a) Dixon's hyperspinor framework C oplus H oplus O acting
on T (which embeds both spacetime and internal sectors into a single
algebra), or (b) Chamseddine--Connes--Marcolli (CCM)-style axiomatic
classification of finite spectral triples.  Both routes invoke additional
structural arguments beyond P1-P3 alone; the current paper records this
explicitly and does not claim a derivation of A_F = C oplus H oplus M_3(C)
solely from the three postulates.

The computation works in R^8 = H^2 (real 8-dim), where:
- Cl(1,3) acts as left-multiplication by M_2(H) (matrix entries in H acting
  via LEFT-H-multiplication on each H component of H^2);
- Right-H acts as right-H-multiplication on each H component;
- LEFT- and RIGHT-multiplication on H commute by associativity, so
  [Cl(1,3), R_q] = 0 by construction.

This is the natural setting (avoids the antilinearity confusion of the
C^4 chiral basis where right-mult by i_H mixes the C-structure).
"""
import numpy as np
import numpy.linalg as la

I4 = np.eye(4, dtype=float)


# ----------- H ~ R^4 with basis (1, i_H, j_H, k_H), LEFT and RIGHT mult ------

def L_H(q):
    """Left-multiplication by q = (q0, q1, q2, q3) on H ~ R^4: returns 4x4 R-matrix."""
    q0, q1, q2, q3 = q
    # Action on basis (1, i, j, k):
    #   q * 1 = q0 + q1 i + q2 j + q3 k             (column: (q0, q1, q2, q3))
    #   q * i = q0 i + q1 i^2 + q2 j i + q3 k i
    #         = q0 i - q1 + q2 (-k) + q3 j
    #         = -q1 + q0 i + q3 j - q2 k            (column: (-q1, q0, q3, -q2))
    #   q * j = q0 j + q1 i j + q2 j^2 + q3 k j
    #         = q0 j + q1 k - q2 + q3 (-i)
    #         = -q2 - q3 i + q0 j + q1 k            (column: (-q2, -q3, q0, q1))
    #   q * k = q0 k + q1 i k + q2 j k + q3 k^2
    #         = q0 k + q1 (-j) + q2 i - q3
    #         = -q3 + q2 i - q1 j + q0 k            (column: (-q3, q2, -q1, q0))
    return np.array([
        [q0, -q1, -q2, -q3],  # row 0 (coeff of 1)
        [q1,  q0, -q3,  q2],  # row 1 (coeff of i)
        [q2,  q3,  q0, -q1],  # row 2 (coeff of j)
        [q3, -q2,  q1,  q0],  # row 3 (coeff of k)
    ], dtype=float)


def R_H(q):
    """Right-multiplication by q on H ~ R^4: returns 4x4 R-matrix.

    Right-mult: psi -> psi * q.  Action on basis (1, i, j, k):
       1 * q = q                                          (column: (q0, q1, q2, q3))
       i * q = i q0 + i q1 i + i q2 j + i q3 k
             = q0 i - q1 + q2 k + q3 (-j)
             = -q1 + q0 i - q3 j + q2 k                   (column: (-q1, q0, -q3, q2))
       j * q = j q0 + j q1 i + j q2 j + j q3 k
             = q0 j + q1 (-k) - q2 + q3 i
             = -q2 + q3 i + q0 j - q1 k                   (column: (-q2, q3, q0, -q1))
       k * q = k q0 + k q1 i + k q2 j + k q3 k
             = q0 k + q1 j + q2 (-i) - q3
             = -q3 - q2 i + q1 j + q0 k                   (column: (-q3, -q2, q1, q0))
    """
    q0, q1, q2, q3 = q
    return np.array([
        [q0, -q1, -q2, -q3],  # row 0
        [q1,  q0,  q3, -q2],  # row 1
        [q2, -q3,  q0,  q1],  # row 2
        [q3,  q2, -q1,  q0],  # row 3
    ], dtype=float)


# ----------- Lift to H^2 ~ R^8 ----------------------------------------------

def block_diag(A, B):
    """Block-diagonal of two 4x4 matrices -> 8x8."""
    out = np.zeros((8, 8), dtype=float)
    out[:4, :4] = A
    out[4:, 4:] = B
    return out


def block_offdiag(A, B):
    """Off-diagonal block [[0, A], [B, 0]] -> 8x8."""
    out = np.zeros((8, 8), dtype=float)
    out[:4, 4:] = A
    out[4:, :4] = B
    return out


def M2H_left_action(a, b, c, d):
    """Left-action of 2x2 H matrix [[a, b], [c, d]] on H^2 ~ R^8.
    Each entry is an H element (4-tuple).
    """
    La = L_H(a); Lb = L_H(b); Lc = L_H(c); Ld = L_H(d)
    out = np.zeros((8, 8), dtype=float)
    out[:4, :4] = La
    out[:4, 4:] = Lb
    out[4:, :4] = Lc
    out[4:, 4:] = Ld
    return out


def right_H_on_H2(q):
    """Right-multiplication by q on H^2 ~ R^8 = R^4 (+) R^4.
    Acts as R_q on each H slot independently.
    """
    Rq = R_H(q)
    return block_diag(Rq, Rq)


# ----------- Dirac gamma matrices in M_2(H) ----------------------------------

ONE = (1.0, 0, 0, 0)
I_H = (0, 1.0, 0, 0)
J_H = (0, 0, 1.0, 0)
K_H = (0, 0, 0, 1.0)
ZERO = (0, 0, 0, 0)


def gamma_matrices():
    """Dirac gamma matrices for Cl(1,3) in the M_2(H) representation.
    eta = diag(+1, -1, -1, -1).

    gamma^0 = [[1, 0], [0, -1]]       (entries = H-identity, -H-identity; gives gamma^0^2 = I)
    gamma^i = [[0, e_i], [e_i, 0]]    (e_i is the H imaginary; gives gamma^i^2 = -I)

    Check: gamma^0 gamma^i = [[0, e_i], [-e_i, 0]] and
           gamma^i gamma^0 = [[0, -e_i], [e_i, 0]], so anticommutator vanishes.
    Check: gamma^i gamma^j = diag(e_i e_j, e_i e_j),
           gamma^j gamma^i = diag(e_j e_i, e_j e_i) = -diag(e_i e_j, e_i e_j),
           so anticommutator vanishes for i != j.
    """
    MINUS_ONE = (-1.0, 0, 0, 0)
    g0 = M2H_left_action(ONE, ZERO, ZERO, MINUS_ONE)
    g1 = M2H_left_action(ZERO, I_H, I_H, ZERO)
    g2 = M2H_left_action(ZERO, J_H, J_H, ZERO)
    g3 = M2H_left_action(ZERO, K_H, K_H, ZERO)
    return [g0, g1, g2, g3]


def gamma5(gammas):
    """gamma^5 = i_H gamma^0 gamma^1 gamma^2 gamma^3 -- but we're real-valued,
    so we use the proper chirality projector.  In M_2(H) Cl(1,3), the volume
    element is gamma^0 gamma^1 gamma^2 gamma^3, which satisfies (vol)^2 = -I.
    For a real chirality projector we need a square-root of +I, so
    gamma_5 := i * vol in the C-complexification; here we use vol directly
    and note that vol acts as the H-volume i_H j_H k_H = -1 (no), or actually
    let me compute.

    gamma^0 gamma^1 gamma^2 gamma^3 = [[0,1],[1,0]][[0,i],[-i,0]][[0,j],[-j,0]][[0,k],[-k,0]]
    """
    return gammas[0] @ gammas[1] @ gammas[2] @ gammas[3]


# ----------- Verification ---------------------------------------------------

def verify_clifford(gammas):
    """{gamma^mu, gamma^nu} = 2 eta^{mu nu} I"""
    eta = np.diag([1, -1, -1, -1])
    I8 = np.eye(8)
    fails = []
    for mu in range(4):
        for nu in range(4):
            anti = gammas[mu] @ gammas[nu] + gammas[nu] @ gammas[mu]
            expected = 2 * eta[mu, nu] * I8
            if la.norm(anti - expected) > 1e-12:
                fails.append((mu, nu, la.norm(anti - expected)))
    return fails


def verify_commute(gammas, R_ops):
    """All gamma^mu commute with all R_q."""
    fails = []
    for mu, g in enumerate(gammas):
        for label, R in R_ops:
            comm = g @ R - R @ g
            n = la.norm(comm)
            if n > 1e-10:
                fails.append((f"gamma^{mu}", label, n))
    return fails


def verify_quaternion(R_iH, R_jH, R_kH):
    I8 = np.eye(8)
    checks = []
    checks.append(("R_iH^2 = -I", la.norm(R_iH @ R_iH + I8)))
    checks.append(("R_jH^2 = -I", la.norm(R_jH @ R_jH + I8)))
    checks.append(("R_kH^2 = -I", la.norm(R_kH @ R_kH + I8)))
    # IMPORTANT: For RIGHT-multiplication on H, the algebra
    # composition reverses.  R_p R_q = R_qp (not R_pq) because:
    #   R_q (R_p psi) = R_q (psi * p) = (psi * p) * q = psi * (p q),
    # which by definition is R_{pq} psi.  Hmm so R_q R_p = R_{pq}.
    # That is, the algebra of right-mults is the OPPOSITE algebra H^{op}.
    # H^{op} has multiplication p *_op q := q p, so the relations are
    # i_op j_op = -k_op (sign flip).
    # Translating: R_iH R_jH psi = R_jH (R_iH psi) = R_jH (psi * i_H)
    #            = (psi i_H) j_H = psi (i_H j_H) = psi k_H = R_kH psi.
    # So R_iH R_jH = R_kH ... wait this is the SAME as standard H.
    # Let me re-derive carefully:
    #   (R_q after R_p)(psi) = R_q(R_p(psi)) = R_q(psi p) = (psi p) q = psi (p q).
    # So composition (R_q . R_p) corresponds to right-mult by (p q).
    # In matrix language: (R_q M_p) (or matrix product R_q @ R_p) corresponds to
    # right-mult by (p q).
    # So R_q @ R_p = R_{pq}.
    # Equivalently: R_p @ R_q = R_{qp}.
    # With i j = k in H: R_j @ R_i = R_{ij} = R_k.
    # So we should expect R_jH @ R_iH = R_kH (NOT R_iH @ R_jH = R_kH).
    checks.append(("R_jH R_iH = R_kH", la.norm(R_jH @ R_iH - R_kH)))
    checks.append(("R_kH R_jH = R_iH", la.norm(R_kH @ R_jH - R_iH)))
    checks.append(("R_iH R_kH = R_jH", la.norm(R_iH @ R_kH - R_jH)))
    return checks


def verify_chirality(gammas, R_ops):
    """The Lorentz pseudoscalar 'volume' gamma^0123 commutes with R_q."""
    g_vol = gammas[0] @ gammas[1] @ gammas[2] @ gammas[3]
    fails = []
    for label, R in R_ops:
        comm = g_vol @ R - R @ g_vol
        n = la.norm(comm)
        if n > 1e-10:
            fails.append((label, n))
    return fails, g_vol


def main():
    print("=" * 90)
    print("  Computation 66 -- Cl(1,3) ~= M_2(H) and the internal SU(2) from right-H")
    print("=" * 90)
    print()
    print("  Verifies the structural starting point of any spacetime-derived SU(2)_L:")
    print("  the 4-d Lorentzian Dirac sector carries an internal SU(2) from")
    print("  right-multiplication by H on its M_2(H) spinor module H^2.")
    print()

    gammas = gamma_matrices()
    R_iH = right_H_on_H2(I_H)
    R_jH = right_H_on_H2(J_H)
    R_kH = right_H_on_H2(K_H)
    R_ops = [("R_iH", R_iH), ("R_jH", R_jH), ("R_kH", R_kH)]

    # 1. Clifford relations
    print("STEP 1: Cl(1,3) Clifford relations on H^2 ~ R^8")
    fails = verify_clifford(gammas)
    if not fails:
        print("    PASS: all 16 anticommutators match diag(+1, -1, -1, -1).")
    else:
        print(f"    FAIL: {len(fails)} relations off")
        for f in fails[:3]:
            print("         ", f)
    print()

    # 2. Quaternion algebra of R_q
    print("STEP 2: Right-H on H^2 forms a quaternion algebra")
    checks = verify_quaternion(R_iH, R_jH, R_kH)
    if all(n < 1e-10 for _, n in checks):
        print("    PASS: R_q satisfy H multiplication table (with R_p R_q = R_{qp}).")
        for label, n in checks:
            print(f"          {label:25s}  residual {n:.2e}")
    else:
        print("    FAIL: quaternion relations off:")
        for label, n in checks:
            print(f"          {label:25s}  residual {n:.2e}")
    print()

    # 3. Internal property: [gamma^mu, R_q] = 0
    print("STEP 3: Right-H commutes with Cl(1,3) -- structurally INTERNAL")
    fails = verify_commute(gammas, R_ops)
    if not fails:
        print("    PASS: all 12 commutators vanish.")
        print("          Right-H is in the centraliser of Cl(1,3) ~ M_2(H) (left action).")
        print("          This is Dixon's structural starting point.")
    else:
        print(f"    FAIL: {len(fails)} commutators non-zero")
        for f in fails[:5]:
            print("         ", f)
    print()

    # 4. Generates SU(2)
    print("STEP 4: {R_q/2} are su(2) generators")
    # [R_p/2, R_q/2] = ? For right-mults: [R_p, R_q] = R_p R_q - R_q R_p = R_{qp} - R_{pq}.
    # With pq = -qp (anticommutation of distinct imaginary quats), this gives:
    # [R_p, R_q] = R_{qp} - R_{pq} = -R_{pq} - R_{pq} = -2 R_{pq}
    # so [R_p, R_q] = -2 R_{pq}.  For i j = k (and so pq = ij = k): [R_i, R_j] = -2 R_k.
    # In Lie algebra normalisation J_a = R_a/2, [J_i, J_j] = -(1/2) R_k = -J_k.
    # The minus sign just reflects that right-mult on H is the opposite algebra
    # H^{op}, and SU(2) ~ Sp(1) acts as the unit quaternions.  Up to a sign
    # convention this is su(2).
    comm_ij = R_iH @ R_jH - R_jH @ R_iH
    expected_ij = -2 * R_kH
    n_ij = la.norm(comm_ij - expected_ij)
    if n_ij < 1e-10:
        print(f"    PASS: [R_iH, R_jH] = -2 R_kH (residual {n_ij:.2e}).")
        print("          Up to sign convention (H^{op} vs H), {R_q/2} is the su(2) Lie algebra,")
        print("          and the exponentiated group is SU(2) ~ Sp(1) ~ unit quaternions.")
    else:
        print(f"    UNEXPECTED:  [R_iH, R_jH] residual from -2 R_kH: {n_ij:.2e}")
    print()

    # 5. Pseudoscalar commutes -> right-H preserves chirality sectors
    print("STEP 5: Lorentz pseudoscalar (gamma^0 gamma^1 gamma^2 gamma^3) commutes with R_q")
    fails, g_vol = verify_chirality(gammas, R_ops)
    if not fails:
        print("    PASS: g_vol commutes with each R_q.")
        print("          Right-H is chirality-preserving (acts within each chirality sector).")
    else:
        print(f"    FAIL: {fails}")
    print()

    # 6. Summary
    print("=" * 90)
    print("  Structural verification complete")
    print("=" * 90)
    print()
    print("  What is established:")
    print("  - The emergent 4-d Lorentzian spacetime M = R x S^3 (derived from")
    print("    P1-P3 via Mosco convergence, Comp 4) has local Lorentz Clifford")
    print("    algebra Cl(1,3) ~= M_2(H), a Cartan-classification fact.")
    print("  - Right-H-multiplication on the H^2 spinor module of M_2(H) is an")
    print("    internal SU(2): commutes with the Cl(1,3) (external) action,")
    print("    satisfies the quaternion algebra, generates Sp(1) ~ SU(2).")
    print("  - This internal SU(2) is chirality-preserving (the Lorentz")
    print("    pseudoscalar commutes with R_q), so it acts within each chirality")
    print("    sector.")
    print()
    print("  What is NOT established by this computation alone:")
    print("  - Identification of this spacetime-Dirac SU(2) with the SU(2)_L factor")
    print("    of A_F = C oplus H oplus M_3(C) in the Connes M x F spectral triple.")
    print("    The A_F factor lives on the F-side (internal), the right-H derived")
    print("    here lives on the M-side (spacetime).  The identification requires")
    print("    additional structural arguments -- either Dixon's hyperspinor")
    print("    framework T = C tensor H tensor O, or CCM-style classification of")
    print("    finite spectral triples.  Both invoke structural inputs beyond")
    print("    P1-P3 alone, and the v20.0 abstract's claim 'No Chamseddine-Connes-")
    print("    Marcolli classification is invoked' is therefore tightened: PST")
    print("    derives the substrate's Cl(0,6) chain algebra structure and the")
    print("    Cl(1,3) spacetime Dirac structure from P1-P3, but the specific")
    print("    decomposition A_F = C oplus H oplus M_3(C) on the F-side rests on")
    print("    CCM (or equivalent) classification arguments.")


if __name__ == "__main__":
    main()
