#!/usr/bin/env python3
"""Audit Chapter 9's Borel--Laplace normalization and lateral signs.

The book uses the shifted transform

    B[a_0 + sum_{n>=1} a_n hbar**n]
        = sum_{n>=1} a_n xi**(n-1) / Gamma(n)

and the inverse

    S_theta f = a_0 + integral exp(-xi/hbar) B[f](xi) dxi.

There is no 1/hbar prefactor in this convention.  Every test raises
RuntimeError explicitly, so Python optimization cannot disable the audit.

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 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(residual)
    if residual != 0:
        raise RuntimeError(f"{label} failed: residual = {residual}")


def require_true(condition: bool, label: str) -> None:
    """Raise a visible failure unless a Boolean condition holds."""
    if not condition:
        raise RuntimeError(f"{label} failed")


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 shifted_borel(
    expression: sp.Expr,
    hbar: sp.Symbol,
    xi: sp.Symbol,
) -> sp.Expr:
    """Apply the shifted Borel transform to a polynomial in hbar."""
    polynomial = sp.Poly(sp.expand(expression), hbar)
    result = sp.Integer(0)
    for (power,), coefficient in polynomial.terms():
        if power == 0:
            continue
        if power < 0:
            raise RuntimeError("ordinary shifted Borel transform received a pole")
        result += coefficient * xi ** (power - 1) / sp.gamma(power)
    return sp.expand(result)


def monomial_and_algebra_audit() -> dict[str, sp.Expr]:
    """Check gamma moments and elementary Borel algebra identities."""
    xi = sp.symbols("xi", nonnegative=True)
    hbar = sp.symbols("hbar", positive=True)

    for n in range(1, 9):
        shifted_moment = sp.integrate(
            sp.exp(-xi / hbar) * xi ** (n - 1) / sp.gamma(n),
            (xi, 0, sp.oo),
        )
        require_zero(shifted_moment - hbar**n, f"shifted moment n={n}")

        unshifted_moment = sp.integrate(
            sp.exp(-xi / hbar) * xi**n / sp.factorial(n) / hbar,
            (xi, 0, sp.oo),
        )
        require_zero(unshifted_moment - hbar**n, f"unshifted moment n={n}")

    a0, a1, a2, a3 = sp.symbols("a0 a1 a2 a3")
    b0, b1, b2 = sp.symbols("b0 b1 b2")
    formal_f = a0 + a1 * hbar + a2 * hbar**2 + a3 * hbar**3
    formal_g = b0 + b1 * hbar + b2 * hbar**2
    borel_f = shifted_borel(formal_f, hbar, xi)
    borel_g = shifted_borel(formal_g, hbar, xi)

    require_zero(
        shifted_borel(hbar**2 * sp.diff(formal_f, hbar), hbar, xi)
        - xi * borel_f,
        "hbar derivative dictionary",
    )

    eta = sp.symbols("eta")
    convolution = sp.integrate(
        borel_f.subs(xi, eta) * borel_g.subs(xi, xi - eta),
        (eta, 0, xi),
    )
    product_dictionary = a0 * borel_g + b0 * borel_f + convolution
    require_zero(
        shifted_borel(formal_f * formal_g, hbar, xi)
        - product_dictionary,
        "product-convolution dictionary",
    )

    return {
        "highest_moment": hbar**8,
        "borel_f": borel_f,
    }


def alternating_euler_audit() -> dict[str, object]:
    """Check the unobstructed Euler sum, ODE, and factorial remainder."""
    xi = sp.symbols("xi")
    order = 10
    alternating_borel = sum((-xi) ** n for n in range(order))
    expected_series = sp.series(1 / (1 + xi), xi, 0, order).removeO()
    require_zero(
        alternating_borel - expected_series,
        "alternating Euler Borel germ",
    )

    mp.mp.dps = 75

    def euler_sum(hbar: mp.mpc) -> mp.mpc:
        return mp.quad(
            lambda value: mp.exp(-value / hbar) / (1 + value),
            [0, 1, mp.inf],
        )

    samples = (
        mp.mpf("0.17"),
        mp.mpf("0.09"),
        mp.mpc("0.12", "0.035"),
    )
    for index, hbar in enumerate(samples):
        require_true(mp.re(1 / hbar) > 0, f"Euler domain sample {index}")
        integral = euler_sum(hbar)
        special = mp.exp(1 / hbar) * mp.e1(1 / hbar)
        require_close(
            integral,
            special,
            f"alternating Euler integral sample {index}",
            mp.mpf("2e-68"),
        )

        derivative = mp.diff(euler_sum, hbar)
        require_close(
            hbar**2 * derivative + integral,
            hbar,
            f"alternating Euler ODE sample {index}",
            mp.mpf("2e-55"),
        )

    positive_hbar = mp.mpf("0.08")
    exact = euler_sum(positive_hbar)
    remainder_ratios = []
    for truncation in range(1, 15):
        partial = mp.fsum(
            (-1) ** n
            * mp.factorial(n)
            * positive_hbar ** (n + 1)
            for n in range(truncation)
        )
        bound = (
            mp.factorial(truncation)
            * positive_hbar ** (truncation + 1)
        )
        ratio = abs(exact - partial) / bound
        require_true(
            ratio <= 1 + mp.mpf("1e-60"),
            f"Euler factorial remainder N={truncation}",
        )
        remainder_ratios.append(ratio)

    return {
        "sample_sum": euler_sum(mp.mpf("0.1")),
        "largest_remainder_ratio": max(remainder_ratios),
    }


def lateral_euler_audit() -> dict[str, object]:
    """Check the simple-pole orientation and the Euler lateral jump."""
    xi, hbar = sp.symbols("xi hbar", positive=True)
    integrand = sp.exp(-xi / hbar) / (1 - xi)
    pole_residue = sp.residue(integrand, xi, 1)
    require_zero(
        pole_residue + sp.exp(-1 / hbar),
        "nonalternating Euler pole residue",
    )

    jump = -2 * sp.pi * sp.I * pole_residue
    require_zero(
        jump - 2 * sp.pi * sp.I * sp.exp(-1 / hbar),
        "nonalternating Euler signed jump",
    )
    require_zero(
        hbar**2 * sp.diff(jump, hbar) - jump,
        "nonalternating Euler homogeneous jump equation",
    )

    mp.mp.dps = 75
    hbar_value = mp.mpf("0.14")
    radius = mp.mpf("0.23")

    def laplace_integrand(value: mp.mpc) -> mp.mpc:
        return mp.exp(-value / hbar_value) / (1 - value)

    def arc_integrand(angle: mp.mpf) -> mp.mpc:
        value = 1 + radius * mp.exp(1j * angle)
        derivative = 1j * radius * mp.exp(1j * angle)
        return laplace_integrand(value) * derivative

    upper_arc = mp.quad(arc_integrand, [mp.pi, 0])
    lower_arc = mp.quad(arc_integrand, [-mp.pi, 0])
    numerical_jump = upper_arc - lower_arc
    expected_jump = 2 * mp.pi * 1j * mp.exp(-1 / hbar_value)
    require_close(
        numerical_jump,
        expected_jump,
        "numerical upper-minus-lower jump",
        mp.mpf("2e-66"),
    )

    pv_expected = mp.exp(-1 / hbar_value) * mp.ei(1 / hbar_value)
    pv_errors = []
    for radius_text in ("0.03", "0.01", "0.003"):
        local_radius = mp.mpf(radius_text)
        pv_cutoff = (
            mp.quad(laplace_integrand, [0, 1 - local_radius])
            + mp.quad(laplace_integrand, [1 + local_radius, mp.inf])
        )
        pv_errors.append(abs(pv_cutoff - pv_expected))
    require_true(
        pv_errors[-1] < pv_errors[0],
        "principal-value convergence",
    )
    require_true(
        pv_errors[-1] < mp.mpf("1e-4"),
        "principal-value numerical accuracy",
    )

    return {
        "pole_residue": pole_residue,
        "jump_at_hbar_0.14": expected_jump,
        "pv_last_error": pv_errors[-1],
    }


def wkb_parity_and_period_audit() -> dict[str, sp.Expr]:
    """Check the odd WKB model and the period/Voros Borel dictionary."""
    xi, hbar, action = sp.symbols("xi hbar action", positive=True)
    order = 7
    wkb_borel_partial = sum(
        xi ** (2 * k - 2) / action ** (2 * k - 1)
        for k in range(1, order)
    )
    wkb_borel_expected = sp.series(
        action / (action**2 - xi**2),
        xi,
        0,
        2 * order - 2,
    ).removeO()
    require_zero(
        wkb_borel_partial - wkb_borel_expected,
        "odd WKB model Borel germ",
    )

    wkb_borel = action / (action**2 - xi**2)
    positive_residue = sp.residue(wkb_borel, xi, action)
    negative_residue = sp.residue(wkb_borel, xi, -action)
    require_zero(
        positive_residue + sp.Rational(1, 2),
        "positive action residue",
    )
    require_zero(
        negative_residue - sp.Rational(1, 2),
        "negative action residue",
    )
    wkb_jump = (
        -2
        * sp.pi
        * sp.I
        * sp.exp(-action / hbar)
        * positive_residue
    )
    require_zero(
        wkb_jump - sp.pi * sp.I * sp.exp(-action / hbar),
        "odd WKB model jump",
    )

    p2, p4, p6, p8 = sp.symbols("Pi2 Pi4 Pi6 Pi8")
    period_tail = (
        p2 * hbar**2
        + p4 * hbar**4
        + p6 * hbar**6
        + p8 * hbar**8
    )
    quantum_voros = sp.expand(period_tail / hbar)
    borel_period = shifted_borel(period_tail, hbar, xi)
    borel_voros = shifted_borel(quantum_voros, hbar, xi)
    eta = sp.symbols("eta")
    integrated_voros = sp.integrate(
        borel_voros.subs(xi, eta),
        (eta, 0, xi),
    )
    require_zero(
        borel_period - integrated_voros,
        "period-Voros Borel integration",
    )

    return {
        "positive_residue": positive_residue,
        "borel_period": borel_period,
        "borel_voros": borel_voros,
    }


def directional_growth_audit() -> dict[str, object]:
    """Check the complex Laplace convergence inequality, not just phases."""
    mp.mp.dps = 75
    theta = mp.mpf("0.63")
    kappa = mp.mpf("2.4")
    hbar = mp.mpf("0.075") * mp.exp(1j * mp.mpf("0.22"))
    decay_rate = mp.re(mp.exp(1j * theta) / hbar)
    require_true(decay_rate > kappa, "compatible complex Laplace domain")

    def ray_integrand(t: mp.mpf) -> mp.mpc:
        xi = mp.exp(1j * theta) * t
        borel_value = mp.exp(kappa * t)
        return (
            mp.exp(-xi / hbar)
            * borel_value
            * mp.exp(1j * theta)
        )

    numerical = mp.quad(ray_integrand, [0, 1, mp.inf])
    expected = (
        mp.exp(1j * theta)
        / (mp.exp(1j * theta) / hbar - kappa)
    )
    require_close(
        numerical,
        expected,
        "complex directional Laplace integral",
        mp.mpf("2e-66"),
    )

    incompatible_hbar = (
        mp.mpf("0.4")
        * mp.exp(1j * (theta + mp.mpf("1.45")))
    )
    incompatible_rate = mp.re(
        mp.exp(1j * theta) / incompatible_hbar
    )
    require_true(
        incompatible_rate <= kappa,
        "incompatible complex Laplace domain",
    )

    return {
        "decay_rate": decay_rate,
        "integral": numerical,
    }


def main() -> None:
    """Run every audit and print a compact deterministic report."""
    moments = monomial_and_algebra_audit()
    alternating = alternating_euler_audit()
    lateral = lateral_euler_audit()
    wkb = wkb_parity_and_period_audit()
    directional = directional_growth_audit()

    print("Advanced ODE Chapter 9 Page 1 audit")
    print(f"Python: {platform.python_version()}")
    print(f"SymPy: {sp.__version__}")
    print(f"mpmath: {mp.__version__}")
    print("PASS monomial moments through n=8")
    print(f"PASS Borel algebra: B[f] = {moments['borel_f']}")
    print(
        "PASS alternating Euler: S_0(0.1) = "
        f"{mp.nstr(alternating['sample_sum'], 18)}"
    )
    print(
        "PASS factorial remainder: largest ratio = "
        f"{mp.nstr(alternating['largest_remainder_ratio'], 12)}"
    )
    print(f"PASS lateral pole residue: {lateral['pole_residue']}")
    print(
        "PASS lateral jump at hbar=0.14: "
        f"{mp.nstr(lateral['jump_at_hbar_0.14'], 18)}"
    )
    print(
        "PASS odd WKB residues: Res(+A) = "
        f"{wkb['positive_residue']}"
    )
    print(
        "PASS period/Voros Borel relation: "
        f"B[Vq] = {wkb['borel_voros']}"
    )
    print(
        "PASS complex convergence rate: "
        f"{mp.nstr(directional['decay_rate'], 16)} > 2.4"
    )
    print("All Borel--Laplace checks passed.")


if __name__ == "__main__":
    main()