# ============================================================================
# ULTIMATE LEECH LATTICE INTERTWINER HESSIAN CERTIFICATE
# ============================================================================
#
# Purpose
# -------
# Reproducible computational study of:
#
#   1. Binary Golay construction from the quadratic-residue seed mod 11.
#   2. Construction of the 196,560-vector minimal Leech shell.
#   3. Exact rational shell geometry.
#   4. Exact 12-dimensional intertwiner action.
#   5. Exact rational spectrum.
#   6. Explicit kernel / rotational-mode certification.
#   7. Parity involution.
#   8. Multi-realization invariance tests.
#   9. Exact E8^3 transverse instability calculation.
#  10. Numerical Niemeier Morse experiments.
#  11. Phenomenological double-slit demonstration.
#
# IMPORTANT STATUS CONVENTION
# ---------------------------
#
# [EXACT]
#   Integer/rational symbolic computations over Q.
#
# [NUMERICAL]
#   Floating-point diagonalization or threshold-based computation.
#
# [MODEL]
#   Illustrative physical/phenomenological interpretation.
#
# The program deliberately does NOT claim that every computational
# observation is a theorem.  In particular, numerical Morse indices and
# the double-slit experiment are explicitly separated from the exact
# algebraic certificate.
#
# Requirements
# ------------
#   Python >= 3.10
#   numpy
#   sympy
#
# ============================================================================

from __future__ import annotations

import hashlib
import json
import math
import platform
import sys
import time

from dataclasses import dataclass
from itertools import combinations
from typing import Dict, List, Tuple, Iterable, Optional

import numpy as np
import sympy as sp


# ============================================================================
# CONFIGURATION
# ============================================================================

@dataclass(frozen=True)
class Config:
    # Mathematical construction parameters
    prime: int = 11
    dimension: int = 24
    shell_norm_squared: int = 32

    # Pair potential:
    #
    #       F(U) = U^(-2)
    #
    # therefore
    #
    #       F'(U)  = -2 U^(-3)
    #       F''(U) =  6 U^(-4)
    #
    potential_power: int = 2

    # Exact invariance realizations
    number_of_realizations: int = 4

    # Deterministic choices used in the invariance experiment.
    #
    # The actual representatives are selected by position in their shell,
    # rather than being encoded geometrically.
    representative_modes: Tuple[str, ...] = (
        "first",
        "last",
        "middle",
        "offset",
    )

    # Numerical tolerances
    numerical_zero_tolerance: float = 1e-8

    # Reproducibility
    random_seed: int = 20260924

    # Run optional expensive sections
    run_morse_experiment: bool = True
    run_double_slit_model: bool = True


CFG = Config()

T0 = time.time()


# ============================================================================
# DISPLAY / CERTIFICATE HELPERS
# ============================================================================

def banner(title: str) -> None:
    print()
    print("=" * 80)
    print(title)
    print("=" * 80)


def section(title: str) -> float:
    print()
    print(f"[{title}]")
    return time.time()


def elapsed(t: float) -> str:
    return f"{time.time() - t:.2f}s"


def exact_assert(condition: bool, message: str) -> None:
    if not condition:
        raise AssertionError("[EXACT CERTIFICATE FAILURE] " + message)


def numerical_assert(condition: bool, message: str) -> None:
    if not condition:
        raise AssertionError("[NUMERICAL CHECK FAILURE] " + message)


def hash_array(X: np.ndarray) -> str:
    h = hashlib.sha256()
    h.update(np.ascontiguousarray(X).tobytes())
    return h.hexdigest()


# ============================================================================
# SECTION 1
# BINARY GOLAY CODE
# ============================================================================

def construct_golay_code(cfg: Config):
    """
    Construct the extended binary Golay code using the quadratic-residue
    circulant construction modulo 11.
    """

    p = cfg.prime

    qr = {
        (x * x) % p
        for x in range(1, p)
    }

    qnr = set(range(1, p)) - qr

    v_circ = [
        1 if (i == 0 or i in qnr) else 0
        for i in range(p)
    ]

    B24 = [[0] * 12 for _ in range(12)]

    for j in range(1, 12):
        B24[0][j] = 1

    for i in range(p):
        B24[i + 1][0] = 1

        for j in range(p):
            B24[i + 1][j + 1] = \
                v_circ[(j - i) % p]

    G = np.zeros((12, 24), dtype=np.uint8)

    for i in range(12):
        G[i, i] = 1

        for j in range(12):
            G[i, 12 + j] = B24[i][j]

    # Enumerate all 2^12 codewords.
    codewords = []

    for mask in range(1 << 12):
        cw = np.zeros(24, dtype=np.uint8)

        for i in range(12):
            if (mask >> i) & 1:
                cw ^= G[i]

        codewords.append(cw)

    codewords = np.asarray(codewords, dtype=np.uint8)

    weights = np.sum(codewords, axis=1)

    unique_weights, counts = np.unique(
        weights,
        return_counts=True
    )

    weight_distribution = dict(
        zip(
            map(int, unique_weights),
            map(int, counts)
        )
    )

    octads = [
        cw
        for cw in codewords
        if int(np.sum(cw)) == 8
    ]

    return {
        "G": G,
        "codewords": codewords,
        "weights": weights,
        "weight_distribution": weight_distribution,
        "octads": octads,
    }


def certify_golay(golay) -> None:
    codewords = golay["codewords"]
    wd = golay["weight_distribution"]

    exact_assert(
        len(codewords) == 4096,
        "Golay code must contain 4096 codewords."
    )

    exact_assert(
        wd.get(8) == 759,
        "Expected 759 octads."
    )

    exact_assert(
        wd.get(12) == 2576,
        "Expected 2576 dodecads."
    )

    expected = {
        0: 1,
        8: 759,
        12: 2576,
        16: 759,
        24: 1,
    }

    exact_assert(
        wd == expected,
        f"Unexpected weight distribution: {wd}"
    )

    # Generator rank over F_2.
    G = golay["G"].copy()

    A = G.astype(np.uint8).copy()
    rank = 0

    rows, cols = A.shape

    for col in range(cols):
        pivot = None

        for r in range(rank, rows):
            if A[r, col]:
                pivot = r
                break

        if pivot is None:
            continue

        A[[rank, pivot]] = A[[pivot, rank]]

        for r in range(rows):
            if r != rank and A[r, col]:
                A[r] ^= A[rank]

        rank += 1

        if rank == rows:
            break

    exact_assert(
        rank == 12,
        f"Golay generator rank is {rank}, expected 12."
    )


# ============================================================================
# SECTION 2
# MINIMAL LEECH SHELL
# ============================================================================

def construct_leech_shell(golay, cfg: Config) -> np.ndarray:

    octads = golay["octads"]
    codewords = golay["codewords"]

    X_list = []

    # ------------------------------------------------------------------------
    # Family A
    #
    # (±4, ±4, 0^22)
    #
    # Count:
    #
    #   C(24,2) * 4 = 1104
    # ------------------------------------------------------------------------

    for i, j in combinations(range(24), 2):

        for s1 in (4, -4):
            for s2 in (4, -4):

                v = [0] * 24
                v[i] = s1
                v[j] = s2

                X_list.append(tuple(v))

    # ------------------------------------------------------------------------
    # Family B
    #
    # (±2^8,0^16)
    #
    # over octads, with even sign parity.
    #
    # 759 * 128 = 97152
    # ------------------------------------------------------------------------

    for octad in octads:

        idx = [
            i
            for i, b in enumerate(octad)
            if b
        ]

        for mask in range(1 << 8):

            if mask.bit_count() & 1:
                continue

            v = [0] * 24

            for k in range(8):
                v[idx[k]] = (
                    -2
                    if ((mask >> k) & 1)
                    else 2
                )

            X_list.append(tuple(v))

    # ------------------------------------------------------------------------
    # Family C
    #
    # (±3, ±1^23)
    #
    # ------------------------------------------------------------------------

    for m in range(24):

        for cw in codewords:

            v = [
                1 - 2 * int(b)
                for b in cw
            ]

            v[m] *= -3

            X_list.append(tuple(v))

    X = np.asarray(
        X_list,
        dtype=np.int64
    )

    return X


def certify_leech_shell(X: np.ndarray, cfg: Config) -> None:

    exact_assert(
        X.shape == (196560, 24),
        f"Unexpected shell shape: {X.shape}"
    )

    norms = np.sum(X * X, axis=1)

    exact_assert(
        np.all(norms == cfg.shell_norm_squared),
        "Not every constructed vector has norm squared 32."
    )

    # Detect accidental duplicate vectors.
    unique_count = len(
        np.unique(X, axis=0)
    )

    exact_assert(
        unique_count == len(X),
        "Leech shell contains duplicate vectors."
    )


# ============================================================================
# SECTION 3
# SHELL GEOMETRY
# ============================================================================

def shell_decomposition(
    X: np.ndarray,
    u: np.ndarray,
    cfg: Config
):
    inner = X @ u

    shell_values = [
        32,
        16,
        8,
        0,
        -8,
        -16,
        -32,
    ]

    shell_indices = {
        s: np.where(inner == s)[0]
        for s in shell_values
    }

    shell_sizes = {
        s: len(shell_indices[s])
        for s in shell_values
    }

    shell_data = {
        s: X[shell_indices[s]]
        for s in shell_values
    }

    return (
        inner,
        shell_values,
        shell_indices,
        shell_sizes,
        shell_data,
    )


def certify_shell_distribution(
    shell_sizes: Dict[int, int]
):

    expected = {
        32: 1,
        16: 4600,
        8: 47104,
        0: 93150,
        -8: 47104,
        -16: 4600,
        -32: 1,
    }

    exact_assert(
        shell_sizes == expected,
        f"Unexpected shell distribution: {shell_sizes}"
    )


# ============================================================================
# SECTION 4
# TANGENT BASIS
# ============================================================================

def tangent_basis(
    u: np.ndarray,
    cfg: Config,
    column_order: Optional[Iterable[int]] = None
):

    R2 = cfg.shell_norm_squared
    D = cfg.dimension

    if column_order is None:
        column_order = range(D)

    candidates = []

    for j in column_order:

        e = np.zeros(D, dtype=np.int64)
        e[j] = 1

        t = (
            R2 * e
            - int(u[j]) * u
        )

        if np.any(t):
            candidates.append(t)

    selected = []

    for r in candidates:

        trial = np.array(
            selected + [r]
        )

        if np.linalg.matrix_rank(trial) > len(selected):
            selected.append(r)

        if len(selected) == D - 1:
            break

    exact_assert(
        len(selected) == D - 1,
        "Could not construct a 23-dimensional tangent basis."
    )

    return selected


# ============================================================================
# SECTION 5
# INTERTWINER DATA
# ============================================================================

FIELD_SPECS = [
    (32, "W1"),

    (16, "W1"),
    (16, "W2"),

    (8, "W1"),
    (8, "W2"),

    (0, "W1"),
    (0, "W2"),

    (-8, "W1"),
    (-8, "W2"),

    (-16, "W1"),
    (-16, "W2"),

    (-32, "W1"),
]


def make_shell_fields():
    result = {}

    for a, (s, kind) in enumerate(FIELD_SPECS):
        result.setdefault(s, []).append(a)

    return result


SHELL_FIELDS = make_shell_fields()


def intertwiner_value(
    field_idx: int,
    z: sp.Matrix,
    v: np.ndarray,
    u: np.ndarray,
    cfg: Config
):

    R2 = cfg.shell_norm_squared

    target_s, kind = FIELD_SPECS[field_idx]

    uv = int(u @ v)

    if uv != target_s:
        return sp.zeros(cfg.dimension, 1)

    v_sp = sp.Matrix(
        v.tolist()
    )

    u_sp = sp.Matrix(
        u.tolist()
    )

    zv = (z.T * v_sp)[0]

    if kind == "W1":

        return (
            z
            - sp.Rational(zv, R2) * v_sp
        )

    if kind == "W2":

        uv_sp = (
            u_sp.T * v_sp
        )[0]

        return (
            zv
            * (
                u_sp
                - sp.Rational(
                    uv_sp,
                    R2
                ) * v_sp
            )
        )

    raise ValueError(
        f"Unknown intertwiner kind: {kind}"
    )


# ============================================================================
# SECTION 6
# DYNAMIC WARD SELF-ENERGY
# ============================================================================

def fp(U):
    return -sp.Rational(2, U**3)


def fpp(U):
    return sp.Rational(6, U**4)


def derive_lambda_S(
    shell_sizes: Dict[int, int],
    cfg: Config
):

    R2 = cfg.shell_norm_squared
    D = cfg.dimension

    K_rot_sum = sp.Rational(0)

    for s in [
        16,
        8,
        0,
        -8,
        -16,
        -32,
    ]:

        count = shell_sizes[s]

        U_val = (
            2 * R2
            - 2 * s
        )

        sum_W_coeff = (
            -s * count
        )

        sum_diff_coeff = (
            -sp.Rational(
                R2 * count,
                D - 1
            )
            * (
                R2
                - sp.Rational(
                    s * s,
                    R2
                )
            )
        )

        K_rot_sum += (
            2 * fp(U_val)
            * sum_W_coeff
            +
            4 * fpp(U_val)
            * sum_diff_coeff
        )

    lambda_S = (
        K_rot_sum / (-R2)
    )

    return sp.factor(lambda_S)


# ============================================================================
# SECTION 7
# ACTION MATRIX CONSTRUCTION
# ============================================================================

def representative_index(
    size: int,
    mode: str
) -> int:

    if size <= 0:
        raise ValueError("Empty shell.")

    if mode == "first":
        return 0

    if mode == "last":
        return size - 1

    if mode == "middle":
        return size // 2

    if mode == "offset":
        return min(
            17,
            size - 1
        )

    raise ValueError(
        f"Unknown representative mode: {mode}"
    )


def construct_action_matrix(
    X: np.ndarray,
    u: np.ndarray,
    representative_mode: str,
    tangent_order: Optional[Iterable[int]] = None,
    measurement_reverse: bool = False,
    cfg: Config = CFG,
):

    D = cfg.dimension
    R2 = cfg.shell_norm_squared

    u_sp = sp.Matrix(
        u.tolist()
    )

    # ------------------------------------------------------------------------
    # Tangent basis
    # ------------------------------------------------------------------------

    T_np_list = tangent_basis(
        u,
        cfg,
        tangent_order
    )

    T = [
        sp.Matrix(v.tolist())
        for v in T_np_list
    ]

    T_np = np.asarray(
        T_np_list,
        dtype=np.int64
    )

    # ------------------------------------------------------------------------
    # Shell decomposition
    # ------------------------------------------------------------------------

    (
        inner,
        shell_values,
        shell_indices,
        shell_sizes,
        shell_data,
    ) = shell_decomposition(
        X,
        u,
        cfg
    )

    rep = {}

    for s in shell_values:

        idx = representative_index(
            shell_sizes[s],
            representative_mode
        )

        if representative_mode == "last":
            idx = shell_sizes[s] - 1

        elif representative_mode == "middle":
            idx = shell_sizes[s] // 2

        elif representative_mode == "offset":
            idx = min(
                17,
                shell_sizes[s] - 1
            )

        rep[s] = X[
            shell_indices[s][idx]
        ]

    # ------------------------------------------------------------------------
    # Evaluation matrix
    # ------------------------------------------------------------------------

    chosen_measurements = []
    exact_rows = []

    for s, f_list in SHELL_FIELDS.items():

        v = rep[s]

        sub_rows = []

        tangent_indices = list(
            range(D - 1)
        )

        if measurement_reverse:
            tangent_indices = tangent_indices[::-1]

        for r in tangent_indices:

            z = T[r]

            q_indices = list(
                range(D)
            )

            if measurement_reverse:
                q_indices = q_indices[::-1]

            for qidx in q_indices:

                sub_row = [
                    intertwiner_value(
                        a,
                        z,
                        v,
                        u,
                        cfg
                    )[qidx]
                    for a in f_list
                ]

                if all(
                    val == 0
                    for val in sub_row
                ):
                    continue

                trial = sp.Matrix(
                    sub_rows + [sub_row]
                )

                if trial.rank() > len(sub_rows):

                    sub_rows.append(
                        sub_row
                    )

                    chosen_measurements.append(
                        (
                            r,
                            s,
                            qidx
                        )
                    )

                    full_row = [
                        intertwiner_value(
                            a,
                            z,
                            v,
                            u,
                            cfg
                        )[qidx]
                        for a in range(12)
                    ]

                    exact_rows.append(
                        full_row
                    )

                    if len(sub_rows) == len(f_list):
                        break

            if len(sub_rows) == len(f_list):
                break

        exact_assert(
            len(sub_rows) == len(f_list),
            f"Could not span shell field group s={s}."
        )

    Eval = sp.Matrix(
        exact_rows
    )

    exact_assert(
        Eval.shape == (12, 12),
        f"Evaluation matrix has shape {Eval.shape}."
    )

    exact_assert(
        Eval.rank() == 12,
        "Evaluation matrix is singular."
    )

    Eval_inv = Eval.inv()

    # ------------------------------------------------------------------------
    # Exact Ward multiplier
    # ------------------------------------------------------------------------

    lambda_S = derive_lambda_S(
        shell_sizes,
        cfg
    )

    # ------------------------------------------------------------------------
    # Group measurement requests
    # ------------------------------------------------------------------------

    grouped = {}

    for m_idx, (r, s, q) in enumerate(
        chosen_measurements
    ):

        grouped.setdefault(
            (r, s),
            []
        ).append(
            (m_idx, q)
        )

    # ------------------------------------------------------------------------
    # Hessian action
    # ------------------------------------------------------------------------

    def fast_hessian_action(
        field_idx,
        r_idx,
        x_shell
    ):

        z_np = T_np[r_idx]
        x_np = rep[x_shell]

        shell_val, kind = \
            FIELD_SPECS[field_idx]

        Y = shell_data[shell_val]

        W_x = intertwiner_value(
            field_idx,
            T[r_idx],
            x_np,
            u,
            cfg
        )

        z_dot_y = Y @ z_np
        u_dot_y = Y @ u

        if kind == "W1":

            W_Y_scaled = (
                R2 * z_np[None, :]
                -
                z_dot_y[:, None] * Y
            )

        else:

            W_Y_scaled = (
                z_dot_y[:, None]
                *
                (
                    R2 * u[None, :]
                    -
                    u_dot_y[:, None] * Y
                )
            )

        xy = Y @ x_np

        U = (
            2 * R2
            - 2 * xy
        )

        valid = (
            U != 0
        )

        U_val = U[valid]
        Y_val = Y[valid]
        W_val = W_Y_scaled[valid]

        x_dot_Wy = (
            W_val @ x_np
        )

        K_total_sp = sp.zeros(
            D,
            1
        )

        for U_lvl in np.unique(U_val):

            mask = (
                U_val == U_lvl
            )

            sum_W = (
                sp.Matrix(
                    np.sum(
                        W_val[mask],
                        axis=0,
                        dtype=np.int64
                    ).tolist()
                )
                / R2
            )

            sum_diff = (
                sp.Matrix(
                    (
                        np.sum(
                            x_dot_Wy[mask],
                            dtype=np.int64
                        )
                        * x_np
                        -
                        x_dot_Wy[mask]
                        @ Y_val[mask]
                    ).tolist()
                )
                / R2
            )

            Uq = sp.Integer(
                int(U_lvl)
            )

            K_total_sp += (
                2 * fp(Uq) * sum_W
                +
                4 * fpp(Uq) * sum_diff
            )

        x_sp = sp.Matrix(
            x_np.tolist()
        )

        Ktan = (
            K_total_sp
            -
            sp.Rational(
                1,
                R2
            )
            *
            (
                x_sp.T
                * K_total_sp
            )[0]
            *
            x_sp
        )

        return (
            lambda_S * W_x
            -
            Ktan
        )

    # ------------------------------------------------------------------------
    # Assemble A
    # ------------------------------------------------------------------------

    hessian_cache = {}

    A = sp.zeros(
        12,
        12
    )

    for b in range(12):

        rhs_values = [
            None
        ] * len(
            chosen_measurements
        )

        for (r, s), requests in grouped.items():

            key = (
                b,
                r,
                s
            )

            if key not in hessian_cache:

                hessian_cache[key] = \
                    fast_hessian_action(
                        b,
                        r,
                        s
                    )

            Hv = hessian_cache[key]

            for m_idx, q in requests:

                rhs_values[m_idx] = \
                    Hv[q]

        coeff = (
            Eval_inv
            * sp.Matrix(rhs_values)
        )

        for a in range(12):

            A[a, b] = sp.factor(
                coeff[a]
            )

    return {
        "A": A,
        "lambda_S": lambda_S,
        "T": T,
        "T_np": T_np,
        "rep": rep,
        "chosen_measurements":
            chosen_measurements,
        "Eval": Eval,
        "Eval_inv": Eval_inv,
        "shell_sizes": shell_sizes,
    }


# ============================================================================
# SECTION 8
# EXACT SPECTRAL CERTIFICATE
# ============================================================================

def exact_spectral_certificate(A):

    eigenvalues_dict = A.eigenvals()

    eigenvalues = sorted(
        list(eigenvalues_dict.keys()),
        key=lambda x: float(x)
    )

    exact_assert(
        len(eigenvalues) == 12,
        "Expected 12 distinct eigenvalues."
    )

    multiplicities = [
        eigenvalues_dict[x]
        for x in eigenvalues
    ]

    exact_assert(
        all(m == 1 for m in multiplicities),
        "Spectrum is not simple."
    )

    exact_assert(
        eigenvalues[0] == 0,
        "Expected a zero eigenvalue."
    )

    exact_assert(
        all(
            x > 0
            for x in eigenvalues[1:]
        ),
        "Not all nonzero eigenvalues are positive."
    )

    return eigenvalues


# ============================================================================
# SECTION 9
# ROTATIONAL KERNEL CERTIFICATE
# ============================================================================

def rotational_kernel_certificate(A):

    # The rotational vector in the chosen intertwiner coordinates.
    #
    # This is a derived coordinate expression for
    #
    #       Phi_rot,z(v) = (u wedge z)v
    #
    # in the particular field ordering used here.
    c_rot = sp.Matrix([
        -32,
        -16,
        1,
        -8,
        1,
        0,
        1,
        8,
        1,
        16,
        1,
        32,
    ])

    exact_assert(
        A * c_rot == sp.zeros(12, 1),
        "Rotational coefficient vector is not in ker(A)."
    )

    nullspace = A.nullspace()

    exact_assert(
        len(nullspace) == 1,
        "Kernel dimension is not one."
    )

    v = nullspace[0]

    pivot = next(
        x
        for x in c_rot
        if x != 0
    )

    scale = sp.simplify(
        v[0] / c_rot[0]
    )

    exact_assert(
        v == scale * c_rot,
        "Kernel vector does not match rotational generator."
    )

    return c_rot


# ============================================================================
# SECTION 10
# PARITY CERTIFICATE
# ============================================================================

def construct_parity_operator(data, u):

    A = data["A"]
    T = data["T"]
    rep = data["rep"]
    chosen = data["chosen_measurements"]
    Eval_inv = data["Eval_inv"]

    P = sp.zeros(
        12,
        12
    )

    for b in range(12):

        rhs = []

        for r, s, q in chosen:

            val = intertwiner_value(
                b,
                T[r],
                -rep[s],
                u,
                CFG
            )[q]

            rhs.append(val)

        coeff = (
            Eval_inv
            * sp.Matrix(rhs)
        )

        for a in range(12):
            P[a, b] = sp.factor(
                coeff[a]
            )

    exact_assert(
        P * P == sp.eye(12),
        "Parity operator does not satisfy P^2=I."
    )

    exact_assert(
        P * A == A * P,
        "Parity does not commute with A."
    )

    return P


def parity_spectrum(P, A):

    eigs = exact_spectral_certificate(A)

    plus = 0
    minus = 0

    for lam in eigs:

        v = (
            A
            - lam * sp.eye(12)
        ).nullspace()[0]

        Pv = P * v

        pivot = next(
            i
            for i in range(12)
            if v[i] != 0
        )

        sign = sp.simplify(
            Pv[pivot]
            / v[pivot]
        )

        exact_assert(
            sign in (1, -1),
            f"Unexpected parity eigenvalue {sign}."
        )

        if sign == 1:
            plus += 1
        else:
            minus += 1

    exact_assert(
        plus == 6 and minus == 6,
        "Expected six even and six odd modes."
    )

    return plus, minus


# ============================================================================
# SECTION 11
# EXACT SIMILARITY BETWEEN TWO REALIZATIONS
# ============================================================================

def exact_similarity(
    A1,
    A2
):

    eigs = exact_spectral_certificate(
        A1
    )

    eigs2 = exact_spectral_certificate(
        A2
    )

    exact_assert(
        eigs == eigs2,
        "Exact spectra differ."
    )

    V_cols = []
    W_cols = []

    for lam in eigs:

        v = (
            A1
            - lam * sp.eye(12)
        ).nullspace()[0]

        w = (
            A2
            - lam * sp.eye(12)
        ).nullspace()[0]

        # Normalize deterministically.
        vp = next(
            x for x in v
            if x != 0
        )

        wp = next(
            x for x in w
            if x != 0
        )

        v = v / vp
        w = w / wp

        V_cols.append(v)
        W_cols.append(w)

    V = sp.Matrix.hstack(
        *V_cols
    )

    W = sp.Matrix.hstack(
        *W_cols
    )

    exact_assert(
        V.det() != 0,
        "Eigenvector matrix V is singular."
    )

    exact_assert(
        W.det() != 0,
        "Eigenvector matrix W is singular."
    )

    S = W * V.inv()

    exact_assert(
        S.det() != 0,
        "Similarity matrix is singular."
    )

    exact_assert(
        A2 * S == S * A1,
        "A2 S != S A1."
    )

    exact_assert(
        A2 == S * A1 * S.inv(),
        "Exact conjugacy failed."
    )

    return S


# ============================================================================
# SECTION 12
# MULTI-REALIZATION INVARIANCE TEST
# ============================================================================

def choose_reference_vectors(
    X: np.ndarray
):

    # These are deliberately geometrically different representatives:
    #
    #   Family A
    #   Family B
    #   Family C
    #
    # We locate them by structural properties rather than relying solely
    # on absolute array positions.

    norms = np.sum(
        X * X,
        axis=1
    )

    candidates = []

    # Family A: exactly two nonzero entries, magnitude 4.
    support2 = np.sum(
        X != 0,
        axis=1
    ) == 2

    A_indices = np.where(
        support2
        & np.all(
            np.abs(X) <= 4,
            axis=1
        )
    )[0]

    candidates.append(
        X[A_indices[0]]
    )

    # Family B: exactly eight nonzero entries, magnitude 2.
    support8 = np.sum(
        X != 0,
        axis=1
    ) == 8

    B_indices = np.where(
        support8
        & np.all(
            np.abs(X) <= 2,
            axis=1
        )
    )[0]

    candidates.append(
        X[B_indices[len(B_indices) // 3]]
    )

    # Family C: exactly 24 nonzero entries, one of magnitude 3.
    support24 = np.sum(
        X != 0,
        axis=1
    ) == 24

    C_indices = np.where(
        support24
    )[0]

    candidates.append(
        X[C_indices[len(C_indices) // 2]]
    )

    return candidates


def run_invariance_experiment(
    X: np.ndarray,
    cfg: Config
):

    references = choose_reference_vectors(
        X
    )

    realizations = []

    for i in range(
        min(
            cfg.number_of_realizations,
            len(references)
        )
    ):

        u = references[i]

        mode = cfg.representative_modes[
            i % len(
                cfg.representative_modes
            )
        ]

        # Reverse measurement order on alternate runs.
        reverse = bool(i % 2)

        tangent_order = (
            list(range(24))
            if i % 2 == 0
            else list(range(23, -1, -1))
        )

        data = construct_action_matrix(
            X,
            u,
            representative_mode=mode,
            tangent_order=tangent_order,
            measurement_reverse=reverse,
            cfg=cfg
        )

        eigs = exact_spectral_certificate(
            data["A"]
        )

        realizations.append(
            {
                "u": u,
                "mode": mode,
                "data": data,
                "eigs": eigs,
            }
        )

    # ------------------------------------------------------------------------
    # Compare every realization to the first.
    # ------------------------------------------------------------------------

    base = realizations[0]

    for i, r in enumerate(
        realizations[1:],
        start=1
    ):

        exact_assert(
            r["eigs"] == base["eigs"],
            f"Exact spectrum differs in realization {i}."
        )

        S = exact_similarity(
            base["data"]["A"],
            r["data"]["A"]
        )

        print(
            f"      realization 0 -> {i}: "
            f"exact similarity det(S) = "
            f"{sp.factor(S.det())}"
        )

    return realizations


# ============================================================================
# SECTION 13
# FUNCTIONAL DETERMINANT
# ============================================================================

def determinant_certificate(eigs):

    nonzero = eigs[1:]

    det_prime = sp.factor(
        sp.prod(nonzero)
    )

    numerator_primes = sp.factorint(
        det_prime.p
    )

    denominator_primes = sp.factorint(
        det_prime.q
    )

    return (
        det_prime,
        numerator_primes,
        denominator_primes,
    )


# ============================================================================
# SECTION 14
# EXACT E8^3 TRANSVERSE INSTABILITY
# ============================================================================

def exact_E8_transverse_certificate():

    # For an E8 root in R^24:
    #
    #   s = 16 : 56
    #   s =  0 : 126
    #   s = -16 : 56
    #   s = -32 : 1
    #
    # Two other E8 blocks contribute:
    #
    #   480 orthogonal roots.
    #
    # This section is an exact rational calculation.

    e8_shells = {
        16: 56,
        0: 126,
        -16: 56,
        -32: 1,
    }

    diag_trans = sp.Rational(0)
    lam_total = sp.Rational(0)

    for s, count in e8_shells.items():

        U = 64 - 2 * s

        fprime = fp(U)

        diag_trans += (
            count
            * 2
            * fprime
        )

        lam_total += (
            count
            * 2
            * fprime
            * sp.Rational(
                32 - s,
                32
            )
        )

    cross_fp = fp(64)

    diag_trans += (
        480
        * 2
        * cross_fp
    )

    lam_total += (
        480
        * 2
        * cross_fp
        * sp.Rational(
            32,
            32
        )
    )

    H = sp.factor(
        diag_trans
        - lam_total
    )

    exact_assert(
        H < 0,
        "E8^3 transverse curvature is not negative."
    )

    return H


# ============================================================================
# SECTION 15
# NUMERICAL NİEMEIER EXPERIMENT
# ============================================================================

def make_Ak_roots(k):

    eye = np.eye(
        k + 1
    )

    raw = [
        4.0 * (
            eye[i] - eye[j]
        )
        for i in range(k + 1)
        for j in range(k + 1)
        if i != j
    ]

    n = (
        np.ones(k + 1)
        / np.sqrt(k + 1)
    )

    basis = []

    for i in range(k):

        v = (
            eye[i]
            - np.dot(
                eye[i],
                n
            ) * n
        )

        for b in basis:

            v -= (
                np.dot(v, b)
                * b
            )

        v /= np.linalg.norm(v)

        basis.append(v)

    return (
        np.array(raw)
        @ np.array(basis).T
    )


def make_Ak_m(k):

    m = 24 // k

    roots = make_Ak_roots(k)

    result = []

    for block in range(m):

        for root in roots:

            result.append(
                np.pad(
                    root,
                    (
                        block * k,
                        24
                        - (block + 1) * k
                    )
                )
            )

    return np.asarray(
        result
    )


def numerical_morse_index(
    X_test,
    cfg: Config
):

    N, D = X_test.shape

    H_test = np.zeros(
        (
            N,
            D,
            N,
            D
        )
    )

    for i in range(N):

        diag = np.zeros(
            (D, D)
        )

        for j in range(N):

            if i == j:
                continue

            diff = (
                X_test[i]
                - X_test[j]
            )

            uv = np.dot(
                X_test[i],
                X_test[j]
            )

            U = (
                64.0
                - 2.0 * uv
            )

            fp_float = (
                -2.0 / U**3
            )

            fpp_float = (
                6.0 / U**4
            )

            block = -(
                2.0 * fp_float
                * np.eye(D)
                +
                4.0 * fpp_float
                * np.outer(
                    diff,
                    diff
                )
            )

            H_test[
                i, :, j, :
            ] = block

            diag -= block

        lam_diag = sum(
            2.0
            * (
                -2.0
                /
                (
                    64.0
                    - 2.0
                    * np.dot(
                        X_test[i],
                        X_test[j]
                    )
                )**3
            )
            * (
                32.0
                -
                np.dot(
                    X_test[i],
                    X_test[j]
                )
            )
            for j in range(N)
            if i != j
        ) / 32.0

        H_test[
            i, :, i, :
        ] = (
            diag
            -
            lam_diag * np.eye(D)
        )

    H = H_test.reshape(
        N * D,
        N * D
    )

    # Tangent projector.
    Ptan = np.zeros_like(H)

    for i in range(N):

        x = X_test[i]

        Pblock = (
            np.eye(D)
            -
            np.outer(x, x)
            / 32.0
        )

        Ptan[
            i * D:(i + 1) * D,
            i * D:(i + 1) * D
        ] = Pblock

    projected = (
        Ptan
        @ H
        @ Ptan
    )

    symmetric = (
        projected
        +
        projected.T
    ) / 2.0

    eigenvalues = np.linalg.eigvalsh(
        symmetric
    )

    negative = int(
        np.sum(
            eigenvalues
            < -cfg.numerical_zero_tolerance
        )
    )

    zero = int(
        np.sum(
            np.abs(eigenvalues)
            < cfg.numerical_zero_tolerance
        )
    )

    rotational_zero_count = (
        zero - N
    )

    return (
        rotational_zero_count,
        negative,
        eigenvalues
    )


# ============================================================================
# SECTION 16
# PHENOMENOLOGICAL DOUBLE-SLIT MODEL
# ============================================================================

def double_slit_model():

    y = np.linspace(
        -30,
        30,
        300
    )

    slit_distance = 8.0
    sigma = 4.0
    wave_number = 5.0
    L = 40.0

    env1 = np.exp(
        -(
            y - slit_distance
        )**2
        /
        (2 * sigma**2)
    )

    env2 = np.exp(
        -(
            y + slit_distance
        )**2
        /
        (2 * sigma**2)
    )

    r1 = np.sqrt(
        L**2
        +
        (
            y - slit_distance
        )**2
    )

    r2 = np.sqrt(
        L**2
        +
        (
            y + slit_distance
        )**2
    )

    psi = (
        np.sqrt(env1)
        * np.exp(
            1j * wave_number * r1
        )
        +
        np.sqrt(env2)
        * np.exp(
            1j * wave_number * r2
        )
    )

    coherent = np.abs(psi)**2

    incoherent = (
        env1 + env2
    )

    coherent_peaks = int(
        np.sum(
            (
                coherent[1:-1]
                >
                coherent[:-2]
            )
            &
            (
                coherent[1:-1]
                >
                coherent[2:]
            )
        )
    )

    incoherent_peaks = int(
        np.sum(
            (
                incoherent[1:-1]
                >
                incoherent[:-2]
            )
            &
            (
                incoherent[1:-1]
                >
                incoherent[2:]
            )
        )
    )

    return (
        coherent_peaks,
        incoherent_peaks
    )


# ============================================================================
# SECTION 17
# MACHINE-READABLE CERTIFICATE
# ============================================================================

def rational_string(x):
    return str(
        sp.factor(x)
    )


def build_certificate(
    golay,
    X,
    lambda_S,
    eigs,
    P,
    determinant,
    H_trans,
    realizations,
):

    return {
        "environment": {
            "python": sys.version,
            "platform": platform.platform(),
            "numpy": np.__version__,
            "sympy": sp.__version__,
        },

        "configuration": {
            "prime": CFG.prime,
            "dimension": CFG.dimension,
            "shell_norm_squared":
                CFG.shell_norm_squared,
            "potential":
                "F(U)=U^(-2)",
            "seed": CFG.random_seed,
        },

        "golay": {
            "codewords":
                len(golay["codewords"]),
            "weight_distribution":
                golay["weight_distribution"],
            "octads":
                len(golay["octads"]),
        },

        "leech_shell": {
            "vectors":
                int(len(X)),
            "dimension":
                int(X.shape[1]),
            "norm_squared":
                CFG.shell_norm_squared,
            "sha256":
                hash_array(X),
        },

        "ward_identity": {
            "lambda_S":
                rational_string(lambda_S),
        },

        "spectrum": {
            "eigenvalues":
                [
                    rational_string(x)
                    for x in eigs
                ],
            "rank":
                int(
                    12
                    -
                    len(
                        (sp.Matrix(
                            np.array(
                                []
                            )
                        )).nullspace()
                    )
                )
                if False else 11,
        },

        "parity": {
            "P_squared_equals_I":
                True,
            "commutes_with_A":
                True,
            "trace":
                int(P.trace()),
        },

        "det_prime": {
            "value":
                rational_string(determinant[0]),
            "numerator_primes":
                {
                    str(k): int(v)
                    for k, v
                    in determinant[1].items()
                },
            "denominator_primes":
                {
                    str(k): int(v)
                    for k, v
                    in determinant[2].items()
                },
        },

        "E8_cubed": {
            "transverse_curvature":
                rational_string(H_trans),
            "negative":
                bool(H_trans < 0),
        },

        "invariance": {
            "realizations":
                len(realizations),
            "exact_spectrum_match":
                True,
            "exact_similarity_verified":
                True,
        },
    }


# ============================================================================
# MAIN CERTIFICATE
# ============================================================================

def main():

    banner(
        "ULTIMATE LEECH LATTICE "
        "INTERTWINER HESSIAN CERTIFICATE"
    )

    print(
        "This program distinguishes EXACT, NUMERICAL, "
        "and MODEL results."
    )

    # ------------------------------------------------------------------------
    # 1. Golay
    # ------------------------------------------------------------------------

    t = section(
        "1. EXACT GOLAY CONSTRUCTION"
    )

    golay = construct_golay_code(
        CFG
    )

    certify_golay(
        golay
    )

    print(
        f"  Codewords       : "
        f"{len(golay['codewords'])}"
    )

    print(
        f"  Octads          : "
        f"{len(golay['octads'])}"
    )

    print(
        f"  Weight spectrum : "
        f"{golay['weight_distribution']}"
    )

    print(
        f"  Completed in {elapsed(t)}"
    )

    # ------------------------------------------------------------------------
    # 2. Leech shell
    # ------------------------------------------------------------------------

    t = section(
        "2. EXACT LEECH SHELL"
    )

    X = construct_leech_shell(
        golay,
        CFG
    )

    certify_leech_shell(
        X,
        CFG
    )

    print(
        f"  |X|             : {len(X)}"
    )

    print(
        f"  Dimension        : {X.shape[1]}"
    )

    print(
        f"  Norm squared     : "
        f"{CFG.shell_norm_squared}"
    )

    print(
        f"  SHA-256          : "
        f"{hash_array(X)}"
    )

    print(
        f"  Completed in {elapsed(t)}"
    )

    # ------------------------------------------------------------------------
    # 3. Reference shell geometry
    # ------------------------------------------------------------------------

    t = section(
        "3. EXACT SHELL DECOMPOSITION"
    )

    u0 = X[0]

    (
        _,
        _,
        _,
        shell_sizes,
        _
    ) = shell_decomposition(
        X,
        u0,
        CFG
    )

    certify_shell_distribution(
        shell_sizes
    )

    print(
        f"  Shell distribution: "
        f"{shell_sizes}"
    )

    print(
        f"  Completed in {elapsed(t)}"
    )

    # ------------------------------------------------------------------------
    # 4. Ward identity
    # ------------------------------------------------------------------------

    t = section(
        "4. EXACT WARD SELF-ENERGY"
    )

    lambda_S = derive_lambda_S(
        shell_sizes,
        CFG
    )

    print(
        f"  lambda_S = {lambda_S}"
    )

    print(
        f"  Completed in {elapsed(t)}"
    )

    # ------------------------------------------------------------------------
    # 5. Primary action matrix
    # ------------------------------------------------------------------------

    t = section(
        "5. EXACT INTERTWINER ACTION"
    )

    primary = construct_action_matrix(
        X,
        u0,
        representative_mode="first",
        tangent_order=list(range(24)),
        measurement_reverse=False,
        cfg=CFG,
    )

    A = primary["A"]

    eigs = exact_spectral_certificate(
        A
    )

    print(
        f"  Matrix dimension : {A.shape}"
    )

    print(
        f"  Rank             : {A.rank()}"
    )

    print(
        f"  Kernel dimension  : "
        f"{12 - A.rank()}"
    )

    print(
        "  Exact spectrum:"
    )

    for i, lam in enumerate(eigs):
        print(
            f"    lambda[{i:2d}] = {lam}"
        )

    print(
        f"  Completed in {elapsed(t)}"
    )

    # ------------------------------------------------------------------------
    # 6. Rotational kernel
    # ------------------------------------------------------------------------

    t = section(
        "6. EXACT ROTATIONAL KERNEL"
    )

    c_rot = rotational_kernel_certificate(
        A
    )

    print(
        f"  c_rot = {list(c_rot)}"
    )

    print(
        "  A*c_rot = 0"
    )

    print(
        "  ker(A) = span_Q{c_rot}"
    )

    print(
        f"  Completed in {elapsed(t)}"
    )

    # ------------------------------------------------------------------------
    # 7. Parity
    # ------------------------------------------------------------------------

    t = section(
        "7. EXACT PARITY INVOLUTION"
    )

    P = construct_parity_operator(
        primary,
        u0
    )

    plus, minus = parity_spectrum(
        P,
        A
    )

    print(
        "  P^2 = I"
    )

    print(
        "  [P,A] = 0"
    )

    print(
        f"  Even modes : {plus}"
    )

    print(
        f"  Odd modes  : {minus}"
    )

    print(
        f"  Tr(P)      : {P.trace()}"
    )

    print(
        f"  Completed in {elapsed(t)}"
    )

    # ------------------------------------------------------------------------
    # 8. Determinant
    # ------------------------------------------------------------------------

    t = section(
        "8. EXACT REGULARIZED DETERMINANT"
    )

    determinant = determinant_certificate(
        eigs
    )

    print(
        f"  det'(A) = {determinant[0]}"
    )

    print(
        "  Numerator prime factors:"
    )

    print(
        f"    {determinant[1]}"
    )

    print(
        "  Denominator prime factors:"
    )

    print(
        f"    {determinant[2]}"
    )

    print(
        f"  Completed in {elapsed(t)}"
    )

    # ------------------------------------------------------------------------
    # 9. Multi-realization invariance
    # ------------------------------------------------------------------------

    t = section(
        "9. MULTI-REALIZATION EXACT EQUIVARIANCE TEST"
    )

    realizations = run_invariance_experiment(
        X,
        CFG
    )

    print(
        f"  Realizations tested : "
        f"{len(realizations)}"
    )

    print(
        "  Exact characteristic-polynomial "
        "agreement verified."
    )

    print(
        "  Exact similarity verified for "
        "every realization against realization 0."
    )

    print(
        f"  Completed in {elapsed(t)}"
    )

    # ------------------------------------------------------------------------
    # 10. E8^3
    # ------------------------------------------------------------------------

    t = section(
        "10. EXACT E8^3 TRANSVERSE CALCULATION"
    )

    H_trans = exact_E8_transverse_certificate()

    print(
        f"  H_transverse = {H_trans}"
    )

    print(
        f"  Negative      : {H_trans < 0}"
    )

    print(
        f"  Completed in {elapsed(t)}"
    )

    # ------------------------------------------------------------------------
    # 11. Numerical Niemeier experiment
    # ------------------------------------------------------------------------

    numerical_results = {}

    if CFG.run_morse_experiment:

        t = section(
            "11. NUMERICAL NİEMEIER EXPERIMENT"
        )

        for k in [1, 2, 3]:

            X_test = make_Ak_m(
                k
            )

            zeros, index, _ = \
                numerical_morse_index(
                    X_test,
                    CFG
                )

            numerical_results[
                f"A_{k}^{24 // k}"
            ] = {
                "rotational_zero_estimate":
                    zeros,
                "morse_index":
                    index,
            }

            print(
                f"  A_{k}^{24 // k}: "
                f"zeros={zeros}, "
                f"Morse index={index}"
            )

        print(
            "  STATUS: NUMERICAL OBSERVATION"
        )

        print(
            f"  Completed in {elapsed(t)}"
        )

    # ------------------------------------------------------------------------
    # 12. Double slit
    # ------------------------------------------------------------------------

    double_slit_results = None

    if CFG.run_double_slit_model:

        t = section(
            "12. PHENOMENOLOGICAL DOUBLE-SLIT MODEL"
        )

        coherent, incoherent = \
            double_slit_model()

        double_slit_results = {
            "coherent_peaks":
                coherent,
            "incoherent_peaks":
                incoherent,
        }

        print(
            f"  Coherent peaks   : "
            f"{coherent}"
        )

        print(
            f"  Incoherent peaks : "
            f"{incoherent}"
        )

        print(
            "  STATUS: MODEL / ILLUSTRATION"
        )

        print(
            f"  Completed in {elapsed(t)}"
        )

    # ------------------------------------------------------------------------
    # 13. Machine-readable certificate
    # ------------------------------------------------------------------------

    certificate = build_certificate(
        golay,
        X,
        lambda_S,
        eigs,
        P,
        determinant,
        H_trans,
        realizations,
    )

    certificate[
        "numerical_observations"
    ] = numerical_results

    certificate[
        "double_slit_model"
    ] = double_slit_results

    certificate[
        "runtime_seconds"
    ] = time.time() - T0

    # ------------------------------------------------------------------------
    # Final report
    # ------------------------------------------------------------------------

    banner(
        "FINAL COMPUTATIONAL CERTIFICATE"
    )

    print(
        """
EXACT RESULTS
-------------
[✓] Binary Golay construction
[✓] 4096 codewords
[✓] 759 octads
[✓] 2576 dodecads
[✓] 196,560 distinct Leech shell vectors
[✓] Every shell vector has norm squared 32
[✓] Seven exact shell inner-product levels
[✓] Exact Ward self-energy
[✓] Exact 12 x 12 rational action matrix
[✓] Rank 11 / kernel dimension 1
[✓] Twelve distinct rational eigenvalues
[✓] Zero eigenvalue
[✓] Eleven strictly positive eigenvalues
[✓] Rotational kernel vector certified
[✓] Exact parity involution
[✓] Six even / six odd modes
[✓] Exact regularized determinant
[✓] Exact E8^3 transverse instability
[✓] Multiple realization similarity tests

NUMERICAL RESULTS
-----------------
[NUMERICAL] Niemeier Morse experiments

MODEL RESULTS
-------------
[MODEL] Double-slit coherence/decoherence demonstration

INTERPRETATION
--------------
The exact algebraic portion establishes a rational finite-dimensional
intertwiner spectrum for the specified Leech-shell construction and
U^(-2) pair potential.

The multi-realization experiment verifies exact similarity for the
tested choices of reference vector, representative, tangent ordering,
and measurement ordering.

This is stronger than merely observing numerical agreement of eigenvalues,
but it is still important to distinguish a finite computational
equivariance test from a separately proved theorem covering every possible
implementation choice.

Likewise, the zero eigenvalue is mathematically certified as a kernel
direction. Calling it a "Goldstone mode" is an additional physical
interpretation.

The first positive eigenvalue is an exact spectral gap of this operator.
Calling it a "Higgs mass" is likewise an interpretive identification,
not a consequence of the algebra alone.
"""
    )

    print(
        "EXACT SPECTRUM"
    )

    for i, lam in enumerate(eigs):

        label = ""

        if i == 0:
            label = "  <-- kernel"
        elif i == 1:
            label = "  <-- spectral gap"

        print(
            f"  {i:2d}: {lam}{label}"
        )

    print()
    print(
        "lambda_S =",
        lambda_S
    )

    print(
        "H_E8^3 =",
        H_trans
    )

    print()
    print(
        "Total runtime:",
        f"{time.time() - T0:.2f}s"
    )

    # ------------------------------------------------------------------------
    # Optional JSON export.
    # ------------------------------------------------------------------------

    with open(
        "leech_certificate.json",
        "w",
        encoding="utf-8"
    ) as f:

        json.dump(
            certificate,
            f,
            indent=2
        )

    print()
    print(
        "Machine-readable certificate written to:"
    )

    print(
        "  leech_certificate.json"
    )

    print()
    print(
        "=" * 80
    )


# ============================================================================
# ENTRY POINT
# ============================================================================

if __name__ == "__main__":
    main()