#!/usr/bin/env python3
"""Audit Chapter 9 Page 6 boundary determinants and exact quantization.

The page uses row frames of solutions and coefficient columns.  If

    F_R = F_L C,

then a right boundary column ``b_R`` is represented in the left frame by
``C b_R``.  The invariant spectral statement is that the left and right
projective lines coincide, equivalently

    det(b_L, C b_R) = 0.

This script checks the convention-sensitive finite-dimensional algebra and
the exactly solvable Weber calibration used on the page:

1. determinant--Wronskian conversion in a common frame;
2. the connection entry selected by each pair of coordinate boundary lines;
3. covariance under endpoint rescaling and independent frame changes;
4. source-facing route order and a two-shear ``1 + V`` factor;
5. reduction of the real-line oscillator to the parabolic-cylinder equation;
6. the exact parabolic-cylinder Wronskian and normalized reciprocal-gamma
   boundary determinant;
7. the reciprocal-gamma/Voros factorization, its physical roots, their
   simplicity, and cancellation of spurious negative-integer cycle roots;
8. the clean half-line Dirichlet/Neumann parity split;
9. the characteristic polynomial of an SL(2) Bloch monodromy matrix.

Every failure raises an explicit exception; no Python ``assert`` statements
are used, so ``python3`` and ``python3 -O`` run the same audit.  Exact checks
use SymPy.  Direct parabolic-cylinder evaluations and the meromorphic
factorization use fixed 80-decimal mpmath arithmetic and deterministic test
points.

These calculations do not prove Borel summability, endpoint admissibility,
self-adjointness of a proposed operator domain, or applicability of an
exact-WKB connection theorem to a new potential.
"""

from __future__ import annotations

import platform
import sys
from collections.abc import Iterable

import mpmath as mp
import sympy as sp


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

    simplified = sp.gammasimp(sp.trigsimp(sp.simplify(expression)))
    return sp.factor(sp.cancel(sp.together(simplified)))


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_nonzero(expression: sp.Expr, message: str) -> None:
    """Require an exact symbolic expression not to vanish."""

    residual = reduced(expression)
    if residual == 0:
        raise RuntimeError(f"{message}: unexpectedly vanished")


def require_matrix_zero(matrix: sp.Matrix, message: str) -> None:
    """Require every entry of a symbolic matrix to vanish."""

    for row in range(matrix.rows):
        for column in range(matrix.cols):
            require_zero(
                matrix[row, column],
                f"{message}, entry ({row}, {column})",
            )


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, 12)} exceeds "
            f"{mp.nstr(tolerance * scale, 12)}"
        )


def two_column_determinant(left: sp.Matrix, right: sp.Matrix) -> sp.Expr:
    """Return the determinant of two two-component coefficient columns."""

    return reduced(sp.Matrix.hstack(left, right).det())


def lower_shear(multiplier: sp.Expr) -> sp.Matrix:
    """Return a lower unitriangular connection matrix."""

    return sp.Matrix([[1, 0], [multiplier, 1]])


def upper_shear(multiplier: sp.Expr) -> sp.Matrix:
    """Return an upper unitriangular connection matrix."""

    return sp.Matrix([[1, multiplier], [0, 1]])


def boundary_determinant_audit() -> None:
    """Check Wronskian conversion and coordinate-line entry passports."""

    f_one, f_two, fp_one, fp_two = sp.symbols(
        "f_1 f_2 fp_1 fp_2"
    )
    c_one, c_two, d_one, d_two = sp.symbols("c_1 c_2 d_1 d_2")
    hbar = sp.symbols("hbar", positive=True)

    frame = sp.Matrix([[f_one, f_two], [fp_one, fp_two]])
    left = sp.Matrix([c_one, c_two])
    right = sp.Matrix([d_one, d_two])
    coefficient_determinant = two_column_determinant(left, right)
    solution_wronskian = two_column_determinant(frame * left, frame * right)

    require_zero(
        solution_wronskian - frame.det() * coefficient_determinant,
        "determinant--Wronskian relation in a common frame",
    )
    require_zero(
        (
            -hbar * solution_wronskian / 2
            - coefficient_determinant
        ).subs(frame.det(), -2 / hbar),
        "unit-WKB determinant--Wronskian normalization",
    )

    c_11, c_12, c_21, c_22 = sp.symbols("C_11 C_12 C_21 C_22")
    connection = sp.Matrix([[c_11, c_12], [c_21, c_22]])
    e_one = sp.Matrix([1, 0])
    e_two = sp.Matrix([0, 1])
    passports = (
        (e_one, e_one, c_21, "first/first boundary entry"),
        (e_one, e_two, c_22, "first/second boundary entry"),
        (e_two, e_one, -c_11, "second/first boundary entry"),
        (e_two, e_two, -c_12, "second/second boundary entry"),
    )
    for left_line, right_line, expected, label in passports:
        require_zero(
            two_column_determinant(left_line, connection * right_line)
            - expected,
            label,
        )


def covariance_audit() -> None:
    """Check projective rescaling, frame covariance, and route reversal."""

    entries = sp.symbols("c_11 c_12 c_21 c_22")
    left_entries = sp.symbols("l_1 l_2")
    right_entries = sp.symbols("r_1 r_2")
    connection = sp.Matrix(2, 2, entries)
    left = sp.Matrix(left_entries)
    right = sp.Matrix(right_entries)
    boundary = two_column_determinant(left, connection * right)

    alpha, beta = sp.symbols("alpha beta", nonzero=True)
    require_zero(
        two_column_determinant(alpha * left, connection * (beta * right))
        - alpha * beta * boundary,
        "projective endpoint rescaling",
    )

    p, q, r, s = sp.symbols("p q r s")
    u, v, w, t = sp.symbols("u v w t")
    left_change = sp.Matrix([[p, q], [r, s]])
    right_change = sp.Matrix([[u, v], [w, t]])
    transformed_connection = (
        left_change.inv() * connection * right_change
    )
    transformed_left = left_change.inv() * left
    transformed_right = right_change.inv() * right
    transformed_boundary = two_column_determinant(
        transformed_left,
        transformed_connection * transformed_right,
    )
    require_zero(
        transformed_boundary - boundary / left_change.det(),
        "independent endpoint-frame covariance",
    )

    reversed_boundary = two_column_determinant(
        right,
        connection.inv() * left,
    )
    require_zero(
        reversed_boundary + boundary / connection.det(),
        "route reversal of the boundary determinant",
    )


def route_and_shear_audit() -> None:
    """Check source-facing order and the finite-well two-shear factor."""

    first = lower_shear(sp.Integer(2))
    second = upper_shear(sp.Integer(3))
    initial = sp.Matrix([0, 1])
    sequential = second * (first * initial)
    ordered_product = second * first
    reversed_product = first * second

    require_matrix_zero(
        sequential - ordered_product * initial,
        "source-facing route product C_2 C_1",
    )
    require_zero(ordered_product.det() - 1, "ordered route determinant")
    require_zero(reversed_product.det() - 1, "reversed route determinant")
    require_nonzero(
        ordered_product - reversed_product,
        "noncommuting route-order regression",
    )
    require_nonzero(
        two_column_determinant(sp.Matrix([1, 0]), ordered_product * initial)
        - two_column_determinant(
            sp.Matrix([1, 0]),
            reversed_product * initial,
        ),
        "route order must change a boundary entry",
    )

    voros = sp.symbols("V")
    imaginary_unit = sp.I
    require_matrix_zero(
        lower_shear(imaginary_unit)
        * lower_shear(imaginary_unit * voros)
        - lower_shear(imaginary_unit * (1 + voros)),
        "lower two-shear 1+V factor",
    )
    require_matrix_zero(
        upper_shear(imaginary_unit)
        * upper_shear(imaginary_unit * voros)
        - upper_shear(imaginary_unit * (1 + voros)),
        "upper two-shear 1+V factor",
    )

    voros_half = sp.symbols("v_half", nonzero=True)
    diagonal_transport = sp.diag(voros_half, 1 / voros_half)
    well_connection = (
        lower_shear(-imaginary_unit)
        * diagonal_transport
        * lower_shear(-imaginary_unit)
    )
    well_connection_expected = sp.Matrix(
        [
            [voros_half, 0],
            [
                -imaginary_unit * (1 + voros_half**2) / voros_half,
                1 / voros_half,
            ],
        ]
    )
    require_matrix_zero(
        well_connection - well_connection_expected,
        "exact word L(-i) diag(v^1/2,v^-1/2) L(-i)",
    )
    require_zero(
        well_connection[1, 0]
        + imaginary_unit
        * (1 + voros_half**2)
        / voros_half,
        "displayed L(-i) D_v L(-i) boundary entry",
    )
    require_zero(
        well_connection.det() - 1,
        "two-turning-point connection determinant",
    )


def oscillator_scaling_audit() -> None:
    """Check the oscillator-to-Weber scaling and its classical cycle."""

    hbar, energy, coordinate = sp.symbols(
        "hbar a_squared X",
        positive=True,
    )
    nu = energy / (2 * hbar) - sp.Rational(1, 2)
    require_zero(
        nu + sp.Rational(1, 2) - energy / (2 * hbar),
        "oscillator parabolic-cylinder parameter",
    )

    original_coefficient = hbar * coordinate**2 / 2 - energy
    factored_coefficient = -2 * hbar * (
        nu + sp.Rational(1, 2) - coordinate**2 / 4
    )
    require_zero(
        original_coefficient - factored_coefficient,
        "oscillator scaling to the Weber equation",
    )

    angle = sp.symbols("angle", real=True)
    half_action = energy * sp.integrate(
        sp.cos(angle) ** 2,
        (angle, -sp.pi / 2, sp.pi / 2),
    )
    require_zero(
        half_action - sp.pi * energy / 2,
        "oscillator real half-cycle action",
    )
    cycle_period = 2 * sp.I * half_action
    require_zero(
        cycle_period - sp.I * sp.pi * energy,
        "oriented Weber cut-cycle period",
    )


def exact_parabolic_wronskian_audit() -> None:
    """Derive the Weber Wronskian exactly from its values at the origin."""

    nu = sp.symbols("nu")
    hbar = sp.symbols("hbar", positive=True)
    value_at_zero = (
        2 ** (nu / 2)
        * sp.sqrt(sp.pi)
        / sp.gamma((1 - nu) / 2)
    )
    derivative_at_zero = (
        -2 ** ((nu + 1) / 2)
        * sp.sqrt(sp.pi)
        / sp.gamma(-nu / 2)
    )

    standard_wronskian = -2 * value_at_zero * derivative_at_zero
    expected_standard = sp.sqrt(2 * sp.pi) / sp.gamma(-nu)
    require_zero(
        standard_wronskian - expected_standard,
        "Wr_X[D_nu(X),D_nu(-X)]",
    )

    left_to_right_wronskian = -standard_wronskian
    physical_wronskian = sp.sqrt(2 / hbar) * left_to_right_wronskian
    normalized_determinant = (
        -sp.sqrt(hbar / (4 * sp.pi)) * physical_wronskian
    )
    require_zero(
        normalized_determinant - 1 / sp.gamma(-nu),
        "normalized Weber boundary determinant",
    )


def parabolic_derivative(nu: complex, coordinate: complex) -> complex:
    """Evaluate D_nu'(coordinate) using a stable three-term identity."""

    return (
        coordinate * mp.pcfd(nu, coordinate) / 2
        - mp.pcfd(nu + 1, coordinate)
    )


def numerical_parabolic_wronskian_audit() -> None:
    """Check Wronskian constancy directly at deterministic regular points."""

    mp.mp.dps = 80
    tolerance = mp.mpf("1e-60")
    parameters: Iterable[complex] = (
        mp.mpf("-0.4"),
        mp.mpf("0.3"),
        mp.mpf("2.25"),
        mp.mpc("0.37", "0.21"),
    )
    coordinates: Iterable[complex] = (
        mp.mpf("-1.1"),
        mp.mpf("0.25"),
        mp.mpf("1.7"),
    )

    for nu in parameters:
        expected = mp.sqrt(2 * mp.pi) / mp.gamma(-nu)
        for coordinate in coordinates:
            right = mp.pcfd(nu, coordinate)
            right_derivative = parabolic_derivative(nu, coordinate)
            left = mp.pcfd(nu, -coordinate)
            left_derivative = -parabolic_derivative(nu, -coordinate)
            wronskian = right * left_derivative - right_derivative * left
            require_close(
                wronskian,
                expected,
                tolerance,
                "parabolic-cylinder Wronskian "
                f"at nu={mp.nstr(nu, 8)}, X={mp.nstr(coordinate, 8)}",
            )


def weber_factorization_and_spectrum_audit() -> None:
    """Check the Voros factorization, roots, slopes, and cancellations."""

    nu, epsilon = sp.symbols("nu epsilon")
    cycle_symbol = -sp.exp(2 * sp.pi * sp.I * nu)
    factorized = (
        sp.gamma(nu + 1)
        * sp.exp(-sp.pi * sp.I * nu)
        / (2 * sp.pi * sp.I)
        * (1 + cycle_symbol)
    )

    mp.mp.dps = 80
    tolerance = mp.mpf("1e-60")
    numeric_parameters: Iterable[complex] = (
        mp.mpf("-0.4"),
        mp.mpf("0.3"),
        mp.mpf("2.25"),
        mp.mpc("0.37", "0.21"),
    )
    for value in numeric_parameters:
        direct_value = 1 / mp.gamma(-value)
        voros_value = -mp.exp(2 * mp.pi * 1j * value)
        factorized_value = (
            mp.gamma(value + 1)
            * mp.exp(-mp.pi * 1j * value)
            / (2 * mp.pi * 1j)
            * (1 + voros_value)
        )
        require_close(
            factorized_value,
            direct_value,
            tolerance,
            f"reciprocal-gamma/Voros factorization at nu={mp.nstr(value, 8)}",
        )

    for level in range(7):
        require_zero(
            1 / sp.gamma(-level),
            f"Weber reciprocal-gamma root n={level}",
        )
        local_series = sp.series(
            1 / sp.gamma(-level - epsilon),
            epsilon,
            0,
            2,
        ).removeO()
        expected_slope = (-1) ** (level + 1) * sp.factorial(level)
        require_zero(
            sp.expand(local_series).coeff(epsilon, 1) - expected_slope,
            f"simple Weber root slope n={level}",
        )
        require_zero(
            (1 + cycle_symbol).subs(nu, level),
            f"physical Weber cycle root n={level}",
        )

    for index in range(6):
        negative_integer = -index - 1
        require_zero(
            (1 + cycle_symbol).subs(nu, negative_integer),
            f"spurious cycle root nu={negative_integer}",
        )
        direct_value = sp.Rational(1, sp.factorial(index))
        require_zero(
            (1 / sp.gamma(-nu)).subs(nu, negative_integer) - direct_value,
            f"finite determinant at nu={negative_integer}",
        )
        require_zero(
            sp.limit(factorized, nu, negative_integer) - direct_value,
            f"gamma-pole cancellation at nu={negative_integer}",
        )


def half_line_boundary_audit() -> None:
    """Check Dirichlet/Neumann spectra and Hermite parity on the half-line."""

    coordinate = sp.symbols("X", real=True)
    for index in range(7):
        even_level = 2 * index
        odd_level = 2 * index + 1

        dirichlet_value = (
            2 ** (sp.Rational(odd_level, 2))
            * sp.sqrt(sp.pi)
            / sp.gamma(sp.Rational(1 - odd_level, 2))
        )
        neumann_value = (
            -2 ** (sp.Rational(even_level + 1, 2))
            * sp.sqrt(sp.pi)
            / sp.gamma(-sp.Rational(even_level, 2))
        )
        require_zero(
            dirichlet_value,
            f"half-line Dirichlet odd level n={odd_level}",
        )
        require_zero(
            neumann_value,
            f"half-line Neumann even level n={even_level}",
        )

    for level in range(7):
        hermite_state = (
            2 ** (-sp.Rational(level, 2))
            * sp.exp(-coordinate**2 / 4)
            * sp.hermite(level, coordinate / sp.sqrt(2))
        )
        require_zero(
            sp.diff(hermite_state, coordinate, 2)
            + (
                level
                + sp.Rational(1, 2)
                - coordinate**2 / 4
            )
            * hermite_state,
            f"Hermite-Gaussian Weber equation n={level}",
        )
        require_zero(
            hermite_state.subs(coordinate, -coordinate)
            - (-1) ** level * hermite_state,
            f"Hermite-Gaussian parity n={level}",
        )


def bloch_audit() -> None:
    """Check the SL(2) Bloch characteristic and trace condition."""

    a, b, c, d, bloch = sp.symbols("a b c d mu_B")
    monodromy = sp.Matrix([[a, b], [c, d]])
    characteristic = reduced((monodromy - bloch * sp.eye(2)).det())
    require_zero(
        characteristic
        - (
            bloch**2
            - bloch * sp.trace(monodromy)
            + monodromy.det()
        ),
        "general two-by-two characteristic polynomial",
    )

    phase = sp.symbols("phase", real=True)
    exponential = sp.exp(sp.I * phase)
    require_zero(
        exponential + exponential**-1 - 2 * sp.cos(phase),
        "Bloch multiplier trace reduction",
    )


def main() -> int:
    """Run every Page-6 exact-quantization audit."""

    print("Advanced ODE Chapter 9 Page 6 exact-quantization audit")
    print(f"Python: {platform.python_version()}")
    print(f"SymPy: {sp.__version__}")
    print(f"mpmath: {mp.__version__}")

    boundary_determinant_audit()
    print("PASS determinant--Wronskian relation and boundary-entry passports")

    covariance_audit()
    print("PASS endpoint rescaling, frame covariance, and route reversal")

    route_and_shear_audit()
    print("PASS source-facing route order and two-shear 1+V factor")

    oscillator_scaling_audit()
    print("PASS oscillator-to-Weber scaling and oriented classical cycle")

    exact_parabolic_wronskian_audit()
    numerical_parabolic_wronskian_audit()
    print("PASS exact and 80-decimal parabolic-cylinder Wronskian checks")

    weber_factorization_and_spectrum_audit()
    print(
        "PASS Weber Gamma/Voros factorization, seven simple roots, "
        "and spurious-root cancellations"
    )

    half_line_boundary_audit()
    print("PASS half-line Dirichlet/Neumann parity split and Hermite states")

    bloch_audit()
    print("PASS SL(2) Bloch characteristic and trace condition")

    print(
        "LIMITATION: these identities do not prove Borel summability, "
        "endpoint admissibility, operator-domain self-adjointness, or "
        "exact-WKB theorem applicability for a new potential."
    )
    print("All exact-quantization checks passed.")
    return 0


if __name__ == "__main__":
    try:
        sys.exit(main())
    except Exception as exc:  # pragma: no cover - command-line failure path
        print(f"FAIL: {exc}", file=sys.stderr)
        sys.exit(1)