#!/usr/bin/env python3
"""Reproduce the formal-WKB formulas for Airy, Weber, and Mathieu.

The convention throughout is

    [-hbar**2 d_x**2 + V(x)] psi = E psi,
    y**2 = R(x, E) = V(x) - E,
    P**2 + hbar P' = R,
    Omega = P_even dx = lambda_0 + hbar**2 lambda_2 + ... .

The audit has four layers:

1. exact local Airy and Weber identities on rational spectral covers;
2. Weber residues, endpoint primitives, and Bernoulli/gamma coefficients;
3. exact Mathieu Picard--Fuchs and quantum-period operator identities; and
4. high-precision Mathieu quadratures and cycle-normalization checks.

Every test raises RuntimeError explicitly.  Python's -O flag therefore cannot
disable any check.  The script verifies formal coefficients and classical
period integrals; it does not claim Borel summability or an exact
quantization 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_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]
    try:
        polynomial = sp.Poly(numerator, sheet_variable)
        curve = sp.Poly(sheet_variable**2 - radicand, sheet_variable)
        residual = sp.rem(polynomial, curve).as_expr()
    except sp.PolynomialError as error:
        raise RuntimeError(
            f"{label} could not be reduced on the spectral curve: {error}"
        ) from error
    require_zero(residual, label)


def require_close(
    actual: mp.mpf,
    expected: mp.mpf,
    label: str,
    tolerance: mp.mpf,
) -> None:
    """Raise a visible failure unless two high-precision numbers 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 y**2 = radicand with respect to coordinate."""
    derivative = sp.diff(expression, coordinate)
    derivative += (
        sp.diff(radicand, coordinate)
        * sp.diff(expression, sheet_variable)
        / (2 * sheet_variable)
    )
    return reduced(derivative)


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


def lambda_two_coefficient(
    radicand: sp.Expr,
    coordinate: sp.Symbol,
    sheet_variable: sp.Symbol,
) -> sp.Expr:
    """Return the coefficient of dx in lambda_2 for y**2 = radicand."""
    first = sp.diff(radicand, coordinate)
    second = sp.diff(radicand, coordinate, 2)
    return reduced(
        second / (8 * sheet_variable**3)
        - 5 * first**2 / (32 * sheet_variable**5)
    )


def riccati_lambda_two_coefficient(
    radicand: sp.Expr,
    coordinate: sp.Symbol,
    sheet_variable: sp.Symbol,
) -> sp.Expr:
    """Derive p_2 directly from p_n recursion on y**2 = radicand."""
    p_zero = sheet_variable
    p_one = reduced(
        -curve_derivative(
            p_zero, coordinate, sheet_variable, radicand
        )
        / (2 * p_zero)
    )
    p_two = reduced(
        -(
            p_one**2
            + curve_derivative(
                p_one, coordinate, sheet_variable, radicand
            )
        )
        / (2 * p_zero)
    )
    return p_two


def generate_branch_on_cover(
    leading_momentum: sp.Expr,
    base_coordinate: sp.Expr,
    cover_coordinate: sp.Symbol,
    order: int,
) -> list[sp.Expr]:
    """Generate one Riccati branch after rationalizing the spectral cover."""
    if order < 0:
        raise RuntimeError(f"formal order must be nonnegative, got {order}")

    dx_dt = sp.diff(base_coordinate, cover_coordinate)
    if dx_dt == 0:
        raise RuntimeError("the cover parametrization has zero derivative")

    def differentiate_in_x(expression: sp.Expr) -> sp.Expr:
        return reduced(sp.diff(expression, cover_coordinate) / dx_dt)

    coefficients = [reduced(leading_momentum)]
    for n in range(1, order + 1):
        convolution = sum(
            coefficients[j] * coefficients[n - j]
            for j in range(1, n)
        )
        next_coefficient = -(
            differentiate_in_x(coefficients[n - 1]) + convolution
        ) / (2 * coefficients[0])
        coefficients.append(reduced(next_coefficient))

    for n in range(1, order + 1):
        convolution = sum(
            coefficients[j] * coefficients[n - j]
            for j in range(1, n)
        )
        residual = (
            2 * coefficients[0] * coefficients[n]
            + convolution
            + differentiate_in_x(coefficients[n - 1])
        )
        require_zero(residual, f"Riccati recursion coefficient hbar^{n}")

    return coefficients


def local_residue_at_infinity(
    one_form_coefficient: sp.Expr,
    coordinate: sp.Symbol,
) -> sp.Expr:
    """Return the residue using the positive local coordinate u=1/t."""
    local_coordinate = sp.symbols(f"u_{coordinate}")
    pulled_back = -one_form_coefficient.subs(
        coordinate, 1 / local_coordinate
    ) / local_coordinate**2
    return reduced(sp.residue(pulled_back, local_coordinate, 0))


def airy_audit() -> dict[str, sp.Expr]:
    """Check the Airy cover, differential operator, and first asymptotic term."""
    x, energy, y = sp.symbols("x E y")
    xi, hbar = sp.symbols("xi hbar", positive=True)
    scaled_coordinate = sp.symbols("X")
    radicand = x - energy
    lambda_two = lambda_two_coefficient(radicand, x, y)

    # In hbar**2 psi_xx = x psi, x=hbar**(2/3) X leaves psi_XX=X psi.
    airy_scaled_coefficient = reduced(
        hbar ** sp.Rational(4, 3)
        * (hbar ** sp.Rational(2, 3) * scaled_coordinate)
        / hbar**2
    )
    require_zero(
        airy_scaled_coefficient - scaled_coordinate,
        "Airy hbar^(2/3) scaling",
    )

    energy_derivatives = [y]
    for _ in range(3):
        energy_derivatives.append(
            energy_derivative(
                energy_derivatives[-1], energy, y, radicand
            )
        )
    require_on_curve_zero(
        lambda_two - sp.Rational(5, 12) * energy_derivatives[3],
        y,
        radicand,
        "Airy pointwise third-energy-derivative operator",
    )

    x_on_cover = energy + xi**2
    y_on_cover = xi
    dx_dxi = sp.diff(x_on_cover, xi)
    cover_coefficients = generate_branch_on_cover(
        y_on_cover, x_on_cover, xi, 2
    )
    pulled_lambda_zero = reduced(y_on_cover * dx_dxi)
    pulled_lambda_two = reduced(
        lambda_two.subs({x: x_on_cover, y: y_on_cover}) * dx_dxi
    )
    require_zero(
        pulled_lambda_two - cover_coefficients[2] * dx_dxi,
        "Airy lambda_2 versus direct Riccati recursion",
    )
    require_zero(
        pulled_lambda_zero - 2 * xi**2,
        "Airy classical differential on the rational cover",
    )
    require_zero(
        pulled_lambda_two + sp.Rational(5, 16) / xi**4,
        "Airy lambda_2 on the rational cover",
    )
    require_zero(
        sp.residue(pulled_lambda_two, xi, 0),
        "Airy lambda_2 turning-point residue",
    )

    action_zero = sp.Rational(2, 3) * y**3
    action_two = sp.Rational(5, 48) / y**3
    require_on_curve_zero(
        curve_derivative(action_zero, x, y, radicand) - y,
        y,
        radicand,
        "Airy classical primitive",
    )
    require_on_curve_zero(
        curve_derivative(action_two, x, y, radicand) - lambda_two,
        y,
        radicand,
        "Airy lambda_2 primitive",
    )

    action_zero_derivative = action_zero
    for _ in range(3):
        action_zero_derivative = energy_derivative(
            action_zero_derivative, energy, y, radicand
        )
    require_on_curve_zero(
        action_two - sp.Rational(5, 12) * action_zero_derivative,
        y,
        radicand,
        "Airy relative-action energy operator",
    )

    action_two_on_cover = action_two.subs(y, xi)
    decaying_first_term = sp.diff(
        sp.exp(-hbar * action_two_on_cover), hbar
    ).subs(hbar, 0)
    require_zero(
        decaying_first_term + sp.Rational(5, 48) / xi**3,
        "first Airy Ai asymptotic coefficient",
    )

    return {
        "lambda_two": lambda_two,
        "cover_lambda_two": pulled_lambda_two,
        "action_two": action_two,
        "ai_first_term": decaying_first_term,
    }


def weber_local_operator_audit() -> dict[str, sp.Expr]:
    """Check the Weber Picard--Fuchs certificate and lambda_2 operators."""
    x, energy, y = sp.symbols("x E y")
    a, hbar = sp.symbols("a hbar", positive=True)
    scaled_coordinate = sp.symbols("X")
    radicand = x**2 - energy
    lambda_two = lambda_two_coefficient(radicand, x, y)
    require_on_curve_zero(
        lambda_two
        - riccati_lambda_two_coefficient(radicand, x, y),
        y,
        radicand,
        "Weber lambda_2 versus direct Riccati recursion",
    )
    expected_lambda_two = -(3 * x**2 + 2 * energy) / (8 * y**5)
    require_on_curve_zero(
        lambda_two - expected_lambda_two,
        y,
        radicand,
        "Weber lambda_2 coefficient",
    )

    energy_derivatives = [y]
    for _ in range(3):
        energy_derivatives.append(
            energy_derivative(
                energy_derivatives[-1], energy, y, radicand
            )
        )

    pf_primitive = x / y
    require_on_curve_zero(
        4 * energy * energy_derivatives[2]
        - curve_derivative(pf_primitive, x, y, radicand),
        y,
        radicand,
        "Weber Picard--Fuchs certificate",
    )

    direct_operator = (
        sp.Rational(3, 2) * energy_derivatives[2]
        + sp.Rational(5, 3) * energy * energy_derivatives[3]
    )
    require_on_curve_zero(
        lambda_two - direct_operator,
        y,
        radicand,
        "Weber pointwise lambda_2 operator",
    )

    primitive_two = (
        5 * x / (24 * y**3) - x / (24 * energy * y)
    )
    require_on_curve_zero(
        curve_derivative(primitive_two, x, y, radicand) - lambda_two,
        y,
        radicand,
        "Weber exact lambda_2 primitive",
    )

    test_function = sp.Function("f")(energy)
    pf_on_test = energy * sp.diff(test_function, energy, 2)
    direct_on_test = (
        sp.Rational(3, 2) * sp.diff(test_function, energy, 2)
        + sp.Rational(5, 3)
        * energy
        * sp.diff(test_function, energy, 3)
    )
    ideal_on_test = (
        sp.Rational(5, 3) * sp.diff(pf_on_test, energy)
        - pf_on_test / (6 * energy)
    )
    require_zero(
        direct_on_test - ideal_on_test,
        "Weber direct operator in the Picard--Fuchs left ideal",
    )

    # X=sqrt(2/hbar) x gives the standard parabolic-cylinder equation.
    parabolic_parameter = a**2 / (2 * hbar) - sp.Rational(1, 2)
    transformed_coefficient = (
        a**2 / (2 * hbar) - scaled_coordinate**2 / 4
    )
    require_zero(
        transformed_coefficient
        - (
            parabolic_parameter
            + sp.Rational(1, 2)
            - scaled_coordinate**2 / 4
        ),
        "Weber parabolic-cylinder parameter",
    )

    return {
        "lambda_two": lambda_two,
        "pf_primitive": pf_primitive,
        "lambda_two_primitive": primitive_two,
        "parabolic_parameter": parabolic_parameter,
    }


def weber_period_audit() -> dict[str, object]:
    """Check Weber residues and relative periods through hbar**6."""
    t = sp.symbols("t")
    a = sp.symbols("a", positive=True)
    energy = sp.symbols("E", positive=True)
    x_on_cover = a * (t + 1 / t) / 2
    y_on_cover = a * (t - 1 / t) / 2
    dx_dt = sp.diff(x_on_cover, t)

    require_zero(
        y_on_cover**2 - (x_on_cover**2 - a**2),
        "Weber rational parametrization",
    )
    coefficients = generate_branch_on_cover(
        y_on_cover, x_on_cover, t, 6
    )
    one_forms = [reduced(coefficient * dx_dt) for coefficient in coefficients]

    expected_lambda_two = -(
        t * (3 * t**4 + 14 * t**2 + 3)
    ) / (2 * a**2 * (t - 1) ** 4 * (t + 1) ** 4)
    require_zero(
        one_forms[2] - expected_lambda_two,
        "Weber lambda_2 in the rational parameter",
    )
    expected_primitive_two = (
        9 * t**4 + 12 * t**2 - 1
    ) / (12 * a**2 * (t - 1) ** 3 * (t + 1) ** 3)
    require_zero(
        sp.diff(expected_primitive_two, t) - one_forms[2],
        "Weber lambda_2 rational primitive",
    )

    residue_plus = local_residue_at_infinity(one_forms[0], t)
    residue_minus = reduced(sp.residue(one_forms[0], t, 0))
    require_zero(
        residue_plus - a**2 / 2,
        "Weber classical residue at infinity_plus",
    )
    require_zero(
        residue_minus + a**2 / 2,
        "Weber classical residue at infinity_minus",
    )
    require_zero(
        residue_plus + residue_minus,
        "sum of Weber classical infinity residues",
    )
    closed_period = reduced(2 * sp.pi * sp.I * residue_plus)
    require_zero(
        closed_period - sp.I * sp.pi * a**2,
        "Weber oriented closed period",
    )

    relative_coefficients: list[sp.Expr] = []
    gamma_coefficients: list[sp.Expr] = []
    for k in range(1, 4):
        order = 2 * k
        one_form = one_forms[order]
        for point, point_label in ((1, "+1"), (-1, "-1"), (0, "0")):
            require_zero(
                sp.residue(one_form, t, point),
                f"Weber lambda_{order} residue at t={point_label}",
            )
        require_zero(
            local_residue_at_infinity(one_form, t),
            f"Weber lambda_{order} residue at t=infinity",
        )

        primitive = reduced(sp.integrate(one_form, t))
        require_zero(
            sp.diff(primitive, t) - one_form,
            f"Weber lambda_{order} rational primitive",
        )
        endpoint_value = reduced(
            sp.limit(primitive, t, sp.oo) - sp.limit(primitive, t, 0)
        )
        bernoulli_value = reduced(
            (1 - 2 ** (2 * k - 1))
            * sp.bernoulli(2 * k)
            / (2 * k * (2 * k - 1) * a ** (4 * k - 2))
        )
        require_zero(
            endpoint_value - bernoulli_value,
            f"Weber Bernoulli relative coefficient hbar^{order}",
        )
        relative_coefficients.append(endpoint_value)

        gamma_value = reduced(
            2 ** (2 * k - 1)
            * sp.bernoulli(2 * k, sp.Rational(1, 2))
            / (
                2
                * k
                * (2 * k - 1)
                * energy ** (2 * k - 1)
            )
        )
        all_order_value = reduced(
            (1 - 2 ** (2 * k - 1))
            * sp.bernoulli(2 * k)
            / (
                2
                * k
                * (2 * k - 1)
                * energy ** (2 * k - 1)
            )
        )
        require_zero(
            gamma_value - all_order_value,
            f"Weber gamma-tail coefficient hbar^{order}",
        )
        require_zero(
            endpoint_value - gamma_value.subs(energy, a**2),
            f"Weber endpoint versus gamma coefficient hbar^{order}",
        )
        gamma_coefficients.append(gamma_value)

    for k in range(4, 7):
        gamma_value = reduced(
            2 ** (2 * k - 1)
            * sp.bernoulli(2 * k, sp.Rational(1, 2))
            / (
                2
                * k
                * (2 * k - 1)
                * energy ** (2 * k - 1)
            )
        )
        all_order_value = reduced(
            (1 - 2 ** (2 * k - 1))
            * sp.bernoulli(2 * k)
            / (
                2
                * k
                * (2 * k - 1)
                * energy ** (2 * k - 1)
            )
        )
        require_zero(
            gamma_value - all_order_value,
            f"Weber all-order gamma/Bernoulli identity k={k}",
        )

    expected_values = [
        -1 / (12 * a**2),
        sp.Rational(7, 360) / a**6,
        -sp.Rational(31, 1260) / a**10,
    ]
    for order, actual, expected in zip(
        (2, 4, 6), relative_coefficients, expected_values
    ):
        require_zero(
            actual - expected,
            f"explicit Weber endpoint coefficient hbar^{order}",
        )

    return {
        "closed_period": closed_period,
        "relative_coefficients": relative_coefficients,
        "gamma_coefficients": gamma_coefficients,
    }


def mathieu_curve_and_operator_audit() -> dict[str, sp.Expr]:
    """Check the Mathieu elliptic curve, PF certificate, and lambda_2."""
    x, energy, y = sp.symbols("x E y")
    coupling = sp.symbols("Lambda", nonzero=True)
    radicand = 2 * coupling**2 * sp.cos(2 * x) - energy
    discriminant_factor = energy**2 - 4 * coupling**4
    lambda_two = lambda_two_coefficient(radicand, x, y)
    require_on_curve_zero(
        lambda_two
        - riccati_lambda_two_coefficient(radicand, x, y),
        y,
        radicand,
        "Mathieu lambda_2 versus direct Riccati recursion",
    )

    trigonometric_form = (
        -coupling**2 * sp.cos(2 * x) / y**3
        - sp.Rational(5, 2)
        * coupling**4
        * sp.sin(2 * x) ** 2
        / y**5
    )
    rational_form = (
        1 / (8 * y)
        + 3 * energy / (4 * y**3)
        + 5 * discriminant_factor / (8 * y**5)
    )
    require_on_curve_zero(
        lambda_two - trigonometric_form,
        y,
        radicand,
        "Mathieu trigonometric lambda_2",
    )
    require_on_curve_zero(
        lambda_two - rational_form,
        y,
        radicand,
        "Mathieu rational lambda_2",
    )

    energy_derivatives = [y]
    for _ in range(3):
        energy_derivatives.append(
            energy_derivative(
                energy_derivatives[-1], energy, y, radicand
            )
        )

    pf_primitive = coupling**2 * sp.sin(2 * x) / (2 * y)
    pf_form = (
        discriminant_factor * energy_derivatives[2] + y / 4
    )
    require_on_curve_zero(
        pf_form
        - curve_derivative(pf_primitive, x, y, radicand),
        y,
        radicand,
        "Mathieu Picard--Fuchs exact-form certificate",
    )

    third_order_form = (
        -sp.Rational(1, 4) * energy_derivatives[1]
        - 3 * energy * energy_derivatives[2]
        - sp.Rational(5, 3)
        * discriminant_factor
        * energy_derivatives[3]
    )
    require_on_curve_zero(
        lambda_two - third_order_form,
        y,
        radicand,
        "Mathieu pointwise third-order lambda_2 operator",
    )

    polynomial_form = (
        sp.Rational(1, 6) * energy_derivatives[1]
        + energy * energy_derivatives[2] / 3
    )
    polynomial_primitive = (
        -5 * coupling**2 * sp.sin(2 * x) / (12 * y**3)
    )
    require_on_curve_zero(
        lambda_two
        - polynomial_form
        - curve_derivative(
            polynomial_primitive, x, y, radicand
        ),
        y,
        radicand,
        "Mathieu reduced polynomial operator and certificate",
    )

    first_order_form = (
        sp.Rational(1, 6) * energy_derivatives[1]
        - energy * y / (12 * discriminant_factor)
    )
    first_order_primitive = (
        polynomial_primitive
        + energy * pf_primitive / (3 * discriminant_factor)
    )
    require_on_curve_zero(
        lambda_two
        - first_order_form
        - curve_derivative(
            first_order_primitive, x, y, radicand
        ),
        y,
        radicand,
        "Mathieu first-order reduced operator and certificate",
    )

    test_function = sp.Function("f")(energy)
    pf_on_test = (
        discriminant_factor * sp.diff(test_function, energy, 2)
        + test_function / 4
    )
    third_on_test = (
        -sp.diff(test_function, energy) / 4
        - 3 * energy * sp.diff(test_function, energy, 2)
        - sp.Rational(5, 3)
        * discriminant_factor
        * sp.diff(test_function, energy, 3)
    )
    polynomial_on_test = (
        sp.diff(test_function, energy) / 6
        + energy * sp.diff(test_function, energy, 2) / 3
    )
    first_order_on_test = (
        sp.diff(test_function, energy) / 6
        - energy * test_function / (12 * discriminant_factor)
    )
    require_zero(
        third_on_test
        - polynomial_on_test
        + sp.Rational(5, 3) * sp.diff(pf_on_test, energy),
        "Mathieu third-order minus polynomial left-ideal identity",
    )
    require_zero(
        polynomial_on_test
        - first_order_on_test
        - energy * pf_on_test / (3 * discriminant_factor),
        "Mathieu polynomial minus first-order PF identity",
    )

    w, capital_y = sp.symbols("w Y")
    cubic = coupling**2 * w**3 - energy * w**2 + coupling**2 * w
    cubic_discriminant = sp.factor(sp.discriminant(cubic, w))
    require_zero(
        cubic_discriminant - coupling**4 * discriminant_factor,
        "Mathieu cubic discriminant",
    )
    cubic_degree = sp.degree(cubic, w)
    generic_hyperelliptic_genus = (cubic_degree - 1) // 2
    if generic_hyperelliptic_genus != 1:
        raise RuntimeError(
            "Mathieu compact curve genus check failed: "
            f"degree={cubic_degree}, genus={generic_hyperelliptic_genus}"
        )
    square_root = sp.sqrt(discriminant_factor)
    root_plus = (energy + square_root) / (2 * coupling**2)
    root_minus = (energy - square_root) / (2 * coupling**2)
    require_zero(
        cubic.subs(w, root_plus),
        "Mathieu turning root w_plus",
    )
    require_zero(
        cubic.subs(w, root_minus),
        "Mathieu turning root w_minus",
    )
    require_zero(
        root_plus * root_minus - 1,
        "Mathieu reciprocal turning roots",
    )

    dw_dx = 2 * sp.I * w
    lambda_zero_in_w = capital_y / (2 * sp.I * w**2)
    require_zero(
        lambda_zero_in_w.subs(capital_y, w * y) * dw_dx - y,
        "Mathieu differential under w=exp(2 i x)",
    )
    capital_y_energy_derivative = -w**2 / (2 * capital_y)
    require_zero(
        capital_y_energy_derivative / (2 * sp.I * w**2)
        + 1 / (4 * sp.I * capital_y),
        "Mathieu energy derivative in exponential coordinates",
    )

    return {
        "cubic_discriminant": cubic_discriminant,
        "genus": sp.Integer(generic_hyperelliptic_genus),
        "lambda_two": lambda_two,
        "pf_primitive": pf_primitive,
        "polynomial_primitive": polynomial_primitive,
        "first_order_primitive": first_order_primitive,
    }


def mathieu_action_audit() -> dict[str, object]:
    """Check the elliptic-integral actions, limits, and cycle basis."""
    parameter = sp.symbols("m", positive=True)
    coupling = sp.symbols("Lambda", positive=True)
    energy = 4 * coupling**2 * parameter - 2 * coupling**2
    discriminant_factor = energy**2 - 4 * coupling**4
    elliptic_k = sp.elliptic_k
    elliptic_e = sp.elliptic_e

    action_minus = 8 * coupling * (
        elliptic_e(parameter)
        - (1 - parameter) * elliptic_k(parameter)
    )
    action_plus = 8 * coupling * (
        elliptic_e(1 - parameter)
        - parameter * elliptic_k(1 - parameter)
    )

    def derivative_in_energy(expression: sp.Expr) -> sp.Expr:
        return reduced(
            sp.diff(expression, parameter) / (4 * coupling**2)
        )

    action_minus_derivative = derivative_in_energy(action_minus)
    action_plus_derivative = derivative_in_energy(action_plus)
    require_zero(
        sp.expand_func(action_minus_derivative)
        - elliptic_k(parameter) / coupling,
        "Mathieu derivative of J_minus",
    )
    require_zero(
        sp.expand_func(action_plus_derivative)
        + elliptic_k(1 - parameter) / coupling,
        "Mathieu derivative of J_plus",
    )

    for action, label in (
        (action_minus, "minus"),
        (action_plus, "plus"),
    ):
        pf_residual = (
            discriminant_factor
            * derivative_in_energy(derivative_in_energy(action))
            + action / 4
        )
        require_zero(
            sp.expand_func(pf_residual),
            f"Mathieu Picard--Fuchs equation for J_{label}",
        )

    def quantum_operator(expression: sp.Expr) -> sp.Expr:
        return reduced(
            derivative_in_energy(expression) / 6
            + energy
            * derivative_in_energy(derivative_in_energy(expression))
            / 3
        )

    quantum_minus = reduced(
        (
            (1 - parameter) * elliptic_k(parameter)
            + (2 * parameter - 1) * elliptic_e(parameter)
        )
        / (12 * coupling * parameter * (1 - parameter))
    )
    quantum_plus = reduced(
        (
            (2 * parameter - 1) * elliptic_e(1 - parameter)
            - parameter * elliptic_k(1 - parameter)
        )
        / (12 * coupling * parameter * (1 - parameter))
    )
    require_zero(
        sp.expand_func(quantum_operator(action_minus) - quantum_minus),
        "Mathieu J_minus order-hbar**2 coefficient",
    )
    require_zero(
        sp.expand_func(quantum_operator(action_plus) - quantum_plus),
        "Mathieu J_plus order-hbar**2 coefficient",
    )

    limit_checks = (
        (
            sp.limit(action_minus / parameter, parameter, 0, dir="+"),
            2 * sp.pi * coupling,
            "Mathieu J_minus harmonic limit",
        ),
        (
            sp.limit(
                action_plus / (1 - parameter),
                parameter,
                1,
                dir="-",
            ),
            2 * sp.pi * coupling,
            "Mathieu J_plus harmonic limit",
        ),
        (
            sp.limit(quantum_minus, parameter, 0, dir="+"),
            sp.pi / (16 * coupling),
            "Mathieu J_minus quantum endpoint limit",
        ),
        (
            sp.limit(quantum_plus, parameter, 1, dir="-"),
            -sp.pi / (16 * coupling),
            "Mathieu J_plus quantum endpoint limit",
        ),
    )
    for actual, expected, label in limit_checks:
        require_zero(actual - expected, label)

    legendre_combination = (
        elliptic_k(parameter) * elliptic_e(1 - parameter)
        + elliptic_e(parameter) * elliptic_k(1 - parameter)
        - elliptic_k(parameter) * elliptic_k(1 - parameter)
    )
    require_zero(
        sp.expand_func(sp.diff(legendre_combination, parameter)),
        "constancy of the Legendre combination",
    )

    period_minus = sp.I * action_minus
    period_plus = action_plus
    period_minus_derivative = sp.I * action_minus_derivative
    period_plus_derivative = action_plus_derivative
    physical_wronskian = reduced(
        period_minus * period_plus_derivative
        - period_minus_derivative * period_plus
    )
    require_zero(
        sp.expand_func(
            physical_wronskian + 8 * sp.I * legendre_combination
        ),
        "Mathieu physical-cycle Wronskian reduction",
    )

    physical_intersection = sp.Matrix([[0, 2], [-2, 0]])
    change_of_basis = sp.Matrix(
        [[1, -sp.Rational(1, 2)], [0, sp.Rational(1, 2)]]
    )
    symplectic_intersection = reduced(
        change_of_basis.T * physical_intersection * change_of_basis
    )
    expected_intersection = sp.Matrix([[0, 1], [-1, 0]])
    if symplectic_intersection != expected_intersection:
        raise RuntimeError(
            "Mathieu symplectic change of basis failed: "
            f"matrix = {symplectic_intersection}"
        )

    period_a = period_minus
    period_b = reduced((period_plus - period_minus) / 2)
    period_b_derivative = reduced(
        (period_plus_derivative - period_minus_derivative) / 2
    )
    require_zero(
        period_plus - period_a - 2 * period_b,
        "Mathieu delta_plus = A + 2 B period relation",
    )
    symplectic_wronskian = reduced(
        period_a * period_b_derivative
        - period_minus_derivative * period_b
    )
    require_zero(
        symplectic_wronskian - physical_wronskian / 2,
        "Mathieu symplectic versus physical Wronskian",
    )

    return {
        "action_minus": action_minus,
        "action_plus": action_plus,
        "quantum_minus": quantum_minus,
        "quantum_plus": quantum_plus,
        "legendre_combination": legendre_combination,
        "physical_wronskian": physical_wronskian,
        "symplectic_wronskian": symplectic_wronskian,
        "symplectic_intersection": symplectic_intersection,
    }


def mathieu_numerical_audit() -> dict[str, mp.mpf]:
    """Independently compare Mathieu action formulas with quadrature."""
    mp.mp.dps = 70
    tolerance = mp.mpf("1e-45")
    coupling = mp.mpf("1.3")
    energy = mp.mpf("0.37")
    parameter = (energy + 2 * coupling**2) / (4 * coupling**2)
    if not (0 < parameter < 1):
        raise RuntimeError(f"Mathieu sample has invalid m={parameter}")

    def nonnegative_square_root(value: mp.mpf) -> mp.mpf:
        return mp.sqrt(max(mp.mpf("0"), value))

    allowed_half_width = mp.asin(mp.sqrt(parameter))
    allowed_left = mp.pi / 2 - allowed_half_width
    allowed_right = mp.pi / 2 + allowed_half_width
    action_minus_quadrature = 2 * (
        mp.quad(
            lambda point: nonnegative_square_root(
                energy - 2 * coupling**2 * mp.cos(2 * point)
            ),
            [allowed_left, mp.pi / 2],
        )
        + mp.quad(
            lambda point: nonnegative_square_root(
                energy - 2 * coupling**2 * mp.cos(2 * point)
            ),
            [mp.pi / 2, allowed_right],
        )
    )

    barrier_half_width = mp.asin(mp.sqrt(1 - parameter))
    action_plus_quadrature = 2 * (
        mp.quad(
            lambda point: nonnegative_square_root(
                2 * coupling**2 * mp.cos(2 * point) - energy
            ),
            [-barrier_half_width, 0],
        )
        + mp.quad(
            lambda point: nonnegative_square_root(
                2 * coupling**2 * mp.cos(2 * point) - energy
            ),
            [0, barrier_half_width],
        )
    )

    elliptic_k = mp.ellipk
    elliptic_e = mp.ellipe
    action_minus_formula = 8 * coupling * (
        elliptic_e(parameter)
        - (1 - parameter) * elliptic_k(parameter)
    )
    action_plus_formula = 8 * coupling * (
        elliptic_e(1 - parameter)
        - parameter * elliptic_k(1 - parameter)
    )
    require_close(
        action_minus_quadrature,
        action_minus_formula,
        "Mathieu allowed-cycle quadrature",
        tolerance,
    )
    require_close(
        action_plus_quadrature,
        action_plus_formula,
        "Mathieu barrier-cycle quadrature",
        tolerance,
    )

    legendre_combination = (
        elliptic_k(parameter) * elliptic_e(1 - parameter)
        + elliptic_e(parameter) * elliptic_k(1 - parameter)
        - elliptic_k(parameter) * elliptic_k(1 - parameter)
    )
    require_close(
        legendre_combination,
        mp.pi / 2,
        "Legendre relation at the Mathieu sample",
        tolerance,
    )

    period_minus = mp.j * action_minus_formula
    period_plus = action_plus_formula
    period_minus_derivative = (
        mp.j * elliptic_k(parameter) / coupling
    )
    period_plus_derivative = -elliptic_k(1 - parameter) / coupling
    physical_wronskian = (
        period_minus * period_plus_derivative
        - period_minus_derivative * period_plus
    )
    require_close(
        physical_wronskian,
        -4 * mp.pi * mp.j,
        "Mathieu physical-cycle Wronskian",
        tolerance,
    )

    period_a = period_minus
    period_b = (period_plus - period_minus) / 2
    period_b_derivative = (
        period_plus_derivative - period_minus_derivative
    ) / 2
    require_close(
        period_plus,
        period_a + 2 * period_b,
        "Mathieu numerical delta_plus = A + 2 B relation",
        tolerance,
    )
    symplectic_wronskian = (
        period_a * period_b_derivative
        - period_minus_derivative * period_b
    )
    require_close(
        symplectic_wronskian,
        -2 * mp.pi * mp.j,
        "Mathieu symplectic-basis Wronskian",
        tolerance,
    )

    return {
        "parameter": parameter,
        "action_minus": action_minus_formula,
        "action_plus": action_plus_formula,
        "physical_wronskian": physical_wronskian,
        "symplectic_wronskian": symplectic_wronskian,
    }


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

    airy = airy_audit()
    print("[PASS] Airy rational cover and formal-WKB normalization")
    print(f"  lambda_2 = ({airy['lambda_two']}) dx")
    print(f"  Ai first correction = {airy['ai_first_term']}")

    weber_local = weber_local_operator_audit()
    print("[PASS] Weber Picard--Fuchs and exact-form operators")
    print(f"  lambda_2 = ({weber_local['lambda_two']}) dx")

    weber_periods = weber_period_audit()
    print("[PASS] Weber residues, Bernoulli coefficients, and gamma tail")
    print(f"  closed period = {weber_periods['closed_period']}")
    print(
        "  relative coefficients = "
        f"{weber_periods['relative_coefficients']}"
    )

    mathieu_local = mathieu_curve_and_operator_audit()
    print("[PASS] Mathieu elliptic curve, PF certificate, and D_2 operators")
    print(f"  cubic discriminant = {mathieu_local['cubic_discriminant']}")

    mathieu_actions = mathieu_action_audit()
    print("[PASS] Mathieu elliptic actions, quantum periods, and limits")
    print(
        "  physical intersection form -> symplectic form: "
        f"{mathieu_actions['symplectic_intersection'].tolist()}"
    )

    mathieu_numeric = mathieu_numerical_audit()
    print("[PASS] Mathieu quadratures and symplectic cycle normalization")
    print(f"  sample m = {mp.nstr(mathieu_numeric['parameter'], 18)}")
    print(f"  J_minus = {mp.nstr(mathieu_numeric['action_minus'], 18)}")
    print(f"  J_plus = {mp.nstr(mathieu_numeric['action_plus'], 18)}")
    print(
        "  W(delta_minus, delta_plus) = "
        f"{mp.nstr(mathieu_numeric['physical_wronskian'], 18)}"
    )
    print(
        "  W(A, B) = "
        f"{mp.nstr(mathieu_numeric['symplectic_wronskian'], 18)}"
    )

    print()
    print("All Airy--Weber--Mathieu checks passed.")
    print("Scope: formal WKB and classical periods; no exact-WKB claim tested.")


if __name__ == "__main__":
    main()