#!/usr/bin/env python3
"""Audit Chapter 8's capstone problems on periods and covariance.

The convention is

    hbar**2 psi_zz = R psi,
    y**2 = R_0,
    P**2 + hbar P_z = R,
    Omega = P_even dz.

Every check raises RuntimeError explicitly, so ``python -O`` cannot disable
the test suite.  The script checks exact formal identities, residue data,
operator reductions, and one high-precision elliptic identity.  It does not
test Borel summability, Stokes automorphisms, or a spectral condition.

Compatibility target: Python 3.9, SymPy 1.14, and mpmath.
"""

from __future__ import annotations

import platform

import mpmath as mp
import sympy as sp


def reduced(expression: sp.Expr) -> sp.Expr:
    """Return a stable exact form for a rational symbolic expression."""
    return sp.factor(sp.cancel(sp.together(expression)))


def require_zero(expression: sp.Expr, label: str) -> None:
    """Raise a visible failure unless an exact expression vanishes."""
    residual = reduced(expression)
    if residual != 0:
        residual = sp.simplify(sp.trigsimp(residual))
    if residual != 0:
        residual = sp.factor(sp.simplify(sp.expand_trig(residual)))
    if residual != 0:
        raise RuntimeError(f"{label} failed: residual = {residual}")


def require_matrix_zero(matrix: sp.Matrix, label: str) -> None:
    """Check every entry of a symbolic matrix."""
    for row in range(matrix.rows):
        for column in range(matrix.cols):
            require_zero(
                matrix[row, column],
                f"{label} entry ({row}, {column})",
            )


def require_close(
    actual: mp.mpf | mp.mpc,
    expected: mp.mpf | mp.mpc,
    label: str,
    tolerance: mp.mpf,
) -> None:
    """Raise a visible failure unless two high-precision values agree."""
    scale = max(mp.mpf("1"), abs(actual), abs(expected))
    error = abs(actual - expected)
    if error > tolerance * scale:
        raise RuntimeError(
            f"{label} failed: actual={actual}, expected={expected}, "
            f"scaled error={error / scale}"
        )


def curve_derivative(
    expression: sp.Expr,
    coordinate: sp.Symbol,
    sheet_variable: sp.Symbol,
    radicand: sp.Expr,
) -> sp.Expr:
    """Differentiate along sheet_variable**2 = radicand."""
    return reduced(
        sp.diff(expression, coordinate)
        + sp.diff(radicand, coordinate)
        * sp.diff(expression, sheet_variable)
        / (2 * sheet_variable)
    )


def energy_derivative(
    expression: sp.Expr,
    energy: sp.Symbol,
    sheet_variable: sp.Symbol,
    radicand: sp.Expr,
) -> sp.Expr:
    """Differentiate at fixed base coordinate along y**2 = R_0."""
    return reduced(
        sp.diff(expression, energy)
        + sp.diff(radicand, energy)
        * sp.diff(expression, sheet_variable)
        / (2 * sheet_variable)
    )


def require_on_curve_zero(
    expression: sp.Expr,
    sheet_variable: sp.Symbol,
    radicand: sp.Expr,
    label: str,
) -> None:
    """Check a rational identity modulo sheet_variable**2 = radicand."""
    numerator = sp.together(expression).as_numer_denom()[0]
    polynomial = sp.Poly(numerator, sheet_variable)
    curve = sp.Poly(sheet_variable**2 - radicand, sheet_variable)
    remainder = sp.rem(polynomial, curve).as_expr()
    require_zero(remainder, label)


def residue_at_infinity(
    coefficient: sp.Expr,
    coordinate: sp.Symbol,
) -> sp.Expr:
    """Compute a one-form residue using the local coordinate u=1/z."""
    local = sp.symbols(f"u_{coordinate}")
    pulled_back = -coefficient.subs(coordinate, 1 / local) / local**2
    return reduced(sp.residue(pulled_back, local, 0))


def schwarzian(function: sp.Expr, coordinate: sp.Symbol) -> sp.Expr:
    """Return the Schwarzian derivative {function, coordinate}."""
    first = sp.diff(function, coordinate)
    return reduced(
        sp.diff(function, coordinate, 3) / first
        - sp.Rational(3, 2)
        * (sp.diff(function, coordinate, 2) / first) ** 2
    )


def geometry_and_relative_audit() -> dict[str, sp.Expr]:
    """Check the quartic, rational-path, and Weber calculations."""
    u, k, sigma = sp.symbols("u k sigma", nonzero=True)
    quartic_form = (
        -sigma
        * u**-4
        * sp.sqrt((1 - u**2) * (1 - k**2 * u**2))
    )
    quartic_series = sp.series(quartic_form, u, 0, 2).removeO().expand()
    require_zero(
        quartic_series.coeff(u, -4) + sigma,
        "quartic leading infinity coefficient",
    )
    require_zero(
        quartic_series.coeff(u, -2)
        - sigma * (1 + k**2) / 2,
        "quartic subleading infinity coefficient",
    )
    require_zero(
        sp.residue(quartic_form, u, 0),
        "quartic infinity residue",
    )

    t, c = sp.symbols("t c")
    omega = c / t + 2 / (t + 1) ** 2
    primitive = c * sp.log(t) + (t - 1) / (t + 1)
    require_zero(sp.diff(primitive, t) - omega, "rational primitive")
    residues = (
        reduced(sp.residue(omega, t, 0)),
        reduced(sp.residue(omega, t, -1)),
        residue_at_infinity(omega, t),
    )
    for actual, expected, label in zip(
        residues,
        (c, sp.Integer(0), -c),
        ("zero", "minus one", "infinity"),
    ):
        require_zero(actual - expected, f"rational residue at {label}")
    endpoint_constant = reduced(
        sp.limit((t - 1) / (t + 1), t, sp.oo)
        - sp.limit((t - 1) / (t + 1), t, 0)
    )
    require_zero(endpoint_constant - 2, "rational endpoint constant")

    a = sp.symbols("a", nonzero=True)
    z_cover = a * (t + 1 / t) / 2
    y_cover = a * (t - 1 / t) / 2
    dz_dt = sp.diff(z_cover, t)
    require_zero(
        y_cover**2 - (z_cover**2 - a**2),
        "Weber rational cover",
    )
    lambda_zero = reduced(y_cover * dz_dt)
    expected_zero = a**2 * (t**2 - 1) ** 2 / (4 * t**3)
    require_zero(lambda_zero - expected_zero, "Weber lambda_0")
    residue_zero = reduced(sp.residue(lambda_zero, t, 0))
    residue_infinity = residue_at_infinity(lambda_zero, t)
    require_zero(residue_zero + a**2 / 2, "Weber residue at zero")
    require_zero(
        residue_infinity - a**2 / 2,
        "Weber residue at infinity",
    )

    lambda_two = -t * (3 * t**4 + 14 * t**2 + 3) / (
        2 * a**2 * (t - 1) ** 4 * (t + 1) ** 4
    )
    primitive_two = (9 * t**4 + 12 * t**2 - 1) / (
        12 * a**2 * (t - 1) ** 3 * (t + 1) ** 3
    )
    require_zero(
        sp.diff(primitive_two, t) - lambda_two,
        "Weber exact lambda_2",
    )
    for point, label in ((0, "zero"), (1, "plus one"), (-1, "minus one")):
        require_zero(
            sp.residue(lambda_two, t, point),
            f"Weber lambda_2 residue at {label}",
        )
    require_zero(
        residue_at_infinity(lambda_two, t),
        "Weber lambda_2 residue at infinity",
    )
    weber_endpoint = reduced(
        sp.limit(primitive_two, t, sp.oo)
        - sp.limit(primitive_two, t, 0)
    )
    require_zero(
        weber_endpoint + 1 / (12 * a**2),
        "Weber lambda_2 endpoint constant",
    )

    return {
        "quartic_series": quartic_series,
        "rational_endpoint": endpoint_constant,
        "weber_endpoint": weber_endpoint,
    }


def mathieu_lattice_and_operator_audit() -> dict[str, object]:
    """Check integral cycles, Picard--Lefschetz matrices, and PF operators."""
    intersection = sp.Matrix([[0, 1], [-1, 0]])
    d_minus = sp.Matrix([1, 0])
    d_plus = sp.Matrix([1, 2])
    change = sp.Matrix.hstack(d_minus, d_plus)
    require_zero(change.det() - 2, "Mathieu physical lattice index")
    require_zero(
        (d_minus.T * intersection * d_plus)[0] - 2,
        "Mathieu physical intersection",
    )

    twist_minus = sp.eye(2) + d_minus * d_minus.T * intersection
    twist_plus = sp.eye(2) + d_plus * d_plus.T * intersection
    require_matrix_zero(
        twist_minus - sp.Matrix([[1, 1], [0, 1]]),
        "twist about delta_minus",
    )
    require_matrix_zero(
        twist_plus - sp.Matrix([[-1, 1], [-4, 3]]),
        "twist about delta_plus",
    )
    for twist, vanishing, label in (
        (twist_minus, d_minus, "minus"),
        (twist_plus, d_plus, "plus"),
    ):
        require_matrix_zero(twist * vanishing - vanishing, f"fixed {label} cycle")
        require_matrix_zero(
            twist.T * intersection * twist - intersection,
            f"symplectic {label} twist",
        )
        require_zero(twist.det() - 1, f"determinant of {label} twist")

    energy, coupling, x, y = sp.symbols("E Lambda x y", nonzero=True)
    radicand = 2 * coupling**2 * sp.cos(2 * x) - energy
    delta = energy**2 - 4 * coupling**4
    derivatives = [y]
    for _ in range(3):
        derivatives.append(
            energy_derivative(derivatives[-1], energy, y, radicand)
        )
    g_primitive = coupling**2 * sp.sin(2 * x) / (2 * y)
    pf_form = delta * derivatives[2] + y / 4
    require_on_curve_zero(
        pf_form - curve_derivative(g_primitive, x, y, radicand),
        y,
        radicand,
        "Mathieu Picard--Fuchs certificate",
    )

    test = sp.Function("f")(energy)
    pf_test = delta * sp.diff(test, energy, 2) + test / 4
    d_high = (
        -sp.diff(test, energy) / 4
        - 3 * energy * sp.diff(test, energy, 2)
        - sp.Rational(5, 3) * delta * sp.diff(test, energy, 3)
    )
    d_poly = (
        sp.diff(test, energy) / 6
        + energy * sp.diff(test, energy, 2) / 3
    )
    d_reduced = sp.diff(test, energy) / 6 - energy * test / (12 * delta)
    require_zero(
        d_high - d_poly + sp.Rational(5, 3) * sp.diff(pf_test, energy),
        "Mathieu high versus polynomial operator",
    )
    require_zero(
        d_poly - d_reduced - energy * pf_test / (3 * delta),
        "Mathieu polynomial versus reduced operator",
    )

    d_high_on_y = (
        -derivatives[1] / 4
        - 3 * energy * derivatives[2]
        - sp.Rational(5, 3) * delta * derivatives[3]
    )
    direct_lambda_two = (
        sp.diff(radicand, x, 2) / (8 * y**3)
        - 5 * sp.diff(radicand, x) ** 2 / (32 * y**5)
    )
    require_on_curve_zero(
        direct_lambda_two - d_high_on_y,
        y,
        radicand,
        "Mathieu pointwise high operator",
    )
    f_poly = -5 * coupling**2 * sp.sin(2 * x) / (12 * y**3)
    d_poly_on_y = derivatives[1] / 6 + energy * derivatives[2] / 3
    require_on_curve_zero(
        direct_lambda_two
        - d_poly_on_y
        - curve_derivative(f_poly, x, y, radicand),
        y,
        radicand,
        "Mathieu polynomial operator primitive",
    )
    f_reduced = f_poly + energy * g_primitive / (3 * delta)
    d_reduced_on_y = derivatives[1] / 6 - energy * y / (12 * delta)
    require_on_curve_zero(
        direct_lambda_two
        - d_reduced_on_y
        - curve_derivative(f_reduced, x, y, radicand),
        y,
        radicand,
        "Mathieu reduced operator primitive",
    )

    # Fixed-w versus fixed-x parameter differentiation: alpha=x dx, x=Ew.
    w = sp.symbols("w")
    coordinate_map = energy * w
    form_coefficient = coordinate_map
    pulled_alpha = form_coefficient * sp.diff(coordinate_map, w)
    endpoint_derivative = sp.integrate(
        sp.diff(pulled_alpha, energy), (w, 0, 1)
    )
    boundary_term = (
        form_coefficient * sp.diff(coordinate_map, energy)
    ).subs(w, 1) - (
        form_coefficient * sp.diff(coordinate_map, energy)
    ).subs(w, 0)
    require_zero(
        endpoint_derivative - boundary_term,
        "moving-coordinate endpoint identity",
    )
    require_zero(boundary_term - energy, "moving-coordinate endpoint value")

    mp.mp.dps = 80
    m_value = mp.mpf("0.371")
    coupling_value = mp.mpf("1.7")
    elliptic_k = mp.ellipk
    elliptic_e = mp.ellipe
    action_minus = 8 * coupling_value * (
        elliptic_e(m_value)
        - (1 - m_value) * elliptic_k(m_value)
    )
    action_plus = 8 * coupling_value * (
        elliptic_e(1 - m_value)
        - m_value * elliptic_k(1 - m_value)
    )
    period_minus = mp.j * action_minus
    period_plus = action_plus
    derivative_minus = mp.j * elliptic_k(m_value) / coupling_value
    derivative_plus = -elliptic_k(1 - m_value) / coupling_value
    physical_wronskian = (
        period_minus * derivative_plus
        - derivative_minus * period_plus
    )
    primitive_wronskian = physical_wronskian / 2
    tolerance = mp.mpf("1e-65")
    require_close(
        physical_wronskian,
        -4 * mp.pi * mp.j,
        "Mathieu physical-period Wronskian",
        tolerance,
    )
    require_close(
        primitive_wronskian,
        -2 * mp.pi * mp.j,
        "Mathieu primitive-period Wronskian",
        tolerance,
    )

    return {
        "twist_minus": twist_minus,
        "twist_plus": twist_plus,
        "physical_wronskian": physical_wronskian,
    }


def projective_covariance_audit() -> dict[str, sp.Expr]:
    """Check half-density, Riccati, Schwarzian, and p_2 covariance."""
    w, hbar = sp.symbols("w hbar", nonzero=True)
    s = sp.Function("s")(w)
    phi = sp.Function("phi")(w)
    potential = sp.Function("R")(w)
    scalar_field = sp.sqrt(s) * phi
    scalar_equation = (
        hbar**2 * sp.diff(scalar_field, w, 2)
        - hbar**2 * sp.diff(s, w) * sp.diff(scalar_field, w) / s
        - s**2 * potential * scalar_field
    ) / sp.sqrt(s)
    sch = (
        sp.diff(s, w, 2) / s
        - sp.Rational(3, 2) * (sp.diff(s, w) / s) ** 2
    )
    expected_normal_form = (
        hbar**2 * sp.diff(phi, w, 2)
        - (s**2 * potential - hbar**2 * sch / 2) * phi
    )
    require_zero(
        scalar_equation - expected_normal_form,
        "inverse-half-density normal form",
    )

    p = sp.Function("p")(w)
    transformed_p = s * p - hbar * sp.diff(s, w) / (2 * s)
    original_r_in_w = p**2 + hbar * sp.diff(p, w) / s
    transformed_r = s**2 * original_r_in_w - hbar**2 * sch / 2
    require_zero(
        transformed_p**2
        + hbar * sp.diff(transformed_p, w)
        - transformed_r,
        "Riccati projective transformation",
    )

    f = sp.Function("f")(w)
    g = sp.Function("g")(w)
    transformed_f = f / sp.sqrt(s)
    transformed_g = g / sp.sqrt(s)
    wronskian_w = (
        transformed_f * sp.diff(transformed_g, w)
        - sp.diff(transformed_f, w) * transformed_g
    )
    wronskian_z = (f * sp.diff(g, w) - sp.diff(f, w) * g) / s
    require_zero(wronskian_w - wronskian_z, "Wronskian covariance")

    u = sp.symbols("u")
    outer = sp.Function("outer")
    inner = sp.Function("inner")(u)
    composite = outer(inner)
    chain_residual = (
        schwarzian(composite, u)
        - sp.diff(inner, u) ** 2
        * schwarzian(outer(w), w).subs(w, inner)
        - schwarzian(inner, u)
    )
    require_zero(chain_residual, "Schwarzian chain rule")

    y = sp.Function("y")(w)
    r_two = sp.Function("r2")(w)
    original_y_z = sp.diff(y, w) / s
    original_y_zz = sp.diff(original_y_z, w) / s
    p_two_original = (
        r_two / (2 * y)
        + original_y_zz / (4 * y**2)
        - 3 * original_y_z**2 / (8 * y**3)
    )
    transformed_y = s * y
    transformed_r_two = s**2 * r_two - sch / 2
    p_two_transformed = (
        transformed_r_two / (2 * transformed_y)
        + sp.diff(transformed_y, w, 2) / (4 * transformed_y**2)
        - 3 * sp.diff(transformed_y, w) ** 2 / (8 * transformed_y**3)
    )
    require_zero(
        p_two_transformed - s * p_two_original,
        "general p_2 covariance",
    )
    p_two_naive = (
        s**2 * r_two / (2 * transformed_y)
        + sp.diff(transformed_y, w, 2) / (4 * transformed_y**2)
        - 3 * sp.diff(transformed_y, w) ** 2 / (8 * transformed_y**3)
    )
    require_zero(
        p_two_naive - s * p_two_original - sch / (4 * s * y),
        "omitted-Schwarzian p_2 error",
    )

    return {
        "schwarzian": sch,
        "transformed_momentum": transformed_p,
        "p_two": p_two_transformed,
    }


def quadratic_map_and_residue_audit() -> dict[str, sp.Expr]:
    """Check z=w**2 for constant and inverse-square equations."""
    w, z_symbol, hbar, kappa = sp.symbols(
        "w z hbar kappa", nonzero=True
    )
    z = w**2
    s = sp.diff(z, w)
    sch = schwarzian(z, w)
    require_zero(sch + 3 / (2 * w**2), "quadratic-map Schwarzian")

    transformed_r = 4 * kappa**2 * w**2 + 3 * hbar**2 / (4 * w**2)
    transformed_momenta: list[sp.Expr] = []
    for sign, label in ((1, "plus"), (-1, "minus")):
        transformed_solution = (2 * w) ** sp.Rational(-1, 2) * sp.exp(
            sign * kappa * w**2 / hbar
        )
        residual = (
            hbar**2 * sp.diff(transformed_solution, w, 2)
            - transformed_r * transformed_solution
        )
        require_zero(residual, f"quadratic-map exact free solution {label}")
        transformed_p = (
            hbar * sp.diff(transformed_solution, w) / transformed_solution
        )
        require_zero(
            transformed_p
            - (2 * sign * kappa * w - hbar / (2 * w)),
            f"quadratic-map exact Riccati momentum {label}",
        )
        transformed_momenta.append(reduced(transformed_p))
    geometric_p_two = -3 / (16 * kappa * w**3)
    source_p_two = 3 / (16 * kappa * w**3)
    require_zero(
        geometric_p_two + source_p_two,
        "quadratic-map p_2 cancellation",
    )

    nu = sp.sqrt(kappa**2 + hbar**2 / 4)
    original_even = nu / z
    pulled_even = s * original_even
    transformed_inverse_square_r = (
        4 * kappa**2 + 3 * hbar**2 / 4
    ) / w**2
    candidate_plus = (hbar / 2 + 2 * nu) / w
    candidate_minus = (hbar / 2 - 2 * nu) / w
    for candidate, label in (
        (candidate_plus, "plus"),
        (candidate_minus, "minus"),
    ):
        require_zero(
            candidate**2
            + hbar * sp.diff(candidate, w)
            - transformed_inverse_square_r,
            f"inverse-square transformed Riccati branch {label}",
        )
    require_zero(
        (candidate_plus - candidate_minus) / 2 - pulled_even,
        "inverse-square phase pullback",
    )
    nu_series = sp.series(nu, hbar, 0, 8).removeO()
    expected_series = (
        kappa
        + hbar**2 / (8 * kappa)
        - hbar**4 / (128 * kappa**3)
        + hbar**6 / (1024 * kappa**5)
    )
    # SymPy's principal square root needs the formal branch sqrt(kappa**2)=kappa.
    nu_series = nu_series.xreplace({sp.sqrt(kappa**2): kappa})
    require_zero(nu_series - expected_series, "inverse-square quantum residue")
    original_form_coefficient = nu / z_symbol
    original_residue = sp.residue(original_form_coefficient, z_symbol, 0)
    pulled_form_coefficient = reduced(original_even * s)
    pulled_residue = sp.residue(pulled_form_coefficient, w, 0)
    require_zero(
        pulled_residue - 2 * original_residue,
        "quadratic-map pullback residue",
    )
    require_zero(
        2 * sp.pi * sp.I * pulled_residue
        - 2 * (2 * sp.pi * sp.I * original_residue),
        "quadratic-map period winding factor",
    )

    naive_even = sp.sqrt(4 * kappa**2 + hbar**2 / 4) / w
    naive_series = sp.series(naive_even * w, hbar, 0, 4).removeO()
    naive_series = naive_series.xreplace({sp.sqrt(kappa**2): kappa})
    require_zero(
        sp.expand(naive_series - (2 * kappa + hbar**2 / (16 * kappa))),
        "omitted-Schwarzian inverse-square series",
    )

    return {
        "schwarzian": sch,
        "free_momentum": transformed_momenta[0],
        "quantum_residue": nu_series,
    }


def endpoint_and_mathieu_chart_audit() -> dict[str, sp.Expr]:
    """Check finite-part chart dependence and the Mathieu chart anomaly."""
    u = sp.symbols("u", positive=True)
    a_two, residue, c_one, d_two, mu_t, mu_u = sp.symbols(
        "a2 r c d mu_t mu_u", nonzero=True, positive=True
    )
    t_of_u = c_one * u + d_two * u**2
    t_primitive = -a_two / t_of_u + residue * sp.log(t_of_u / mu_t)
    u_singular = (
        -a_two / (c_one * u)
        + residue * sp.log(u / mu_u)
    )
    constant_shift = sp.limit(t_primitive - u_singular, u, 0, dir="+")
    expected_shift = (
        a_two * d_two / c_one**2
        + residue * sp.log(c_one * mu_u / mu_t)
    )
    require_zero(
        sp.expand_log(constant_shift - expected_shift, force=True),
        "finite-part coordinate shift",
    )

    epsilon, endpoint, nonlinear = sp.symbols(
        "epsilon b nonlinear", positive=True
    )
    integral_w = -1 / endpoint + 1 / (epsilon + nonlinear * epsilon**2)
    naive_finite_part = sp.limit(integral_w - 1 / epsilon, epsilon, 0)
    transported_finite_part = sp.limit(
        integral_w - (1 / epsilon - nonlinear),
        epsilon,
        0,
    )
    require_zero(
        naive_finite_part + 1 / endpoint + nonlinear,
        "naive nonlinear finite part",
    )
    require_zero(
        transported_finite_part + 1 / endpoint,
        "transported nonlinear finite part",
    )

    w, energy, coupling, sheet = sp.symbols(
        "w E Lambda Y", nonzero=True
    )
    x_of_w = sp.log(w) / (2 * sp.I)
    sch = schwarzian(x_of_w, w)
    require_zero(sch - 1 / (2 * w**2), "Mathieu logarithmic Schwarzian")
    polynomial = coupling**2 * w**3 - energy * w**2 + coupling**2 * w
    source_r_zero = coupling**2 * (w + 1 / w) - energy
    coordinate_jacobian = sp.diff(x_of_w, w)
    transformed_r_zero = coordinate_jacobian**2 * source_r_zero
    expected_r_zero = -source_r_zero / (4 * w**2)
    require_zero(
        transformed_r_zero - expected_r_zero,
        "Mathieu transformed leading coefficient",
    )
    lambda_zero = sheet / (2 * sp.I * w**2)
    require_on_curve_zero(
        lambda_zero**2 - transformed_r_zero,
        sheet,
        polynomial,
        "Mathieu transformed classical one-form coefficient",
    )
    lambda_zero_energy = energy_derivative(
        lambda_zero,
        energy,
        sheet,
        polynomial,
    )
    correct_r_two = -1 / (4 * w**2)
    leading_y = sheet / (2 * sp.I * w**2)
    omitted_schwarzian_error = -correct_r_two / (2 * leading_y)
    require_on_curve_zero(
        omitted_schwarzian_error - lambda_zero_energy,
        sheet,
        polynomial,
        "Mathieu omitted-Schwarzian period density",
    )
    require_on_curve_zero(
        lambda_zero_energy - sp.I / (4 * sheet),
        sheet,
        polynomial,
        "Mathieu classical energy derivative",
    )

    return {
        "finite_part_shift": expected_shift,
        "mathieu_schwarzian": sch,
        "mathieu_error": omitted_schwarzian_error,
    }


def main() -> None:
    """Run every audit and print a compact reproducibility ledger."""
    print("Advanced Linear ODEs -- cycles, residues, covariance audit")
    print(f"Python: {platform.python_version()}")
    print(f"SymPy: {sp.__version__}")
    print(f"mpmath: {mp.__version__}")
    print()

    geometry = geometry_and_relative_audit()
    print("[PASS] Quartic punctures, relative paths, and Weber endpoint data")
    print(f"  quartic infinity series = {geometry['quartic_series']}")
    print(f"  Weber lambda_2 endpoint = {geometry['weber_endpoint']}")

    mathieu = mathieu_lattice_and_operator_audit()
    print("[PASS] Mathieu lattice, Picard--Lefschetz, and PF operators")
    print(f"  T_minus = {mathieu['twist_minus'].tolist()}")
    print(f"  T_plus = {mathieu['twist_plus'].tolist()}")
    print(
        "  W(delta_minus, delta_plus) = "
        f"{mp.nstr(mathieu['physical_wronskian'], 18)}"
    )

    covariance = projective_covariance_audit()
    print("[PASS] Half-density, Schwarzian cocycle, Wronskian, and p_2")
    print(f"  transformed P = {covariance['transformed_momentum']}")

    quadratic = quadratic_map_and_residue_audit()
    print("[PASS] Quadratic-map solutions and inverse-square quantum residue")
    print(f"  Schwarzian = {quadratic['schwarzian']}")
    print(f"  nu(hbar) = {quadratic['quantum_residue']} + O(hbar**8)")

    endpoints = endpoint_and_mathieu_chart_audit()
    print("[PASS] Endpoint finite parts and Mathieu chart anomaly")
    print(f"  finite-part shift = {endpoints['finite_part_shift']}")
    print(f"  Mathieu lambda_2 error = {endpoints['mathieu_error']} dw")

    print()
    print("All cycles--residues--covariance checks passed.")
    print("Scope: formal identities and period data; no exact-WKB claim tested.")


if __name__ == "__main__":
    main()