#!/usr/bin/env python3
"""Bounded WKB order comparison and independently derived exact checks.

The exact oscillator calculation derives the two local correction polynomials
without using the target Iyer--Will formulas. A separate RW continued fraction
compares global resonances at ell=2,10 only, fundamental positive-Re branch.
There is no live solver in a website audit, no Kerr spectrum, and no assertion
that a local WKB or finite continued-fraction residual certifies a continuum
resonance. Run with Python 3, SymPy 1.14.0, mpmath 1.3.0; no installation.
"""
import argparse
import hashlib
import json
from pathlib import Path

import mpmath as mp
import sympy as s


EXACT = {}
NEGATIVE = {}


def zero(name, value):
    value = s.factor(s.simplify(value))
    if value != 0:
        raise AssertionError((name, value))
    EXACT[name] = str(value)


def nonzero(name, value):
    value = s.factor(s.simplify(value))
    if value == 0:
        raise AssertionError((name, value))
    NEGATIVE[name] = {"residual": str(value), "rejected": True}


def add(target, key, value):
    target[key] = s.expand(target.get(key, 0) + value)


def oscillator_energy():
    """Formal recursion with every intermediate ladder path retained.

    Uses |m>=(a_dagger)^m|0>, so q|m>=(m|m-1>+|m+1>)/sqrt(2).
    Dictionary keys are shifts m-n. For integer n>=0, every path crossing
    below the vacuum automatically gets a zero lowering factor. No finite
    projection of q is multiplied to approximate a projected power of q.
    """
    n = s.symbols("n", integer=True, nonnegative=True)
    c3, c4, c5, c6 = s.symbols("c3 c4 c5 c6")
    coefficients = {1: c3, 2: c4, 3: c5, 4: c6}

    def q_action(state):
        result = {}
        for shift, value in state.items():
            add(result, shift - 1, (n + shift) * value / s.sqrt(2))
            add(result, shift + 1, value / s.sqrt(2))
        return result

    def interaction(order, state):
        for _ in range(order + 2):
            state = q_action(state)
        return {j: s.expand(coefficients[order] * v) for j, v in state.items()}

    states = [{0: s.Integer(1)}]
    energies = [2*n + 1]
    for order in range(1, 5):
        current = {}
        for j in range(1, order + 1):
            for shift, value in interaction(j, states[order-j]).items():
                add(current, shift, value)
        energies.append(s.factor(current.get(0, 0)))
        if order == 4:
            break
        state = {}
        for shift in current:
            if shift == 0:
                continue
            numerator = -current[shift]
            for j in range(1, order):
                numerator += energies[j] * states[order-j].get(shift, 0)
            state[shift] = s.expand(numerator/(2*shift))
        states.append(state)
        # Full vector recursion: not merely the eigenvalue projection.
        shifts = set(current) | set(state)
        for earlier in states:
            shifts |= set(earlier)
        for shift in shifts:
            residual = 2*shift*state.get(shift, 0) + current.get(shift, 0)
            for j in range(1, order+1):
                residual -= energies[j] * states[order-j].get(shift, 0)
            zero("oscillator_order_%d_shift_%+d" % (order, shift), residual)
    zero("oscillator_first_order_parity", energies[1])
    zero("oscillator_third_order_parity", energies[3])
    return n, (c3, c4, c5, c6), energies, states


def corrections(v2, v3, v4, v5, v6, nu):
    """Iyer 1987 II p.3633 Eqs.(1.2--1.3), expanded-E convention.

    V_j is the j-th TORTOISE derivative at the real maximum; v2<0.
    B2 has the units of V; B3 has the units of v2, so Omega3=B3/(-2v2)
    is dimensionless. These are not the implicit quantization Lambdas.
    """
    B2 = (v4/v2)*(s.Rational(1, 4)+nu**2)/8
    B2 -= (v3/v2)**2*(7+60*nu**2)/288
    B3 = 5*(v3/v2)**4*(77+188*nu**2)/6912
    B3 -= v3**2*v4*(51+100*nu**2)/(384*v2**3)
    B3 += (v4/v2)**2*(67+68*nu**2)/2304
    B3 += v3*v5*(19+28*nu**2)/(288*v2**2)
    B3 -= (v6/v2)*(5+4*nu**2)/288
    return B2, B3


def exact_checks():
    n, cs, oscillator, states = oscillator_energy()
    nu = n + s.Rational(1, 2)
    A = s.symbols("A", positive=True)
    v2 = -2*A
    v3, v4, v5, v6 = s.symbols("V3 V4 V5 V6", real=True)
    B2, B3 = corrections(v2, v3, v4, v5, v6, nu)
    rotation = {
        cs[0]: -(1+s.I)*v3/(6*s.sqrt(2)*A**s.Rational(5, 4)),
        cs[1]: -s.I*v4/(24*A**s.Rational(3, 2)),
        cs[2]: (1-s.I)*v5/(120*s.sqrt(2)*A**s.Rational(7, 4)),
        cs[3]: v6/(720*A**2),
    }
    independent_B2 = -s.I*s.sqrt(A)*oscillator[2].subs(rotation)
    independent_E3 = -s.I*s.sqrt(A)*oscillator[4].subs(rotation)
    zero("independent_oscillator_second_WKB_correction", independent_B2-B2)
    zero("independent_oscillator_third_WKB_correction",
         independent_E3 + s.I*nu*2*s.sqrt(A)*B3/(4*A))
    nonzero("wrong_second_correction_sign", independent_B2+B2)
    nonzero("wrong_third_correction_sign",
            independent_E3-s.I*nu*2*s.sqrt(A)*B3/(4*A))
    nonzero("wrong_cubic_term_coefficient", independent_B2-B2-v3**2/(288*v2**2))
    # P=|0><0| gives PqP=0, whereas Pq^4P=3P/4. Obtain the latter
    # from the independently derived first-order quartic coefficient.
    projected_q = s.zeros(1)
    full_q4_ground = s.diff(oscillator[2], cs[1]).subs(n, 0)
    nonzero("project_q_before_fourth_power", (projected_q**4)[0, 0]-full_q4_ground)

    # The full radial potential is differentiated BEFORE any peak restriction.
    R = s.symbols("R", positive=True)
    J = s.symbols("J", positive=True)
    f = 1-2/R
    potential = f*(J/R**2-6/R**3)
    derivatives = [potential]
    for order in range(1, 7):
        derivatives.append(s.factor(f*s.diff(derivatives[-1], R)))
    polynomial = 2*J*R**2-(6*J+18)*R+48
    discriminant = 9*J**2-42*J+81
    peak = (3*(J+3)+s.sqrt(discriminant))/(2*J)
    zero("radial_peak_polynomial", s.diff(potential, R)+polynomial/R**5)
    zero("exterior_peak_is_stationary", polynomial.subs(R, peak))
    zero("peak_polynomial_at_horizon", polynomial.subs(R, 2)-(12-4*J))
    zero("peak_polynomial_at_photon_orbit", polynomial.subs(R, 3)+6)
    zero("peak_polynomial_derivative", s.diff(polynomial, R).subs(R, peak)-2*s.sqrt(discriminant))
    # These modulo-P identities avoid using the stationary equation too early.
    curvature_identity = s.together(derivatives[2]+2*f**2*s.sqrt(discriminant)/R**5)
    zero("peak_tortoise_curvature", curvature_identity.subs(R, peak))
    jinv = s.symbols("jinv", positive=True)
    zero("large_J_peak_expansion",
         s.series(peak.subs(J, 1/jinv), jinv, 0, 3).removeO()-3-jinv-s.Rational(8, 3)*jinv**2)
    # At a stationary point D_*^3 V=f^3 V_rrr+3 f^2 f_r V_rr.
    premature_third = f*s.diff(f**2*s.diff(potential, R, 2), R)
    legitimate_third = derivatives[3]
    p2 = peak.subs(J, 6)
    nonzero("substitute_peak_condition_before_third_derivative",
            (legitimate_third-premature_third).subs({J: 6, R: p2}))
    nonzero("use_radial_instead_of_tortoise_curvature",
            (derivatives[2]-s.diff(potential, R, 2)).subs({J: 6, R: p2}))
    nonzero("finite_ell_peak_is_exactly_photon_orbit", p2-3)
    zero("ell2_peak_radical", p2-(9+s.sqrt(17))/4)

    # Dimension homogeneity, with b=GM and V_j=b^(-j-2) W_j.
    b = s.symbols("b", positive=True)
    scaled = [v2/b**4, v3/b**5, v4/b**6, v5/b**7, v6/b**8]
    scaled_B2, scaled_B3 = corrections(*scaled, nu)
    zero("B2_has_inverse_length_squared", scaled_B2-B2/b**2)
    zero("B3_has_inverse_length_to_fourth", scaled_B3-B3/b**4)
    zero("Omega3_is_dimensionless", scaled_B3/(-2*scaled[0])-B3/(-2*v2))
    zero("leading_curvature_scale_has_inverse_length_squared",
         s.sqrt(-2*scaled[0])-s.sqrt(-2*v2)/b**2)

    # Derive this computation's CF coefficients from the full nonconstant RW
    # equation, using a different representation from the WKB derivatives.
    z = s.symbols("z")
    w = s.symbols("w")
    j = s.symbols("j", integer=True, nonnegative=True)
    z_of_R = 1-2/R
    log_prefactor = s.I*w-2*s.I*w/(R-2)+4*s.I*w/R
    # Exterior R>2 fixes the real logarithms; the single exponential avoids
    # assuming generally invalid identities between complex power products.
    prefactor = s.exp(s.I*w*(R-2)-2*s.I*w*s.log(R-2)+4*s.I*w*s.log(R))
    zero("CF_outgoing_prefactor_log_derivative", s.diff(prefactor, R)/prefactor-log_prefactor)
    multiplication = R**3/(R-2)
    change = {R: 2/(1-z)}
    op_a = s.factor((multiplication*f**2*s.diff(z_of_R, R)**2).subs(change))
    op_b = s.factor((multiplication*(f**2*(s.diff(z_of_R, R, 2)+2*log_prefactor*s.diff(z_of_R, R))
                       +f*s.diff(f, R)*s.diff(z_of_R, R))).subs(change))
    op_c = s.factor((multiplication*(f**2*(s.diff(log_prefactor, R)+log_prefactor**2)
                       +f*s.diff(f, R)*log_prefactor+w*w-potential)).subs(change))
    alpha = s.expand(op_a).coeff(z, 1)*j*(j+1)+s.expand(op_b).coeff(z, 0)*(j+1)
    beta = s.expand(op_a).coeff(z, 2)*j*(j-1)+s.expand(op_b).coeff(z, 1)*j+s.expand(op_c).coeff(z, 0)
    gamma = s.expand(op_a).coeff(z, 3)*(j-1)*(j-2)+s.expand(op_b).coeff(z, 2)*(j-1)+s.expand(op_c).coeff(z, 1)
    zero("independent_CF_alpha_from_RW_ODE", alpha-(j+1)*(j+1-4*s.I*w))
    zero("independent_CF_beta_from_RW_ODE", beta-(-2*j*j+(16*s.I*w-2)*j+32*w*w+8*s.I*w+3-J))
    zero("independent_CF_gamma_from_RW_ODE", gamma-((j-4*s.I*w)**2-4))
    nonzero("CF_ell10_cannot_reuse_ell2_beta", beta.subs(J, 110)-beta.subs(J, 6))

    # Independent exact Pöschl--Teller family, away from V0=alpha^2/4.
    x, h = s.symbols("x h", real=True)
    height, width = s.symbols("V0 alpha", positive=True)
    pt = height/s.cosh(width*x)**2
    pt_derivatives = [s.diff(pt, x, k).subs(x, 0) for k in range(7)]
    expected_pt = [height, 0, -2*height*width**2, 0, 16*height*width**4, 0, -272*height*width**6]
    for k, (actual, expected) in enumerate(zip(pt_derivatives, expected_pt)):
        zero("Poschl_Teller_derivative_%d" % k, actual-expected)
    pt_B2, pt_B3 = corrections(*pt_derivatives[2:], nu)
    zero("Poschl_Teller_second_correction", pt_B2+width**2*(nu**2+s.Rational(1, 4)))
    zero("Poschl_Teller_third_factor", pt_B3/(-2*pt_derivatives[2])+width**2/(8*height))
    pt_exact_E = height-h*h*width**2*(nu**2+s.Rational(1, 4))-2*s.I*h*width*nu*s.sqrt(height-h*h*width**2/4)
    pt_wkb_E = height-2*s.I*h*width*nu*s.sqrt(height)+h*h*pt_B2
    pt_wkb_E -= 2*s.I*h**3*width*nu*s.sqrt(height)*pt_B3/(-2*pt_derivatives[2])
    zero("Poschl_Teller_exact_energy_series_through_h3",
         s.series(pt_exact_E-pt_wkb_E, h, 0, 4).removeO())
    nonzero("Poschl_Teller_first_WKB_is_exact",
            s.expand((s.sqrt(s.Rational(99, 100))-s.I/s.Integer(10))**2)-(1-s.I/s.Integer(5)))

    # A damped positive-frequency root and its conjugate-frequency partner.
    zero("damped_positive_frequency_square_root", (2-s.I)**2-(3-4*s.I))
    zero("negative_frequency_retarded_partner", (-2-s.I)**2-s.conjugate(3-4*s.I))
    nonzero("other_root_of_same_energy_is_damped", s.im(-2+s.I))
    nonzero("negative_frequency_partner_has_same_energy", (-2-s.I)**2-(3-4*s.I))

    # General static spherical, areal radius and coordinate time t.
    r = s.symbols("r", positive=True)
    F, H = s.Function("F", positive=True)(r), s.Function("H", positive=True)(r)
    U = F/r**2
    impact2 = s.symbols("impact2", positive=True)
    speed2 = F*H*(1-impact2*U)
    def at_orbit(expression):
        # Replace the first derivative as an atomic value only: the orbit
        # relation must not be differentiated to constrain F'' off the orbit.
        return expression.subs(impact2, r*r/F).xreplace({s.diff(F, r): 2*F/r})
    lyapunov2 = -F*H*s.diff(U, r, 2)/(2*U)
    zero("null_orbit_radial_speed_vanishes", at_orbit(speed2))
    zero("null_orbit_radial_acceleration_vanishes", at_orbit(s.diff(speed2, r)))
    zero("coordinate_time_Lyapunov_from_radial_geodesic",
         at_orbit(s.diff(speed2, r, 2)/2-lyapunov2))
    U0, Up, Upp, F0, H0, Fp, Hp = s.symbols("U0 Up Upp F0 H0 Fp Hp")
    dstar2_U = F0*H0*Upp+(Fp*H0+F0*Hp)*Up/2
    zero("Lyapunov_tortoise_radial_equivalence", (dstar2_U-F0*H0*Upp).subs(Up, 0))
    scale = s.symbols("C", positive=True)
    zero("time_rescaling_orbital_frequency_squared", (F/scale**2)/r**2-U/scale**2)
    zero("time_rescaling_Lyapunov_squared",
         -(F/scale**2)*H*s.diff(U/scale**2, r, 2)/(2*U/scale**2)-lyapunov2/scale**2)
    eta = s.symbols("eta", nonnegative=True)
    Ftoy = 1-2*b/r
    Htoy = Ftoy*(1+eta*(3*b/r)**2)
    Utoy = Ftoy/r**2
    lamtoy = -Ftoy*Htoy*s.diff(Utoy, r, 2)/(2*Utoy)
    zero("static_counterexample_circular_radius", s.diff(Utoy, r).subs(r, 3*b))
    zero("static_counterexample_Omega_squared", Utoy.subs(r, 3*b)-1/(27*b*b))
    zero("static_counterexample_Lyapunov_squared", lamtoy.subs(r, 3*b)-(1+eta)/(27*b*b))
    zero("Schwarzschild_Lyapunov_equals_Omega", (lamtoy-Utoy).subs({eta: 0, r: 3*b}))
    nonzero("all_static_spherical_Lyapunov_equals_Omega", (lamtoy-Utoy).subs({eta: 1, r: 3*b}))
    zero("Schwarzschild_eikonal_U_tortoise_curvature",
         (f*s.diff(f*s.diff(f/R**2, R), R)).subs(R, 3)+s.Rational(2, 729))

    return {
        "R": R, "J": J, "potential": potential, "derivatives": derivatives, "peak": peak,
        "oscillator": {
            "energies": list(map(str, oscillator)),
            "intermediate_shifts": [sorted(state) for state in states],
            "rotation_coefficients": {str(k): str(v) for k, v in rotation.items()},
        },
        "B2": B2, "B3": B3,
    }


def dec(value, digits=110):
    # Enough decimal places to round-trip every retained 100-dps value.
    # Printed places are NOT a claim of certified spectral accuracy.
    return mp.nstr(value, digits, strip_zeros=False)


def complex_json(value):
    return {"real": dec(mp.re(value)), "imag": dec(mp.im(value))}


def cf_coefficients(j, w, ell):
    return ((j+1)*(j+1-4*mp.j*w),
            -2*j*j+(16*mp.j*w-2)*j+32*w*w+8*mp.j*w+3-ell*(ell+1),
            (j-4*mp.j*w)**2-4)


def cf(w, ell, depth):
    """Zero terminal ratio a_(N+1)/a_N=0, evaluating rows N down to zero."""
    denominator = cf_coefficients(depth, w, ell)[1]
    for j in range(depth-1, -1, -1):
        alpha, beta, _ = cf_coefficients(j, w, ell)
        gamma_next = cf_coefficients(j+1, w, ell)[2]
        denominator = beta-alpha*gamma_next/denominator
    return denominator


def cf_ratio_independent(w, ell, depth):
    """Independent top-level ratio recurrence; same continuum ODE and tail.

    This checks indexing/normalization of the finite CF, not an independent
    physical spectrum. Ratios R_j=a_(j+1)/a_j use row j+1, not row j.
    """
    ratio = mp.mpc(0)
    for j in range(depth, 0, -1):
        alpha, beta, gamma = cf_coefficients(j, w, ell)
        ratio = -gamma/(beta+alpha*ratio)
    alpha0, beta0, _ = cf_coefficients(0, w, ell)
    return beta0+alpha0*ratio


def solve_cf(ell, depth, dps):
    seed_strings = {2: ("0.374", "-0.089"), 10: ("1.997", "-0.096")}
    with mp.workdps(dps):
        seed = mp.mpc(*seed_strings[ell])
        tolerance = mp.power(10, -(dps-8))
        root = mp.findroot(lambda w: cf(w, ell, depth),
                           (seed, seed+mp.mpc("0.001", "-0.001")),
                           tol=tolerance, maxsteps=80, verify=True)
        residual = abs(cf(root, ell, depth))
        real_window = {2: ("0.30", "0.45"), 10: ("1.90", "2.10")}[ell]
        if not (mp.mpf(real_window[0])<mp.re(root)<mp.mpf(real_window[1])
                and mp.mpf("-0.12")<mp.im(root)<mp.mpf("-0.05")
                and residual<mp.power(10, -(dps-10))):
            raise AssertionError(("CF branch/residual", ell, depth, dps, root, residual))
        return +root


def numerical_checks(symbolic):
    output = {"continued_fraction": [], "wkb": [], "derivative_comparison": [],
              "numeric_checks": {}, "negative_controls": {}}
    numeric_checks = output["numeric_checks"]
    negatives = output["negative_controls"]

    def below(name, residual, tolerance):
        residual, tolerance = mp.mpf(residual), mp.mpf(tolerance)
        if not (mp.isfinite(residual) and residual>=0 and residual<tolerance):
            raise AssertionError((name, residual, tolerance))
        numeric_checks[name] = {"residual": dec(residual), "tolerance": dec(tolerance)}

    def rejected(name, residual, minimum):
        residual, minimum = mp.mpf(residual), mp.mpf(minimum)
        if not (mp.isfinite(residual) and residual>minimum):
            raise AssertionError((name, residual, minimum))
        negatives[name] = {"residual": dec(residual), "minimum": dec(minimum), "rejected": True}

    roots = {}
    for ell in (2, 10):
        for dps in (50, 80):
            for depth in (512, 1024, 2048):
                root = solve_cf(ell, depth, dps)
                roots[(ell, depth, dps)] = root
                with mp.workdps(100):
                    saved = complex_json(root)
                    saved_root = mp.mpc(saved["real"], saved["imag"])
                    solved = abs(cf(saved_root, ell, depth))
                    deep = abs(cf(saved_root, ell, 4096))
                    top = cf_coefficients(0, saved_root, ell)
                    normalized_deep = deep/(1+abs(top[0])+abs(top[1]))
                    alternate = abs(cf(root, ell, depth)-cf_ratio_independent(root, ell, depth))
                    below("CF_ratio_agreement_ell%d_N%d_dps%d" % (ell, depth, dps), alternate, "1e-90")
                    below("CF_saved_root_residual_ell%d_N%d_dps%d" % (ell, depth, dps), solved,
                          mp.power(10, -(dps-10)))
                    output["continued_fraction"].append({
                        "ell": ell, "n": 0, "inversion": 0, "depth": depth, "dps": dps,
                        "terminal_ratio": "a_(N+1)/a_N=0", "root": saved,
                        "saved_root_residual_at_100dps": dec(solved),
                        "deep_check_depth": 4096, "deep_check_dps": 100,
                        "deep_normalized_residual": dec(normalized_deep),
                    })

    with mp.workdps(100):
        output["reference_checks"] = []
        for ell in (2, 10):
            ref = roots[(ell, 2048, 80)]
            precision_drift = abs(ref-roots[(ell, 2048, 50)])
            depth_drift = abs(ref-roots[(ell, 1024, 80)])
            deeper_ref = solve_cf(ell, 4096, 100)
            reference_drift = abs(ref-deeper_ref)
            below("CF_precision_drift_ell%d" % ell, precision_drift, "1e-45")
            below("CF_depth_1024_to_2048_ell%d" % ell, depth_drift, "1e-40")
            below("CF_reference_to_4096_ell%d" % ell, reference_drift, "1e-58")
            output["reference_checks"].append({"ell": ell,
                "reference_depth": 2048, "reference_dps": 80,
                "reference_root": complex_json(ref),
                "precision_50_to_80_drift_at_2048": dec(precision_drift),
                "depth_1024_to_2048_drift_at_80dps": dec(depth_drift),
                "depth_4096_100dps_root": complex_json(deeper_ref),
                "reference_to_deeper_root_drift": dec(reference_drift)})

    R, J = symbolic["R"], symbolic["J"]
    derivative_functions = [s.lambdify((R, J), expr, "mpmath") for expr in symbolic["derivatives"]]
    with mp.workdps(80):
        for ell in (2, 10):
            angular = ell*(ell+1)
            peak = (3*(angular+3)+mp.sqrt(9*angular**2-42*angular+81))/(2*angular)
            xpeak = peak+2*mp.log(peak/2-1)
            # Independent coordinate inversion, not repeated f*d/dR.
            def exact_R_of_x(x):
                return 2*(1+mp.lambertw(mp.exp(x/2-1)))
            def potential_of_x(x):
                r = exact_R_of_x(x)
                return (1-2/r)*(angular/r**2-6/r**3)
            derivatives = []
            for order, function in enumerate(derivative_functions):
                symbolic_value = function(peak, angular)
                inverse_value = mp.diff(potential_of_x, xpeak, order)
                difference = abs(symbolic_value-inverse_value)/(1+abs(symbolic_value))
                below("independent_tortoise_derivative_ell%d_order%d" % (ell, order), difference, "1e-65")
                derivatives.append(symbolic_value)
                output["derivative_comparison"].append({"ell": ell, "order": order,
                    "off_shell_chain_then_peak": dec(symbolic_value),
                    "Lambert_W_coordinate_inverse": dec(inverse_value),
                    "normalized_difference": dec(difference)})
            if not (peak>3 and derivatives[0]>0 and derivatives[2]<0):
                raise AssertionError(("peak regime", ell, peak, derivatives))
            D = mp.sqrt(-2*derivatives[2])
            for overtone in (0,):
                nu = mp.mpf(overtone)+mp.mpf("0.5")
                # mpmath arithmetic: do not let Python rational divisions
                # create double-precision constants in high-precision data.
                v2, v3, v4, v5, v6 = derivatives[2:]
                B2 = (v4/v2)*(mp.mpf("0.25")+nu*nu)/8-(v3/v2)**2*(7+60*nu*nu)/288
                B3 = 5*(v3/v2)**4*(77+188*nu*nu)/6912
                B3 -= v3*v3*v4*(51+100*nu*nu)/(384*v2**3)
                B3 += (v4/v2)**2*(67+68*nu*nu)/2304
                B3 += v3*v5*(19+28*nu*nu)/(288*v2*v2)
                B3 -= (v6/v2)*(5+4*nu*nu)/288
                Omega3 = B3/(D*D)
                energies = [derivatives[0]-1j*nu*D,
                            derivatives[0]+B2-1j*nu*D,
                            derivatives[0]+B2-1j*nu*D*(1+Omega3)]
                for order, energy in enumerate(energies, 1):
                    frequency = mp.sqrt(energy)
                    if not (mp.re(frequency)>0 and mp.im(frequency)<0):
                        raise AssertionError(("WKB branch", ell, overtone, order, frequency))
                    below("WKB_square_root_ell%d_n%d_order%d" % (ell, overtone, order),
                          abs(frequency*frequency-energy), "1e-70")
                    row = {"ell": ell, "n": overtone, "order": order, "dps": 80,
                           "peak_R": dec(peak), "D": dec(D), "B2": dec(B2),
                           "Omega3": dec(Omega3), "omega_squared": complex_json(energy),
                           "w": complex_json(frequency)}
                    if overtone == 0:
                        reference = roots[(ell, 2048, 80)]
                        row["absolute_complex_error_to_CF"] = dec(abs(frequency-reference))
                        row["relative_complex_error_to_CF"] = dec(abs(frequency-reference)/abs(reference))
                    output["wkb"].append(row)
            rejected("ell%d_naive_radial_curvature_changes_WKB" % ell,
                     abs(mp.sqrt(-2*derivatives[2])-
                         mp.sqrt(-2*mp.diff(lambda r: (1-2/r)*(angular/r**2-6/r**3), peak, 2))),
                     "1e-3")
            rejected("ell%d_zero_WKB_root_residual_is_not_spectral_accuracy" % ell,
                     abs(mp.sqrt(derivatives[0]-0.5j*D)-roots[(ell, 2048, 80)]), "1e-4")
        # The first CF truncation can solve to tiny residual without giving
        # the deeply continued root. This is separate from local WKB error.
        shallow = solve_cf(2, 32, 80)
        below("shallow_CF_solved_residual", abs(cf(shallow, 2, 32)), "1e-65")
        rejected("shallow_CF_tiny_residual_is_not_deep_accuracy",
                 abs(shallow-roots[(2, 2048, 80)]), "1e-10")

    return output


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--output", type=Path, required=True)
    args = parser.parse_args()
    requirements = Path(__file__).with_name("wkb-order-comparison-requirements.txt")
    requirements_bytes = requirements.read_bytes()
    symbolic = exact_checks()
    print("Exact derivation checks passed; evaluating finite CF and WKB comparisons.", flush=True)
    numeric = numerical_checks(symbolic)
    result = {
        "status": "finite reproducible comparison; not a certified resonance-error bound",
        "scope": {
            "model": "Schwarzschild vacuum axial gravitational RW potential, ell=2 and ell=10",
            "dimensionless_variables": "b=GM, R=r/b, x=r_star/b, w=b*omega, W=b^2*V, c=1",
            "Fourier_convention": "exp(-i*omega*t)",
            "derivative_rule": "W_j=(f*d/dR)^j W before peak substitution, j=0,...,6",
            "WKB_orders": [1, 2, 3], "WKB_overtone_inputs": [0],
            "CF_mode_inputs": [
                {"ell": 2, "n": 0, "seed": ["0.374", "-0.089"],
                 "real_selection_window": ["0.30", "0.45"], "imag_selection_window": ["-0.12", "-0.05"]},
                {"ell": 10, "n": 0, "seed": ["1.997", "-0.096"],
                 "real_selection_window": ["1.90", "2.10"], "imag_selection_window": ["-0.12", "-0.05"]}],
            "CF_depths": [512, 1024, 2048], "CF_precisions_dps": [50, 80],
            "CF_deep_check_depth": 4096, "CF_deep_check_dps": 100,
            "CF_boundary_condition": "Leaver RW outgoing prefactor and minimal series, zero ratio terminal approximation",
            "limitations": [
                "Finite WKB order is a local asymptotic approximation, not a certified error bound.",
                "The general-n exact oscillator identity is not a numerical overtone survey.",
                "CF depth/precision and ratio checks share one RW ODE; they are not an independent continuum or Kerr calculation.",
                "The general static F,H counterexample is kinematic, not asserted to solve vacuum Einstein equations.",
                "Pöschl--Teller check assumes V0>alpha^2/4 and is not a Schwarzschild spectrum.",
                "The semiclassical h expansion is formal bookkeeping; h=1 is not automatically a small parameter.",
                "Decreasing WKB errors in these two cases do not establish general monotonic convergence.",
                "110-place strings round-trip retained values; they are not 110 certified spectral digits. A zero depth drift can be a precision plateau.",
            ],
        },
        "source_sha256": hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
        "requirements_sha256": hashlib.sha256(requirements_bytes).hexdigest(),
        "source_conventions": {
            "Iyer_1987_II": "Phys.Rev.D35,3632; p.3633 Eqs.(1.2),(1.3a),(1.3b), expanded omega_squared with exp(-i omega t)",
            "corrections": "D=sqrt(-2 V2)>0; E1=V0-i nu D; E2=E1+B2; E3=E2-i nu D Omega3; Omega3=B3/D^2",
            "general_static_metric": "ds^2=-F dt^2+dr^2/H+r^2 dOmega^2; U=F/r^2; lambda_t^2=-Fc Hc Urr/(2 Uc)",
            "counterexample": "F=1-2b/r, H=F[1+eta(3b/r)^2], eta>=0; rc=3b, Omega_c^2=1/(27b^2), lambda_t^2=(1+eta)/(27b^2)",
        },
        "versions": {"mpmath": mp.__version__, "sympy": s.__version__},
        "exact_residuals": EXACT, "exact_negative_controls": NEGATIVE,
        "independent_oscillator": symbolic["oscillator"],
        "off_shell_tortoise_derivatives": list(map(str, symbolic["derivatives"])),
        "exterior_peak_R": str(symbolic["peak"]),
        "numerical": numeric,
    }
    result["counts"] = {"exact_zero_residuals": len(EXACT), "exact_negative_controls": len(NEGATIVE),
                        "numeric_residual_checks": len(numeric["numeric_checks"]),
                        "numeric_negative_controls": len(numeric["negative_controls"]),
                        "CF_convergence_runs": len(numeric["continued_fraction"]),
                        "CF_additional_deep_reference_solves": 2,
                        "CF_shallow_negative_control_solves": 1,
                        "WKB_rows": len(numeric["wkb"])}
    args.output.write_text(json.dumps(result, indent=2, sort_keys=True, ensure_ascii=False)+"\n")
    print(json.dumps(result["counts"], sort_keys=True))


if __name__ == "__main__":
    main()