#!/usr/bin/env python3
"""Reproduce the local connection checks for Chapter 9, Page 4.

The book uses

    [-hbar**2 d_z**2 + V(z)] psi(z) = E psi(z),
    hbar**2 psi'' = R psi,
    Wr[f, g] = f*g' - f'*g,

and the unit-normalized formal WKB pair

    psi_+/- = P_even**(-1/2)
              * exp(+/- integral(P_even dz) / hbar).

Consequently Wr[psi_+, psi_-] = -2/hbar.  Sectorial row frames act on
the right.  If region 2 lies counterclockwise from region 1, the
source-facing theorem is AC(Phi_1) = Phi_2 V_12.  Replacing the
continued source frame by the canonical destination frame uses the
inverse relation Phi_2 = AC(Phi_1) V_12**(-1).

This program verifies:

1. the invariant Airy scale at a simple zero and its exact linear model;
2. the unit WKB amplitude, source-normalization translation, and Wronskian;
3. exact rotated-Airy identities and the three canonical Airy frames;
4. triangular connection and frame-replacement matrices and three-edge closure;
5. conjugation of multipliers under a diagonal normalization change;
6. reduction of the simple-pole model to the modified Bessel equation;
7. exact Bessel Wronskians and continuation formulas giving Koike's
   multiplier 2*i*cos(pi*nu); and
8. the vanishing multiplier at half-integer nu.

Every check raises an explicit exception, so Python's -O flag does not
disable the audit.  The special-function identities calibrate local
connection formulas.  They do not prove Borel summability, path
admissibility, or any global connection or monodromy formula.
"""

from __future__ import annotations

import platform
import sys

import mpmath as mp
import sympy as sp


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

    if not condition:
        raise RuntimeError(message)


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

    residual = sp.simplify(sp.trigsimp(expression))
    if residual != 0:
        residual = sp.simplify(residual.rewrite(sp.sin))
    if residual != 0:
        raise RuntimeError(f"{message}: residual = {residual}")


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


def wronskian(
    first: tuple[complex, complex],
    second: tuple[complex, complex],
) -> complex:
    """Return Wr[first, second] from (value, derivative) pairs."""

    return first[0] * second[1] - first[1] * second[0]


def scaled_solution(
    solution: tuple[complex, complex],
    scalar: complex,
) -> tuple[complex, complex]:
    """Multiply a value-and-derivative pair by a constant."""

    return scalar * solution[0], scalar * solution[1]


def lower(multiplier: sp.Expr) -> sp.Matrix:
    """Lower-triangular unipotent right-action matrix."""

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


def upper(multiplier: sp.Expr) -> sp.Matrix:
    """Upper-triangular unipotent right-action matrix."""

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


def airy_scaling_and_formal_normalization_audit() -> None:
    """Check the Airy scale and the book/source normalization ledger."""

    local_x, c, hbar = sp.symbols(
        "local_x c hbar",
        positive=True,
        real=True,
    )
    zeta = c ** sp.Rational(1, 3) * local_x
    zeta_prime = sp.diff(zeta, local_x)
    require_zero(
        zeta * zeta_prime**2 - c * local_x,
        "simple-zero natural Airy coordinate",
    )

    airy_x = zeta / hbar ** sp.Rational(2, 3)
    airy_jacobian = sp.diff(airy_x, local_x)
    require_zero(
        hbar**2 * airy_jacobian**2 * airy_x - c * local_x,
        "linear turning-point reduction to the Airy equation",
    )

    airy_action = sp.Rational(2, 3) * airy_x ** sp.Rational(3, 2)
    wkb_action = (
        sp.Rational(2, 3)
        * sp.sqrt(c)
        * local_x ** sp.Rational(3, 2)
        / hbar
    )
    require_zero(
        airy_action - wkb_action,
        "Airy exponential versus WKB action",
    )

    airy_common = (
        2
        * sp.sqrt(sp.pi)
        * hbar ** -sp.Rational(1, 6)
        * c ** -sp.Rational(1, 6)
    )
    normalized_airy_amplitude = (
        airy_common
        * airy_x ** -sp.Rational(1, 4)
        / (2 * sp.sqrt(sp.pi))
    )
    require_zero(
        normalized_airy_amplitude
        - (c * local_x) ** -sp.Rational(1, 4),
        "Airy recessive unit WKB amplitude",
    )

    normalized_bi_amplitude = (
        airy_common
        * airy_x ** -sp.Rational(1, 4)
        / 2
        / sp.sqrt(sp.pi)
    )
    require_zero(
        normalized_bi_amplitude
        - (c * local_x) ** -sp.Rational(1, 4),
        "Airy dominant unit WKB amplitude",
    )

    momentum = sp.symbols("P", positive=True, real=True)
    source_prefactor = (momentum / hbar) ** -sp.Rational(1, 2)
    book_prefactor = momentum ** -sp.Rational(1, 2)
    require_zero(
        source_prefactor - sp.sqrt(hbar) * book_prefactor,
        "traditional odd-part versus book normalization",
    )

    amplitude_log_derivative = sp.symbols("amplitude_log_derivative")
    branch_product = 1 / momentum
    formal_wronskian = branch_product * (
        amplitude_log_derivative
        - momentum / hbar
        - amplitude_log_derivative
        - momentum / hbar
    )
    require_zero(
        formal_wronskian + 2 / hbar,
        "book unit-normalized formal Wronskian",
    )
    require_zero(
        hbar * formal_wronskian + 2,
        "traditional source-normalized formal Wronskian",
    )


def airy_solutions(
    airy_x: complex,
    hbar: mp.mpf,
    slope: mp.mpf,
) -> tuple[
    tuple[complex, complex],
    tuple[complex, complex],
    tuple[complex, complex],
]:
    """Return the three normalized rotated Airy solutions and z derivatives."""

    omega = mp.exp(2j * mp.pi / 3)
    common = (
        2
        * mp.sqrt(mp.pi)
        * hbar ** (-mp.mpf(1) / 6)
        * slope ** (-mp.mpf(1) / 6)
    )
    d_x_d_z = slope ** (mp.mpf(1) / 3) / hbar ** (mp.mpf(2) / 3)

    value_0 = common * mp.airyai(airy_x)
    derivative_0 = common * d_x_d_z * mp.airyai(airy_x, 1)

    value_1 = (
        common
        * mp.exp(1j * mp.pi / 6)
        * mp.airyai(omega * airy_x)
    )
    derivative_1 = (
        common
        * mp.exp(1j * mp.pi / 6)
        * omega
        * d_x_d_z
        * mp.airyai(omega * airy_x, 1)
    )

    value_2 = (
        common
        * mp.exp(-1j * mp.pi / 6)
        * mp.airyai(omega**2 * airy_x)
    )
    derivative_2 = (
        common
        * mp.exp(-1j * mp.pi / 6)
        * omega**2
        * d_x_d_z
        * mp.airyai(omega**2 * airy_x, 1)
    )

    return (
        (value_0, derivative_0),
        (value_1, derivative_1),
        (value_2, derivative_2),
    )


def airy_exact_audit() -> None:
    """Check rotated Airy identities, asymptotic constants, and Wronskians."""

    hbar = mp.mpf("0.73")
    slope = mp.mpf("1.41")
    tolerance = mp.mpf("1e-55")
    points = [
        mp.mpc("0.43", "0.27"),
        mp.mpc("-0.61", "0.38"),
        mp.mpc("1.17", "-0.29"),
    ]

    for point_index, airy_x in enumerate(points):
        airy_0, airy_1, airy_2 = airy_solutions(airy_x, hbar, slope)

        require_close(
            airy_1[0] - airy_2[0],
            1j * airy_0[0],
            tolerance,
            f"rotated Airy value identity at point {point_index}",
        )
        require_close(
            airy_1[1] - airy_2[1],
            1j * airy_0[1],
            tolerance,
            f"rotated Airy derivative identity at point {point_index}",
        )

        bi_scale = (
            mp.sqrt(mp.pi)
            * hbar ** (-mp.mpf(1) / 6)
            * slope ** (-mp.mpf(1) / 6)
        )
        d_x_d_z = (
            slope ** (mp.mpf(1) / 3)
            / hbar ** (mp.mpf(2) / 3)
        )
        bi_value = mp.airybi(airy_x)
        ai_value = mp.airyai(airy_x)
        bi_derivative = d_x_d_z * mp.airybi(airy_x, 1)
        ai_derivative = d_x_d_z * mp.airyai(airy_x, 1)
        plus_below = (
            bi_scale * (bi_value + 1j * ai_value),
            bi_scale * (bi_derivative + 1j * ai_derivative),
        )
        plus_above = (
            bi_scale * (bi_value - 1j * ai_value),
            bi_scale * (bi_derivative - 1j * ai_derivative),
        )
        for component in range(2):
            require_close(
                airy_1[component],
                plus_below[component],
                tolerance,
                f"lower-bank Airy/Bi identity at point {point_index}",
            )
            require_close(
                airy_2[component],
                plus_above[component],
                tolerance,
                f"upper-bank Airy/Bi identity at point {point_index}",
            )

        expected_wronskian = -2 / hbar
        frame_2_second = scaled_solution(airy_1, -1j)
        require_close(
            wronskian(airy_1, airy_0),
            expected_wronskian,
            tolerance,
            f"Airy frame 0 Wronskian at point {point_index}",
        )
        require_close(
            wronskian(airy_2, airy_0),
            expected_wronskian,
            tolerance,
            f"Airy frame 1 Wronskian at point {point_index}",
        )
        require_close(
            wronskian(airy_2, frame_2_second),
            expected_wronskian,
            tolerance,
            f"Airy frame 2 Wronskian at point {point_index}",
        )

    # A large positive point independently detects missing factors of two,
    # pi, or hbar in the unit WKB normalization.
    large_x = mp.mpf("25")
    airy_0, airy_1, airy_2 = airy_solutions(large_x, 1, 1)
    action = mp.mpf(2) * large_x ** (mp.mpf(3) / 2) / 3
    recessive_wkb = large_x ** (-mp.mpf(1) / 4) * mp.exp(-action)
    dominant_wkb = large_x ** (-mp.mpf(1) / 4) * mp.exp(action)
    asymptotic_tolerance = mp.mpf("0.002")
    require_close(
        airy_0[0] / recessive_wkb,
        1,
        asymptotic_tolerance,
        "Airy recessive leading normalization",
    )
    require_close(
        airy_1[0] / dominant_wkb,
        1,
        asymptotic_tolerance,
        "Airy lower-bank dominant leading normalization",
    )
    require_close(
        airy_2[0] / dominant_wkb,
        1,
        asymptotic_tolerance,
        "Airy upper-bank dominant leading normalization",
    )


def connection_matrix_audit() -> None:
    """Check oriented triangular matrices, inverses, and Airy closure."""

    identity = sp.eye(2)
    symplectic_form = sp.Matrix([[0, 1], [-1, 0]])
    sheet_swap = sp.Matrix([[0, -sp.I], [-sp.I, 0]])

    stokes_0 = lower(-sp.I)
    stokes_1 = upper(-sp.I)
    stokes_2 = lower(-sp.I)
    matrices = [stokes_0, stokes_1, stokes_2]

    for index, matrix in enumerate(matrices):
        require_zero(matrix.det() - 1, f"Airy matrix {index} determinant")
        require_matrix_zero(
            matrix.T * symplectic_form * matrix - symplectic_form,
            f"Airy matrix {index} Wronskian preservation",
        )

    require_matrix_zero(
        stokes_0 * lower(sp.I) - identity,
        "positive-action connection/replacement inverse pair",
    )
    require_matrix_zero(
        stokes_1 * upper(sp.I) - identity,
        "negative-action connection/replacement inverse pair",
    )

    require_matrix_zero(
        stokes_0 * stokes_1 * stokes_2 - sheet_swap,
        "three-edge Airy product",
    )
    require_matrix_zero(
        stokes_0 * stokes_1 * stokes_2 * sheet_swap.inv() - identity,
        "Airy local analytic-monodromy closure",
    )
    require_matrix_zero(
        sheet_swap**2 + identity,
        "two-turn WKB half-density sign",
    )

    # Express every sector frame in the initial exact basis (A_1, A_0).
    frame_0 = identity
    frame_1 = sp.Matrix([[1, 0], [-sp.I, 1]])
    frame_2 = sp.Matrix([[1, -sp.I], [-sp.I, 0]])
    frame_3 = sheet_swap
    require_matrix_zero(
        frame_1 - frame_0 * stokes_0,
        "Airy frame 0-to-1 relation",
    )
    require_matrix_zero(
        frame_2 - frame_1 * stokes_1,
        "Airy frame 1-to-2 relation",
    )
    require_matrix_zero(
        frame_3 - frame_2 * stokes_2,
        "Airy frame 2-to-3 relation",
    )

    # Swapping the +/- order exchanges lower and upper triangular matrices.
    permutation = sp.Matrix([[0, 1], [1, 0]])
    multiplier = sp.symbols("multiplier")
    require_matrix_zero(
        permutation.inv()
        * lower(multiplier)
        * permutation
        - upper(multiplier),
        "branch reordering exchanges triangular type",
    )


def normalization_conjugation_audit() -> None:
    """Check how independent column rescalings change a multiplier."""

    multiplier = sp.symbols("multiplier")
    d_plus, d_minus = sp.symbols(
        "d_plus d_minus",
        nonzero=True,
    )
    diagonal = sp.diag(d_plus, d_minus)

    require_matrix_zero(
        diagonal.inv()
        * lower(multiplier)
        * diagonal
        - lower(multiplier * d_plus / d_minus),
        "lower multiplier under diagonal normalization",
    )
    require_matrix_zero(
        diagonal.inv()
        * upper(multiplier)
        * diagonal
        - upper(multiplier * d_minus / d_plus),
        "upper multiplier under diagonal normalization",
    )

    common_scale = sp.symbols("common_scale", nonzero=True)
    common_diagonal = common_scale * sp.eye(2)
    require_matrix_zero(
        common_diagonal.inv()
        * lower(multiplier)
        * common_diagonal
        - lower(multiplier),
        "common normalization leaves multiplier fixed",
    )

    hbar = sp.symbols("hbar", nonzero=True)
    original_wronskian = -2 / hbar
    require_zero(
        d_plus * d_minus * original_wronskian
        + 2 * d_plus * d_minus / hbar,
        "Wronskian under independent column rescaling",
    )


def simple_pole_symbolic_audit() -> None:
    """Reduce the canonical simple-pole equation to modified Bessel form."""

    x, residue, hbar = sp.symbols(
        "x residue hbar",
        positive=True,
        real=True,
    )
    nu = sp.symbols("nu", real=True)
    b = (nu**2 - 1) / 4
    bessel_x = 2 * sp.sqrt(residue * x) / hbar
    require_zero(
        sp.diff(bessel_x, x) - bessel_x / (2 * x),
        "simple-pole Bessel coordinate derivative",
    )
    require_zero(
        bessel_x**2 - 4 * residue * x / hbar**2,
        "simple-pole Bessel coordinate scale",
    )

    bessel_function = sp.Function("u")
    dummy = sp.symbols("s")
    composed_solution = sp.sqrt(x) * bessel_function(bessel_x)
    differentiated_solution = (
        4
        * x ** sp.Rational(3, 2)
        * sp.diff(composed_solution, x, 2)
    )
    first_composed_derivative = sp.Subs(
        sp.diff(bessel_function(dummy), dummy),
        dummy,
        bessel_x,
    )
    second_composed_derivative = sp.Subs(
        sp.diff(bessel_function(dummy), dummy, 2),
        dummy,
        bessel_x,
    )
    expected_second_derivative = (
        bessel_x**2 * second_composed_derivative
        + bessel_x * first_composed_derivative
        - bessel_function(bessel_x)
    )
    require_zero(
        differentiated_solution - expected_second_derivative,
        "simple-pole square-root gauge chain rule",
    )

    function_value, first_derivative, second_derivative = sp.symbols(
        "u u_s u_ss"
    )
    chain_rule_equation = (
        bessel_x**2 * second_derivative
        + bessel_x * first_derivative
        - function_value
        - (bessel_x**2 + 4 * b) * function_value
    )
    modified_bessel_equation = (
        bessel_x**2 * second_derivative
        + bessel_x * first_derivative
        - (bessel_x**2 + nu**2) * function_value
    )
    require_zero(
        chain_rule_equation - modified_bessel_equation,
        "simple-pole reduction to modified Bessel equation",
    )

    rho_plus = (1 + nu) / 2
    rho_minus = (1 - nu) / 2
    require_zero(
        rho_plus * (rho_plus - 1) - b,
        "simple-pole positive Frobenius exponent",
    )
    require_zero(
        rho_minus * (rho_minus - 1) - b,
        "simple-pole negative Frobenius exponent",
    )
    require_zero(
        rho_plus - rho_minus - nu,
        "simple-pole exponent difference",
    )

    expected_amplitude = (residue / x) ** -sp.Rational(1, 4)
    i_leading_amplitude = (
        2
        * sp.sqrt(sp.pi)
        / sp.sqrt(hbar)
        * sp.sqrt(x)
        / sp.sqrt(2 * sp.pi * bessel_x)
    )
    k_leading_amplitude = (
        2
        / sp.sqrt(sp.pi * hbar)
        * sp.sqrt(x)
        * sp.sqrt(sp.pi / (2 * bessel_x))
    )
    require_zero(
        i_leading_amplitude - expected_amplitude,
        "simple-pole growing unit WKB amplitude",
    )
    require_zero(
        k_leading_amplitude - expected_amplitude,
        "simple-pole recessive unit WKB amplitude",
    )

    # Wr_s[I_nu, K_nu] = -1/s.  The square-root gauge and
    # ds/dx=s/(2x) turn the raw x-Wronskian into -1/2.
    raw_wronskian = x * (bessel_x / (2 * x)) * (-1 / bessel_x)
    normalization_product = 4 / hbar
    require_zero(
        raw_wronskian + sp.Rational(1, 2),
        "raw simple-pole Bessel Wronskian",
    )
    require_zero(
        normalization_product * raw_wronskian + 2 / hbar,
        "unit-normalized simple-pole Wronskian",
    )

    t = sp.exp(sp.pi * sp.I * nu)
    koike_multiplier = 2 * sp.I * sp.cos(sp.pi * nu)
    require_zero(
        sp.I * (t + 1 / t) - koike_multiplier,
        "Koike t-plus-inverse multiplier",
    )
    require_zero(
        lower(-koike_multiplier).det() - 1,
        "simple-pole positive-action matrix determinant",
    )
    require_zero(
        upper(-koike_multiplier).det() - 1,
        "simple-pole negative-action matrix determinant",
    )
    require_matrix_zero(
        lower(-koike_multiplier) * lower(koike_multiplier) - sp.eye(2),
        "simple-pole positive-action connection/replacement inverse pair",
    )
    require_matrix_zero(
        upper(-koike_multiplier) * upper(koike_multiplier) - sp.eye(2),
        "simple-pole negative-action connection/replacement inverse pair",
    )

    monodromy_plus = -sp.exp(sp.pi * sp.I * nu)
    monodromy_minus = -sp.exp(-sp.pi * sp.I * nu)
    monodromy_trace = monodromy_plus + monodromy_minus
    require_zero(
        monodromy_plus * monodromy_minus - 1,
        "simple-pole Frobenius monodromy determinant",
    )
    require_zero(
        monodromy_trace + 2 * sp.cos(sp.pi * nu),
        "simple-pole Frobenius monodromy trace",
    )
    require_zero(
        koike_multiplier + sp.I * monodromy_trace,
        "Koike multiplier versus local monodromy trace",
    )


def besseli_prime(order: complex, argument: complex) -> complex:
    """Derivative of I_order by its argument."""

    return (
        mp.besseli(order - 1, argument)
        + mp.besseli(order + 1, argument)
    ) / 2


def besselk_prime(order: complex, argument: complex) -> complex:
    """Derivative of K_order by its argument."""

    return -(
        mp.besselk(order - 1, argument)
        + mp.besselk(order + 1, argument)
    ) / 2


def square_root_bessel_solution(
    coefficient: complex,
    x: mp.mpf,
    bessel_x: mp.mpf,
    function_value: complex,
    function_derivative: complex,
) -> tuple[complex, complex]:
    """Return coefficient*sqrt(x)*F(s) and its x derivative."""

    value = coefficient * mp.sqrt(x) * function_value
    derivative = (
        coefficient
        / (2 * mp.sqrt(x))
        * (function_value + bessel_x * function_derivative)
    )
    return value, derivative


def simple_pole_solutions(
    nu: mp.mpf,
    x: mp.mpf,
    residue: mp.mpf,
    hbar: mp.mpf,
) -> tuple[
    tuple[complex, complex],
    tuple[complex, complex],
    tuple[complex, complex],
    tuple[complex, complex],
]:
    """Return I, recessive K, and the lower/upper lateral growing pairs."""

    bessel_x = 2 * mp.sqrt(residue * x) / hbar
    positive_rotation = mp.exp(1j * mp.pi)
    negative_rotation = mp.exp(-1j * mp.pi)

    i_value = mp.besseli(nu, bessel_x)
    i_derivative = besseli_prime(nu, bessel_x)
    k_value = mp.besselk(nu, bessel_x)
    k_derivative = besselk_prime(nu, bessel_x)

    i_solution = square_root_bessel_solution(
        2 * mp.sqrt(mp.pi) / mp.sqrt(hbar),
        x,
        bessel_x,
        i_value,
        i_derivative,
    )
    k_solution = square_root_bessel_solution(
        2 / mp.sqrt(mp.pi * hbar),
        x,
        bessel_x,
        k_value,
        k_derivative,
    )

    lower_argument = positive_rotation * bessel_x
    upper_argument = negative_rotation * bessel_x
    lower_value = mp.besselk(nu, lower_argument)
    upper_value = mp.besselk(nu, upper_argument)
    lower_derivative = (
        positive_rotation * besselk_prime(nu, lower_argument)
    )
    upper_derivative = (
        negative_rotation * besselk_prime(nu, upper_argument)
    )
    lower_solution = square_root_bessel_solution(
        2j / mp.sqrt(mp.pi * hbar),
        x,
        bessel_x,
        lower_value,
        lower_derivative,
    )
    upper_solution = square_root_bessel_solution(
        -2j / mp.sqrt(mp.pi * hbar),
        x,
        bessel_x,
        upper_value,
        upper_derivative,
    )

    return i_solution, k_solution, lower_solution, upper_solution


def simple_pole_bessel_audit() -> None:
    """Check Bessel continuation, Wronskians, and Koike's multiplier."""

    x = mp.mpf("0.82")
    residue = mp.mpf("1.31")
    hbar = mp.mpf("0.71")
    bessel_x = 2 * mp.sqrt(residue * x) / hbar
    tolerance = mp.mpf("1e-55")
    positive_rotation = mp.exp(1j * mp.pi)
    negative_rotation = mp.exp(-1j * mp.pi)

    for index, nu in enumerate(
        [mp.mpf("0.37"), mp.mpf("0.73"), mp.mpf("1.21")]
    ):
        i_value = mp.besseli(nu, bessel_x)
        k_value = mp.besselk(nu, bessel_x)
        continued_above = mp.besselk(nu, positive_rotation * bessel_x)
        continued_below = mp.besselk(nu, negative_rotation * bessel_x)
        expected_above = (
            mp.exp(-1j * mp.pi * nu) * k_value
            - mp.pi * 1j * i_value
        )
        expected_below = (
            mp.exp(1j * mp.pi * nu) * k_value
            + mp.pi * 1j * i_value
        )
        require_close(
            continued_above,
            expected_above,
            tolerance,
            f"upper K continuation for nu sample {index}",
        )
        require_close(
            continued_below,
            expected_below,
            tolerance,
            f"lower K continuation for nu sample {index}",
        )

        i_solution, k_solution, lower_solution, upper_solution = (
            simple_pole_solutions(nu, x, residue, hbar)
        )
        expected_wronskian = -2 / hbar
        require_close(
            wronskian(i_solution, k_solution),
            expected_wronskian,
            tolerance,
            f"I/K unit Wronskian for nu sample {index}",
        )
        require_close(
            wronskian(lower_solution, k_solution),
            expected_wronskian,
            tolerance,
            f"lower lateral Wronskian for nu sample {index}",
        )
        require_close(
            wronskian(upper_solution, k_solution),
            expected_wronskian,
            tolerance,
            f"upper lateral Wronskian for nu sample {index}",
        )

        multiplier = 2j * mp.cos(mp.pi * nu)
        for component in range(2):
            require_close(
                lower_solution[component] - upper_solution[component],
                multiplier * k_solution[component],
                tolerance,
                f"Koike multiplier for nu sample {index}, component {component}",
            )


def half_integer_multiplier_audit() -> None:
    """Check that the simple-pole multiplier can vanish geometrically."""

    x = mp.mpf("0.82")
    residue = mp.mpf("1.31")
    hbar = mp.mpf("0.71")
    tolerance = mp.mpf("1e-55")

    for index, nu in enumerate(
        [mp.mpf("0.5"), mp.mpf("1.5"), mp.mpf("2.5")]
    ):
        _, k_solution, lower_solution, upper_solution = simple_pole_solutions(
            nu,
            x,
            residue,
            hbar,
        )
        multiplier = 2j * mp.cos(mp.pi * nu)
        require_close(
            multiplier,
            0,
            tolerance,
            f"half-integer Koike multiplier {index}",
        )
        for component in range(2):
            require_close(
                lower_solution[component] - upper_solution[component],
                multiplier * k_solution[component],
                tolerance,
                f"half-integer lateral equality {index}, component {component}",
            )


def main() -> int:
    """Run every local connection audit."""

    mp.mp.dps = 80

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

    airy_scaling_and_formal_normalization_audit()
    print("PASS Airy scaling and formal WKB normalization")

    airy_exact_audit()
    print("PASS exact Airy identities, asymptotics, and Wronskians")

    connection_matrix_audit()
    print("PASS oriented triangular matrices, reversal, and closure")

    normalization_conjugation_audit()
    print("PASS normalization conjugation and multiplier covariance")

    simple_pole_symbolic_audit()
    print("PASS simple-pole Bessel reduction and exact normalization")

    simple_pole_bessel_audit()
    print("PASS Bessel continuation and Koike multiplier")

    half_integer_multiplier_audit()
    print("PASS half-integer simple-pole zero multiplier")

    print(
        "LIMITATION: these local identities do not prove Borel summability, "
        "path admissibility, or global connection data."
    )
    print("All canonical sectorial connection checks passed.")
    return 0


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