"""RUN-010 - Round 09 (v0.9): the 270-degree cusp in arbitrary precision, its monotone Darboux enclosure, the branch
audit, the cusps' one-sided slopes and the peak fit.

Round 09 takes Round 08's best candidate (eps = 0.037) as a curve defined by decimal constants, with tangent angle
theta(u) = beta + delta (g(u) - g(1)) / (g(0) - g(1)), g(u) = (1 - tanh((u - c)/eps)) / 2, and reduces its scale at
phi = 3pi/2 to s270 = y1/A + (x1 + l0)/B, two one-dimensional integrals. It reports s270 to 120 digits from mpmath's
tanh-sinh and Gauss-Legendre quadratures at 40 to 120 digits; the event five-link's four branches and their minimum
s0; the difference s270 - s0 = 1.0581952124e-5; an enclosure of s270 by monotone left and right sums (lower bound
minus s0 > 7.8e-6 at 4,096 parts); the candidate's four branches in high precision; difference quotients for the
slopes at the two cusps; and low-order fits of Round 08's peak scan. It leaves the full phase-circle exclusion open.

Checks (each returns a list of failures; empty = pass)
  P1 definitions  the candidate's and the event control's decimal constants are Round 08's best parameters and
                  Round 05's event root, as shortest decimals; the verifier script and the report's section 2 use
                  them; the g-form of the layer is RUN-009's F
  P2 cusp value   this run's s270, by two quadratures of its own at 200 digits, holds the summary's, the manifest's
                  and the audit note's 120 digits, and each convergence row to its working precision; the two
                  methods' rows are identical
  P3 event        the event five-link's four branches, from its vertices at 200 digits: the summary's values, its
                  minimum and the difference, the manifest's and the audit note's, each to its printed digits (the
                  difference of two 120-digit values to 120 decimals)
  P4 Darboux      theta decreases inside (0, pi/2), so the left and right sums enclose s270; each row's lower,
                  upper, width and lower minus the event minimum, and the audit note's, to 99 decimals (the sums
                  carry about 100 digits, the package verifier's working precision)
  P5 branches     each branch-audit value is the exact scale at its phase (a low branch's printed phase, a cusp's
                  exact phase, against which the printed one is checked), to 1e-68; its labels name supports
                  active there; each low phase lies within 1e-7 rad of a stationary point; the margins
  P6 slopes       each slope is the one-sided difference quotient at its h, to 1e-45, and within h of the exact
                  one-sided derivative; at both cusps D- < 0 < D+
  P7 peak fit     each fit is the least-squares polynomial of its degree through Round 08's six peak-scan scales:
                  its coefficients, eps*, scale and second derivative
  P8 consistency  the report's tables and short figures, the README's lists and Darboux figure, the draft, the
                  manifest, and the summary's copies of the data files
  P9 digits       every long figure of the README and the report, the branch audit's values and the report's
                  derivatives hold to their printed digits (fixed-text failures, one per figure)

Drill: each plant names the checks that must trip - trip meaning a failure the unplanted package does not already
have; controls must trip nothing.

Usage:
  python code/run010_round09_precision_cusp.py [--drill [--log PATH] | --blind P#]
"""
import argparse
import copy
import csv
import io
import json
import math
import pathlib
import re
import warnings

import mpmath as mp
import numpy as np

from corpus import TREE, sha256, source_bytes
import gate
from run001_corpus_identity import restore_escapes
from run004_round03_symmetric_links import load

RID = "round-v0.9"
LOG = TREE / "data" / "gate-logs" / "run010-round09-precision-cusp.json"
iv = mp.iv


# ---------------------------------------------------------------- the world

def cusp_270(x1, y1, l0, A, B):
    """the scale at 3pi/2, where the two endpoints tie on the leg: y1/A + (x1 + l0)/B"""
    return y1 / A + (x1 + l0) / B


def build_world():
    f, f8, f5 = load(RID), load("round-v0.8.3"), load("round-v0.5")
    md = lambda n: restore_escapes(f[n].decode("utf-8"))
    rows = lambda n: list(csv.DictReader(io.StringIO(f[n].decode("utf-8"))))
    return {
        "readme": md("README.md"), "report": md("reports/ROUND_09_REPORT.md"),
        "audit": md("reports/ROUND_09_ARBITRARY_PRECISION_AUDIT.md"),
        "draft": md("reports/ROUND_09_INTERVAL_CERTIFICATE_DRAFT.md"),
        "verifier": f["src/high_precision_verify.py"].decode("utf-8"),
        "manifest": json.loads(f["MANIFEST.json"]),
        "members": sorted(f),
        "summary": json.loads(f["data/round09_summary.json"]),
        "branch_audit": json.loads(f["data/branch_audit.json"]),
        "darboux": rows("data/darboux_cusp_enclosure.csv"),
        "convergence": rows("data/precision_convergence.csv"),
        "peak_fit": json.loads(f["data/epsilon_peak_fit.json"]),
        "round08_best": json.loads(f8["data/round08_summary.json"])["best_verified_candidate"],
        "round08_peak": json.loads(f8["data/round08_peak_scan.json"])["peak_scan"],
        "round05": json.loads(f5["data/round05_summary.json"])["event_solution"]["parameters"],
        "legs": ("1", "1"),
        "cusp_formula": cusp_270,
    }


def ctx(w):
    """what the high-precision values depend on: the definitions, the triangle and the cusp formula"""
    return json.dumps(w["summary"]["definitions"], sort_keys=True), tuple(w["legs"]), id(w["cusp_formula"])


def legs_mp(w):
    fa, fb = (mp.mpf(x) for x in w["legs"])
    return (3 + 4 * mp.sqrt(3)) / 18 / mp.mpf("1.0048") * fa, (4 + mp.sqrt(3)) / 6 / mp.mpf("1.0048") * fb


def holds(s, true, cap=None):
    """a printed decimal holds when it is within one unit of its last digit of the value (rounded or cut); with a
    cap, a figure printed past its working precision is asked to hold to that many decimals only"""
    s = s.strip().strip("`").rstrip("…")
    n = len(s.split(".")[1])
    if cap is not None:
        n = min(n, cap)
    with mp.workdps(n + 60):
        true = mp.mpf(true) if isinstance(true, str) else true
        return abs(mp.mpf(s) - true) <= mp.mpf(10) ** (-n)


def digits_held(s, true):
    """how many decimals of a printed figure are the value's: all of them, or those before its error's first digit"""
    s = s.strip().strip("`").rstrip("…")
    n = len(s.split(".")[1])
    with mp.workdps(n + 40):
        err = abs(mp.mpf(s) - true)
        return n if err <= mp.mpf(10) ** (-n) else int(mp.floor(-mp.log10(err)))


def rounds_to(printed, value):
    """a short printed figure (fixed or e-notation) is the value rounded to its significant digits"""
    m = re.fullmatch(r"(-?)(\d+)\.?(\d*)(?:e([+-]?\d+))?", printed.strip())
    if not m:
        return False
    mant = (m.group(2) + m.group(3)).lstrip("0") or "0"
    p = float(printed)
    if p == 0:
        return abs(value) < 1e-300 or value == 0
    unit = 10 ** (math.floor(math.log10(abs(p))) - len(mant) + 1)
    return abs(p - value) <= 0.501 * unit


# ---------------------------------------------------------------- the candidate in high precision

FUNS = ("max_d", "min_x", "min_y")


class Cand:
    """the candidate at a working precision: the tangent angle, wing positions by quadrature split at the layer, and
    the exact support of a linear functional - the four corners, and on each wing the one point whose tangent is
    perpendicular to it (theta is monotone, so that point comes from inverting the tanh)"""

    def __init__(self, w, dps):
        self.dps = dps
        self.memo = {}
        with mp.workdps(dps):
            D = w["summary"]["definitions"]["candidate"]
            self.eps, self.w, self.b, self.d, self.c = (mp.mpf(D[k]) for k in
                                                        ("epsilon", "wing_length", "beta", "delta", "transition_center"))
            self.T0, self.T1 = mp.tanh(self.c / self.eps), mp.tanh((self.c - 1) / self.eps)
            self.l0 = 1 - 2 * self.w
            self.A, self.B = legs_mp(w)
            self.a = {"max_d": (1 / self.A, 1 / self.B), "min_x": (-1 / self.A, mp.mpf(0)), "min_y": (mp.mpf(0), -1 / self.B)}
            x1, y1 = self.pos(mp.mpf(1))
            self.P = {"p0": (mp.mpf(0), mp.mpf(0)), "p1": (x1, y1), "p2": (x1 + self.l0, y1), "p3": (2 * x1 + self.l0, mp.mpf(0))}

    def theta(self, u):
        return self.b + self.d * (mp.tanh((self.c - u) / self.eps) - self.T1) / (self.T0 - self.T1)

    def pos(self, u, method="quad"):
        with mp.workdps(self.dps):
            sp = sorted({mp.mpf(0), u} | {self.c + k * self.eps for k in (-12, -6, -2, 0, 2, 6, 12) if 0 < self.c + k * self.eps < u})
            q = getattr(mp, method)
            return (self.w * q(lambda t: mp.cos(self.theta(t)), sp), self.w * q(lambda t: mp.sin(self.theta(t)), sp))

    def u_of(self, th):
        f = (th - self.b) / self.d
        if not 0 < f < 1:
            return None
        return self.c - self.eps * mp.atanh(self.T1 + f * (self.T0 - self.T1))

    def support(self, name, ph):
        """the functional's maximum on the curve rotated by ph, the points attaining it, and their phase derivatives"""
        a = self.a[name]
        cs, sn = mp.cos(ph), mp.sin(ph)
        n = (a[0] * cs + a[1] * sn, -a[0] * sn + a[1] * cs)
        dn = (-a[0] * sn + a[1] * cs, -a[0] * cs - a[1] * sn)
        pts = dict(self.P)
        x1 = self.P["p1"][0]
        for side, ang in (("L", mp.atan2(n[1], n[0])), ("R", mp.atan2(n[1], -n[0]))):
            for k in (-1, 1):
                t = ang + k * mp.pi / 2
                t -= 2 * mp.pi * mp.floor((t + mp.pi) / (2 * mp.pi))
                u = self.u_of(t)
                if u is not None:
                    px, py = self.pos(u)
                    pts[side] = (px, py) if side == "L" else (2 * x1 + self.l0 - px, py)
        vals = {k: n[0] * p[0] + n[1] * p[1] for k, p in pts.items()}
        top = max(vals.values())
        act = sorted(k for k, v in vals.items() if top - v <= mp.mpf(10) ** (30 - self.dps))
        return top, act, {k: dn[0] * pts[k][0] + dn[1] * pts[k][1] for k in act}

    def terms(self, ph):
        key = mp.nstr(ph, self.dps)
        if key not in self.memo:
            with mp.workdps(self.dps):
                self.memo[key] = {nm: self.support(nm, ph) for nm in FUNS}
        return self.memo[key]

    def s(self, ph):
        with mp.workdps(self.dps):
            return mp.fsum(v[0] for v in self.terms(ph).values())

    def dsides(self, ph):
        """the exact one-sided phase derivatives: the sum over the three maxima of the least (left) and the largest
        (right) derivative among their active points"""
        with mp.workdps(self.dps):
            T = self.terms(ph)
            return mp.fsum(min(v[2].values()) for v in T.values()), mp.fsum(max(v[2].values()) for v in T.values())


_cache = {}


def memo(tag, w, *extra):
    def deco(fn):
        key = (tag, ctx(w)) + extra
        if key not in _cache:
            _cache[key] = fn()
        return _cache[key]
    return deco


def cand(w, dps):
    return memo("cand", w, dps)(lambda: Cand(w, dps))


def s270(w, method="quad", dps=200):
    def go():
        C = cand(w, dps)
        with mp.workdps(dps):
            x1, y1 = C.pos(mp.mpf(1), method)
            return w["cusp_formula"](x1, y1, C.l0, C.A, C.B)
    return memo("s270", w, method, dps)(go)


def event_branches(w, dps=200):
    """the event five-link's scale at its four branch phases, from its six vertices"""
    def go():
        with mp.workdps(dps):
            E = w["summary"]["definitions"]["event_control"]
            l1, l2, b, d = (mp.mpf(E[k]) for k in ("l1", "l2", "beta", "delta"))
            al = b + d
            V = [(mp.mpf(0), mp.mpf(0))]
            for L, ang in ((l1, al), (l2, b), (1 - 2 * (l1 + l2), mp.mpf(0)), (l2, -b), (l1, -al)):
                V.append((V[-1][0] + L * mp.cos(ang), V[-1][1] + L * mp.sin(ang)))
            A, B = legs_mp(w)

            def s_at(ph):
                c, s = mp.cos(ph), mp.sin(ph)
                X = [x * c - y * s for x, y in V]
                Y = [x * s + y * c for x, y in V]
                return max(x / A + y / B for x, y in zip(X, Y)) - min(X) / A - min(Y) / B
            return [s_at(ph) for ph in (mp.pi / 2 - al, mp.pi / 2 - b, 2 * mp.pi / 3, 3 * mp.pi / 2)]
    return memo("event", w, dps)(go)


def darboux(w, N, dps=120):
    """the left and right sums of N parts: theta decreases, so sin(theta) decreases and cos(theta) increases, and
    each sum bounds its integral from one side"""
    def go():
        C = cand(w, dps)
        with mp.workdps(dps):
            th = [C.theta(mp.mpf(k) / N) for k in range(N + 1)]
            sv, cv = [mp.sin(t) for t in th], [mp.cos(t) for t in th]
            sinL, sinR = mp.fsum(sv[:-1]) / N, mp.fsum(sv[1:]) / N
            cosL, cosR = mp.fsum(cv[:-1]) / N, mp.fsum(cv[1:]) / N
            f = w["cusp_formula"]
            return f(C.w * cosL, C.w * sinR, C.l0, C.A, C.B), f(C.w * cosR, C.w * sinL, C.l0, C.A, C.B)
    return memo("darboux", w, N, dps)(go)


def audit_blocks(w):
    """the audit note's code blocks: s270, the event minimum, the difference, and the 4,096-part enclosure"""
    return re.findall(r"```text\n(.*?)\n```", w["audit"].replace("\r\n", "\n"), re.S)


# ---------------------------------------------------------------- the checks

def p1_definitions(w):
    failures = []
    D = w["summary"]["definitions"]
    B8 = w["round08_best"]
    want = {"epsilon": B8["epsilon"], "wing_length": B8["parameters"][0], "beta": B8["parameters"][1],
            "delta": B8["parameters"][2], "transition_center": B8["parameters"][3]}
    for k, v in want.items():
        if D["candidate"][k] != repr(float(v)):
            failures.append(f"candidate {k}: {D['candidate'][k]} is not Round 08's {v!r}")
    for k, src in (("l1", "l1"), ("l2", "l2"), ("beta", "beta_radians"), ("delta", "delta_radians")):
        if D["event_control"][k] != repr(float(w["round05"][src])):
            failures.append(f"event {k}: {D['event_control'][k]} is not Round 05's {w['round05'][src]!r}")
    V = w["verifier"]
    for name, k in (("eps", "epsilon"), ("w", "wing_length"), ("beta", "beta"), ("delta", "delta"), ("center", "transition_center")):
        m = re.search(rf'\b{name}=mp\.mpf\("([\d.]+)"\)', V)
        if not m or m.group(1) != D["candidate"][k]:
            failures.append(f"the verifier's {name} is not the definition's")
    for frag in ('(1-mp.tanh((u-center)/eps))', "(raw(u)-r1)/(r0-r1)", "w*iy/A+(w*ix+l0)/B", 'mp.mpf("1.0048")'):
        if frag not in V:
            failures.append(f"the verifier lacks {frag}")
    rp = gate.squash(w["report"])
    for tex, k in ((r"\\varepsilon=([\d.]+),", "epsilon"), (r"\$\$w=([\d.]+),", "wing_length"), (r"\\beta=([\d.]+),", "beta"),
                   (r"\\delta=([\d.]+),", "delta"), (r"\$\$c=([\d.]+)\.", "transition_center")):
        m = re.search(tex, rp)
        if not m or m.group(1) != D["candidate"][k]:
            failures.append(f"the report's {k} is not the definition's")
    # the g-form of the layer is RUN-009's F
    C = cand(w, 50)
    with mp.workdps(50):
        g = lambda u: (1 - mp.tanh((u - C.c) / C.eps)) / 2
        for u in (mp.mpf(0), mp.mpf("0.25"), C.c, mp.mpf("0.75"), mp.mpf(1)):
            Fg = (g(u) - g(1)) / (g(0) - g(1))
            F9 = (mp.tanh((C.c - u) / C.eps) - C.T1) / (C.T0 - C.T1)
            if abs(Fg - F9) > mp.mpf(10) ** -40:
                failures.append(f"the g-form is not RUN-009's F at u = {mp.nstr(u, 5)}")
    return failures


def p2_cusp(w):
    failures = []
    s1, s2 = s270(w, "quad"), s270(w, "quadgl")
    with mp.workdps(200):
        if abs(s1 - s2) > mp.mpf(10) ** -180:
            failures.append(f"the two quadratures differ by {mp.nstr(s1 - s2, 3)}")
    blocks = audit_blocks(w)
    for label, v in (("summary", w["summary"]["arbitrary_precision"]["candidate_270_scale"]),
                     ("manifest", w["manifest"]["candidate_270_scale"]), ("audit note", blocks[0] if blocks else None)):
        if v is None or not holds(v, s1):
            failures.append(f"the {label}'s 270-degree value does not hold to its digits")
    by = {}
    for r in w["convergence"]:
        with mp.workdps(200):
            if abs(mp.mpf(r["scale"]) - s1) > mp.mpf(10) ** (2 - int(r["dps"])):
                failures.append(f"convergence row {r['dps']} {r['method']} is off at its working precision")
        by.setdefault(r["dps"], set()).add(r["scale"])
    if any(len(v) != 1 for v in by.values()):
        failures.append("the two methods' rows are not identical")
    return failures


def p3_event(w):
    failures = []
    br = event_branches(w)
    S = w["summary"]["arbitrary_precision"]
    for i, (v, pub) in enumerate(zip(br, S["event_branch_values"])):
        if not holds(pub, v):
            failures.append(f"event branch B{i + 1} does not hold to its digits")
    with mp.workdps(200):
        m, d = min(br), s270(w) - min(br)
    blocks = audit_blocks(w)
    for label, v, true in (("summary minimum", S["event_minimum"], m), ("summary difference", S["candidate_minus_event"], d),
                           ("manifest minimum", w["manifest"]["event_minimum"], m),
                           ("manifest difference", w["manifest"]["candidate_minus_event"], d),
                           ("audit note minimum", blocks[1] if len(blocks) > 1 else None, m),
                           ("audit note difference", blocks[2] if len(blocks) > 2 else None, d)):
        # a difference of two 120-digit values is printed to 124 decimals: it is asked to hold to 120
        if v is None or not holds(v, true, 120 if "difference" in label else None):
            failures.append(f"the {label} does not hold to its digits")
    return failures


def p4_darboux(w):
    failures = []
    C = cand(w, 120)
    with mp.workdps(120):
        if not (C.d > 0 and C.T0 > C.T1 and 0 < C.theta(mp.mpf(1)) and C.theta(mp.mpf(0)) < mp.pi / 2):
            failures.append("theta is not decreasing inside (0, pi/2): the sums are not bounds")
        m = min(event_branches(w))
    for r in w["darboux"]:
        lo, hi = darboux(w, int(r["partitions"]))
        with mp.workdps(120):
            for key, true in (("lower", lo), ("upper", hi), ("width", hi - lo), ("lower_minus_event", lo - m)):
                # the sums are printed to about 105 decimals from about 100-digit arithmetic: 99 are asked for
                if not holds(r[key], true, 99):
                    failures.append(f"{r['partitions']} parts: {key} does not hold to its digits")
    blocks = audit_blocks(w)
    enc = dict(re.findall(r"^(lower|upper|lower - event) = ([\d.]+)$", blocks[3], re.M)) if len(blocks) > 3 else {}
    lo, hi = darboux(w, 4096)
    with mp.workdps(120):
        for key, true in (("lower", lo), ("upper", hi), ("lower - event", lo - m)):
            if key not in enc or not holds(enc[key], true, 99):
                failures.append(f"the audit note's {key} does not hold to its digits")
    return failures


def branch_records(w):
    BA = w["branch_audit"]
    return [(f"low {i + 1}", r) for i, r in enumerate(BA["low_branches"])] + [(r["name"], r) for r in BA["special_branches"]]


def phase_of(name, r):
    """a low branch's printed phase; a cusp's exact phase"""
    return {"120deg": 2 * mp.pi / 3, "270deg": 3 * mp.pi / 2}.get(name, None) or mp.mpf(r["phase"])


def p5_branches(w):
    failures = []
    C = cand(w, 100)
    s270_pub = w["summary"]["arbitrary_precision"]["candidate_270_scale"]
    for name, r in branch_records(w):
        with mp.workdps(100):
            ph = phase_of(name, r)
            if abs(C.s(ph) - mp.mpf(r["scale"])) > mp.mpf(10) ** -68:
                failures.append(f"{name}: the value is not the scale at its phase")
            T = C.terms(ph)
            for k, lab in r["labels"].items():
                if lab not in T[k][1]:
                    failures.append(f"{name}: {lab} is not an active support of {k}")
        with mp.workdps(160):
            if abs(mp.mpf(r["margin_above_270"]) - (mp.mpf(r["scale"]) - mp.mpf(s270_pub))) > mp.mpf(10) ** -66:
                failures.append(f"{name}: the margin is not the value minus s270")
        if name.startswith("low"):
            with mp.workdps(100):
                h = mp.mpf("1e-7")
                if not (C.dsides(ph - h)[1] < 0 < C.dsides(ph + h)[1]):
                    failures.append(f"{name}: no stationary minimum within 1e-7 rad")
        else:
            with mp.workdps(100):
                exact = 2 * mp.pi / 3 if name == "120deg" else 3 * mp.pi / 2
                if not holds(r["phase"], exact):
                    failures.append(f"{name}: the phase is not the cusp's")
    return failures


def p6_slopes(w):
    failures = []
    C = cand(w, 100)
    for name, r in branch_records(w):
        if "slopes" not in r:
            continue
        with mp.workdps(100):
            ph = phase_of(name, r)
            Dm, Dp = C.dsides(ph)
            if not Dm < 0 < Dp:
                failures.append(f"{name}: not D- < 0 < D+")
            s0 = C.s(ph)
            for h, sl in r["slopes"].items():
                H = mp.mpf(h)
                qr, ql = (C.s(ph + H) - s0) / H, (s0 - C.s(ph - H)) / H
                if abs(qr - mp.mpf(sl["right"])) > mp.mpf(10) ** -45 or abs(ql - mp.mpf(sl["left"])) > mp.mpf(10) ** -45:
                    failures.append(f"{name} h = {h}: the slopes are not the difference quotients")
                if abs(qr - Dp) > H or abs(ql - Dm) > H:
                    failures.append(f"{name} h = {h}: a quotient is not within h of the derivative")
    return failures


def peak_fits(w):
    pk = w["round08_peak"]
    x = np.array([float(k) for k in pk])
    y = np.array([R["best_scale"] for R in pk.values()])
    out = {}
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        for deg in (2, 3, 4):
            co = np.polyfit(x, y, deg)
            p = np.poly1d(co)
            roots = [z.real for z in p.deriv().roots if abs(z.imag) < 1e-12 and x.min() <= z.real <= x.max()]
            e = max(roots, key=p)
            out[deg] = {"coefficients": co, "eps": float(e), "scale": float(p(e)), "second": float(p.deriv(2)(e))}
    return out


def p7_fit(w):
    failures = []
    mine = peak_fits(w)
    for F in w["peak_fit"]:
        m = mine[F["degree"]]
        if np.abs(m["coefficients"] / np.array(F["coefficients"]) - 1).max() > 1e-6:
            failures.append(f"degree {F['degree']}: the coefficients are not the fit's")
        for key, mk, tol in (("estimated_epsilon", "eps", 1e-10), ("estimated_scale", "scale", 1e-14), ("second_derivative", "second", 1e-9)):
            if abs(F[key] - m[mk]) > tol:
                failures.append(f"degree {F['degree']}: {key} {F[key]!r}, the fit's {m[mk]!r}")
    return failures


STATUS = {"通過": "passed", "**未完成**": "not completed", "未開始": "not attempted"}


def p8_consistency(w):
    failures = []
    rp, S = w["report"].replace("\r\n", "\n"), w["summary"]
    rs = gate.squash(rp)
    # section 4: each row prints the start of its convergence row, and its seconds
    conv = {(r["dps"], r["method"]): r for r in w["convergence"]}
    rows = re.findall(r"^\|\s*(\d+)\s*\|\s*(quad|quadgl)\s*\|\s*`([\d.]+)…`\s*\|\s*([\d.]+)\s*\|", rp, re.M)
    if len(rows) != 8:
        failures.append(f"section 4: {len(rows)} rows")
    for dps, meth, val, sec in rows:
        r = conv.get((dps, meth))
        if not r or not r["scale"].startswith(val) or f"{float(r['elapsed_seconds']):.3f}" != sec:
            failures.append(f"section 4 row {dps} {meth} is not the data's")
    # section 5: the branch values' first digits
    rows = re.findall(r"^\|\s*B(\d)\s*\|\s*`([\d.]+)…`\s*\|", rp, re.M)
    if len(rows) != 4 or any(not S["arbitrary_precision"]["event_branch_values"][int(i) - 1].startswith(v) for i, v in rows):
        failures.append("section 5's table is not the data's")
    # section 6: the enclosure table
    dar = {r["partitions"]: r for r in w["darboux"]}
    rows = re.findall(r"^\|\s*(\d+)\s*\|\s*`([\d.]+)…`\s*\|\s*`([\d.]+)…`\s*\|\s*([\d.e-]+)\s*\|\s*([\d.e-]+)\s*\|", rp, re.M)
    if len(rows) != 4:
        failures.append(f"section 6: {len(rows)} rows")
    for n, lo, hi, wd, le in rows:
        r = dar.get(n)
        if not r or not r["lower"].startswith(lo) or not r["upper"].startswith(hi) \
                or not rounds_to(wd, float(r["width"])) or not rounds_to(le, float(r["lower_minus_event"])):
            failures.append(f"section 6 row {n} is not the data's")
    m = re.search(r"s_\{270\}\^\{\\mathrm\{lower\}\}-s_0>([\d.]+)\\times10\^\{-(\d+)\}", rs)
    if not m or not float(dar["4096"]["lower_minus_event"]) > float(m.group(1)) * 10 ** -int(m.group(2)):
        failures.append("section 6's bound is not the 4,096-part margin")
    # section 7: the branch table
    recs = dict(branch_records(w))
    names = {"low-1": "low 1", "low-2": "low 2", "120°": "120deg", "270°": "270deg"}
    rows = re.findall(r"^\|\s*(low-1|low-2|120°|270°)\s*\|\s*([^|]+?)\s*\|\s*`([\d.]+)…`\s*\|\s*([\d.e-]+)\s*\|", rp, re.M)
    if len(rows) != 4:
        failures.append(f"section 7: {len(rows)} rows")
    for nm, ph, val, mar in rows:
        r = recs[names[nm]]
        ok = r["scale"].startswith(val) and rounds_to(mar, float(r["margin_above_270"]) if nm != "270°" else 0.0)
        ok = ok and (holds(ph, r["phase"]) if nm.startswith("low") else ph == {"120°": r"$2\pi/3$", "270°": r"$3\pi/2$"}[nm])
        if not ok:
            failures.append(f"section 7 row {nm} is not the data's")
    m = re.search(r"s_\{120\}-s_\{270\}\\approx([\d.]+)\\times10\^\{-(\d+)\}", rs)
    if not m or not rounds_to(f"{m.group(1)}e-{m.group(2)}", float(recs["120deg"]["margin_above_270"])):
        failures.append("section 7's gap is not the 120-degree margin")
    # section 8: the h = 1e-9 quotients
    got = re.findall(r"D_([-+])s\\approx(-?[\d.]+),?", rs)
    want = [(sd, recs[nm]["slopes"]["1e-9"][k]) for nm in ("120deg", "270deg") for sd, k in (("-", "left"), ("+", "right"))]
    if [(a, b.rstrip(".")) for a, b in got] != want:
        failures.append("section 8's slopes are not the data's h = 1e-9 quotients")
    # section 9: the fits
    fits = {str(F["degree"]): F for F in w["peak_fit"]}
    rows = re.findall(r"^\|\s*([234])\s*\|\s*([\d.]+)\s*\|\s*([\d.]+)\s*\|\s*(-?[\d.e+-]+)\s*\|", rp, re.M)
    if len(rows) != 3:
        failures.append(f"section 9: {len(rows)} rows")
    for d, e, s, sec in rows:
        F = fits[d]
        if not (rounds_to(e, F["estimated_epsilon"]) and rounds_to(s, F["estimated_scale"]) and rounds_to(sec, F["second_derivative"])):
            failures.append(f"section 9 row {d} is not the fit's")
    m = re.search(r"([\d.]+)\\lesssim\\varepsilon_\*\\lesssim(\d+\.\d+)", rs)
    eps = [F["estimated_epsilon"] for F in w["peak_fit"]]
    if not m or not (rounds_to(m.group(1), min(eps)) or float(m.group(1)) <= min(eps)) \
            or not (rounds_to(m.group(2), max(eps)) or float(m.group(2)) >= max(eps)):
        failures.append("section 9's range is not the fits'")
    # section 10 and the README's lists: the certificate statuses
    rows = re.findall(r"^\|\s*([^|]+?)\s*\|\s*(通過|\*\*未完成\*\*|未開始)\s*\|", rp, re.M)
    counts = {v: sum(1 for _, s in rows if STATUS[s] == v) for v in STATUS.values()}
    cs = list(S["certificate_status"].values())
    if counts["not completed"] != cs.count("not completed") or counts["not attempted"] != cs.count("not attempted") \
            or counts["passed"] < cs.count("passed"):
        failures.append("section 10's statuses are not the summary's")
    rd = "\n".join(l.rstrip() for l in w["readme"].replace("\r\n", "\n").split("\n"))
    lists = re.findall(r"## (已通過|尚未通過)\n\n((?:- .*\n)+)", rd)
    got = {k: len(v.strip().split("\n")) for k, v in lists}
    if got.get("已通過") != counts["passed"] or got.get("尚未通過") != counts["not completed"] + counts["not attempted"]:
        failures.append("the README's lists are not section 10's")
    m = re.search(r"4096 分割 Darboux 下界減事件控制：`([\d.]+)`", rd)
    if not m or not holds(m.group(1), dar["4096"]["lower_minus_event"]):
        failures.append("the README's Darboux margin is not the data's")
    # the draft's target and the manifest
    if "\\eta<s_{270}-s_0" not in gate.squash(w["draft"]) or "目前尚未完成" not in w["draft"]:
        failures.append("the draft does not state its target as open")
    M = w["manifest"]
    if M["version"] != "0.9" or M["status"] != S["status"] or M["date"] != "2026-07-26" or "2026年7月26日" not in rs:
        failures.append("the manifest's version, status or date is not the package's")
    if sorted(M["files"]) != sorted(set(w["members"]) - {"MANIFEST.json"}):
        failures.append("the manifest's file list is not the package's")
    # the summary's copies of the data files
    if S["branch_audit"] != w["branch_audit"] or S["epsilon_peak_fit"] != w["peak_fit"]:
        failures.append("the summary's branch audit or fit is not the data file")
    norm = lambda rows: [{k: str(v) for k, v in r.items()} for r in rows]
    if norm(S["darboux_cusp_enclosure"]) != norm(w["darboux"]) or norm(S["precision_convergence"]) != norm(w["convergence"]):
        failures.append("the summary's enclosure or convergence rows are not the data files")
    return failures


def long_figures(w):
    """the README's and the report's long figures: (label, printed, value)"""
    s = s270(w)
    br = event_branches(w)
    with mp.workdps(200):
        m = min(br)
        d = s - m
    rd, rs = w["readme"], gate.squash(w["report"])
    one = lambda pat, text: (re.search(pat, text) or [None, None])[1]
    v270 = re.findall(r"\\boxed\{s_\{270\}=([\d.]+)\}", rs)
    v0 = re.findall(r"\\boxed\{s_0=([\d.]+)\}", rs)
    vd = re.findall(r"\\boxed\{s_\{270\}-s_0=([\d.]+)\}", rs)
    two = lambda L, i: L[i] if len(L) == 2 else None
    return [("the README's 270-degree value", one(r"平滑候選 270 度：`([\d.]+)`", rd), s),
            ("the README's event minimum", one(r"事件控制最低值：`([\d.]+)`", rd), m),
            ("the README's difference", one(r"高精度差值：`([\d.]+)`", rd), d),
            ("the report's section 4 value", two(v270, 0), s),
            ("the report's section 5 event minimum", two(v0, 0), m),
            ("the report's section 5 difference", two(vd, 0), d),
            ("the report's section 13 value", two(v270, 1), s),
            ("the report's section 13 event minimum", two(v0, 1), m),
            ("the report's section 13 difference", two(vd, 1), d)]


def p9_digits(w):
    failures = []
    for label, printed, true in long_figures(w):
        if printed is None or not holds(printed, true):
            failures.append(f"{label} does not hold to its digits")
    C = cand(w, 100)
    with mp.workdps(100):
        if any(not holds(r["scale"], C.s(phase_of(name, r))) for name, r in branch_records(w)):
            failures.append("the branch audit's values do not hold to their printed digits")
        got = re.findall(r"D_[-+]s\\approx(-?[\d.]+)", gate.squash(w["report"]))
        exact = [x for ph in (2 * mp.pi / 3, 3 * mp.pi / 2) for x in C.dsides(ph)]
        if len(got) != 4 or any(not holds(g.rstrip("."), e) for g, e in zip(got, exact)):
            failures.append("the report's section 8 derivatives do not hold to their printed digits")
    return failures


CHECKS = {"P1": p1_definitions, "P2": p2_cusp, "P3": p3_event, "P4": p4_darboux, "P5": p5_branches, "P6": p6_slopes,
          "P7": p7_fit, "P8": p8_consistency, "P9": p9_digits}
HP = {"summary", "legs", "cusp_formula"}
DEPENDS = {
    "P1": {"summary", "round08_best", "round05", "verifier", "report"},
    "P2": HP | {"manifest", "audit", "convergence"},
    "P3": HP | {"manifest", "audit"},
    "P4": HP | {"darboux", "audit"},
    "P5": HP | {"branch_audit"},
    "P6": HP | {"branch_audit"},
    "P7": {"peak_fit", "round08_peak"},
    "P8": {"readme", "report", "draft", "manifest", "members", "summary", "branch_audit", "darboux", "convergence", "peak_fit"},
    "P9": HP | {"readme", "report", "branch_audit"},
}


# ---------------------------------------------------------------- the interval certificate

def dec(x, n, up=False):
    """x to n decimals, cut toward the safe side: down for a lower bound, up for an upper one"""
    with mp.workdps(n + 60):
        q = int(mp.ceil(x * mp.mpf(10) ** n) if up else mp.floor(x * mp.mpf(10) ** n))
    s = str(abs(q)).rjust(n + 1, "0")
    return ("-" if q < 0 else "") + s[:-n] + "." + s[-n:]


def lo(x):
    """an interval's lower endpoint as a plain number (exact: 40 digits hold its 30-digit endpoint)"""
    with mp.workdps(40):
        return mp.mpf(x._mpi_[0])


def hi(x):
    with mp.workdps(40):
        return mp.mpf(x._mpi_[1])


def tanh_iv(x):
    return 1 - 2 / (iv.exp(2 * x) + 1)


def certificate_boxes(w, K):
    """interval boxes for the candidate's points at u = j/K on both wings, in mpmath's interval arithmetic: each
    step's integral of cos(theta) and of sin(theta) is its trapezoid value within h^3 M2 / 12, where M2 bounds the
    second derivative of both - |f''| <= theta'^2 + |theta''|, with |F'| <= 1/(eps (T0 - T1)) and
    |F''| <= (4/(3 sqrt 3)) / (eps^2 (T0 - T1)), sech^2 z |tanh z| being at most 2/(3 sqrt 3)"""
    D = w["summary"]["definitions"]["candidate"]
    eps, wl, b, d, c = (iv.mpf(D[k]) for k in ("epsilon", "wing_length", "beta", "delta", "transition_center"))
    T0, T1 = tanh_iv(c / eps), tanh_iv((c - 1) / eps)
    den = T0 - T1
    M2 = iv.mpf(hi(d * (4 / (3 * iv.sqrt(3))) / (eps ** 2 * den) + (d / (eps * den)) ** 2))
    l0 = 1 - 2 * wl
    cs, sn = [], []
    for j in range(K + 1):
        th = b + d * (tanh_iv((c - iv.mpf(j) / K) / eps) - T1) / den
        cs.append(iv.cos(th))
        sn.append(iv.sin(th))
    h = iv.mpf(1) / K
    step_err = iv.mpf(hi(h ** 3 * M2 / 12))
    X, Y = [iv.mpf(0)], [iv.mpf(0)]
    sx, sy = iv.mpf(0), iv.mpf(0)
    for j in range(K):
        sx += h * (cs[j] + cs[j + 1]) / 2
        sy += h * (sn[j] + sn[j + 1]) / 2
        e = hi(step_err * (j + 1))
        X.append(wl * (sx + iv.mpf([-e, e])))
        Y.append(wl * (sy + iv.mpf([-e, e])))
    x1, y1 = X[-1], Y[-1]
    # the right wing mirrors the left about x = x1 + l0/2; all points are moved by a fixed (exact) centre, which
    # leaves the scale unchanged (the three functionals sum to zero)
    tx, ty = iv.mpf(float((x1 + l0 / 2).mid)), iv.mpf(float((y1 / 2).mid))
    PX = [v - tx for v in X] + [2 * x1 + l0 - v - tx for v in X[::-1]]
    PY = [v - ty for v in Y] + [v - ty for v in Y[::-1]]
    return {"PX": PX, "PY": PY, "x1": x1, "y1": y1, "l0": l0, "M2": hi(M2),
            "box_x": max(float(v.delta) for v in X), "box_y": max(float(v.delta) for v in Y)}


def legs_iv(w):
    iv.dps = 30
    r3 = iv.sqrt(3)
    fa, fb = (iv.mpf(x) for x in w["legs"])
    return (3 + 4 * r3) / 18 / iv.mpf("1.0048") * fa, (4 + r3) / 6 / iv.mpf("1.0048") * fb


def event_least_upper(w):
    """the least of the event five-link's scales at its four branch phases (enclosing intervals), from its vertices
    in interval arithmetic: an upper bound for its least scale"""
    A, B = legs_iv(w)
    funs = [(1 / A, 1 / B), (-1 / A, iv.mpf(0)), (iv.mpf(0), -1 / B)]
    E = w["summary"]["definitions"]["event_control"]
    l1, l2, be, de = (iv.mpf(E[k]) for k in ("l1", "l2", "beta", "delta"))
    al = be + de
    V = [(iv.mpf(0), iv.mpf(0))]
    for L, ang in ((l1, al), (l2, be), (1 - 2 * (l1 + l2), iv.mpf(0)), (l2, -be), (l1, -al)):
        V.append((V[-1][0] + L * iv.cos(ang), V[-1][1] + L * iv.sin(ang)))

    def s_poly_upper(ph):
        cph, sph = iv.cos(ph), iv.sin(ph)
        tot = iv.mpf(0)
        for a in funs:
            n = (a[0] * cph + a[1] * sph, -a[0] * sph + a[1] * cph)
            tot += iv.mpf(max(hi(n[0] * x + n[1] * y) for x, y in V))
        return hi(tot)
    return min(s_poly_upper(ph) for ph in (iv.pi / 2 - al, iv.pi / 2 - be, 2 * iv.pi / 3, 3 * iv.pi / 2))


def certificate(w, K=65536, target=1.05e-5):
    """a lower bound for the candidate's least scale over every phase, against an upper bound for the event
    five-link's: for each phase interval [a, b] with midpoint m, s(phi) >= sum over the three maxima of the
    functional at m on one chosen point's box (interval arithmetic) - Lip |phi - m|, where Lip bounds the phase
    derivative of every functional on every point. A floating-point search bisects the phases until each
    interval's estimate clears the target; each interval's bound is then taken in interval arithmetic"""
    A, B = legs_iv(w)
    funs = [(1 / A, 1 / B), (-1 / A, iv.mpf(0)), (iv.mpf(0), -1 / B)]
    ev_up = event_least_upper(w)
    bx = certificate_boxes(w, K)
    PX, PY = bx["PX"], bx["PY"]
    n_pts = len(PX)
    cx = np.array([float(v.mid) for v in PX])
    cy = np.array([float(v.mid) for v in PY])
    # the radius and the Lipschitz constant, upper bounds in interval arithmetic: every point's box lies within
    # max|x| and max|y| of the centre
    with mp.workdps(40):
        mxe = max(max(abs(mp.mpf(v._mpi_[0])), abs(mp.mpf(v._mpi_[1]))) for v in PX)
        mye = max(max(abs(mp.mpf(v._mpi_[0])), abs(mp.mpf(v._mpi_[1]))) for v in PY)
    R_up = hi(iv.sqrt(iv.mpf(mxe) ** 2 + iv.mpf(mye) ** 2))
    Lip_up = hi(iv.mpf(R_up) * sum((iv.mpf(hi(iv.sqrt(a[0] ** 2 + a[1] ** 2))) for a in funs), iv.mpf(0)))
    anorm = [float(hi(iv.sqrt(a[0] ** 2 + a[1] ** 2))) for a in funs]
    Lip = float(Lip_up) * (1 + 1e-12)            # the float search's copy; the bounds use Lip_up
    box_x = max(float(v.delta) for v in PX)
    box_y = max(float(v.delta) for v in PY)
    af = [(float(a[0].mid), float(a[1].mid)) for a in funs]
    neg_tang = -np.unwrap(np.arctan2(np.diff(cy), np.diff(cx)))

    def best_idx(ph):
        """for each functional the point of largest float value among the ends and the points beside the one place
        the tangent turns past perpendicular to it (the curve turns clockwise through less than 180 degrees)"""
        c, s = math.cos(ph), math.sin(ph)
        out = []
        for a0, a1 in af:
            n0, n1 = a0 * c + a1 * s, -a0 * s + a1 * c
            psi = math.atan2(n1, n0)
            cand_i = [0, n_pts - 1]
            for bnd in (psi - math.pi / 2, psi + math.pi / 2):
                k = int(np.searchsorted(neg_tang, -math.remainder(bnd, 2 * math.pi)))
                cand_i += [i for i in range(k - 2, k + 3) if 0 <= i < n_pts]
            ci = np.array(cand_i)
            vals = cx[ci] * n0 + cy[ci] * n1
            out.append((int(ci[int(np.argmax(vals))]), float(vals.max())))
        return out

    box_pen = sum(anorm) * math.hypot(box_x, box_y)
    ev_f = float(ev_up)

    def cover(tgt, start=0.0, end=6.2832, parts=4096):
        """the search, in floating point: bisect [start, end] until every interval's float estimate clears the
        event plus the target; None when an interval narrower than 1e-9 cannot. Its leaves are then bounded in
        interval arithmetic"""
        todo = [(start + k * (end - start) / parts, start + (k + 1) * (end - start) / parts) for k in range(parts)][::-1]
        leaves = []
        while todo:
            a_, b_ = todo.pop()
            m = (a_ + b_) / 2
            hw = (b_ - a_) / 2
            if sum(v for _, v in best_idx(m)) - Lip * hw - box_pen - 1e-12 > ev_f + tgt:
                leaves.append((a_, b_, m))
            elif hw < 1e-9:
                return None, m
            else:
                todo += [(m, b_), (a_, m)]
        return leaves, None

    leaves, _ = cover(target)
    # the rigorous pass: every leaf in interval arithmetic, with its chosen points
    LipI = iv.mpf(Lip_up)
    worst = None
    for a_, b_, m in leaves:
        M = iv.mpf(m)
        hw = max((M - iv.mpf(a_)).b, (iv.mpf(b_) - M).b)
        cph, sph = iv.cos(M), iv.sin(M)
        tot = iv.mpf(0)
        for (i, _), a in zip(best_idx(m), funs):
            tot += (a[0] * cph + a[1] * sph) * PX[i] + (-a[0] * sph + a[1] * cph) * PY[i]
        lb = lo(tot - LipI * iv.mpf(hw))
        if worst is None or lb < worst[0]:
            worst = (lb, m)
    # an upper bound for the candidate's least: its scale at 3pi/2, from the end boxes
    up270 = hi(cusp_270(bx["x1"], bx["y1"], bx["l0"], A, B))
    # the refusing control, beside 3pi/2: a target above the candidate's own certified value there must not be met
    gain_up = float(up270 - ev_up)
    refused, near = cover(gain_up + 1e-9, 3 * math.pi / 2 - 0.01, 3 * math.pi / 2 + 0.01, 64)
    sl = sorted(leaves)
    covered = sl[0][0] == 0 and sl[-1][1] == 4096 * 6.2832 / 4096 > 2 * math.pi and all(
        x[1] == y[0] for x, y in zip(sl, sl[1:]))
    mp.mp.dps = 40
    out = {"wing_steps": K, "interval_digits": 30, "second_derivative_bound": dec(bx["M2"], 8, up=True),
            "largest_box_left_wing": [bx["box_x"], bx["box_y"]], "largest_box": [box_x, box_y],
            "radius_at_most": dec(R_up, 8, up=True), "lipschitz_at_most": dec(Lip_up, 8, up=True),
            "phase_intervals": len(leaves), "intervals_cover_0_to_6.2832": covered,
            "certified": bool(covered and worst[0] > ev_up),
            "candidate_least_at_least": dec(worst[0], 20), "worst_interval_midpoint": worst[1],
            "event_least_at_most": dec(ev_up, 30, up=True),
            "certified_margin_at_least": dec(worst[0] - ev_up, 25),
            "candidate_least_at_most": dec(up270, 20, up=True),
            "control_beside_3pi_2_target_above_s270_refused": refused is None, "control_refused_near": near}
    mp.mp.dps = 15
    return out


def darboux_interval(w, N=4096):
    """Round 09's 4,096-part lower bound, in interval arithmetic"""
    iv.dps = 30
    D = w["summary"]["definitions"]["candidate"]
    eps, wl, b, d, c = (iv.mpf(D[k]) for k in ("epsilon", "wing_length", "beta", "delta", "transition_center"))
    T0, T1 = tanh_iv(c / eps), tanh_iv((c - 1) / eps)
    r3 = iv.sqrt(3)
    fa, fb = (iv.mpf(x) for x in w["legs"])
    A = (3 + 4 * r3) / 18 / iv.mpf("1.0048") * fa
    B = (4 + r3) / 6 / iv.mpf("1.0048") * fb
    th = [b + d * (tanh_iv((c - iv.mpf(j) / N) / eps) - T1) / (T0 - T1) for j in range(N + 1)]
    sinR = sum((iv.sin(t) for t in th[1:]), iv.mpf(0)) / N
    cosL = sum((iv.cos(t) for t in th[:-1]), iv.mpf(0)) / N
    return lo(cusp_270(wl * cosL, wl * sinR, 1 - 2 * wl, A, B))


# ---------------------------------------------------------------- findings

def findings(w):
    s = s270(w)
    br = event_branches(w)
    C = cand(w, 100)
    out = {}
    with mp.workdps(200):
        m = min(br)
        out["s270_120_decimals_cut"] = dec(s, 120)
        defs = w["summary"]["definitions"]
        out["decimal_minus_double_constants_largest"] = mp.nstr(max(abs(mp.mpf(v) - mp.mpf(float(v))) for part in defs.values()
                                                                    for v in part.values()), 3)
        out["quadratures_differ_by"] = mp.nstr(abs(s270(w, "quadgl") - s), 3)
        out["event_branches_minus_minimum"] = [mp.nstr(v - m, 5) for v in br]
        out["difference_60_digits"] = mp.nstr(s - m, 60)
    out["long_figures_digits_held"] = {label: (digits_held(p, t), len(p.split(".")[1])) for label, p, t in long_figures(w)}
    out["convergence_rows_digits_held"] = {f"{r['dps']} {r['method']}": digits_held(r["scale"], s) for r in w["convergence"]}
    rows = {}
    with mp.workdps(100):
        for name, r in branch_records(w):
            ph = phase_of(name, r)
            v = C.s(ph)
            rec = {"active": {k: C.terms(ph)[k][1] for k in FUNS}, "labels": r["labels"],
                   "digits_held": (digits_held(r["scale"], v), len(r["scale"].split(".")[1]))}
            if name.startswith("low"):
                f = lambda q: C.dsides(q)[1]
                q = mp.findroot(f, ph)
                rec.update({"derivative_at_printed_phase": mp.nstr(f(ph), 5),
                            "stationary_phase": mp.nstr(q, 30), "printed_minus_stationary": mp.nstr(ph - q, 5),
                            "minimum": mp.nstr(C.s(q), 40), "printed_value_minus_minimum": mp.nstr(mp.mpf(r["scale"]) - C.s(q), 5)})
            else:
                Dm, Dp = C.dsides(ph)
                rec.update({"exact_left": mp.nstr(Dm, 30), "exact_right": mp.nstr(Dp, 30),
                            "quotient_minus_exact": {h: [mp.nstr(mp.mpf(sl["left"]) - Dm, 5), mp.nstr(mp.mpf(sl["right"]) - Dp, 5)]
                                                     for h, sl in r["slopes"].items()}})
            rows[name] = rec
    out["branches"] = rows
    fits = peak_fits(w)
    out["peak_fit_minus_stated"] = {str(F["degree"]): [fits[F["degree"]]["eps"] - F["estimated_epsilon"],
                                                         fits[F["degree"]]["scale"] - F["estimated_scale"],
                                                         fits[F["degree"]]["second"] - F["second_derivative"]] for F in w["peak_fit"]}
    with mp.workdps(40):
        lo = darboux_interval(w)
        ev = min(br)
        out["darboux_4096_interval_lower"] = dec(lo, 25)
        out["darboux_4096_interval_lower_minus_event_upper"] = dec(lo - event_least_upper(w), 25)
        out["darboux_4096_interval_lower_minus_stated"] = mp.nstr(lo - mp.mpf(next(r["lower"] for r in w["darboux"] if r["partitions"] == "4096")), 5)
    # what the long data figures hold, and each tolerance's margin on the package
    with mp.workdps(200):
        d = s - m
    S = w["summary"]["arbitrary_precision"]
    blocks = audit_blocks(w)
    out["differences_digits_held"] = {label: (digits_held(v, d), len(v.split(".")[1])) for label, v in
                                      (("summary", S["candidate_minus_event"]), ("manifest", w["manifest"]["candidate_minus_event"]),
                                       ("audit note", blocks[2]))}
    out["convergence_rows_printed_decimals"] = {f"{r['dps']} {r['method']}": len(r["scale"].split(".")[1]) for r in w["convergence"]}
    dh = {}
    for r in w["darboux"]:
        lo_, hi_ = darboux(w, int(r["partitions"]))
        with mp.workdps(120):
            dh[r["partitions"]] = {key: (digits_held(r[key], v), len(r[key].split(".")[1])) for key, v in
                                   (("lower", lo_), ("upper", hi_), ("width", hi_ - lo_), ("lower_minus_event", lo_ - m))}
    out["darboux_rows_digits_held"] = dh
    with mp.workdps(100):
        val_err = max(abs(C.s(phase_of(name, r)) - mp.mpf(r["scale"])) for name, r in branch_records(w))
        q_err, q_ratio = mp.mpf(0), mp.mpf(0)
        for name, r in branch_records(w):
            if "slopes" not in r:
                continue
            ph = phase_of(name, r)
            Dm, Dp = C.dsides(ph)
            for h, sl in r["slopes"].items():
                H = mp.mpf(h)
                qr, ql = (C.s(ph + H) - C.s(ph)) / H, (C.s(ph) - C.s(ph - H)) / H
                q_err = max(q_err, abs(qr - mp.mpf(sl["right"])), abs(ql - mp.mpf(sl["left"])))
                q_ratio = max(q_ratio, abs(qr - Dp) / H, abs(ql - Dm) / H)
        out["check_margins"] = {
            "P3_difference_decimals_over_120": min(v[0] for v in out["differences_digits_held"].values()) - 120,
            "P4_decimals_over_99": min(v[0] for row in dh.values() for v in row.values()) - 99,
            "P5_value_tolerance_over_error": mp.nstr(mp.mpf(10) ** -68 / val_err, 3),
            "P6_quotient_tolerance_over_error": mp.nstr(mp.mpf(10) ** -45 / q_err, 3) if q_err else None,
            "P6_largest_distance_over_h": mp.nstr(q_ratio, 6),
            "P7_tolerance_over_deviation": [1e-10 / max(abs(v[0]) for v in out["peak_fit_minus_stated"].values()),
                                            1e-14 / max(abs(v[1]) for v in out["peak_fit_minus_stated"].values()),
                                            1e-9 / max(abs(v[2]) for v in out["peak_fit_minus_stated"].values())]}
    out["certificate"] = certificate(w)
    return out


# ---------------------------------------------------------------- the drill

def bump(s, pos):
    """the decimal string with its pos-th decimal changed"""
    i = s.index(".") + pos
    return s[:i] + str((int(s[i]) + 1) % 10) + s[i + 1:]


def drill_cases(base):
    def edit(key, fn):
        v = copy.deepcopy(base[key])
        fn(v)
        return v
    S = base["summary"]
    b8 = edit("round08_best", lambda b: b["parameters"].__setitem__(3, b["parameters"][3] + 1e-16))
    ver = base["verifier"].replace('eps=mp.mpf("0.037")', 'eps=mp.mpf("0.038")')
    s270p = edit("summary", lambda s: s["arbitrary_precision"].__setitem__("candidate_270_scale", bump(S["arbitrary_precision"]["candidate_270_scale"], 100)))
    conv = edit("convergence", lambda rows: rows[2].__setitem__("scale", bump(rows[2]["scale"], 50)))
    s_conv = edit("summary", lambda s: s["precision_convergence"][2].__setitem__("scale", conv[2]["scale"]))
    evp = edit("summary", lambda s: s["arbitrary_precision"]["event_branch_values"].__setitem__(2, bump(S["arbitrary_precision"]["event_branch_values"][2], 90)))
    dar = edit("darboux", lambda rows: rows[2].__setitem__("upper", bump(rows[2]["upper"], 80)))
    s_dar = edit("summary", lambda s: s["darboux_cusp_enclosure"][2].__setitem__("upper", dar[2]["upper"]))
    # branch-audit plants: the file and the summary's copy alike
    def ba(fn):
        return edit("branch_audit", fn), edit("summary", lambda s: fn(s["branch_audit"]))
    mar = ba(lambda b: b["low_branches"][0].__setitem__("margin_above_270", bump(b["low_branches"][0]["margin_above_270"], 20)))
    lab = ba(lambda b: b["special_branches"][0]["labels"].__setitem__("min_x", "p0"))
    slo = ba(lambda b: b["special_branches"][1]["slopes"]["1e-7"].__setitem__("right", bump(b["special_branches"][1]["slopes"]["1e-7"]["right"], 30)))
    coef = lambda f: f[0]["coefficients"].__setitem__(0, f[0]["coefficients"][0] * (1 + 1e-5))
    fit, s_fit = edit("peak_fit", coef), edit("summary", lambda s: coef(s["epsilon_peak_fit"]))
    pk8 = edit("round08_peak", lambda p: p["0.039"].__setitem__("best_scale", p["0.039"]["best_scale"] + 1e-12))
    rep7 = base["report"].replace("| 4.5788994e-8 |", "| 4.5788984e-8 |")
    lme = re.search(r"4096 分割 Darboux 下界減事件控制：`([\d.]+)`", base["readme"]).group(1)
    readme = base["readme"].replace(lme, bump(lme, 20))
    v0 = re.findall(r"s_0\s*=\s*(0\.9989037571325086749551298249005564802172682267121921037894386799512563)\b", base["report"])
    i13 = base["report"].rindex(v0[-1]) if v0 else 0
    rep13 = base["report"][:i13] + bump(v0[-1], 60) + base["report"][i13 + len(v0[-1]):] if v0 else base["report"]
    man = edit("manifest", lambda m: m.__setitem__("candidate_minus_event", bump(m["candidate_minus_event"], 50)))
    s_copy = edit("summary", lambda s: s["epsilon_peak_fit"][1].__setitem__("estimated_scale", s["epsilon_peak_fit"][1]["estimated_scale"] + 1e-9))
    d8 = "0.43268675794937853584820971062851505492461178023733"
    rep8 = base["report"].replace(d8, bump(d8, 40))
    for a, b in ((ver, base["verifier"]), (rep7, base["report"]), (readme, base["readme"]), (rep13, base["report"]), (rep8, base["report"])):
        assert a != b, "a plant must hit the raw text"
    return [
        ("K0-unmodified", "control", {}, set()),
        ("K1-report-crlf", "control", {"report": base["report"].replace("\n", "\r\n")}, set()),
        ("K2-readme-trailing-spaces", "control", {"readme": "\n".join(l + " " for l in base["readme"].split("\n"))}, set()),
        ("K3-branch-audit-copied", "control", {"branch_audit": copy.deepcopy(base["branch_audit"])}, set()),
        ("M1-triangle-leg", "mutation", {"legs": ("1.0005", "1")}, {"P2", "P3", "P4", "P5", "P6", "P9"}),
        ("M2-round08-centre", "mutation", {"round08_best": b8}, {"P1"}),
        ("M3-verifier-constant", "mutation", {"verifier": ver}, {"P1"}),
        ("M4-summary-s270-digit", "mutation", {"summary": s270p}, {"P2"}),
        ("M5-convergence-row", "mutation", {"convergence": conv, "summary": s_conv}, {"P2"}),
        ("M6-event-branch-digit", "mutation", {"summary": evp}, {"P3"}),
        ("M7-darboux-upper", "mutation", {"darboux": dar, "summary": s_dar}, {"P4"}),
        ("M8-branch-margin", "mutation", {"branch_audit": mar[0], "summary": mar[1]}, {"P5"}),
        ("M9-branch-label", "mutation", {"branch_audit": lab[0], "summary": lab[1]}, {"P5"}),
        ("M10-cusp-slope", "mutation", {"branch_audit": slo[0], "summary": slo[1]}, {"P6"}),
        ("M11-fit-coefficient", "mutation", {"peak_fit": fit, "summary": s_fit}, {"P7"}),
        ("M12-round08-peak-scale", "mutation", {"round08_peak": pk8}, {"P7"}),
        ("M13-report-margin", "mutation", {"report": rep7}, {"P8"}),
        ("M14-readme-darboux", "mutation", {"readme": readme}, {"P8"}),
        ("M15-section13-digit", "mutation", {"report": rep13}, {"P9"}),
        ("M16-formula-term", "mutation", {"cusp_formula": lambda x1, y1, l0, A, B: y1 / A + (x1 - l0) / B}, {"P2", "P3", "P4", "P9"}),
        ("M17-manifest-difference", "mutation", {"manifest": man}, {"P3"}),
        ("M18-summary-copy", "mutation", {"summary": s_copy}, {"P8"}),
        ("M19-report-derivative", "mutation", {"report": rep8}, {"P8"}),
    ]


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--drill", action="store_true")
    ap.add_argument("--blind", metavar="CID")
    ap.add_argument("--log", default=None)
    args = ap.parse_args()
    if args.blind:
        CHECKS[args.blind] = lambda w: []
        args.drill = True
    w = build_world()
    results = {cid: fn(w) for cid, fn in CHECKS.items()}
    for cid, f in results.items():
        print(f"{cid} {'PASS' if not f else 'FAIL'}  ({len(f)} failures)")
        for x in f[:8]:
            print("    ", x)
    log = {"run": "RUN-010",
           "subject": "Round 09 (v0.9): the 270-degree cusp in arbitrary precision, the monotone Darboux enclosure, the branch audit",
           "package": {"id": RID, "sha256": sha256(source_bytes(RID))},
           "checks": {cid: {"pass": not f, "failures": f} for cid, f in results.items()},
           "findings": None if args.blind else findings(w)}
    if log["findings"]:
        print(json.dumps(log["findings"], indent=1, ensure_ascii=False)[:8000])
    if args.drill:
        drill = gate.run_drill(build_world, CHECKS, DEPENDS, drill_cases(build_world()))
        log["drill"] = gate.summarize(drill, CHECKS)
        print(gate.report_drill(log["drill"]))
        for d in drill:
            if not d["ok"]:
                print("   DRILL FAIL", d)
    if args.blind:
        seen = sum(1 for d in log["drill"]["cases"] if args.blind in d["expected"] and args.blind not in d["tripped"])
        print(f"blind {args.blind}: reported blind in {seen} cases (log not written)")
        return
    if args.drill:
        path = pathlib.Path(args.log) if args.log else LOG
        path.parent.mkdir(parents=True, exist_ok=True)
        path.write_text(json.dumps(log, ensure_ascii=False, indent=1) + "\n", encoding="utf-8", newline="\n")


if __name__ == "__main__":
    main()