#!/usr/bin/env python3
"""Exact audit of Picard--Fuchs and WKB differential-operator formulas.

The checks use the book convention

    hbar**2 psi'' = R psi,
    lambda_0 = y dz,
    y**2 = R_0,
    Omega = sum_{k >= 0} hbar**(2*k) lambda_(2*k).

The script verifies:

1. the energy-derivative ladder for lambda_0;
2. the Legendre-family Picard--Fuchs exact-form certificate;
3. the universal source-free reduction of lambda_2;
4. the discriminant and classical Picard--Fuchs certificate for
   R_0 = z**2 + g*z**4 - E;
5. the exact order-hbar**2 quantum differential operator;
6. its reduced representative modulo the Picard--Fuchs operator; and
7. the closed Weber limit at g = 0.

Every check raises ``RuntimeError`` explicitly, so optimized Python executes
the same audit.  These are algebraic identities of meromorphic forms and
formal periods; no cycle continuation, endpoint regularization, Borel
summation, or quantization condition is tested.

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

from __future__ import annotations

import platform

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


def require_on_curve_zero(
    expression: sp.Expr,
    cover_coordinate: sp.Symbol,
    radicand: sp.Expr,
    label: str,
) -> None:
    """Raise unless a rational expression vanishes modulo y**2-radicand."""
    numerator, _ = sp.cancel(sp.together(expression)).as_numer_denom()
    numerator_polynomial = sp.Poly(
        sp.expand(numerator),
        cover_coordinate,
    )
    curve_relation = sp.Poly(
        cover_coordinate**2 - radicand,
        cover_coordinate,
    )
    remainder = numerator_polynomial.rem(curve_relation).as_expr()
    require_zero(remainder, label)


def total_curve_derivative(
    expression: sp.Expr,
    base_coordinate: sp.Symbol,
    cover_coordinate: sp.Symbol,
    radicand: sp.Expr,
) -> sp.Expr:
    """Differentiate along y**2=radicand with respect to the base variable."""
    dy_dz = sp.diff(radicand, base_coordinate) / (2 * cover_coordinate)
    return reduced(
        sp.diff(expression, base_coordinate)
        + sp.diff(expression, cover_coordinate) * dy_dz
    )


def total_parameter_derivative(
    expression: sp.Expr,
    parameter: sp.Symbol,
    cover_coordinate: sp.Symbol,
    radicand: sp.Expr,
) -> sp.Expr:
    """Differentiate along y**2=radicand with respect to a parameter."""
    dy_du = sp.diff(radicand, parameter) / (2 * cover_coordinate)
    return reduced(
        sp.diff(expression, parameter)
        + sp.diff(expression, cover_coordinate) * dy_du
    )


def energy_derivative_ladder_audit() -> list[sp.Expr]:
    """Verify partial_E**n(y dz) through n=6 for y**2=V-E."""
    y = sp.symbols("y", nonzero=True)
    current = y
    coefficients: list[sp.Expr] = []

    # At fixed z, dy/dE=-1/(2y).  No explicit E dependence is needed for
    # this local differential-algebra calculation.
    for n in range(1, 7):
        current = reduced(-sp.diff(current, y) / (2 * y))
        expected = -sp.factorial2(2 * n - 3) / (
            2**n * y ** (2 * n - 1)
        )
        require_zero(
            current - expected,
            f"energy derivative ladder at n={n}",
        )
        coefficients.append(current)

    return coefficients


def legendre_picard_fuchs_audit() -> dict[str, sp.Expr]:
    """Verify the Legendre derivative ladder and exact PF certificate."""
    x, t = sp.symbols("x t")
    capital_y = sp.symbols("Y", nonzero=True)
    legendre_polynomial = x * (x - 1) * (x - t)
    omega = 1 / capital_y

    omega_t = total_parameter_derivative(
        omega,
        t,
        capital_y,
        legendre_polynomial,
    )
    omega_tt = total_parameter_derivative(
        omega_t,
        t,
        capital_y,
        legendre_polynomial,
    )
    require_on_curve_zero(
        omega_t - omega / (2 * (x - t)),
        capital_y,
        legendre_polynomial,
        "Legendre first parameter derivative",
    )
    require_on_curve_zero(
        omega_tt - 3 * omega / (4 * (x - t) ** 2),
        capital_y,
        legendre_polynomial,
        "Legendre second parameter derivative",
    )

    pf_left = reduced(
        t * (1 - t) * omega_tt
        + (1 - 2 * t) * omega_t
        - omega / 4
    )
    certificate = x * (x - 1) / (2 * (x - t) * capital_y)
    certificate_derivative = total_curve_derivative(
        certificate,
        x,
        capital_y,
        legendre_polynomial,
    )
    require_on_curve_zero(
        pf_left - certificate_derivative,
        capital_y,
        legendre_polynomial,
        "Legendre Picard--Fuchs exact-form certificate",
    )

    discriminant = reduced(sp.discriminant(legendre_polynomial, x))
    require_zero(
        discriminant - t**2 * (t - 1) ** 2,
        "Legendre discriminant",
    )

    return {
        "omega_t": omega_t,
        "omega_tt": omega_tt,
        "discriminant": discriminant,
    }


def universal_lambda_two_audit() -> dict[str, sp.Expr]:
    """Verify two equivalent source-free reductions of lambda_2."""
    z = sp.symbols("z")
    y = sp.symbols("y", nonzero=True)
    radicand = sp.Function("R")(z)
    first = sp.diff(radicand, z)
    second = sp.diff(radicand, z, 2)

    lambda_two = second / (8 * y**3) - 5 * first**2 / (32 * y**5)
    derivative = reduced(
        second / y**3 - 3 * first**2 / (2 * y**5)
    )

    # First form: isolate the square of R'.
    require_zero(
        lambda_two - derivative / 8 - first**2 / (32 * y**5),
        "universal lambda_2 reduction to (R')**2",
    )

    # Second form: the normalization used in the quartic calculation.
    require_zero(
        lambda_two - second / (48 * y**3) - 5 * derivative / 48,
        "universal lambda_2 reduction to R''",
    )

    return {
        "lambda_two": lambda_two,
        "exact_derivative": derivative,
    }


def anharmonic_quartic_audit() -> dict[str, object]:
    """Verify the quartic PF and order-hbar**2 operator certificates."""
    z, energy, coupling = sp.symbols("z E g")
    y = sp.symbols("y", nonzero=True)
    radicand = z**2 + coupling * z**4 - energy
    radicand_z = sp.diff(radicand, z)

    discriminant = reduced(sp.discriminant(radicand, z))
    expected_discriminant = -16 * energy * coupling * (
        1 + 4 * coupling * energy
    ) ** 2
    require_zero(
        discriminant - expected_discriminant,
        "anharmonic quartic discriminant",
    )

    lambda_zero = y
    energy_two = -1 / (4 * y**3)
    energy_three = -3 / (8 * y**5)

    # [4E(1+4gE) partial_E**2 + 3g] lambda_0 = d_z F_PF.
    pf_left = reduced(
        4
        * energy
        * (1 + 4 * coupling * energy)
        * energy_two
        + 3 * coupling * lambda_zero
    )
    pf_certificate = (
        z
        * ((1 + coupling * z**2) ** 2 + coupling * energy)
        / y
    )
    pf_certificate_derivative = total_curve_derivative(
        pf_certificate,
        z,
        y,
        radicand,
    )
    require_on_curve_zero(
        pf_left - pf_certificate_derivative,
        y,
        radicand,
        "anharmonic quartic Picard--Fuchs certificate",
    )

    # Generate p_2 independently from the source-free Riccati recurrence.
    square_root = sp.sqrt(radicand)
    p_one = reduced(-sp.diff(square_root, z) / (2 * square_root))
    p_two = reduced(
        -(sp.diff(p_one, z) + p_one**2) / (2 * square_root)
    )
    printed_lambda_two = reduced(
        sp.diff(radicand, z, 2) / (8 * radicand ** sp.Rational(3, 2))
        - 5
        * radicand_z**2
        / (32 * radicand ** sp.Rational(5, 2))
    )
    require_zero(
        p_two - printed_lambda_two,
        "quartic Riccati generation of lambda_2",
    )

    lambda_two = (
        sp.diff(radicand, z, 2) / (8 * y**3)
        - 5 * radicand_z**2 / (32 * y**5)
    )
    operator_two = (3 + 20 * coupling * energy) / 6
    operator_three = 2 * energy * (1 + 4 * coupling * energy) / 3
    operator_form = reduced(
        operator_two * energy_two + operator_three * energy_three
    )
    quantum_certificate = (
        z * (3 - 4 * coupling * energy + 8 * coupling * z**2)
        / (24 * y**3)
    )
    quantum_certificate_derivative = total_curve_derivative(
        quantum_certificate,
        z,
        y,
        radicand,
    )
    require_on_curve_zero(
        lambda_two - operator_form - quantum_certificate_derivative,
        y,
        radicand,
        "anharmonic quartic lambda_2 differential-operator certificate",
    )

    return {
        "z": z,
        "E": energy,
        "g": coupling,
        "y": y,
        "R": radicand,
        "discriminant": discriminant,
        "operator_two": reduced(operator_two),
        "operator_three": reduced(operator_three),
        "pf_certificate": pf_certificate,
        "quantum_certificate": quantum_certificate,
    }


def reduced_operator_audit(data: dict[str, object]) -> dict[str, sp.Expr]:
    """Check D_2 modulo the noncommutative PF left ideal on periods."""
    energy = data["E"]
    coupling = data["g"]
    operator_two = data["operator_two"]
    operator_three = data["operator_three"]
    if not isinstance(energy, sp.Symbol) or not isinstance(coupling, sp.Symbol):
        raise RuntimeError("invalid quartic parameter data")
    if not isinstance(operator_two, sp.Expr) or not isinstance(
        operator_three, sp.Expr
    ):
        raise RuntimeError("invalid quartic operator data")

    period_zero, period_one = sp.symbols("Pi_0 Pi_1")
    pf_leading = 4 * energy * (1 + 4 * coupling * energy)
    period_two = reduced(-3 * coupling * period_zero / pf_leading)
    period_three = reduced(
        -(
            4 * (1 + 8 * coupling * energy) * period_two
            + 3 * coupling * period_one
        )
        / pf_leading
    )

    high_order_value = reduced(
        operator_two * period_two + operator_three * period_three
    )
    reduced_value = reduced(
        -coupling * period_one / 2
        + coupling
        * (1 + 12 * coupling * energy)
        * period_zero
        / (8 * energy * (1 + 4 * coupling * energy))
    )
    require_zero(
        high_order_value - reduced_value,
        "quartic high- and low-order D_2 representatives",
    )

    # Explicit left-ideal certificate.  Here dL means partial_E composed
    # with L, so differentiation acts on both its coefficients and period.
    period_two_symbol, period_three_symbol = sp.symbols("Pi_2 Pi_3")
    pf_value = reduced(
        pf_leading * period_two_symbol + 3 * coupling * period_zero
    )
    derivative_pf_value = reduced(
        pf_leading * period_three_symbol
        + 4 * (1 + 8 * coupling * energy) * period_two_symbol
        + 3 * coupling * period_one
    )
    high_symbolic = reduced(
        operator_two * period_two_symbol
        + operator_three * period_three_symbol
    )
    low_symbolic = reduced(
        -coupling * period_one / 2
        + coupling
        * (1 + 12 * coupling * energy)
        * period_zero
        / (8 * energy * (1 + 4 * coupling * energy))
    )
    left_ideal_value = reduced(
        derivative_pf_value / 6
        - (1 + 12 * coupling * energy)
        * pf_value
        / (24 * energy * (1 + 4 * coupling * energy))
    )
    require_zero(
        high_symbolic - low_symbolic - left_ideal_value,
        "quartic D_2 left-ideal certificate",
    )

    return {
        "high_order_value": high_order_value,
        "reduced_value": reduced_value,
    }


def weber_limit_audit(data: dict[str, object]) -> dict[str, sp.Expr]:
    """Check the g=0 closed Weber period and exact lambda_2 primitive."""
    z = data["z"]
    energy = data["E"]
    y = data["y"]
    if not isinstance(z, sp.Symbol) or not isinstance(energy, sp.Symbol):
        raise RuntimeError("invalid Weber coordinate data")
    if not isinstance(y, sp.Symbol):
        raise RuntimeError("invalid Weber cover data")

    weber_radicand = z**2 - energy
    weber_lambda_two = (
        sp.diff(weber_radicand, z, 2) / (8 * y**3)
        - 5 * sp.diff(weber_radicand, z) ** 2 / (32 * y**5)
    )
    weber_primitive = 5 * z / (24 * y**3) - z / (24 * energy * y)
    weber_primitive_derivative = total_curve_derivative(
        weber_primitive,
        z,
        y,
        weber_radicand,
    )
    require_on_curve_zero(
        weber_lambda_two - weber_primitive_derivative,
        y,
        weber_radicand,
        "Weber lambda_2 exact primitive",
    )

    closed_period = sp.I * sp.pi * energy
    pf_period_value = reduced(4 * energy * sp.diff(closed_period, energy, 2))
    require_zero(pf_period_value, "Weber closed classical PF equation")

    high_operator_value = reduced(
        sp.Rational(1, 2) * sp.diff(closed_period, energy, 2)
        + 2 * energy * sp.diff(closed_period, energy, 3) / 3
    )
    require_zero(
        high_operator_value,
        "Weber closed order-hbar**2 period from unreduced D_2",
    )

    return {
        "closed_period": closed_period,
        "quantum_period_two": high_operator_value,
        "primitive": weber_primitive,
    }


def main() -> None:
    """Run every exact audit and print a concise reproducibility ledger."""
    print("Advanced Linear ODEs -- Picard--Fuchs operator audit")
    print(f"Python: {platform.python_version()}")
    print(f"SymPy: {sp.__version__}")
    print()

    ladder = energy_derivative_ladder_audit()
    print("[PASS] Energy-derivative ladder through n=6")
    print(f"  partial_E^3 lambda_0 coefficient = {ladder[2]}")

    legendre = legendre_picard_fuchs_audit()
    print("[PASS] Legendre Picard--Fuchs exact-form certificate")
    print(f"  discriminant = {legendre['discriminant']}")

    universal = universal_lambda_two_audit()
    print("[PASS] Universal source-free lambda_2 reductions")
    print(f"  lambda_2 coefficient = {universal['lambda_two']}")

    quartic = anharmonic_quartic_audit()
    print("[PASS] Anharmonic quartic discriminant and exact certificates")
    print(f"  discriminant = {quartic['discriminant']}")
    print(
        "  D_2 = "
        f"({quartic['operator_two']}) partial_E^2 + "
        f"({quartic['operator_three']}) partial_E^3"
    )

    reduced_operator_audit(quartic)
    print("[PASS] Quartic reduced operator and left-ideal equivalence")
    print(
        "  D_2^red = -(g/2) partial_E "
        "+ g(1+12gE)/(8E(1+4gE))"
    )

    weber = weber_limit_audit(quartic)
    print("[PASS] Closed Weber limit at g=0")
    print(f"  Pi_0 = {weber['closed_period']}")
    print(f"  Pi_2 = {weber['quantum_period_two']}")

    print()
    print("All Picard--Fuchs and quantum-operator checks passed.")
    print("Scope: formal closed-period identities; no analytic summation tested.")


if __name__ == "__main__":
    main()