#!/usr/bin/env python3
"""Audit Chapter 9 Page 7 transseries and ambiguity-cancellation formulas.

The checks are deliberately split across independent representations:

1. formal implicit inversion versus an exact Lambert-W root;
2. moving-action terms in the lateral symmetric-double-well displacement;
3. the signed Euler Stokes jump and balanced parameter shift;
4. quartic saddle geometry, Gaussian moments, an ODE recurrence, and
   hypergeometric Borel germs;
5. the quartic Borel cut, companion-sector large order, and the stable
   Bessel/quadrature evaluation;
6. extraction of a farther signed Borel action after subtracting the
   leading pole;
7. the instanton/anti-instanton event lattice behind the double-well
   neutral-sector selection rule.

Every failure raises an explicit exception.  No Python assert statements
are used, so normal and optimized execution perform the same audit.
These calculations do not prove resurgence, uniform Borel summability, or
exact-WKB applicability for a new spectral problem.
"""

from __future__ import annotations

import platform
import sys

import mpmath as mp
import sympy as sp


def reduced(expression: sp.Expr) -> sp.Expr:
    """Return a stable exact simplification."""

    return sp.factor(
        sp.cancel(
            sp.together(
                sp.gammasimp(
                    sp.powsimp(sp.simplify(expression), force=True)
                )
            )
        )
    )


def require(condition: bool, message: str) -> None:
    """Raise an optimization-safe failure unless condition is true."""

    if not condition:
        raise RuntimeError(message)


def require_zero(expression: sp.Expr, message: str) -> None:
    """Require an exact symbolic expression to vanish."""

    residual = reduced(expression)
    if residual != 0:
        residual = sp.simplify(sp.expand_complex(residual))
    if residual != 0:
        raise RuntimeError(f"{message}: residual = {residual}")


def require_close(
    actual: complex,
    expected: complex,
    tolerance: mp.mpf,
    message: str,
) -> None:
    """Compare high-precision numbers with a scale-aware tolerance."""

    error = abs(actual - expected)
    scale = max(mp.mpf(1), abs(actual), abs(expected))
    if error > tolerance * scale:
        raise RuntimeError(
            f"{message}: error {mp.nstr(error, 14)} exceeds "
            f"{mp.nstr(tolerance * scale, 14)}"
        )


def implicit_root_audit() -> None:
    """Compare recursive coefficients with the exact Lambert-W root."""

    q, lam, e_star = sp.symbols("q lambda E_star", nonzero=True)
    order = 7
    coefficients = [
        -((-n) ** (n - 1))
        * lam ** (n - 1)
        * sp.exp(n * lam * e_star)
        / sp.factorial(n)
        for n in range(1, order)
    ]
    candidate = e_star + sum(
        coefficients[n - 1] * q**n for n in range(1, order)
    )
    determinant = candidate - e_star + q * sp.exp(lam * candidate)
    residual_series = sp.series(determinant, q, 0, order).removeO()
    require_zero(
        residual_series,
        "Lambert-W coefficients solve the implicit determinant",
    )

    exact = e_star - sp.LambertW(
        lam * q * sp.exp(lam * e_star)
    ) / lam
    exact_series = sp.series(exact, q, 0, order).removeO()
    require_zero(
        exact_series - candidate,
        "recursive and exact Lambert-W series",
    )

    mp.mp.dps = 80
    lam_mp = mp.mpf("0.7")
    e_star_mp = mp.mpf("0.4")
    q_mp = mp.mpf("0.001")
    exact_mp = e_star_mp - mp.lambertw(
        lam_mp * q_mp * mp.exp(lam_mp * e_star_mp)
    ) / lam_mp
    require_close(
        exact_mp - e_star_mp
        + q_mp * mp.exp(lam_mp * exact_mp),
        0,
        mp.mpf("1e-72"),
        "80-decimal Lambert-W root residual",
    )


def moving_action_audit() -> None:
    """Check the first two exact double-well displacement coefficients."""

    q, hbar, epsilon, td_prime = sp.symbols(
        "q hbar epsilon td_prime", nonzero=True
    )
    sign = sp.symbols("s", nonzero=True)
    d_one = -epsilon * hbar / (2 * sp.pi)
    d_two = (
        -sign * sp.I * hbar / (4 * sp.pi)
        - hbar * td_prime / (8 * sp.pi**2)
    )
    displacement = d_one * q + d_two * q**2
    log_derivative = -td_prime / (2 * hbar)
    shifted_q = q * (1 + log_derivative * displacement)
    implicit_right = (
        -sign
        * sp.I
        * hbar
        / (2 * sp.pi)
        * sp.log(1 - sign * sp.I * epsilon * shifted_q)
    )
    for lateral_sign in (sp.Integer(1), sp.Integer(-1)):
        residual = sp.series(
            displacement.subs(sign, lateral_sign)
            - implicit_right.subs(sign, lateral_sign),
            q,
            0,
            3,
        ).removeO()
        residual = residual.subs(epsilon**2, 1)
        require_zero(
            residual,
            f"double-well displacement for sign {lateral_sign}",
        )


def euler_median_audit() -> None:
    """Verify the signed Euler jump and balanced parameter values."""

    mp.mp.dps = 80
    sigma = mp.mpf("0.375")
    for hbar in (mp.mpf("0.07"), mp.mpf("0.11"), mp.mpf("0.19")):
        flat = mp.exp(-1 / hbar)
        principal = flat * mp.ei(1 / hbar)
        upper = principal + mp.pi * 1j * flat
        lower = principal - mp.pi * 1j * flat
        require_close(
            upper - lower,
            2 * mp.pi * 1j * flat,
            mp.mpf("1e-72"),
            "Euler upper-minus-lower jump",
        )
        balanced_upper = upper + (sigma - mp.pi * 1j) * flat
        balanced_lower = lower + (sigma + mp.pi * 1j) * flat
        target = principal + sigma * flat
        require_close(
            balanced_upper,
            target,
            mp.mpf("1e-72"),
            "Euler upper balanced value",
        )
        require_close(
            balanced_lower,
            target,
            mp.mpf("1e-72"),
            "Euler lower balanced value",
        )
        require_close(
            mp.im(target),
            0,
            mp.mpf("1e-72"),
            "Euler median reality",
        )


def quartic_coefficient(n: int) -> sp.Rational:
    """Return the positive unstable quartic coefficient a_n."""

    return sp.Rational(
        sp.factorial(4 * n),
        96**n * sp.factorial(2 * n) * sp.factorial(n),
    )


def quartic_geometry_and_series_audit() -> None:
    """Check the saddle action and three independent coefficient routes."""

    x = sp.symbols("x")
    unstable_action = x**2 / 2 - x**4 / 24
    require_zero(
        sp.diff(unstable_action, x)
        - x * (1 - x**2 / 6),
        "quartic critical-point factorization",
    )
    saddles = (sp.Integer(0), sp.sqrt(6), -sp.sqrt(6))
    expected_hessians = (sp.Integer(1), sp.Integer(-2), sp.Integer(-2))
    for saddle, expected_hessian in zip(saddles, expected_hessians):
        require_zero(
            sp.diff(unstable_action, x).subs(x, saddle),
            f"quartic saddle {saddle}",
        )
        require_zero(
            sp.diff(unstable_action, x, 2).subs(x, saddle)
            - expected_hessian,
            f"quartic Hessian at {saddle}",
        )
    require_zero(
        unstable_action.subs(x, sp.sqrt(6))
        - unstable_action.subs(x, 0)
        - sp.Rational(3, 2),
        "quartic saddle action A=3/2",
    )

    a_recurrence = sp.Integer(1)
    b_ode = sp.Integer(1)
    xi = sp.symbols("xi")
    action = sp.Rational(3, 2)
    borel_unstable = sp.hyper(
        [sp.Rational(1, 4), sp.Rational(3, 4)],
        [sp.Integer(1)],
        xi / action,
    )
    borel_stable = sp.hyper(
        [sp.Rational(1, 4), sp.Rational(3, 4)],
        [sp.Integer(1)],
        -xi / action,
    )
    unstable_series = sp.series(borel_unstable, xi, 0, 14).removeO()
    stable_series = sp.series(borel_stable, xi, 0, 14).removeO()

    for n in range(14):
        direct = quartic_coefficient(n)
        require_zero(
            a_recurrence - direct,
            f"quartic Gaussian-moment recurrence at n={n}",
        )
        require_zero(
            b_ode - (-1) ** n * direct,
            f"stable ODE recurrence at n={n}",
        )
        require_zero(
            unstable_series.coeff(xi, n)
            - direct / sp.factorial(n),
            f"unstable hypergeometric Borel coefficient n={n}",
        )
        require_zero(
            stable_series.coeff(xi, n)
            - (-1) ** n * direct / sp.factorial(n),
            f"stable hypergeometric Borel coefficient n={n}",
        )
        ratio = sp.Rational(
            (4 * n + 1) * (4 * n + 3),
            24 * (n + 1),
        )
        a_recurrence *= ratio
        b_ode *= -ratio


def quartic_cut_audit() -> None:
    """Check the exact hypergeometric cut relation at high precision."""

    mp.mp.dps = 80
    action = mp.mpf(3) / 2
    epsilon = mp.mpf("1e-45")
    tolerance = mp.mpf("2e-38")
    for point in (mp.mpf("1.8"), mp.mpf("2.4"), mp.mpf("4.0")):
        upper = mp.hyp2f1(
            mp.mpf(1) / 4,
            mp.mpf(3) / 4,
            1,
            (point + 1j * epsilon) / action,
        )
        lower = mp.hyp2f1(
            mp.mpf(1) / 4,
            mp.mpf(3) / 4,
            1,
            (point - 1j * epsilon) / action,
        )
        companion = mp.hyp2f1(
            mp.mpf(1) / 4,
            mp.mpf(3) / 4,
            1,
            -(point - action) / action,
        )
        require_close(
            upper - lower,
            1j * mp.sqrt(2) * companion,
            tolerance,
            f"quartic Borel cut at xi={point}",
        )


def quartic_large_order_audit() -> None:
    """Recover the saddle action and test companion-loop corrections."""

    mp.mp.dps = 80
    action = mp.mpf(3) / 2
    amplitude = 1 / (mp.pi * mp.sqrt(2))

    def a_mp(n: int) -> mp.mpf:
        return (
            mp.factorial(4 * n)
            / (
                mp.mpf(96) ** n
                * mp.factorial(2 * n)
                * mp.factorial(n)
            )
        )

    n = 100
    action_estimate = n * a_mp(n) / a_mp(n + 1)
    amplitude_estimate = action**n * a_mp(n) / mp.gamma(n)
    require_close(
        action_estimate,
        action,
        mp.mpf("3e-5"),
        "quartic action ratio estimate",
    )
    require_close(
        amplitude_estimate,
        amplitude,
        mp.mpf("5e-4"),
        "quartic leading-amplitude estimate",
    )

    for n in (40, 70, 100):
        predicted = mp.mpf(0)
        for k in range(6):
            b_k = (-1) ** k * a_mp(k)
            predicted += (
                amplitude
                * b_k
                * mp.gamma(n - k)
                / action ** (n - k)
            )
        require_close(
            a_mp(n),
            predicted,
            mp.mpf("2e-7"),
            f"quartic six-loop large order at n={n}",
        )


def stable_quartic_audit() -> None:
    """Compare direct stable-integral quadrature with its Bessel form."""

    mp.mp.dps = 80
    for coupling in (
        mp.mpf("0.08"),
        mp.mpf("0.15"),
        mp.mpf("0.30"),
    ):
        integrand = lambda x: mp.exp(
            -(x**2 / 2 + x**4 / 24) / coupling
        )
        quadrature = (
            2
            * mp.quad(integrand, [0, mp.inf])
            / mp.sqrt(2 * mp.pi * coupling)
        )
        bessel = (
            mp.sqrt(3 / (2 * mp.pi * coupling))
            * mp.exp(3 / (4 * coupling))
            * mp.besselk(mp.mpf(1) / 4, 3 / (4 * coupling))
        )
        require_close(
            quadrature,
            bessel,
            mp.mpf("1e-68"),
            f"stable quartic integral at g={coupling}",
        )


def farther_action_audit() -> None:
    """Extract a signed farther action after exact leading subtraction."""

    action_a = sp.Rational(3, 2)
    action_b = sp.Rational(-7, 2)
    residue_a = sp.Rational(5, 7)
    residue_b = sp.Rational(-2, 5)

    def normalized_coefficient(n: int) -> sp.Expr:
        return residue_a / action_a**n + residue_b / action_b**n

    for n in (8, 17, 31):
        residual_n = (
            normalized_coefficient(n) - residue_a / action_a**n
        )
        residual_next = (
            normalized_coefficient(n + 1)
            - residue_a / action_a ** (n + 1)
        )
        require_zero(
            residual_n / residual_next - action_b,
            f"farther signed action at n={n}",
        )
        require_zero(
            residual_n * action_b**n - residue_b,
            f"farther pole residue at n={n}",
        )

    mp.mp.dps = 80
    n = 60
    a = mp.mpf(3) / 2
    b = -mp.mpf(7) / 2
    ra = mp.mpf(5) / 7
    rb = -mp.mpf(2) / 5
    coefficient_n = mp.gamma(n) * (ra / a**n + rb / b**n)
    coefficient_next = mp.gamma(n + 1) * (
        ra / a ** (n + 1) + rb / b ** (n + 1)
    )
    estimate = n * coefficient_n / coefficient_next
    require_close(
        estimate,
        a,
        mp.mpf("1e-20"),
        "dominant two-pole action estimate",
    )


def event_lattice_audit() -> None:
    """Check the double-well event/charge lattice and its neutral column."""

    neutral_nonzero = []
    for total in range(8):
        charges = list(range(-total, total + 1, 2))
        require(
            len(charges) == total + 1,
            f"unexpected number of lattice nodes at total={total}",
        )
        for charge in charges:
            instantons = (total + charge) // 2
            anti_instantons = (total - charge) // 2
            require(
                instantons >= 0 and anti_instantons >= 0,
                f"negative event count at ({total}, {charge})",
            )
            require(
                instantons + anti_instantons == total,
                f"total-event mismatch at ({total}, {charge})",
            )
            require(
                instantons - anti_instantons == charge,
                f"charge mismatch at ({total}, {charge})",
            )
            if charge == 0 and total > 0:
                neutral_nonzero.append(total)
            for step in (-1, 1):
                next_charge = charge + step
                require(
                    next_charge in range(
                        -(total + 1), total + 2, 2
                    ),
                    f"missing event-lattice edge from ({total}, {charge})",
                )
    require(
        neutral_nonzero[0] == 2,
        "the first nonzero neutral sector must contain two events",
    )
    require(
        all(total % 2 == 0 for total in neutral_nonzero),
        "neutral sectors must have even event number",
    )


def main() -> None:
    """Run every audit and print a compact status ledger."""

    print("Advanced ODE Chapter 9 Page 7 transseries audit")
    print(f"Python: {platform.python_version()}")
    print(f"SymPy: {sp.__version__}")
    print(f"mpmath: {mp.__version__}")

    audits = (
        (
            "implicit Lambert-W root and formal coefficients",
            implicit_root_audit,
        ),
        (
            "moving-action double-well displacement",
            moving_action_audit,
        ),
        (
            "Euler lateral jump and balanced median",
            euler_median_audit,
        ),
        (
            "quartic saddles, moments, ODE, and Borel germs",
            quartic_geometry_and_series_audit,
        ),
        (
            "80-decimal quartic hypergeometric cut",
            quartic_cut_audit,
        ),
        (
            "quartic companion-sector large order",
            quartic_large_order_audit,
        ),
        (
            "stable quartic quadrature and Bessel form",
            stable_quartic_audit,
        ),
        (
            "dominant and farther signed Borel actions",
            farther_action_audit,
        ),
        (
            "double-well neutral event lattice",
            event_lattice_audit,
        ),
    )

    for label, audit in audits:
        audit()
        print(f"PASS {label}")

    print(
        "LIMITATION: these identities do not prove resurgence, "
        "summability, thimble completeness, or exact-WKB applicability "
        "for a new spectral problem."
    )
    print("All instanton-transseries checks passed.")


if __name__ == "__main__":
    try:
        main()
    except Exception as error:
        print(f"FAIL {error}", file=sys.stderr)
        raise