#!/usr/bin/env python3
"""Audit the Voros-symbol and wall-crossing formulas for Chapter 9, Page 5.

The book uses the small parameter ``hbar`` and lateral Borel contours with
``+`` above the oriented ray and ``-`` below it.  For an oriented saddle
class ``gamma_0`` whose cycle symbol is exponentially small, the ordinary
Delabaere--Dillinger--Pham map is

    X_beta  -> X_beta (1 + Y_gamma0)**(-<gamma_0, beta>),
    Y_gamma -> Y_gamma (1 + Y_gamma0)**(-(gamma_0, gamma)).

The operational convention is

    S_(theta+) = S_(theta-) o Stokes_(gamma_0).

This script checks the convention-sensitive algebra behind that statement:

1. cycle-character addition and inversion, followed by a rank-two DDP map,
   active-symbol invariance, multiplicativity, and inverse;
2. preservation of the log-canonical Poisson bracket;
3. the pentagon identity, with the rightmost operator acting first;
4. the distinction between an analytic Stokes map and a geometric basis flip;
5. diagonal conjugation of Page 4's local triangular connection matrices;
6. the type-II simple-pole wall polynomial and its Frobenius parameter;
7. invariance of all cycle symbols under a degenerate-saddle pop; and
8. the Weber Bernoulli series, shifted-Borel germ and pole lattice, the
   finite-grade Log(1+q) jump, its log-Gamma asymptotics, and saddle orientation.

Every test raises an explicit exception.  Python's ``-O`` flag therefore does
not disable any check.  Exact identities use SymPy; the one asymptotic
log-Gamma calibration uses fixed, 80-decimal working-precision mpmath
arithmetic.  The runtime target is Python 3.9+, SymPy 1.14, and mpmath 1.3.

These calculations verify algebraic consequences of the cited exact-WKB
theorems.  They do not prove Borel summability, admissibility of paths,
existence of a saddle trajectory, or applicability of a wall-crossing theorem
to a particular differential equation.
"""

from __future__ import annotations

import platform
import sys
from collections.abc import Iterable, Mapping

import mpmath as mp
import sympy as sp


LatticeVector = tuple[int, int]
BirationalMap = dict[sp.Symbol, sp.Expr]


def reduced(expression: sp.Expr) -> sp.Expr:
    """Return a stable exact form of a rational or trigonometric expression."""

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


def intersection(left: LatticeVector, right: LatticeVector) -> int:
    """Return the rank-two signed pairing with (e_1, e_2) = +1."""

    return left[0] * right[1] - left[1] * right[0]


def character(
    charge: LatticeVector,
    first_symbol: sp.Symbol,
    second_symbol: sp.Symbol,
) -> sp.Expr:
    """Return the untwisted torus character associated with ``charge``."""

    return first_symbol ** charge[0] * second_symbol ** charge[1]


def apply_map(expression: sp.Expr, mapping: Mapping[sp.Symbol, sp.Expr]) -> sp.Expr:
    """Apply one algebra automorphism by simultaneous symbol replacement."""

    return reduced(expression.xreplace(dict(mapping)))


def apply_operator_word(
    expression: sp.Expr,
    factors_left_to_right: Iterable[Mapping[sp.Symbol, sp.Expr]],
) -> sp.Expr:
    """Apply a displayed operator product, with its rightmost factor first."""

    value = expression
    factors = list(factors_left_to_right)
    for mapping in reversed(factors):
        value = apply_map(value, mapping)
    return reduced(value)


def cycle_ddp_map(
    active: LatticeVector,
    first_symbol: sp.Symbol,
    second_symbol: sp.Symbol,
    *,
    inverse: bool = False,
) -> BirationalMap:
    """Return the ordinary DDP map on a rank-two cycle torus."""

    active_symbol = character(active, first_symbol, second_symbol)
    wall_factor = 1 + active_symbol
    sign = 1 if inverse else -1
    basis = ((1, 0), (0, 1))
    symbols = (first_symbol, second_symbol)
    return {
        symbol: reduced(
            symbol * wall_factor ** (sign * intersection(active, basis_vector))
        )
        for symbol, basis_vector in zip(symbols, basis)
    }


def log_canonical_bracket(
    first: sp.Expr,
    second: sp.Expr,
    first_symbol: sp.Symbol,
    second_symbol: sp.Symbol,
) -> sp.Expr:
    """Return the bracket determined by {Y_1, Y_2} = Y_1 Y_2."""

    return reduced(
        first_symbol
        * second_symbol
        * (
            sp.diff(first, first_symbol) * sp.diff(second, second_symbol)
            - sp.diff(first, second_symbol) * sp.diff(second, first_symbol)
        )
    )


def rank_two_ddp_audit() -> None:
    """Check the fixed-homology DDP map and its elementary invariants."""

    x_one, x_two, y_one, y_two = sp.symbols(
        "X_1 X_2 Y_1 Y_2",
        nonzero=True,
    )
    wall_factor = 1 + y_one

    first_charge = (2, -3)
    second_charge = (-5, 7)
    summed_charge = (
        first_charge[0] + second_charge[0],
        first_charge[1] + second_charge[1],
    )
    reversed_charge = (-first_charge[0], -first_charge[1])
    first_character = character(first_charge, y_one, y_two)
    second_character = character(second_charge, y_one, y_two)
    require_zero(
        character(summed_charge, y_one, y_two)
        - first_character * second_character,
        "addition law for cycle characters",
    )
    require_zero(
        character(reversed_charge, y_one, y_two) - first_character**-1,
        "orientation reversal for cycle characters",
    )

    forward: BirationalMap = {
        x_one: x_one / wall_factor,
        x_two: x_two,
        y_one: y_one,
        y_two: y_two / wall_factor,
    }
    backward: BirationalMap = {
        x_one: x_one * wall_factor,
        x_two: x_two,
        y_one: y_one,
        y_two: y_two * wall_factor,
    }

    require_zero(
        apply_map(y_one, forward) - y_one,
        "active cycle symbol under the DDP map",
    )
    require_zero(
        apply_map(wall_factor, forward) - wall_factor,
        "active wall factor under the DDP map",
    )

    cycle_charge = (2, -3)
    cycle_monomial = character(cycle_charge, y_one, y_two)
    expected_cycle_image = (
        cycle_monomial
        * wall_factor ** (-intersection((1, 0), cycle_charge))
    )
    require_zero(
        apply_map(cycle_monomial, forward) - expected_cycle_image,
        "rank-two DDP character formula",
    )

    path_monomial = x_one**2 / x_two**3
    require_zero(
        apply_map(path_monomial, forward)
        - path_monomial / wall_factor**2,
        "rank-two relative-path pairing",
    )

    first_factor = x_one * y_two**2
    second_factor = x_two**-1 * y_one**3
    require_zero(
        apply_map(first_factor * second_factor, forward)
        - apply_map(first_factor, forward) * apply_map(second_factor, forward),
        "multiplicativity of the DDP map",
    )

    for generator in (x_one, x_two, y_one, y_two):
        require_zero(
            apply_operator_word(generator, [backward, forward]) - generator,
            f"DDP inverse after forward map on {generator}",
        )
        require_zero(
            apply_operator_word(generator, [forward, backward]) - generator,
            f"DDP forward after inverse map on {generator}",
        )


def poisson_audit() -> None:
    """Check that ordinary DDP transformations preserve the Poisson torus."""

    y_one, y_two = sp.symbols("Y_1 Y_2", nonzero=True)
    for active in ((1, 0), (0, 1), (1, 1), (2, -1)):
        mapping = cycle_ddp_map(active, y_one, y_two)
        transformed_one = mapping[y_one]
        transformed_two = mapping[y_two]
        require_zero(
            log_canonical_bracket(
                transformed_one,
                transformed_two,
                y_one,
                y_two,
            )
            - transformed_one * transformed_two,
            f"log-canonical Poisson bracket for active charge {active}",
        )


def pentagon_audit() -> None:
    """Verify the rank-two pentagon with the rightmost map acting first."""

    y_one, y_two = sp.symbols("Y_1 Y_2", nonzero=True)
    stokes_one = cycle_ddp_map((1, 0), y_one, y_two)
    stokes_two = cycle_ddp_map((0, 1), y_one, y_two)
    stokes_sum = cycle_ddp_map((1, 1), y_one, y_two)

    # S_1 o S_(1+2) o S_2 = S_2 o S_1.  The helper deliberately
    # evaluates the rightmost operator first, as ordinary composition does.
    for generator in (y_one, y_two):
        left = apply_operator_word(
            generator,
            [stokes_one, stokes_sum, stokes_two],
        )
        right = apply_operator_word(generator, [stokes_two, stokes_one])
        require_zero(left - right, f"pentagon identity on {generator}")

    # A leftmost-first reading is genuinely different and must not pass by
    # accidental commutativity.
    wrong_value = y_one
    for mapping in (stokes_one, stokes_sum, stokes_two):
        wrong_value = apply_map(wrong_value, mapping)
    correct_value = apply_operator_word(
        y_one,
        [stokes_one, stokes_sum, stokes_two],
    )
    require_nonzero(
        wrong_value - correct_value,
        "pentagon order-sensitivity check",
    )


def geometric_flip_audit() -> None:
    """Separate the geometric basis mutation from the analytic DDP map."""

    x_one, x_two, y_one, y_two = sp.symbols(
        "X_1 X_2 Y_1 Y_2",
        nonzero=True,
    )
    adjacency = sp.Matrix([[0, 1], [-1, 0]])

    # Columns give the primed basis vectors in the unprimed basis:
    # gamma'_1=-gamma_1, gamma'_2=gamma_1+gamma_2;
    # beta'_1=-beta_1+beta_2, beta'_2=beta_2.
    cycle_change = sp.Matrix([[-1, 1], [0, 1]])
    path_change = sp.Matrix([[-1, 0], [1, 1]])
    require_matrix_zero(
        cycle_change.T * path_change - sp.eye(2),
        "duality after the geometric flip",
    )
    require_matrix_zero(
        cycle_change.T * adjacency * cycle_change + adjacency,
        "intersection matrix after the geometric flip",
    )

    analytic_map: BirationalMap = {
        x_one: x_one / (1 + y_one),
        x_two: x_two,
        y_one: y_one,
        y_two: y_two / (1 + y_one),
    }
    geometric_cycles = (y_one**-1, y_one * y_two)
    geometric_paths = (x_one**-1 * x_two, x_two)

    require_nonzero(
        analytic_map[y_one] - geometric_cycles[0],
        "analytic map versus geometric active-cycle inversion",
    )
    require_nonzero(
        analytic_map[y_two] - geometric_cycles[1],
        "analytic map versus geometric neighboring-cycle mutation",
    )

    combined_cycles = tuple(
        apply_map(symbol, analytic_map) for symbol in geometric_cycles
    )
    combined_paths = tuple(
        apply_map(symbol, analytic_map) for symbol in geometric_paths
    )
    require_zero(
        combined_cycles[0] - y_one**-1,
        "combined active cycle mutation",
    )
    require_zero(
        combined_cycles[1] - y_one * y_two / (1 + y_one),
        "combined neighboring cycle mutation",
    )
    require_zero(
        combined_paths[0] - x_two * (1 + y_one) / x_one,
        "combined active path mutation",
    )
    require_zero(
        combined_paths[1] - x_two,
        "combined spectator path mutation",
    )


def diagonal_shear_audit() -> None:
    """Check Page 4's diagonal conjugation of lower and upper shears."""

    d_plus, d_minus = sp.symbols("d_plus d_minus", nonzero=True)
    multiplier = sp.symbols("m")
    diagonal = sp.diag(d_plus, d_minus)
    lower = sp.Matrix([[1, 0], [multiplier, 1]])
    upper = sp.Matrix([[1, multiplier], [0, 1]])

    expected_lower = sp.Matrix(
        [[1, 0], [multiplier * d_plus / d_minus, 1]]
    )
    expected_upper = sp.Matrix(
        [[1, multiplier * d_minus / d_plus], [0, 1]]
    )
    require_matrix_zero(
        diagonal.inv() * lower * diagonal - expected_lower,
        "diagonal conjugation of a lower shear",
    )
    require_matrix_zero(
        diagonal.inv() * upper * diagonal - expected_upper,
        "diagonal conjugation of an upper shear",
    )
    require_zero(expected_lower.det() - 1, "lower conjugated determinant")
    require_zero(expected_upper.det() - 1, "upper conjugated determinant")

    transport = sp.symbols("transport", nonzero=True)
    voros_diagonal = sp.diag(transport, 1 / transport)
    airy_lower = sp.Matrix([[1, 0], [sp.I, 1]])
    airy_upper = sp.Matrix([[1, sp.I], [0, 1]])
    require_zero(
        (voros_diagonal.inv() * airy_lower * voros_diagonal)[1, 0]
        - sp.I * transport**2,
        "Voros transport dressing of the lower Airy multiplier",
    )
    require_zero(
        (voros_diagonal.inv() * airy_upper * voros_diagonal)[0, 1]
        - sp.I / transport**2,
        "Voros transport dressing of the upper Airy multiplier",
    )


def type_ii_simple_pole_audit() -> None:
    """Check the simple-pole polynomial, Frobenius data, and wall map."""

    nu = sp.symbols("nu", real=True)
    y_one, y_two, x_one = sp.symbols("Y_1 Y_2 X_1", nonzero=True)
    t_value = sp.exp(sp.pi * sp.I * nu)
    kappa = 2 * sp.cos(sp.pi * nu)
    local_multiplier = sp.I * kappa
    polynomial = 1 + kappa * y_one + y_one**2

    require_zero(
        t_value + t_value**-1 - kappa,
        "simple-pole trace parameter",
    )
    require_zero(
        polynomial - (1 + t_value * y_one) * (1 + y_one / t_value),
        "factorization of the type-II wall polynomial",
    )
    require_zero(
        polynomial.subs(nu, -nu) - polynomial,
        "simple-pole branch invariance",
    )
    require_zero(
        polynomial - (1 - sp.I * local_multiplier * y_one + y_one**2),
        "local multiplier versus type-II polynomial coefficient",
    )

    b_value = (nu**2 - 1) / 4
    rho_plus = (1 + nu) / 2
    rho_minus = (1 - nu) / 2
    require_zero(
        rho_plus * (rho_plus - 1) - b_value,
        "positive Frobenius indicial root",
    )
    require_zero(
        rho_minus * (rho_minus - 1) - b_value,
        "negative Frobenius indicial root",
    )

    eigenvalue_plus = sp.exp(2 * sp.pi * sp.I * rho_plus)
    eigenvalue_minus = sp.exp(2 * sp.pi * sp.I * rho_minus)
    require_zero(
        eigenvalue_plus + t_value,
        "positive Frobenius monodromy eigenvalue",
    )
    require_zero(
        eigenvalue_minus + t_value**-1,
        "negative Frobenius monodromy eigenvalue",
    )
    monodromy_trace = -t_value - t_value**-1
    require_zero(
        monodromy_trace + kappa,
        "Frobenius monodromy trace",
    )
    require_zero(
        -sp.I * monodromy_trace - local_multiplier,
        "Koike multiplier from Frobenius trace",
    )

    type_ii_forward: BirationalMap = {
        x_one: x_one / polynomial,
        y_one: y_one,
        y_two: y_two / polynomial,
    }
    type_ii_backward: BirationalMap = {
        x_one: x_one * polynomial,
        y_one: y_one,
        y_two: y_two * polynomial,
    }
    require_zero(
        apply_map(y_one, type_ii_forward) - y_one,
        "active symbol under the type-II map",
    )
    for generator in (x_one, y_one, y_two):
        require_zero(
            apply_operator_word(
                generator,
                [type_ii_backward, type_ii_forward],
            )
            - generator,
            f"type-II inverse on {generator}",
        )
    require_zero(
        log_canonical_bracket(
            type_ii_forward[y_one],
            type_ii_forward[y_two],
            y_one,
            y_two,
        )
        - type_ii_forward[y_one] * type_ii_forward[y_two],
        "type-II log-canonical Poisson bracket",
    )

    for half_integer in (
        sp.Rational(1, 2),
        sp.Rational(3, 2),
        sp.Rational(5, 2),
    ):
        require_zero(
            local_multiplier.subs(nu, half_integer),
            f"half-integer local multiplier at nu={half_integer}",
        )
        require_zero(
            polynomial.subs(nu, half_integer) - (1 + y_one**2),
            f"half-integer type-II polynomial at nu={half_integer}",
        )

    require_nonzero(
        (1 + y_one**2).subs(y_one, sp.Rational(1, 3)) - 1,
        "nontrivial half-integer wall despite a zero local shear",
    )


def pop_audit() -> None:
    """Check that a degenerate-saddle pop leaves every cycle symbol fixed."""

    x_one, x_two, y_zero, y_one, y_two = sp.symbols(
        "X_1 X_2 Y_0 Y_1 Y_2",
        nonzero=True,
    )
    pop_factor = 1 - y_zero
    forward: BirationalMap = {
        x_one: x_one * pop_factor,
        x_two: x_two / pop_factor**2,
        y_zero: y_zero,
        y_one: y_one,
        y_two: y_two,
    }
    backward: BirationalMap = {
        x_one: x_one / pop_factor,
        x_two: x_two * pop_factor**2,
        y_zero: y_zero,
        y_one: y_one,
        y_two: y_two,
    }

    for cycle_symbol in (y_zero, y_one, y_two, y_one**2 / y_two):
        require_zero(
            apply_map(cycle_symbol, forward) - cycle_symbol,
            f"cycle invariance under a pop for {cycle_symbol}",
        )
    require_zero(
        apply_map(x_one * x_two, forward)
        - apply_map(x_one, forward) * apply_map(x_two, forward),
        "multiplicativity of the pop map",
    )
    for generator in (x_one, x_two, y_zero, y_one, y_two):
        require_zero(
            apply_operator_word(generator, [backward, forward]) - generator,
            f"pop inverse on {generator}",
        )


def mp_rational(value: sp.Rational) -> mp.mpf:
    """Convert an exact SymPy rational to an mpmath number."""

    rational = sp.Rational(value)
    return mp.mpf(int(rational.p)) / int(rational.q)


def weber_audit() -> None:
    """Audit the Weber series, Borel poles, jump, and saddle orientation."""

    hbar, mu, scale = sp.symbols("hbar mu a", positive=True)
    xi, q_value = sp.symbols("xi q")
    expected_mu = (
        -hbar / (24 * mu),
        sp.Rational(7, 2880) * hbar**3 / mu**3,
        -sp.Rational(31, 40320) * hbar**5 / mu**5,
    )
    expected_scale = (
        -hbar / (12 * scale**2),
        sp.Rational(7, 360) * hbar**3 / scale**6,
        -sp.Rational(31, 1260) * hbar**5 / scale**10,
    )

    for k in range(1, 4):
        bernoulli_at_half = sp.bernoulli(2 * k, sp.Rational(1, 2))
        require_zero(
            bernoulli_at_half
            - (sp.Integer(2) ** (1 - 2 * k) - 1)
            * sp.bernoulli(2 * k),
            f"Bernoulli polynomial identity at k={k}",
        )
        term = reduced(
            bernoulli_at_half
            / ((2 * k - 1) * (2 * k))
            * (hbar / mu) ** (2 * k - 1)
        )
        require_zero(
            term - expected_mu[k - 1],
            f"Weber coefficient in mu-normalization at k={k}",
        )
        require_zero(
            term.subs(mu, scale**2 / 2) - expected_scale[k - 1],
            f"Weber coefficient after mu=a^2/2 at k={k}",
        )

    # In the shifted convention B[hbar**n] = xi**(n-1)/Gamma(n),
    # the Bernoulli series sums to this germ.  Its apparent pole at the
    # origin is removable.  The exponential form is the one displayed
    # on Page 5; the hyperbolic form makes the pole lattice transparent.
    borel_germ = (
        1
        / (2 * xi)
        * (
            1 / (sp.exp(xi / (2 * mu)) - 1)
            + 1 / (sp.exp(xi / (2 * mu)) + 1)
            - 2 * mu / xi
        )
    )
    borel_germ_hyperbolic = (
        1 / (2 * xi * sp.sinh(xi / (2 * mu))) - mu / xi**2
    )
    require_zero(
        borel_germ - borel_germ_hyperbolic.rewrite(sp.exp),
        "exponential and hyperbolic Weber Borel germs",
    )

    number_of_taylor_terms = 6
    borel_taylor = sp.series(
        borel_germ,
        xi,
        0,
        2 * number_of_taylor_terms,
    ).removeO()
    for k in range(1, number_of_taylor_terms + 1):
        expected_borel_coefficient = (
            sp.bernoulli(2 * k, sp.Rational(1, 2))
            / (sp.factorial(2 * k) * mu ** (2 * k - 1))
        )
        require_zero(
            sp.expand(borel_taylor).coeff(xi, 2 * k - 2)
            - expected_borel_coefficient,
            f"Weber shifted-Borel Taylor coefficient k={k}",
        )

    # For beta_W oriented so that <delta,beta_W>=-1, the positive
    # imaginary poles produce +Log(1+q).  Reversing beta_W changes the
    # pairing to +1 and the discontinuity to -Log(1+q), as in the DDP
    # example on Page 5.  A finite-grade identity does not by itself
    # prove convergence of the complete residue sum.
    maximum_action_grade = 6
    residues: dict[int, sp.Expr] = {}
    for k in range(1, maximum_action_grade + 1):
        singularity = 2 * sp.pi * sp.I * mu * k
        residue = reduced(
            sp.limit(
                (xi - singularity) * borel_germ,
                xi,
                singularity,
            )
        )
        expected_residue = sp.Rational((-1) ** k, k) / (2 * sp.pi * sp.I)
        require_zero(
            residue - expected_residue,
            f"Weber Borel residue at positive lattice pole k={k}",
        )
        residues[k] = residue

    jump_from_residues = reduced(
        sum(
            -2 * sp.pi * sp.I * residues[k] * q_value**k
            for k in range(1, maximum_action_grade + 1)
        )
    )
    logarithmic_jump = sp.series(
        sp.log(1 + q_value),
        q_value,
        0,
        maximum_action_grade + 1,
    ).removeO()
    require_zero(
        jump_from_residues - logarithmic_jump,
        "Weber residue sum versus +Log(1+q) through action grade 6",
    )
    require_zero(
        -jump_from_residues + logarithmic_jump,
        "reversed Weber path versus -Log(1+q) through action grade 6",
    )

    # The same coefficients are the large-t expansion of
    # log Gamma(t+1/2)-t log t+t-(1/2)log(2*pi), with t=mu/hbar.
    # This is an 80-decimal working-precision computation; the requested
    # comparison tolerance is about 1e-24, not a claim of 80 correct digits.
    mp.mp.dps = 80
    mu_value = mp.mpf("1.3")
    hbar_value = mp.mpf("0.04")
    t_value = mu_value / hbar_value
    number_of_terms = 8
    exact_remainder_function = (
        mp.loggamma(t_value + mp.mpf("0.5"))
        - t_value * mp.log(t_value)
        + t_value
        - mp.log(2 * mp.pi) / 2
    )
    partial_sum = mp.mpf(0)
    for k in range(1, number_of_terms + 1):
        coefficient = sp.bernoulli(2 * k, sp.Rational(1, 2))
        partial_sum += (
            mp_rational(coefficient)
            / (2 * k * (2 * k - 1))
            * (hbar_value / mu_value) ** (2 * k - 1)
        )
    next_index = number_of_terms + 1
    next_coefficient = sp.bernoulli(
        2 * next_index,
        sp.Rational(1, 2),
    )
    next_term = (
        mp_rational(next_coefficient)
        / (2 * next_index * (2 * next_index - 1))
        * (hbar_value / mu_value) ** (2 * next_index - 1)
    )
    asymptotic_error = abs(exact_remainder_function - partial_sum)
    require(
        asymptotic_error < mp.mpf("1.1") * abs(next_term),
        "Weber log-Gamma remainder is not controlled by the next term: "
        f"error={mp.nstr(asymptotic_error, 10)}, "
        f"next={mp.nstr(abs(next_term), 10)}",
    )
    require_close(
        partial_sum,
        exact_remainder_function,
        mp.mpf("1e-24"),
        "Weber log-Gamma calibration at 80-decimal working precision",
    )

    # Page 3 uses Z_delta=i*pi*a^2 at theta=pi/2.  The DDP active class
    # is -delta, so its phased action is negative and its symbol is small.
    z_delta = sp.I * sp.pi * scale**2
    z_active = -z_delta
    phased_action = sp.exp(-sp.I * sp.pi / 2) * z_active
    require_zero(
        phased_action + sp.pi * scale**2,
        "Weber active saddle orientation",
    )
    numerical_symbol = mp.exp(
        -mp.pi * mp.mpf("1.3") ** 2 / mp.mpf("0.2")
    )
    require(
        mp.mpf(0) < numerical_symbol < mp.mpf(1),
        "Weber active cycle symbol is not exponentially small",
    )


def main() -> int:
    """Run every Page-5 wall-crossing audit."""

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

    rank_two_ddp_audit()
    print(
        "PASS cycle-character laws and rank-two DDP active-symbol, product, "
        "and inverse checks"
    )

    poisson_audit()
    print("PASS log-canonical Poisson preservation")

    pentagon_audit()
    print("PASS pentagon identity with the rightmost map acting first")

    geometric_flip_audit()
    print("PASS geometric basis flip separated from the analytic Stokes map")

    diagonal_shear_audit()
    print("PASS Page 4 diagonal conjugation of local shears")

    type_ii_simple_pole_audit()
    print("PASS type-II polynomial, Frobenius trace, and half-integer case")

    pop_audit()
    print("PASS degenerate-saddle pop and cycle-symbol invariance")

    weber_audit()
    print(
        "PASS Weber Bernoulli series, Borel residues k=1..6, Log(1+q) "
        "through q^6, and saddle orientation"
    )
    print(
        "PASS Weber-Gamma calibration at mu=1.3, hbar=0.04 with 8 terms "
        "and 80-decimal working precision"
    )

    print(
        "LIMITATION: these identities do not prove Borel summability, "
        "convergence of the infinite pole-lattice sum, path admissibility, "
        "saddle existence, or theorem applicability for a given ODE."
    )
    print("All Voros wall-crossing 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)