"""Bounded Schwarzschild l=2, s=-2 continued-fraction benchmark.

Run: python3 frobenius-schwarzschild.py --output results.json
Requires mpmath 1.3.0 and SymPy 1.14.0; installs nothing.
The case list, depths, precisions, tolerances and independent checks are
deliberately explicit. This is not a Kerr solver or a generic mode-search API.
"""

import argparse
import hashlib
import json
from pathlib import Path
import platform

import mpmath as mp
import sympy as sy


CASES = (
    {"name": "fundamental", "n": 0, "inversion": 0,
     "seeds": (("0.37", "-0.09"), ("0.38", "-0.08")),
     "real_window": ("0.36", "0.39"), "imag_window": ("-0.11", "-0.07")},
    {"name": "first_overtone", "n": 1, "inversion": 1,
     "seeds": (("0.34", "-0.27"), ("0.35", "-0.28")),
     "real_window": ("0.33", "0.36"), "imag_window": ("-0.30", "-0.25")},
)
DEPTHS = (64, 128, 256, 512, 1024)
PRECISIONS = (30, 50, 80)
TAILS = ("zero", "nollert_C1_C2")
CHECK_DEPTHS = (2048, 4096)
REFERENCE_DEPTHS = (4096, 8192)
CHECK_DPS = 100
STORED_DIGITS = 90
MAX_STEPS = 50
SEAM = 12
PREFIX_LAST = 65
ACCURACY_TOL = "1e-25"
SELF_CONSISTENCY_TOL = "1e-85"


def scalar(value, digits=STORED_DIGITS):
    return mp.nstr(value, n=digits, strip_zeros=False)


def complex_data(value, digits=STORED_DIGITS):
    return {"real": scalar(mp.re(value), digits), "imag": scalar(mp.im(value), digits)}


def complex_value(value):
    return mp.mpc(value["real"], value["imag"])


def seed(pair):
    return mp.mpc(*pair)


def finite(value):
    return bool(mp.isfinite(mp.re(value)) and mp.isfinite(mp.im(value)))


def rw_coefficients(w):
    """RW recurrence, w=GM omega, s=-2 and ell=2."""
    iw = 1j * w
    def coefficient(j):
        return ((j + 1) * (j + 1 - 4 * iw),
                -2 * j * j + (16 * iw - 2) * j + 32 * w * w + 8 * iw - 3,
                (j - 4 * iw) ** 2 - 4)
    return coefficient


def teukolsky_coefficients(w):
    """Independent unsimplified D-coefficient specialization at a=0.

    Do not substitute the simplified RW alpha/gamma here. The off-diagonal
    coefficients differ, though adjacent alpha_j gamma_(j+1) products agree.
    """
    spin, ell = -2, 2
    angular = ell * (ell + 1) - spin * (spin + 1)
    zeta, xi, eta, p = 1j * w, -spin - 2j * w, 0, 1j * w
    aa = 1 + spin + xi + eta - 2 * zeta + spin
    gg, dd = 1 + spin + 2 * eta, 1 + spin + 2 * xi
    ss = (angular - 8 * w * w + p * (2 * aa + gg - dd)
          + (1 + spin - (gg + dd) / 2) * (spin + (gg + dd) / 2))
    d0, d1, d2 = dd, 4 * p - 2 * aa + gg - dd - 2, 2 * aa - gg + 2
    d3, d4 = aa * (4 * p - dd) - ss, aa * (aa - gg + 1)
    def coefficient(j):
        return (j * j + (d0 + 1) * j + d0,
                -2 * j * j + (d1 + 2) * j + d3,
                j * j + (d2 - 3) * j + d4 - d2 + 2)
    return coefficient


def forward_tail(w, depth, tail, formulation="RW"):
    if tail == "zero":
        return mp.mpc(0)
    if tail != "nollert_C1_C2":
        raise ValueError("Unsupported tail")
    c1 = mp.sqrt(-4j * w)
    if mp.re(c1) > 0:
        c1 = -c1
    if not mp.re(c1) < 0:
        raise ValueError("This finite benchmark excludes the equal-real-part branch case")
    c2 = -mp.mpf(3) / 4 - 4j * w
    if formulation == "Teukolsky":
        c2 -= 2
    elif formulation != "RW":
        raise ValueError("Unsupported recurrence")
    return 1 + c1 / mp.sqrt(depth) + c2 / depth


def rw_fraction(w, depth, inversion, tail):
    """Forward ratio R_j=a_(j+1)/a_j, propagated through row j+1."""
    coefficient = rw_coefficients(w)
    right_ratio = forward_tail(w, depth, tail)
    for j in range(depth, inversion, -1):
        alpha, beta, gamma = coefficient(j)
        right_ratio = -gamma / (beta + alpha * right_ratio)
    alpha, beta, gamma = coefficient(inversion)
    right = alpha * right_ratio
    left = mp.mpc(0)
    if inversion:
        a0, b0, _ = coefficient(0)
        left_ratio = -b0 / a0
        for j in range(1, inversion):
            aj, bj, gj = coefficient(j)
            left_ratio = -(bj + gj / left_ratio) / aj
        left = gamma / left_ratio
    terms = (beta, right, left)
    return sum(terms), sum(abs(term) for term in terms)


def teukolsky_fraction(w, depth, inversion, keep_denominators=False):
    """Independent denominator elimination with a zero tail, not ratio code."""
    coefficient = teukolsky_coefficients(w)
    _, denominator, _ = coefficient(depth)
    denominators = {depth: denominator} if keep_denominators else None
    for j in range(depth - 1, inversion, -1):
        aj, bj, _ = coefficient(j)
        _, _, next_gamma = coefficient(j + 1)
        denominator = bj - aj * next_gamma / denominator
        if keep_denominators:
            denominators[j] = denominator
    alpha, beta, _ = coefficient(inversion)
    _, _, next_gamma = coefficient(inversion + 1)
    right = -alpha * next_gamma / denominator
    left = mp.mpc(0)
    if inversion:
        _, left_denominator, _ = coefficient(0)
        for j in range(1, inversion):
            previous_alpha, _, _ = coefficient(j - 1)
            _, bj, gj = coefficient(j)
            left_denominator = bj - previous_alpha * gj / left_denominator
        previous_alpha, _, _ = coefficient(inversion - 1)
        _, _, gamma = coefficient(inversion)
        left = -previous_alpha * gamma / left_denominator
    terms = (beta, right, left)
    return sum(terms), sum(abs(term) for term in terms), denominators


def normalized(value, scale):
    if not finite(value) or not finite(scale) or not scale > 0:
        raise ArithmeticError("Nonfinite or zero residual scale")
    return abs(value) / scale


def in_window(w, case):
    return bool(mp.mpf(case["real_window"][0]) < mp.re(w) < mp.mpf(case["real_window"][1])
                and mp.mpf(case["imag_window"][0]) < mp.im(w) < mp.mpf(case["imag_window"][1]))


def equation_residual(coefficient, j, previous, current, following):
    alpha, beta, gamma = coefficient(j)
    terms = [alpha * following, beta * current]
    if j:
        terms.append(gamma * previous)
    return normalized(sum(terms), sum(abs(term) for term in terms))


def reconstruct_sequence(w, depth):
    """Left initial solution and independently selected right solution meet at J.

    Only the seam equation is not enforced while constructing the sequence.
    Non-seam row residuals measure arithmetic/serialization self-consistency,
    not convergence of the infinite minimal solution or of the physical ODE.
    """
    rw = rw_coefficients(w)
    tk = teukolsky_coefficients(w)
    _, _, denominators = teukolsky_fraction(w, depth, 0, True)
    a = [mp.mpc(1)]
    alpha, beta, _ = rw(0)
    a.append(-beta / alpha)
    for j in range(1, SEAM):
        alpha, beta, gamma = rw(j)
        a.append(-(beta * a[j] + gamma * a[j - 1]) / alpha)
    for j in range(SEAM, depth):
        _, _, next_gamma = tk(j + 1)
        teuk_ratio = -next_gamma / denominators[j + 1]
        dnext = (j + 1 - 4j * w) / (j + 3 - 4j * w)
        a.append(a[-1] * teuk_ratio / dnext)
    a.append(mp.mpc(0))
    indices = list(range(PREFIX_LAST + 1)) + [depth - 1, depth, depth + 1]
    stored = {str(j): complex_data(a[j]) for j in indices}
    # Compute all quoted residuals AFTER parsing exactly the stored values.
    parsed = {int(j): complex_value(value) for j, value in stored.items()}
    row_residuals = {}
    for j in list(range(PREFIX_LAST)) + [depth]:
        row_residuals[str(j)] = scalar(equation_residual(
            rw, j, parsed.get(j - 1, mp.mpc(0)), parsed[j], parsed[j + 1]))
    nonseam = {j: mp.mpf(value) for j, value in row_residuals.items() if int(j) != SEAM}
    allrows = {j: mp.mpf(value) for j, value in row_residuals.items()}
    worst = max(allrows, key=allrows.get)
    return {
        "normalization": "a_0=1", "seam_index": SEAM,
        "right_selection": "deep raw-Teukolsky zero-tail ratios mapped to RW",
        "right_depth": depth,
        "stored_digits": STORED_DIGITS,
        "stored_coefficients": stored,
        "normalized_row_residuals": row_residuals,
        "initial_row_residual": row_residuals["0"],
        "seam_row_residual": row_residuals[str(SEAM)],
        "far_truncation_row_residual": row_residuals[str(depth)],
        "imposed_far_boundary_abs": scalar(abs(parsed[depth + 1])),
        "max_nonseam_residual": scalar(max(nonseam.values())),
        "max_all_stored_row_residual": scalar(allrows[worst]),
        "worst_stored_row_index": int(worst),
        "nonseam_self_consistency_pass": bool(max(nonseam.values()) <= mp.mpf(SELF_CONSISTENCY_TOL)),
        "seam_accuracy_pass": bool(mp.mpf(row_residuals[str(SEAM)]) <= mp.mpf(ACCURACY_TOL)),
    }


def exact_checks():
    n, w, spin, ell = sy.symbols("n w spin ell")
    checks, controls = {}, {}
    def zero(name, value):
        result = sy.simplify(value)
        if result != 0:
            raise AssertionError((name, result))
        checks[name] = str(result)
    def reject(name, value):
        result = sy.simplify(value)
        if result == 0 or result.has(sy.nan, sy.zoo, sy.oo):
            raise AssertionError((name, result))
        controls[name] = {"rejected": True, "residual": str(result)}
    ar = (n + 1) * (n + 1 - 4 * sy.I * w)
    gr = (n - 4 * sy.I * w) ** 2 - spin ** 2
    at = (n + 1) * (n + 1 - spin - 4 * sy.I * w)
    gt = (n - 4 * sy.I * w) * (n + spin - 4 * sy.I * w)
    beta = -2*n*n + (16*sy.I*w-2)*n + 32*w*w + 8*sy.I*w + spin*spin-1-ell*(ell+1)
    dnext = (n + 1 - 4 * sy.I * w) / (n + 1 - spin - 4 * sy.I * w)
    zero("RW_Teukolsky_adjacent_products", ar * gr.subs(n, n+1) - at * gt.subs(n, n+1))
    zero("coefficient_rescaling_alpha", at*dnext-ar)
    zero("coefficient_rescaling_gamma", gt/dnext.subs(n,n-1)-gr)
    angular = ell*(ell+1)-spin*(spin+1)
    zeta, xi, p = sy.I*w, -spin-2*sy.I*w, sy.I*w
    aa, gg, dd = 1+2*spin+xi-2*zeta, 1+spin, 1+spin+2*xi
    ss=angular-8*w*w+p*(2*aa+gg-dd)+(1+spin-(gg+dd)/2)*(spin+(gg+dd)/2)
    D=[dd,4*p-2*aa+gg-dd-2,2*aa-gg+2,aa*(4*p-dd)-ss,aa*(aa-gg+1)]
    zero("raw_D_alpha_specialization",n*n+(D[0]+1)*n+D[0]-at)
    zero("raw_D_beta_specialization",-2*n*n+(D[1]+2)*n+D[3]-beta)
    zero("raw_D_gamma_specialization",n*n+(D[2]-3)*n+D[4]-D[2]+2-gt)
    h, c1, c2, c3 = sy.symbols("h c1 c2 c3")
    rn=1+c1*h+c2*h*h+c3*h**3
    rm=1+c1*h/sy.sqrt(1-h*h)+c2*h*h/(1-h*h)+c3*h**3/(1-h*h)**sy.Rational(3,2)
    for name, alpha, gamma, correct_c2 in (
        ("RW",ar,gr,-sy.Rational(3,4)-4*sy.I*w),
        ("Teukolsky",at,gt,spin-sy.Rational(3,4)-4*sy.I*w)):
        expansion=sy.series(h**4*(alpha.subs(n,h**-2)*rn+beta.subs(n,h**-2)+gamma.subs(n,h**-2)/rm),h,0,4).removeO().expand()
        zero(name+"_tail_constant",expansion.coeff(h,0))
        zero(name+"_tail_half_power",expansion.coeff(h,1))
        zero(name+"_tail_C1_equation",expansion.coeff(h,2)-c1*c1-4*sy.I*w)
        reduced=sy.rem(expansion.coeff(h,3),c1*c1+4*sy.I*w,c1)
        zero(name+"_tail_C2_equation",reduced.subs(c2,correct_c2))
    r=lambda j:(j+1)/(j+2)
    correct=-(n+1)/(-2*(n+2)+(n+3)*r(n+1))
    wrong=-n/(-2*(n+1)+(n+2)*r(n+1))
    zero("forward_ratio_next_row_index",correct-r(n))
    reject("forward_ratio_unshifted_row_index",(wrong-r(n)).subs(n,1))
    reject("RW_C2_used_for_spin_minus_two_Teukolsky",sy.Integer(-2))
    reject("forward_tail_used_as_negative_backward_ratio",sy.Integer(2))
    return checks, controls


def references():
    result = {}
    with mp.workdps(CHECK_DPS):
        for case in CASES:
            rows=[]
            for depth in REFERENCE_DEPTHS:
                calls=0
                def function(w):
                    nonlocal calls
                    calls+=1
                    return teukolsky_fraction(w,depth,case["inversion"])[0]
                w=mp.findroot(function,tuple(seed(value) for value in case["seeds"]),
                              tol=mp.mpf("1e-90"),maxsteps=MAX_STEPS)
                solved_residual=abs(function(w))
                if not solved_residual<=mp.mpf("1e-90"):
                    raise AssertionError("Reference solver absolute tolerance failed")
                saved=complex_data(w)
                w=complex_value(saved)
                f,scale,_=teukolsky_fraction(w,depth,case["inversion"])
                rows.append({"depth":depth,"dps":CHECK_DPS,"tail":"zero",
                             "omega_GM":saved,"function_evaluations":calls,
                             "pre_serialization_absolute_residual":scalar(solved_residual),
                             "solver_absolute_tolerance_pass":True,
                             "stored_digits":STORED_DIGITS,
                             "normalized_residual":scalar(normalized(f,scale)),
                             "in_mode_window":in_window(w,case)})
            drift=abs(complex_value(rows[1]["omega_GM"])-complex_value(rows[0]["omega_GM"]))
            if not all(row["in_mode_window"] for row in rows) or not drift < mp.mpf("1e-55"):
                raise AssertionError("Reference depth comparison failed")
            result[case["name"]]={"runs":rows,"reference_depth_drift":scalar(drift),
                                  "adopted_omega_GM":rows[-1]["omega_GM"],
                                  "reference_uncertainty_indicator":scalar(max(10*drift,mp.mpf("1e-88")))}
            print("Reference ready: "+case["name"],flush=True)
    return result


def solve_case(case, depth, dps, tail, reference):
    record={"case":case["name"],"n":case["n"],"inversion_index":case["inversion"],
            "depth":depth,"dps":dps,"tail":tail,"maxsteps":MAX_STEPS,
            "seeds":[{"real":a,"imag":b} for a,b in case["seeds"]],
            "left_elimination_steps":case["inversion"],
            "right_elimination_steps":depth-case["inversion"]}
    try:
        with mp.workdps(dps):
            tolerance=mp.power(10,-(dps-8))
            calls=0
            def function(w):
                nonlocal calls
                calls+=1
                return rw_fraction(w,depth,case["inversion"],tail)[0]
            w=mp.findroot(function,tuple(seed(value) for value in case["seeds"]),
                          tol=tolerance,maxsteps=MAX_STEPS)
            record["omega_GM"]=complex_data(w,dps)
            record["function_evaluations"]=calls
            record["solver_absolute_tolerance"]=scalar(tolerance,dps)
        with mp.workdps(CHECK_DPS):
            w=complex_value(record["omega_GM"])
            f,scale=rw_fraction(w,depth,case["inversion"],tail)
            record["serialized_root_solve_absolute_residual"]=scalar(abs(f))
            record["serialized_root_solve_normalized_residual"]=scalar(normalized(f,scale))
            record["solver_tolerance_pass"]=bool(abs(f)<=mp.mpf(record["solver_absolute_tolerance"]))
            record["in_mode_window"]=in_window(w,case)
            record["independent_checks"]=[]
            for check_depth in CHECK_DEPTHS:
                df,ds,_=teukolsky_fraction(w,check_depth,case["inversion"])
                record["independent_checks"].append({"formulation":"raw_Teukolsky_D",
                    "tail":"zero","depth":check_depth,"dps":CHECK_DPS,
                    "absolute_residual":scalar(abs(df)),"normalized_residual":scalar(normalized(df,ds))})
            record["reference_absolute_error"]=scalar(abs(w-complex_value(reference["adopted_omega_GM"])))
            other_inversion=1-case["inversion"]
            xf,xs,_=teukolsky_fraction(w,CHECK_DEPTHS[-1],other_inversion)
            record["cross_inversion_check"]={"inversion_index":other_inversion,
                "formulation":"raw_Teukolsky_D","depth":CHECK_DEPTHS[-1],"dps":CHECK_DPS,
                "tail":"zero","normalized_residual":scalar(normalized(xf,xs))}
            record["coefficient_sequence"]=reconstruct_sequence(w,CHECK_DEPTHS[-1])
            independent_ok=all(mp.mpf(check["normalized_residual"])<=mp.mpf(ACCURACY_TOL)
                               for check in record["independent_checks"])
            record["independent_accuracy_pass"]=bool(independent_ok
                and normalized(xf,xs)<=mp.mpf(ACCURACY_TOL)
                and mp.mpf(record["reference_absolute_error"])<=mp.mpf(ACCURACY_TOL))
            record["accepted_at_target_tolerance"]=bool(record["solver_tolerance_pass"]
                and record["in_mode_window"] and record["independent_accuracy_pass"]
                and record["coefficient_sequence"]["nonseam_self_consistency_pass"]
                and record["coefficient_sequence"]["seam_accuracy_pass"])
            record["status"]="converged_root"
    except (ValueError, ZeroDivisionError, ArithmeticError) as error:
        record["status"]="failed"
        record["error_type"]=type(error).__name__
        record["error"]=str(error)
        record["accepted_at_target_tolerance"]=False
    return record


def numerical_controls(records, reference):
    controls={}
    with mp.workdps(CHECK_DPS):
        def reject(name,value,tolerance):
            if not finite(value) or not abs(value)>tolerance:
                raise AssertionError((name,value,tolerance))
            controls[name]={"rejected":True,"nonzero_witness":scalar(abs(value)),
                            "rejection_threshold":scalar(tolerance)}
        probe=mp.mpc("0.3","-0.2")
        f0=rw_fraction(probe,128,0,"zero")[0]
        f1=rw_fraction(probe,128,1,"zero")[0]
        reject("inversion_index_ignored_at_generic_probe",f0-f1,mp.mpf("1e-3"))
        fundamental=complex_value(reference["fundamental"]["adopted_omega_GM"])
        if in_window(fundamental,CASES[1]):
            raise AssertionError("Wrong mode not rejected")
        reject("fundamental_substituted_for_first_overtone",fundamental-complex_value(
            reference["first_overtone"]["adopted_omega_GM"]),mp.mpf("0.1"))
        f,scale,_=teukolsky_fraction(fundamental+mp.mpf("1e-5"),4096,0)
        reject("perturbed_root_fails_deep_fraction",normalized(f,scale),mp.mpf(ACCURACY_TOL))
        row=next(item for item in records if item["case"]=="first_overtone"
                 and item["depth"]==64 and item["dps"]==80 and item["tail"]=="zero")
        if not row["solver_tolerance_pass"] or row["accepted_at_target_tolerance"]:
            raise AssertionError("Need an honestly retained solver-only false positive")
        reject("solver_residual_alone_false_positive",mp.mpf(row["independent_checks"][-1]["normalized_residual"]),
               mp.mpf(ACCURACY_TOL))
        good=next(item for item in records if item["case"]=="fundamental"
                  and item["depth"]==1024 and item["dps"]==80 and item["tail"]=="nollert_C1_C2")
        w=complex_value(good["omega_GM"])
        stored=good["coefficient_sequence"]["stored_coefficients"]
        a={int(j):complex_value(value) for j,value in stored.items()}
        coeff=rw_coefficients(w)
        reject("poisoned_initial_coefficient",equation_residual(coeff,0,0,a[0],a[1]*(1+mp.mpf("1e-3"))),mp.mpf("1e-5"))
        reject("poisoned_seam_coefficient",equation_residual(coeff,SEAM,a[SEAM-1],a[SEAM],
               a[SEAM+1]*(1+mp.mpf("1e-3"))),mp.mpf("1e-5"))
        far=CHECK_DEPTHS[-1]
        reject("poisoned_far_boundary",equation_residual(coeff,far,a[far-1],a[far],a[far]),mp.mpf("0.1"))
        reject("Teukolsky_coefficient_used_without_RW_rescaling",equation_residual(teukolsky_coefficients(w),
               5,a[4],a[5],a[6]),mp.mpf("1e-3"))
        # Evaluate the alternate asymptotic branch at a finite depth. This is a
        # tail initialization defect, not a claim that every wrong-tail solve fails.
        depth=64
        c1=mp.sqrt(-4j*w)
        if mp.re(c1)<0:
            c1=-c1
        wrong=1+c1/mp.sqrt(depth)+(-mp.mpf(3)/4-4j*w)/depth
        correct=forward_tail(w,depth,"nollert_C1_C2")
        reject("dominant_tail_branch_selected",wrong-correct,mp.mpf("0.1"))
        wrong_c2=forward_tail(w,depth,"nollert_C1_C2","Teukolsky")
        reject("Teukolsky_tail_inserted_into_RW",wrong_c2-correct,mp.mpf("0.01"))
    return controls


def main():
    parser=argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--output",type=Path,required=True)
    args=parser.parse_args()
    exact,exact_controls=exact_checks()
    reference=references()
    records=[]
    for case in CASES:
        for dps in PRECISIONS:
            for tail in TAILS:
                for depth in DEPTHS:
                    records.append(solve_case(case,depth,dps,tail,reference[case["name"]]))
                print(f"Sweep ready: {case['name']} dps={dps} tail={tail}",flush=True)
    verification_failure=None
    try:
        controls=numerical_controls(records,reference)
    except (AssertionError,StopIteration,KeyError) as error:
        # Preserve failed solves and their inputs even if a dependent negative
        # control cannot run; do not present an incomplete check set as passed.
        controls={}
        verification_failure={"error_type":type(error).__name__,"error":str(error)}
    with mp.workdps(CHECK_DPS):
        for case in CASES:
            for dps in PRECISIONS:
                for tail in TAILS:
                    subset=[r for r in records if r["case"]==case["name"] and r["dps"]==dps and r["tail"]==tail]
                    for previous,current in zip(subset,subset[1:]):
                        if previous["status"]==current["status"]=="converged_root":
                            current["previous_depth_root_drift"]=scalar(abs(complex_value(current["omega_GM"])
                                -complex_value(previous["omega_GM"])))
        plots=[]
        for row in records:
            if row["status"]!="converged_root":
                continue
            floor=max(mp.power(10,-(row["dps"]-3)),mp.mpf(reference[row["case"]]["reference_uncertainty_indicator"]))
            plots.append({"case":row["case"],"n":row["n"],"depth":row["depth"],"dps":row["dps"],"tail":row["tail"],
                "absolute_frequency_error":row["reference_absolute_error"],"display_floor":scalar(floor),
                "at_or_below_display_floor":bool(mp.mpf(row["reference_absolute_error"])<=floor),
                "deep_normalized_residual":row["independent_checks"][-1]["normalized_residual"],
                "seam_normalized_residual":row["coefficient_sequence"]["seam_row_residual"]})
    summary={"exact_zero_count":len(exact),"exact_negative_control_count":len(exact_controls),
        "numerical_negative_control_count":len(controls),"scheduled_solve_count":len(records),
        "converged_root_count":sum(r["status"]=="converged_root" for r in records),
        "failed_solve_count":sum(r["status"]=="failed" for r in records),
        "accepted_at_target_tolerance_count":sum(r["accepted_at_target_tolerance"] for r in records),
        "retained_below_target_accuracy_count":sum(r["status"]=="converged_root" and not r["accepted_at_target_tolerance"] for r in records)}
    output={"schema_version":1,"scope":"Schwarzschild axial ell=2, s=-2, explicit fundamental and first-overtone cases only; no Kerr implementation or independent continuum ODE certificate.",
        "conventions":{"frequency":"w=GM omega; exp(-i omega t); c=1 with G retained on restoring units",
            "forward_ratio":"R_j=a_(j+1)/a_j",
            "backward_relation":"R_(j-1)=-gamma_j/(beta_j+alpha_j R_j)",
            "inverted_fraction":"F_v=beta_v+alpha_v R_v+gamma_v/R_(v-1,left), with no left term at v=0",
            "zero_tail":"R_N=0, equivalently a_(N+1)=0",
            "RW_tail":"R_N=1+C1/sqrt(N)+C2/N; C1^2=-4 i w; Re C1<0; C2=-3/4-4 i w",
            "Teukolsky_tail":"C2=s-3/4-4 i w; not used in the zero-tail independent checks",
            "residual_normalization":"absolute sum divided by sum of absolute term magnitudes; zero/nonfinite scale rejected"},
        "parameters":{"spin":-2,"ell":2,"cases":CASES,"depths":DEPTHS,"precisions":PRECISIONS,
            "tails":TAILS,"solver_tolerance_rule":"10^(-(dps-8)) absolute; explicitly reevaluated at serialized root",
            "maxsteps":MAX_STEPS,"independent_check_depths":CHECK_DEPTHS,"independent_check_dps":CHECK_DPS,
            "reference_depths":REFERENCE_DEPTHS,"reference_dps":CHECK_DPS,"reference_solver_tolerance":"1e-90",
            "target_accuracy_tolerance":ACCURACY_TOL,"nonseam_self_consistency_tolerance":SELF_CONSISTENCY_TOL,
            "coefficient_stored_digits":STORED_DIGITS,"coefficient_seam_index":SEAM,"coefficient_prefix_last_index":PREFIX_LAST},
        "limitations":[
            "The two recurrences encode the same Schwarzschild spectral representation. Separate formulas, ratio-versus-denominator evaluation, greater depth and precision test implementation and truncation, not an independent continuum method.",
            "An inversion index rearranges the same spectral equation where its denominators are nonzero; it is not a theorem assigning the nth root. Explicit seeds and disjoint windows select these two cases.",
            "All numerical divisions require nonzero denominators. Exceptional frequencies, algebraically special modes, pure-imaginary/Stokes branch cases and root searches outside these windows are excluded.",
            "A two-term asymptotic tail is not exact. Differences between tails and depths are empirical convergence evidence for these cases, not a universal error bound.",
            "The reference uncertainty indicator is ten times a finite depth drift, floored at 1e-88; it is not a rigorous error bound.",
            "Reference roots are solved at 100 digits and their absolute residual is checked before serialization. Stored reference roots and their recomputed normalized residuals use the declared 90-digit decimals.",
            "Stored non-seam recurrence residuals verify construction and decimal serialization. The seam is independently matched; the far boundary is an explicitly imposed finite zero-tail condition, not the exact infinite endpoint.",
            "A small root-finder residual alone is insufficient. Coarse converged roots failing independent checks remain present and labeled below target accuracy.",
            "The displayed error floor is a conservative plotting convention based on decimal precision and reference drift, not measured physical uncertainty.",
            "No table of externally tabulated frequencies, qnm package, angular solver, Kerr root, or time-domain evolution is used."],
        "exact_checks":exact,"exact_negative_controls":exact_controls,"numerical_negative_controls":controls,
        "references":reference,"runs":records,"plot_data":plots,"summary":summary,
        "source_sha256":hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),
        "runtime":{"python":platform.python_version(),"mpmath":mp.__version__,"sympy":sy.__version__}}
    output["verification_failure"]=verification_failure
    args.output.write_text(json.dumps(output,indent=2,sort_keys=True)+"\n")
    print(json.dumps(summary,sort_keys=True),flush=True)
    if verification_failure is not None or summary["failed_solve_count"]:
        raise SystemExit(1)


if __name__=="__main__":
    main()