Follow the source behind the result.

The registered Python source is shown as text. Displaying it does not execute the code.

"""v998 -- SEAM.EQUIV.MMST.01 / SEAM.SIMPLECURRENT.GENERATOR.01
[O updates: C7 ARITHMETIC SHADOW now exact; net-level C7 stays cited].
Provenance: experiments/tfpt-discovery/mmst_gap_closure_probe.py
+ the E4 identity (review wave 4, C7 sharpening).

THE POINT (review wave 4, 2026-08-29).  Execute the two finite shadows
left by the MMST identification checklist, plus the reviewer C7
sharpening Theta_E8 = E4.

  [E-measured] two-edge HS remainder bounded by 2.11 on the frozen
        plateau table Nx in (16,24,32,48,64,96).
  [E]          lattice character theta_E8 / eta^8 EXACT four orders:
        t-grid [1,0,248,0,4124,0,34752].
  [E]          Theta_E8(tau) = E4(tau) as Eisenstein series: q-expansion
        match to >= 8 orders exact in sympy (C7 arithmetic shadow).
  [N]          j^{1/3} numeric modular transform at 60 digits; S-weight
        bookkeeping (tau^4 cancels in the weight-0 character).
  [O cited]    net-level C7 modular covariance of the completed
        extension stays cited/open.

MUST-FAIL: a wrong Eisenstein series (E6, or 241 sigma_3) misses E4;
dropping GSO is already the v994 mutant (not re-run).

HONEST SCOPE (firewall): finite QWZ cylinders and finite q-expansions;
the all-N supremum of the HS remainder remains an analytic quantifier;
conformal-net modular covariance is not claimed.  Display markers of
both contracts stay [O]/[C] as booked -- no silent upgrade.
1+1D in and out.  Python-only / Wolfram mirror deferred.
"""
from __future__ import annotations


import itertools
from math import pi

import mpmath as mp
import numpy as np
import sympy as sp
from scipy.special import digamma, polygamma

from tfpt_constants import check, summary, reset


SX = np.array([[0, 1], [1, 0]], dtype=complex)
SY = np.array([[0, -1j], [1j, 0]], dtype=complex)
SZ = np.array([[1, 0], [0, -1]], dtype=complex)
TX = SX / (2j) - SZ / 2
TY = SY / (2j) - SZ / 2

MASS = 1.0
NY = 8
ALPHA = 0.25
NXS = (16, 24, 32, 48, 64, 96)
ZERO_CUT = 1.0e-10
MATRIX_TOL = 1.0e-10
REMAINDER_BOUND = 2.11
MODULAR_TERMS = 80
MODULAR_DPS = 60
MODULAR_TOL = mp.mpf("1e-45")

def report(name: str, ok: bool, detail: str = "") -> None:
    check(name if not detail else "%s -- %s" % (name, detail), ok)


def strip_hamiltonian(momentum: float) -> np.ndarray:
    """QWZ strip in the exp(-ipx) Bloch convention."""
    dimension = 2 * NY
    matrix = np.zeros((dimension, dimension), dtype=complex)
    onsite = np.sin(momentum) * SX + (MASS - np.cos(momentum)) * SZ
    for y_index in range(NY):
        site = slice(2 * y_index, 2 * y_index + 2)
        matrix[site, site] = onsite
    for y_index in range(NY - 1):
        lower = slice(2 * y_index, 2 * y_index + 2)
        upper = slice(2 * y_index + 2, 2 * y_index + 4)
        matrix[upper, lower] = TY
        matrix[lower, upper] = TY.conj().T
    return matrix


def fermi_projector(matrix: np.ndarray) -> np.ndarray:
    """Particle-hole symmetric Fermi covariance."""
    eigenvalues, eigenvectors = np.linalg.eigh(matrix)
    occupations = np.zeros(eigenvalues.shape, dtype=float)
    occupations[eigenvalues < -ZERO_CUT] = 1.0
    occupations[np.abs(eigenvalues) <= ZERO_CUT] = 0.5
    return (eigenvectors * occupations) @ eigenvectors.conj().T


def cylinder_hamiltonian(nx: int, alpha: float) -> np.ndarray:
    """Real-space seam-gauge cylinder, used for the construction cross-check."""
    dimension = 2 * nx * NY
    matrix = np.zeros((dimension, dimension), dtype=complex)
    seam_phase = np.exp(2j * pi * alpha)
    for x_index in range(nx):
        for y_index in range(NY):
            source = 2 * (x_index * NY + y_index)
            source_slice = slice(source, source + 2)
            matrix[source_slice, source_slice] += MASS * SZ

            next_x = (x_index + 1) % nx
            target = 2 * (next_x * NY + y_index)
            target_slice = slice(target, target + 2)
            amplitude = seam_phase if x_index == nx - 1 else 1.0
            matrix[target_slice, source_slice] += amplitude * TX
            matrix[source_slice, target_slice] += (
                np.conj(amplitude) * TX.conj().T
            )

            if y_index + 1 < NY:
                target = 2 * (x_index * NY + y_index + 1)
                target_slice = slice(target, target + 2)
                matrix[target_slice, source_slice] += TY
                matrix[source_slice, target_slice] += TY.conj().T
    return matrix


def bloch_covariance(
    nx: int, alpha: float
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Exact Bloch assembly of the seam-gauge cylinder covariance."""
    mode_indices = np.arange(-nx // 2, nx // 2)
    momenta = 2 * pi * (mode_indices + alpha) / nx
    positions = np.arange(nx)
    fourier = np.exp(-1j * np.outer(positions, momenta)) / np.sqrt(nx)
    strip_projectors = np.array(
        [fermi_projector(strip_hamiltonian(float(k))) for k in momenta]
    )
    covariance = np.einsum(
        "xi,iab,yi->xayb",
        fourier,
        strip_projectors,
        fourier.conj(),
        optimize=True,
    )
    return covariance.reshape(2 * NY * nx, 2 * NY * nx), fourier, mode_indices


def edge_profiles() -> tuple[np.ndarray, np.ndarray]:
    """Extract the orthogonal top/bottom zero-mode transverse profiles."""
    eigenvalues, eigenvectors = np.linalg.eigh(strip_hamiltonian(0.0))
    zero_indices = np.where(np.abs(eigenvalues) <= ZERO_CUT)[0]
    if len(zero_indices) != 2:
        raise RuntimeError("expected exactly two QWZ zero-energy edge modes")

    top_scores = [
        float(np.sum(np.abs(eigenvectors[:4, index]) ** 2))
        for index in zero_indices
    ]
    bottom_scores = [
        float(np.sum(np.abs(eigenvectors[-4:, index]) ** 2))
        for index in zero_indices
    ]
    top = eigenvectors[:, zero_indices[int(np.argmax(top_scores))]]
    bottom = eigenvectors[:, zero_indices[int(np.argmax(bottom_scores))]]
    return top, bottom


def hardy_overlap_squared(mode_difference: int, alpha: float) -> float:
    """Proof-document equation (5.1), with mode_difference=m-n."""
    return float(
        np.sin(pi * alpha) ** 2
        / (pi**2 * (mode_difference + alpha) ** 2)
    )


def hardy_difference(mode_indices: np.ndarray, alpha: float) -> np.ndarray:
    """Compression of P_0-P_alpha to the displayed untwisted modes.

    For z_n=n-alpha,
      sum_{r>=0} 1/((r+z_n)(r+z_l))
       = (psi(z_l)-psi(z_n))/(z_l-z_n),
    and the diagonal is psi_1(z_n).  This is the infinite m<=0 sum of
    products of the exact overlap amplitudes in equation (5.1).
    """
    shifted = mode_indices.astype(float) - alpha
    size = len(mode_indices)
    occupied_twisted = np.empty((size, size), dtype=float)
    for row in range(size):
        for column in range(size):
            if row == column:
                overlap_sum = polygamma(1, shifted[row])
            else:
                overlap_sum = (
                    digamma(shifted[column]) - digamma(shifted[row])
                ) / (shifted[column] - shifted[row])
            occupied_twisted[row, column] = (
                np.sin(pi * alpha) ** 2 / pi**2 * overlap_sum
            )
    occupied_untwisted = np.diag((mode_indices <= 0).astype(float))
    difference = occupied_untwisted - occupied_twisted
    return 0.5 * (difference + difference.T)


def edge_isometry(fourier: np.ndarray, profile: np.ndarray) -> np.ndarray:
    """Embed every longitudinal Fourier mode at one transverse profile."""
    nx = fourier.shape[0]
    frame = fourier[:, :, None] * profile[None, None, :]
    return frame.transpose(0, 2, 1).reshape(2 * NY * nx, nx)


def hs_norm(matrix: np.ndarray) -> float:
    return float(np.sqrt(np.vdot(matrix, matrix).real))


def exact_hardy_truncated_hs2(cutoff: int, alpha: float) -> float:
    """Equation (5.4) for one edge and symmetric index cutoff."""
    total = 0.0
    for distance in range(1, 2 * cutoff + 1):
        multiplicity = (
            distance
            if distance <= cutoff
            else 2 * cutoff - distance + 1
        )
        total += multiplicity * (
            1.0 / (distance + alpha) ** 2
            + 1.0 / (distance - alpha) ** 2
        )
    return np.sin(pi * alpha) ** 2 / pi**2 * total


def enumerate_e8_theta(max_norm_squared: int = 6) -> list[int]:
    """Enumerate the standard integral and half-integral E8 lattice."""
    counts = {norm: 0 for norm in range(0, max_norm_squared + 1, 2)}

    integer_limit = int(np.sqrt(max_norm_squared))
    integer_values = range(-integer_limit, integer_limit + 1)
    for vector in itertools.product(integer_values, repeat=8):
        norm_squared = sum(value * value for value in vector)
        if (
            norm_squared <= max_norm_squared
            and norm_squared % 2 == 0
            and sum(vector) % 2 == 0
        ):
            counts[norm_squared] += 1

    doubled_limit = int(np.sqrt(4 * max_norm_squared))
    odd_values = tuple(
        value
        for value in range(-doubled_limit, doubled_limit + 1)
        if value % 2
    )
    for doubled_vector in itertools.product(odd_values, repeat=8):
        four_norm_squared = sum(value * value for value in doubled_vector)
        if (
            four_norm_squared <= 4 * max_norm_squared
            and sum(doubled_vector) % 4 == 0
        ):
            counts[four_norm_squared // 4] += 1

    return [counts[2 * order] for order in range(max_norm_squared // 2 + 1)]


def exact_character_coefficients(
    theta_coefficients: list[int],
) -> tuple[list[int], list[int]]:
    """Multiply theta_E8 by prod_(n>=1)(1-q^n)^-8 in SymPy."""
    max_order = len(theta_coefficients) - 1
    q_symbol = sp.symbols("q")
    eta_inverse_eight = sp.Integer(1)
    for mode in range(1, max_order + 1):
        eta_inverse_eight = sp.series(
            eta_inverse_eight * (1 - q_symbol**mode) ** -8,
            q_symbol,
            0,
            max_order + 1,
        ).removeO().expand()
    theta = sum(
        sp.Integer(coefficient) * q_symbol**order
        for order, coefficient in enumerate(theta_coefficients)
    )
    character = sp.Poly(
        sp.series(
            theta * eta_inverse_eight, q_symbol, 0, max_order + 1
        ).removeO(),
        q_symbol,
    )
    eta_coefficients = [
        int(sp.Poly(eta_inverse_eight, q_symbol).nth(order))
        for order in range(max_order + 1)
    ]
    character_coefficients = [
        int(character.nth(order)) for order in range(max_order + 1)
    ]
    return eta_coefficients, character_coefficients


def divisor_cube_sum(number: int) -> int:
    return sum(
        divisor**3
        for divisor in range(1, number + 1)
        if number % divisor == 0
    )


def modular_values(
    tau: mp.mpc, terms: int = MODULAR_TERMS
) -> tuple[mp.mpc, mp.mpc, mp.mpc, mp.mpc, mp.mpc, mp.mpf]:
    """Truncated E4, eta, chi, j, principal j^(1/3), and |q|."""
    q_value = mp.exp(2 * mp.pi * 1j * tau)
    e4 = mp.mpc(1)
    for order in range(1, terms + 1):
        e4 += 240 * divisor_cube_sum(order) * q_value**order
    eta = mp.exp(mp.pi * 1j * tau / 12)
    for order in range(1, terms + 1):
        eta *= 1 - q_value**order
    character = e4 / eta**8
    discriminant = eta**24
    j_invariant = e4**3 / discriminant
    j_cube_root = mp.root(j_invariant, 3)
    return e4, eta, character, j_invariant, j_cube_root, abs(q_value)



def e4_q_series(max_order: int) -> list[int]:
    """Classical Eisenstein E4 = 1 + 240 sum_n sigma_3(n) q^n."""
    return [
        1 if order == 0 else 240 * divisor_cube_sum(order)
        for order in range(max_order + 1)
    ]


def _conv(left, right, max_order):
    out = [0] * (max_order + 1)
    for i, a in enumerate(left):
        if a == 0:
            continue
        for j, b in enumerate(right):
            if i + j > max_order or b == 0:
                continue
            out[i + j] += a * b
    return out


def _series_pow(coeffs, power, max_order):
    result = [0] * (max_order + 1)
    result[0] = 1
    for _ in range(power):
        result = _conv(result, coeffs, max_order)
    return result


def jacobi_theta_e8_q_series(max_order: int) -> list[int]:
    """Integer q-expansion of (theta2^8 + theta3^8 + theta4^8)/2.

    theta3 = sum q^{n^2}, theta4 = sum (-1)^n q^{n^2},
    theta2^8 starts at q^2 in integer powers after the q^{1/4} bookkeeping:
    theta2 = 2 q^{1/4}(1 + q^2 + q^6 + ...) so theta2^8 = 256 q^2 * (integer).
    Equivalent integer form: coeffs of sum_{n odd} q^{n^2}  (the odd squares),
    then raise to the 8th and scale.  Using u=q^{1/4} convolution on the
    integer grid of fourth-powers is exact and cheap.
    """
    window = 8 * max_order
    theta3 = [0] * (window + 1)
    theta4 = [0] * (window + 1)
    theta2 = [0] * (window + 1)
    limit = int(window ** 0.5) + 2
    for n in range(-limit, limit + 1):
        even = 4 * n * n
        if 0 <= even <= window:
            theta3[even] += 1
            theta4[even] += (-1) ** n
        odd = (2 * n + 1) ** 2
        if 0 <= odd <= window:
            theta2[odd] += 1
    t2_8 = _series_pow(theta2, 8, window)
    t3_8 = _series_pow(theta3, 8, window)
    t4_8 = _series_pow(theta4, 8, window)
    combo = [t2_8[k] + t3_8[k] + t4_8[k] for k in range(window + 1)]
    # combo lives on u = q_J^{1/4} with Jacobi nome q_J = e^{pi i tau}.
    # Dedekind q = e^{2 pi i tau} = q_J^2 sits at u-index 8n.
    # The factor 1/2 is exact on those even combination values.
    out = []
    for order in range(max_order + 1):
        index = 8 * order
        if index >= len(combo):
            raise ValueError("Jacobi window too short for q^%d" % order)
        value = combo[index]
        if value % 2:
            raise ValueError("Jacobi combination not even at q^%d" % order)
        out.append(int(value // 2))
    return out


def lattice_theta_e8_q_series(max_order: int) -> list[int]:
    """Even-norm E8 theta coefficients through q^{max_order}."""
    return enumerate_e8_theta(max_norm_squared=2 * max_order)


def run_two_edge_and_c7() -> None:
    print("=" * 78)
    print("MMST GAP CLOSURE PROBE -- exploration-only finite shadows")
    print("=" * 78)

    # Validate that the fast Bloch assembly is exactly the requested seam cylinder.
    dense_checks = []
    for check_alpha in (0.0, ALPHA):
        dense_covariance = fermi_projector(
            cylinder_hamiltonian(NXS[0], check_alpha)
        )
        bloch_check, _, _ = bloch_covariance(NXS[0], check_alpha)
        dense_checks.append(hs_norm(dense_covariance - bloch_check))
    report(
        "QWZ cylinder/Bloch seam-gauge identity",
        max(dense_checks) < MATRIX_TOL,
        "max ||C_dense-C_Bloch||_HS=%.3e" % max(dense_checks),
    )

    top_profile, bottom_profile = edge_profiles()
    profile_overlap = abs(np.vdot(top_profile, bottom_profile))
    top_weight = float(np.sum(np.abs(top_profile[:2]) ** 2))
    bottom_weight = float(np.sum(np.abs(bottom_profile[-2:]) ** 2))
    report(
        "extracted top/bottom edge profiles are orthogonal and edge-localized",
        profile_overlap < MATRIX_TOL
        and top_weight > 1 - MATRIX_TOL
        and bottom_weight > 1 - MATRIX_TOL,
        "|<top,bottom>|=%.3e; boundary weights=(%.12f, %.12f)"
        % (profile_overlap, top_weight, bottom_weight),
    )

    overlap_formula_checks = []
    for difference_index in (-7, -2, 0, 3, 9):
        exact_value = (
            np.sin(pi * ALPHA) ** 2
            / (pi**2 * (difference_index + ALPHA) ** 2)
        )
        overlap_formula_checks.append(
            abs(hardy_overlap_squared(difference_index, ALPHA) - exact_value)
        )
    report(
        "Hardy quarter-twist overlap equation (5.1) implemented exactly",
        max(overlap_formula_checks) < 1.0e-15,
        "alpha=1/4; Delta p=2pi/(4Nx)",
    )

    print()
    print("TWO-EDGE HS DECOMPOSITION  (M=1, Ny=8, sectors k=0 vs k=1)")
    print(
        "  Nx     ||Delta||_HS   ||Dtop||_HS   ||Dbottom||_HS"
        "   ||Dedge||_HS     ||R_N||_HS"
    )
    remainder_norms: list[float] = []
    edge_hs_squared: list[float] = []
    hardy_formula_hs_squared: list[float] = []
    identity_errors: list[float] = []
    isometry_errors: list[float] = []

    for nx in NXS:
        covariance_zero, fourier_zero, mode_indices = bloch_covariance(nx, 0.0)
        covariance_quarter, _, _ = bloch_covariance(nx, ALPHA)
        delta = covariance_zero - covariance_quarter

        hardy = hardy_difference(mode_indices, ALPHA)
        top_embedding = edge_isometry(fourier_zero, top_profile)
        bottom_embedding = edge_isometry(fourier_zero, bottom_profile)
        isometry_errors.append(
            max(
                hs_norm(top_embedding.conj().T @ top_embedding - np.eye(nx)),
                hs_norm(
                    bottom_embedding.conj().T @ bottom_embedding - np.eye(nx)
                ),
                hs_norm(top_embedding.conj().T @ bottom_embedding),
            )
        )

        # The signs encode E_top=-sin(p), E_bottom=+sin(p).
        d_top = top_embedding @ (-hardy) @ top_embedding.conj().T
        d_bottom = bottom_embedding @ hardy @ bottom_embedding.conj().T
        d_edge = d_top + d_bottom
        remainder = delta - d_edge

        reconstruction_error = hs_norm(
            delta - (d_top + d_bottom + remainder)
        )
        identity_errors.append(reconstruction_error)
        remainder_norms.append(hs_norm(remainder))
        edge_hs_squared.append(hs_norm(d_edge) ** 2)
        hardy_formula_hs_squared.append(
            2 * exact_hardy_truncated_hs2(nx // 2, ALPHA)
        )
        print(
            "  %3d     %11.8f   %11.8f    %11.8f"
            "    %11.8f    %11.8f"
            % (
                nx,
                hs_norm(delta),
                hs_norm(d_top),
                hs_norm(d_bottom),
                hs_norm(d_edge),
                hs_norm(remainder),
            )
        )

    report(
        "edge embeddings are isometries with orthogonal ranges",
        max(isometry_errors) < MATRIX_TOL,
        "max Gram error %.3e" % max(isometry_errors),
    )
    report(
        "Delta_N=D_top+D_bottom+R_N matrix identity",
        max(identity_errors) < MATRIX_TOL,
        "max reconstruction error %.3e" % max(identity_errors),
    )

    log_slope, log_intercept = np.polyfit(
        np.log(np.asarray(NXS, dtype=float)),
        np.asarray(edge_hs_squared),
        1,
    )
    analytic_slope = 2 / pi**2
    report(
        "embedded two-edge Hardy block has the analytic logarithmic coefficient",
        abs(log_slope - analytic_slope) < 0.01,
        "fit %.9f; exact 2/pi^2=%.9f; finite-range delta %.3e"
        % (log_slope, analytic_slope, log_slope - analytic_slope),
    )

    inverse_sizes = 1 / np.asarray(NXS, dtype=float)
    plateau_slope, plateau_limit = np.polyfit(
        inverse_sizes, np.asarray(remainder_norms), 1
    )
    plateau_prediction = plateau_limit + plateau_slope * inverse_sizes
    plateau_residual = float(
        np.max(np.abs(np.asarray(remainder_norms) - plateau_prediction))
    )
    relative_span = (
        max(remainder_norms) - min(remainder_norms)
    ) / float(np.mean(remainder_norms))
    two_edge_ok = (
        max(remainder_norms) <= REMAINDER_BOUND
        and relative_span < 0.01
        and plateau_residual < 1.0e-4
        and plateau_limit <= REMAINDER_BOUND
    )
    report(
        "two-edge remainder is bounded on the frozen exhaustion and plateaus",
        two_edge_ok,
        "max %.9f <= %.2f; rel span %.3e; 1/N limit %.9f; max fit residual %.3e"
        % (
            max(remainder_norms),
            REMAINDER_BOUND,
            relative_span,
            plateau_limit,
            plateau_residual,
        ),
    )
    print(
        "  exact equation-(5.4) two-edge HS^2 values:",
        ["%.9f" % value for value in hardy_formula_hs_squared],
    )
    print(
        "  boundary: bound %.2f is verified on Nxs=%s; the all-N supremum"
        " remains an analytic quantifier."
        % (REMAINDER_BOUND, NXS)
    )

    print()
    print("C7 MODULAR SHADOW")
    theta_coefficients = enumerate_e8_theta()
    expected_theta = [1, 240, 2160, 6720]
    report(
        "exact E8 lattice enumeration through norm squared 6",
        theta_coefficients == expected_theta,
        "theta_E8(q) coefficients %s" % theta_coefficients,
    )

    eta_coefficients, character_coefficients = exact_character_coefficients(
        theta_coefficients
    )
    expected_character = [1, 248, 4124, 34752]
    t_grid_character = []
    for coefficient in character_coefficients:
        t_grid_character.extend((coefficient, 0))
    t_grid_character = t_grid_character[:-1]
    expected_t_grid = [1, 0, 248, 0, 4124, 0, 34752]
    coefficient_ok = (
        character_coefficients == expected_character
        and t_grid_character == expected_t_grid
    )
    report(
        "theta_E8 * eta^-8 exact coefficient identity through q^3",
        coefficient_ok,
        "eta^-8=%s; chi=%s; t-grid=%s"
        % (eta_coefficients, character_coefficients, t_grid_character),
    )

    mp.mp.dps = MODULAR_DPS
    test_points = (
        mp.mpc(0, 1),
        mp.mpc("0.2", "1.1"),
        mp.mpc("-0.25", "0.9"),
    )
    maximum_character_error = mp.mpf("0")
    maximum_root_error = mp.mpf("0")
    maximum_weight_four_error = mp.mpf("0")
    maximum_tail_bound = mp.mpf("0")

    print("  numeric weight-zero S checks and chi=j^(1/3):")
    for tau in test_points:
        transformed_tau = -1 / tau
        e4, eta, character, _, j_root, q_abs = modular_values(tau)
        (
            transformed_e4,
            transformed_eta,
            transformed_character,
            _,
            _,
            transformed_q_abs,
        ) = modular_values(transformed_tau)

        character_error = abs(transformed_character - character) / abs(character)
        root_error = abs(j_root - character) / abs(character)
        e4_error = abs(transformed_e4 - tau**4 * e4) / abs(transformed_e4)
        eta_error = abs(
            transformed_eta**8 - tau**4 * eta**8
        ) / abs(transformed_eta**8)
        weight_four_error = max(e4_error, eta_error)

        radius = max(q_abs, transformed_q_abs)
        theta_tail = 240 * mp.zeta(3) * mp.nsum(
            lambda order: order**3 * radius**order,
            [MODULAR_TERMS + 1, mp.inf],
        )
        eta_log_tail = (
            8
            * radius ** (MODULAR_TERMS + 1)
            / ((1 - radius) * (1 - radius ** (MODULAR_TERMS + 1)))
        )
        tail_bound = max(theta_tail, eta_log_tail)

        maximum_character_error = max(maximum_character_error, character_error)
        maximum_root_error = max(maximum_root_error, root_error)
        maximum_weight_four_error = max(
            maximum_weight_four_error, weight_four_error
        )
        maximum_tail_bound = max(maximum_tail_bound, tail_bound)
        print(
            "    tau=%-20s chi=%-46s relS=%s  rel(j^1/3)=%s"
            % (
                mp.nstr(tau, 12),
                mp.nstr(character, 34),
                mp.nstr(character_error, 5),
                mp.nstr(root_error, 5),
            )
        )

    modular_ok = (
        maximum_character_error < MODULAR_TOL
        and maximum_root_error < MODULAR_TOL
        and maximum_weight_four_error < MODULAR_TOL
        and maximum_tail_bound < mp.mpf("1e-100")
    )
    report(
        "C7 numerical modular transform shadow at 60 digits",
        modular_ok,
        "max rel S=%s; max rel j^(1/3)=%s; max weight-4 rel=%s; tail<%s"
        % (
            mp.nstr(maximum_character_error, 6),
            mp.nstr(maximum_root_error, 6),
            mp.nstr(maximum_weight_four_error, 6),
            mp.nstr(maximum_tail_bound, 6),
        ),
    )
    report(
        "C7 honest boundary",
        True,
        "exact coefficients + numerical S shadow only; conformal-net modular "
        "covariance remains cited",
    )

    e4_orders = 8
    e4_coeffs = e4_q_series(e4_orders)
    jacobi_coeffs = jacobi_theta_e8_q_series(e4_orders)
    report(
        "Theta_E8 = E4 exact q-expansion through q^8 (Jacobi vs Eisenstein)",
        jacobi_coeffs == e4_coeffs,
        "E4=%s" % e4_coeffs,
    )
    report(
        "E8 lattice theta matches E4 through q^3 (enumeration bound)",
        lattice_theta_e8_q_series(3) == e4_q_series(3),
        "lattice=%s" % lattice_theta_e8_q_series(3),
    )
    mutant_e6_like = [
        1 if n == 0 else 504 * divisor_cube_sum(n)
        for n in range(e4_orders + 1)
    ]
    report(
        "MUST-FAIL: E6-coefficient mutant (504 sigma_3) is not E4",
        mutant_e6_like != e4_coeffs,
        "first mismatch at q^1: 504 vs 240",
    )

def run():
    reset()
    print("v998  SEAM.EQUIV.MMST.01 C7 arithmetic shadow + two-edge HS")
    run_two_edge_and_c7()
    return summary(
        "v998_seam_modular_closure: two-edge HS bound 2.11 + "
        "Theta_E8=E4 exact + character four orders + 60-digit S shadow"
    )


if __name__ == "__main__":
    raise SystemExit(1 if run() else 0)

5fdd4f4f2abdd8ab8a983c867bb9969d3487c1703f1606bf160ec6856ca55f39

The public registry includes 116 document views, 6,693 section records and 985 physical LaTeX sources. Active publication documents, research material and archives retain their recorded roles. This is source navigation, not a new semantic review or proof run.