Fix: Stop lens-map vertex decoding on read failure

Generate plain long-double and experimental FMA variants from the production entry kernel. Compare exact-rational references, geometric validation, fallback counts and repeated timings with self-contained fixtures.

Validate cached build dependencies and reject stale timing output. Count all unconfirmed entry outcomes independently of reference classification.

Makefile

@@ -108,6 +108,8 @@ TEST_OUT_DIR := $(OBJECT_DIR)/$(HDR_BUILD_TAG)
 TEST_TARGET := $(TEST_OUT_DIR)/test_geodesic
 ADAPTIVE_GEODESIC_TEST_TARGET := $(TEST_OUT_DIR)/test_geodesic_adaptive
 ASYMPTOTIC_TEST_TARGET := $(TEST_OUT_DIR)/test_asymptotic
+ASYMPTOTIC_ENTRY_TEST_TARGET := $(TEST_OUT_DIR)/test_asymptotic_entry
+ASYMPTOTIC_QUADRATIC_TEST_TARGET := $(TEST_OUT_DIR)/test_asymptotic_quadratic
 ASYMPTOTIC_SCHWARZSCHILD_TEST_TARGET := $(TEST_OUT_DIR)/test_asymptotic_schwarzschild
 TERMINATION_ORACLE_TEST_TARGET := $(TEST_OUT_DIR)/test_termination_oracle
 FRAME_TEST_TARGET := $(TEST_OUT_DIR)/test_frame
@@ -217,9 +219,20 @@ $(ADAPTIVE_GEODESIC_TEST_TARGET): tests/test_geodesic_adaptive.c $(COMMON_SOURCE
 $(ASYMPTOTIC_TEST_TARGET): tests/test_asymptotic.c $(CORE_MINKOWSKI_SOURCES) $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR)
 	$(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@
 
+# The test includes asymptotic.c to cover its private floating-point kernel.
+$(ASYMPTOTIC_QUADRATIC_TEST_TARGET): tests/test_asymptotic_quadratic.c $(CORE_MINKOWSKI_SOURCES) $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR)
+	$(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $(filter-out src/asymptotic.c,$^) $(LDLIBS) -o $@
+
 $(ASYMPTOTIC_SCHWARZSCHILD_TEST_TARGET): tests/test_asymptotic_schwarzschild.c $(COMMON_SOURCES) src/spacetime_schwarzschild.c $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR)
 	$(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -DSPACETIME_SCHWARZSCHILD -Isrc $^ $(LDLIBS) -o $@
 
+# Independent core regression for the backend-free numerical entry localizer.
+# It links only the new module and the shared spacetime dispatch wrapper: no
+# analytic backend, no geodesic integrator and no asymptotic.c are required,
+# so it stays exercisable independently of the route integration.
+$(ASYMPTOTIC_ENTRY_TEST_TARGET): tests/test_asymptotic_entry.c src/asymptotic_entry.c src/spacetime_common.c src/asymptotic_entry.h src/asymptotic.h src/geodesic.h src/spacetime.h src/observer.h | $(TEST_OUT_DIR)
+	$(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc tests/test_asymptotic_entry.c src/asymptotic_entry.c src/spacetime_common.c $(LDLIBS) -o $@
+
 $(FRAME_TEST_TARGET): tests/test_frame.c $(CORE_MINKOWSKI_SOURCES) $(CPU_FFTW_SOURCES) | $(TEST_OUT_DIR)
 	$(CC) $(CPPFLAGS) $(BUILD_CPPFLAGS) $(CFLAGS) $(BUILD_CFLAGS) $(OPENMP_FLAGS) -Isrc $^ $(LDLIBS) -o $@
 
@@ -280,12 +293,14 @@ FAST_PSF_FFTW_TEST_DEP :=
 FAST_PSF_FFTW_TEST_RUN :=
 endif
 
-test: $(CAMERA_TEST_TARGETS) $(TEST_TARGET) $(ADAPTIVE_GEODESIC_TEST_TARGET) $(ASYMPTOTIC_TEST_TARGET) $(ASYMPTOTIC_SCHWARZSCHILD_TEST_TARGET) $(TERMINATION_ORACLE_TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD_TEST_TARGET) $(ALCUBIERRE_TEST_TARGET) $(OBSERVER_TRACK_TEST_TARGET) $(CATALOG_PREFETCH_TEST_TARGET) $(FAST_PSF_FFTW_TEST_DEP) $(TONE_MAP_TEST_TARGET) $(MOVIE_OUTPUT_TEST_TARGET) $(SENSOR_BLOOM_TEST_TARGET)
+test: $(CAMERA_TEST_TARGETS) $(TEST_TARGET) $(ADAPTIVE_GEODESIC_TEST_TARGET) $(ASYMPTOTIC_TEST_TARGET) $(ASYMPTOTIC_ENTRY_TEST_TARGET) $(ASYMPTOTIC_QUADRATIC_TEST_TARGET) $(ASYMPTOTIC_SCHWARZSCHILD_TEST_TARGET) $(TERMINATION_ORACLE_TEST_TARGET) $(FRAME_TEST_TARGET) $(SCHWARZSCHILD_TEST_TARGET) $(ALCUBIERRE_TEST_TARGET) $(OBSERVER_TRACK_TEST_TARGET) $(CATALOG_PREFETCH_TEST_TARGET) $(FAST_PSF_FFTW_TEST_DEP) $(TONE_MAP_TEST_TARGET) $(MOVIE_OUTPUT_TEST_TARGET) $(SENSOR_BLOOM_TEST_TARGET)
 	$(TEST_OUT_DIR)/test_observer_minkowski
 	$(TEST_OUT_DIR)/test_observer_schwarzschild
 	$(TEST_TARGET)
 	$(ADAPTIVE_GEODESIC_TEST_TARGET)
 	$(ASYMPTOTIC_TEST_TARGET)
+	$(ASYMPTOTIC_ENTRY_TEST_TARGET)
+	$(ASYMPTOTIC_QUADRATIC_TEST_TARGET)
 	$(ASYMPTOTIC_SCHWARZSCHILD_TEST_TARGET)
 	$(TERMINATION_ORACLE_TEST_TARGET)
 	$(FRAME_TEST_TARGET)

benchmarks/quadratic_precision/README.md

@@ -0,0 +1,232 @@
+# quadratic_precision — stable entry quadratic: precision vs. cost
+
+Self-contained, reproducible benchmark comparing four arithmetic realisations of
+the **same** asymptotic-entry algebra.  Nothing here modifies production code,
+the `Makefile`, or git state.
+
+## Question
+
+The production camera pre-route solves the relative-distance quadratic
+
+```
+F(s) = |d + q s|^2 - (R0 - rr s)^2 = a s^2 + b s + c
+d = x_cur - c_frame,  q = w_frame + v_frame
+```
+
+in `long double`, with a stable root formula, an ordinary discriminant
+`b*b - 4*a*c` and entry slope, a common power-of-two scaling, an
+uncertainty band, geometric validation and a numerical fallback.  Which matters
+for accuracy and which for speed?
+
+Four variants share the *entire* rest of the module (scaling, uncertainty band,
+validation, fallback, dispatch) and are generated from the current
+`src/asymptotic.c` by rewriting exactly one lexical region:
+
+| variant             | coefficients | type        | discriminant / slope |
+|---------------------|--------------|-------------|----------------------|
+| `ld_plain`          | as production| `long double` | plain `b*b - 4*a*c`, plain slope |
+| `ld_fma`            | as production| `long double` | compensated `fmal` (experimental arm) |
+| `double_fma`        | ordinary `+=`| `double`    | compensated `fma` |
+| `double_fma_coeff`  | `fma` dots, `fma` products | `double` | compensated `fma` |
+
+`double_fma` is the control that isolates *type* precision from *FMA*; the two
+double variants isolate *coefficient accumulation*.
+
+## Generation and build
+
+`build.py` locates the kernel region between two lexical anchors in the current
+working-tree `src/asymptotic.c`:
+
+* start: the comment `/* Long-double coefficients of the relative-distance
+  quadratic`
+* end:   the forward declaration `static void minkowski_route_entry(...)`
+
+It asserts each anchor and the presence of `entry_quadratic_coeffs`,
+`entry_solve`, `entry_discriminant` exactly once, then emits four full-module
+copies under `<output>/generated/`.  Every other line (dispatch, fallback,
+`asymptotic_route_camera`, Schwarzschild path, …) is copied verbatim.  If a
+future edit moves an anchor or changes the region, generation **fails** rather
+than silently measuring the wrong code.
+
+`ld_plain` retains the production kernel verbatim; `ld_fma` rewrites only the
+discriminant body and the entry slope. Double
+variants replace `long double`→`double`, `fmal`/`fabsl`/`fmaxl`/`frexpl`/
+`scalbnl`/`sqrtl`/`copysignl`→their double forms, `LDBL_*`→`DBL_*`, and strip
+`L` suffixes from literals.  A guard rejects any remaining `long double`
+promotion, `L` literal or `*l`/`fmal` call in the double kernel section, so no
+double expression can be re-promoted to `long double`.  `double_fma_coeff`
+additionally changes the coefficient accumulation to `fma`.
+
+The probe `#include`s the generated module, so the private static
+`entry_quadratic_coeffs`/`entry_solve` and the public `asymptotic_route_camera`
+that are exercised are the **actual generated code**, not an independent copy.
+
+Compile flags, identical for all four variants:
+
+```
+-std=c11 -march=native -O2 -DNDEBUG -ffp-contract=off -fopenmp
+```
+
+`-ffp-contract=off` prevents implicit FMA contraction from silently changing
+the `ld_plain` arm; explicit `fma`/`fmal` still lower as written.
+
+The probe links a minimal production subset (`geodesic.c`,
+`asymptotic_entry.c`, `asymptotic_schwarzschild.c`, `spacetime_common.c`,
+`spacetime_minkowski.c`) plus libm.  No FFTW, PNG, catalog or observer-track
+data are needed.
+
+## Inputs
+
+All fixtures are frozen as IEEE-754 hex literals (`cases.json` and the
+generated C header are bit-identical), including the two production grazing
+rows reproduced from the Alcubierre pre-route (`x=(0,-24,0)`, `centre=2t`,
+`v=2`, `R=5`, canonical `w` hex values).  Families:
+
+* **curated / adversarial**: head-on hit/miss, clear miss, grazing
+  `y=nextafter(R,0)`, exact tangent, near-boundary, `1e10 + 0.5, R=1`
+  cancellation, large-`t` `centre=0.1t`, exact linear `a==0` (inward and
+  initially outward), translated-origin cancellation.
+* **fixed**: axes and three oblique frames, radii
+  `1e-100…1e100`, `D/R ∈ {1+ulp,2,10,100,512,1024,1e4,1e8,1e10}`, impact
+  ratios `{0,.5,.99,1-1e-6,nextafter(1,0),1,nextafter(1,∞),1+1e-6,1.1}`.
+  The `512`/`1024` ratios sample the crossover where long-double coefficient
+  rounding meets the `128 ε_D R²` geometry tolerance (`≈512`) and the double
+  crossover (`≈11`), separating fallback rates from the all-UNCERTAIN `1e10`
+  regime.
+* **moving**: `v ∈ {0,.1,2,10}` along a transverse direction.
+* **growing**: `rr = -1` and `nextafter(-1,∓∞)`, **inward and
+  initially-outward** photons (`w = ∓tow`), axis and oblique.
+* **shrinking**: `rr>0`, labelled/domain-checked algebra (radius would go
+  negative on the open past; roots outside `R0 - rr·s > 0` are non-physical).
+* **route**: `rr = 0` plus **`rr<0` growing inward/outward** families
+  (positive radius on the whole past, `valid_t_min=-1e300`) across all three
+  direction frames, driven through the public `asymptotic_route_camera`.  The
+  oblique-2 growing-outward family reproduces the double kernel false-MISS
+  cases, where a kernel MISS bypasses the fallback and the route escapes
+  directly; the route reference uses the actual normalised canonical `w`, so
+  normalisation effects are explicit.
+
+Exact counts are emitted by the run (and in `summary.json`); they are not
+hard-coded here.
+
+## Reference oracle
+
+Inputs are exact, so they are represented as `fractions.Fraction` (Python
+stdlib).  `d`, `q`, `a`, `b`, `c`, the discriminant and polynomial residuals
+are **exact rationals**; only `sqrt` uses `decimal` (`--precision`, default
+160).  A custom hex parser handles both `%a` (double) and `%La` (long double)
+without `float.fromhex`, which would silently drop a 64-bit long-double mantissa
+to 53 bits.  Reference classification mirrors production semantics: `c<0` is
+INSIDE (not an entry failure), `c==0` uses the boundary slope, `a==0` is the
+exact linear branch, a real double root/tangent is MISS.
+
+Route references are built from the **actual canonical state the probe printed**
+(`canonx*`, `canonw*`) joined with the callback samples at `t0`, so the 1-ulp
+normalisation of `w` and both callback models (`centre = c0 + v(t-t0)` and the
+production `centre = v t`) are handled explicitly.
+
+## Metrics
+
+Kernel:
+
+* status counts (ENTRY/MISS/UNCERTAIN) vs the reference, **false MISS**
+  (reference ENTER, kernel MISS) and **false candidate** (reference MISS,
+  kernel ENTRY);
+* root error in ULP of the returned double root and relative error, for the
+  common reference-ENTER set and the intersection where *all* variants returned
+  ENTRY (so more UNCERTAIN cannot look better by selection);
+* polynomial residual at the candidate root using the actual printed
+  coefficients, the exact ideal polynomial, and the nearest-double-rounded
+  reference root as a baseline;
+* `a` exact-zero vs kernel-zero, `a` sign mismatches, `a` collapse/spurious
+  non-zero, and conditioning labels (`near_linear`, `late`, `ill_conditioned`).
+  Outward/growing and near-linear cases are ill-conditioned; a rounded-zero `a`
+  can turn a very late true entry into a false MISS, so these are reported
+  separately and are not treated as a universal precision claim.
+
+Route:
+
+* kind/status vs reference, false MISS, false candidate, camera-INSIDE
+  mismatches, total `INVALID/ENTRY_UNCONFIRMED` outcomes independent of reference
+  classification (a conservative refusal to confirm is not a false escape);
+* the **public kernel classification at the actual canonical inputs**
+  (`kern_status` in the route CSV), so a reference ENTER / kernel MISS / route
+  ESCAPED (`kernel_false_miss_escaped`) is visible and is not confused with
+  canonical normalisation;
+* fast-path vs fallback (`entry_fallback_evaluations`), `F/tol` for accepted
+  entries, and Pi/`L_camera` preservation.
+
+The reference enforces `R0 - rr·s > 0` for quadratic roots; a positive root
+outside the physical radius domain is reported as `domain_clipped` and is not
+counted as a physical ENTER/false-MISS.  A 160-vs-240-digit consistency check is
+run for every curated and near-linear/late/domain-clipped id and must report the
+same status, discriminant sign and nearest-double root.
+
+Timing:
+
+* kernel: coefficient assembly + `entry_solve`, serial, `noinline` +
+  `volatile` sink, and an inline-asm opaque loop index so GCC cannot hoist the
+  pure kernel out of the repetition loop (identical for all variants).  A
+  linearity check runs the kernel at 1× and 2× calls and reports the ratio
+  (expect ≈2) to prove per-call execution.
+* route: public pre-route, OpenMP `static` schedule + reductions, with
+  route-kind and fallback counts captured so a timing gain cannot come from
+  silently escaping more rays.  `--threads` defaults to 4.
+
+Codegen evidence is collected with `objdump`/`nm`: for the double variants
+`entry_solve` contains hardware `vfmadd` and zero x87; for `ld_fma` it contains
+x87 plus `fmal` PLT relocations (extended-precision FMA is a libm software
+routine on x86-64).
+
+## Run
+
+```sh
+python3 benchmarks/quadratic_precision/run.py \
+    --output-dir /tmp/opencode/quadratic-comparison
+```
+
+Useful overrides: `--rounds N`, `--threads T`, `--precision P`,
+`--kernel-target-calls N`, `--route-target-calls N`, `--skip-build` (reuse the
+last build), `--repo PATH`, `--cc CC`.
+
+Artifacts (all under `--output-dir`):
+
+```
+generated/{ld_plain,ld_fma,double_fma,double_fma_coeff}.c
+cases/quadratic_cases.h, cases/cases.json
+environment.json, build_manifest.json
+raw/kernel_<v>.csv, raw/route_<v>.csv
+raw/kernel_metrics_<v>.csv, raw/route_metrics_<v>.csv
+raw/curated_kernel.csv, raw/curated_route.csv
+raw/microbench_<v>_r*.json, raw/routebench_<v>_r*.json
+raw/fma_verify.json, raw/{linearity,curated}*
+summary.json
+logs/{build_commands,accuracy_*,microbench_*,routebench_*}.log
+```
+
+`summary.json` is the machine-readable aggregate; the script also prints the
+coverage counts, per-variant ULP distributions, the curated case table, codegen
+evidence and median ns/call.
+
+## Interpretation caveats
+
+* FMA improves the rounding of one product only.  It **cannot** restore the
+  information already lost in the double inputs and in the coefficient sums
+  (`d = x - c`, `q = w + v`, and the `qq`/`dq`/`dd` accumulations).  The data
+  should be read stage by stage.
+* A kernel UNCERTAIN is not a failure: the shared route fallback may confirm the
+  entry.  A kernel **MISS bypasses the fallback entirely**, so growing/outward
+  false-MISS families are exercised through dedicated `rr<0` route fixtures and
+  reported both as the raw kernel decision (`kern_status`) and the route
+  outcome.  This benchmark does **not** claim that the fallback rescues a raw
+  kernel MISS on the same case; it only observes that more fallback is used to
+  confirm inaccurate candidates.
+* The double uncertainty band differs from the long-double one by the precision
+  term: production uses `64·(DBL_EPSILON + LDBL_EPSILON)·scale` while the
+  double variants use `64·(DBL_EPSILON + DBL_EPSILON)·scale`; since
+  `LDBL_EPSILON ≪ DBL_EPSILON` this is roughly a factor of two in the band
+  width, which only matters exactly at the UNCERTAIN/MISS border.
+* Number/scale results are machine- and compiler-specific; `environment.json`
+  records the source hash, compiler, flags and platform.
+* This benchmark does not establish a universal guarantee, only the sampled
+  behaviour on the frozen fixture set.

benchmarks/quadratic_precision/build.py

@@ -0,0 +1,489 @@
+#!/usr/bin/env python3
+"""Generate the four precision variants of the asymptotic entry kernel and build
+one self-contained probe executable per variant.
+
+This module never edits the production tree.  It reads the *current working
+tree* ``src/asymptotic.c``, extracts the "entry kernel" region delimited by two
+lexical anchors, emits four generated full-module copies under the scratch
+directory, and compiles a probe that textually includes exactly one copy so the
+private static ``entry_solve``/``entry_quadratic_coeffs`` are exercised as the
+real generated code (not an independent toy copy).
+
+Variants (all share the same stable-entry algebra, scaling, uncertainty policy,
+validation and fallback; only the arithmetic under test changes):
+
+  ld_plain         long double coefficients/solver, ordinary b*b-4ac and slope
+  ld_fma           reconstructed experimental compensated fma disc/slope
+  double_fma       double type/solver with compensated fma disc/slope,
+                   ordinary coefficient accumulation (type-precision control)
+  double_fma_coeff double like double_fma plus fma coefficient accumulation
+
+Only the kernel region is rewritten.  Every other production dispatch/fallback
+path is copied verbatim.
+"""
+from __future__ import annotations
+
+import hashlib
+import json
+import os
+import re
+import shutil
+import subprocess
+import sys
+from pathlib import Path
+
+# ---------------------------------------------------------------------------
+# Fixed compiler configuration.  Identical for all four variants so the
+# comparison is not confounded by flags.  ``-ffp-contract=off`` disables the
+# implicit FMA contraction GCC enables at -O2 for this target, so the ld_plain
+# arm cannot silently gain extra precision; explicit fma()/fmal() calls still
+# lower to hardware FMA.
+# ---------------------------------------------------------------------------
+BASE_FLAGS = [
+    "-std=c11",
+    "-march=native",
+    "-O2",
+    "-DNDEBUG",
+    "-ffp-contract=off",
+    "-fopenmp",
+]
+
+VARIANTS = ("ld_plain", "ld_fma", "double_fma", "double_fma_coeff")
+
+# Common production translation units needed by the probe.  Chosen as the
+# minimal dependency closure of the entry kernel plus the public pre-route:
+# geodesic initialisation, the numerical entry localizer, the Schwarzschild
+# exterior referenced by the full module, the dispatch wrapper and the
+# Minkowski provider.  main/asymptotic/frame/catalog/PSF/FFTW are excluded.
+COMMON_SOURCES = (
+    "src/geodesic.c",
+    "src/asymptotic_entry.c",
+    "src/asymptotic_schwarzschild.c",
+    "src/spacetime_common.c",
+    "src/spacetime_minkowski.c",
+)
+
+START_ANCHOR = "/* Long-double coefficients of the relative-distance quadratic"
+FWD_PREFIX = "static void minkowski_route_entry("
+END_ANCHOR = (
+    "static void minkowski_route_entry(const SpacetimeAsymptoticEnd *end, "
+    "double t0,"
+)
+
+
+class BuildError(RuntimeError):
+    pass
+
+
+def _replace_once(text: str, old: str, new: str, what: str) -> str:
+    count = text.count(old)
+    if count != 1:
+        raise BuildError(f"expected exactly one occurrence of {what}, found {count}")
+    return text.replace(old, new, 1)
+
+
+def locate_kernel(lines: list[str]) -> tuple[int, int]:
+    """Return [start, end) line indices of the rewritten kernel region.
+
+    start is the 'Long-double coefficients' comment; end is the forward
+    declaration of ``minkowski_route_entry`` (preserved verbatim).
+    """
+    starts = [i for i, l in enumerate(lines) if l.startswith(START_ANCHOR)]
+    if len(starts) != 1:
+        raise BuildError(
+            f"kernel start anchor must appear exactly once, found {len(starts)}"
+        )
+    start = starts[0]
+    ends = [
+        i
+        for i, l in enumerate(lines)
+        if i > start and l.startswith(FWD_PREFIX)
+    ]
+    if not ends:
+        raise BuildError("kernel end anchor (minkowski_route_entry) not found")
+    end = ends[0]
+    if not lines[end].startswith(END_ANCHOR):
+        raise BuildError(
+            "kernel end anchor does not match the expected signature; "
+            "production source changed"
+        )
+    region = "".join(lines[start:end])
+    for needle in (
+        "typedef struct {",
+        "static EntryQuadratic entry_quadratic_coeffs(",
+        "static EntrySolveResult entry_solve(",
+        "static long double entry_discriminant(",
+        "EntryQuadratic;",
+    ):
+        if region.count(needle) != 1:
+            raise BuildError(
+                f"kernel region must contain exactly one '{needle}', "
+                f"found {region.count(needle)}"
+            )
+    # The region must not contain the forward declaration or any later route
+    # code: those must be byte-for-byte preserved.
+    if "minkowski_preroute" in region or "AsymptoticStatus asymptotic_route" in region:
+        raise BuildError("kernel region overran into preserved route code")
+    return start, end
+
+
+def transform_ld_plain(kernel: str) -> str:
+    """Convert a compensated experimental kernel back to plain long double."""
+    old_disc = (
+        "  const long double four_a = 4.0L * k->a;\n"
+        "  const long double ac = four_a * k->c;\n"
+        "  const long double ac_error = fmal(four_a, k->c, -ac);\n"
+        "  *scale = k->b * k->b + fabsl(ac);\n"
+        "  return fmal(k->b, k->b, -ac) - ac_error;\n"
+    )
+    new_disc = (
+        "  const long double four_a = 4.0L * k->a;\n"
+        "  const long double ac = four_a * k->c;\n"
+        "  *scale = k->b * k->b + fabsl(ac);\n"
+        "  return k->b * k->b - ac;\n"
+    )
+    out = _replace_once(kernel, old_disc, new_disc, "ld_plain discriminant body")
+    old_slope = "const long double slope = fmal(2.0L * k->a, s, k->b);"
+    new_slope = "const long double slope = 2.0L * k->a * s + k->b;"
+    out = _replace_once(out, old_slope, new_slope, "ld_plain slope")
+    return out
+
+
+def _assert_double_clean(t: str) -> None:
+    """The rewritten double kernel must contain no long-double promotion."""
+    forbidden = ("long double", "fmal(", "fabsl(", "fmaxl(", "frexpl(",
+                 "scalbnl(", "sqrtl(", "copysignl(", "LDBL_")
+    for tok in forbidden:
+        if tok in t:
+            raise BuildError(f"double kernel still contains '{tok}'")
+    if re.search(r"(?<=[0-9])L(?![A-Za-z0-9_])", t):
+        raise BuildError("double kernel still contains an L-suffixed literal")
+
+
+def transform_to_double(kernel: str) -> str:
+    """Convert the kernel from long double to double arithmetic.
+
+    Explicit fmal->fma, *l->*, LDBL_*->DBL_*, long double->double and strips the
+    L suffix from the (few) decimal literals.  Explicit fma calls are retained
+    (double_fma keeps compensated discriminant/slope).
+    """
+    t = kernel
+    for old, new in (
+        ("fmal(", "fma("),
+        ("fabsl(", "fabs("),
+        ("fmaxl(", "fmax("),
+        ("frexpl(", "frexp("),
+        ("scalbnl(", "scalbn("),
+        ("sqrtl(", "sqrt("),
+        ("copysignl(", "copysign("),
+    ):
+        t = t.replace(old, new)
+    t = t.replace("LDBL_MIN", "DBL_MIN")
+    t = t.replace("LDBL_EPSILON", "DBL_EPSILON")
+    t = t.replace("long double", "double")
+    # Strip L suffixes from decimal literals; otherwise 4.0L/2.0L would
+    # re-promote the double expression and confound type accuracy/timing.
+    t = re.sub(r"(?<=[0-9])L(?![A-Za-z0-9_])", "", t)
+    _assert_double_clean(t)
+    return t
+
+
+def transform_double_fma_coeff(kernel: str) -> str:
+    """double variants plus fma coefficient accumulation and products."""
+    dbl = transform_to_double(kernel)
+    old_block = (
+        "  double qq = 0.0, dot_dq = 0.0, dot_dd = 0.0;\n"
+        "  for (int i = 0; i < 3; ++i) {\n"
+        "    qq += q[i] * q[i];\n"
+        "    dot_dq += d[i] * q[i];\n"
+        "    dot_dd += d[i] * d[i];\n"
+        "  }\n"
+        "  const double R0_ld = (double)R0;\n"
+        "  const double rr_ld = (double)rr;\n"
+        "  EntryQuadratic k;\n"
+        "  k.a = qq - rr_ld * rr_ld;\n"
+        "  k.b = 2.0 * (dot_dq + R0_ld * rr_ld);\n"
+        "  k.c = dot_dd - R0_ld * R0_ld;\n"
+    )
+    new_block = (
+        "  double qq = 0.0, dot_dq = 0.0, dot_dd = 0.0;\n"
+        "  for (int i = 0; i < 3; ++i) {\n"
+        "    qq = fma(q[i], q[i], qq);\n"
+        "    dot_dq = fma(d[i], q[i], dot_dq);\n"
+        "    dot_dd = fma(d[i], d[i], dot_dd);\n"
+        "  }\n"
+        "  const double R0_ld = (double)R0;\n"
+        "  const double rr_ld = (double)rr;\n"
+        "  EntryQuadratic k;\n"
+        "  k.a = fma(-rr_ld, rr_ld, qq);\n"
+        "  k.b = 2.0 * fma(R0_ld, rr_ld, dot_dq);\n"
+        "  k.c = fma(-R0_ld, R0_ld, dot_dd);\n"
+    )
+    return _replace_once(dbl, old_block, new_block, "double_fma_coeff block")
+
+
+def generate_variants(repo: Path, out_dir: Path) -> dict:
+    """Read current src/asymptotic.c and emit the four variant modules."""
+    src_path = repo / "src" / "asymptotic.c"
+    text = src_path.read_text()
+    lines = text.splitlines(keepends=True)
+    start, end = locate_kernel(lines)
+    kernel = "".join(lines[start:end])
+    prefix = "".join(lines[:start])
+    suffix = "".join(lines[end:])
+    sha = hashlib.sha256(text.encode()).hexdigest()
+
+    # Production now uses plain long-double arithmetic. Reconstruct the former
+    # compensated arm so the original four-way comparison remains reproducible.
+    plain_disc = (
+        "  const long double four_a = 4.0L * k->a;\n"
+        "  const long double ac = four_a * k->c;\n"
+        "  *scale = k->b * k->b + fabsl(ac);\n"
+        "  return k->b * k->b - ac;\n"
+    )
+    fused_disc = plain_disc.replace(
+        "  *scale =", "  const long double ac_error = fmal(four_a, k->c, -ac);\n  *scale ="
+    ).replace("return k->b * k->b - ac;", "return fmal(k->b, k->b, -ac) - ac_error;")
+    fused = _replace_once(kernel, plain_disc, fused_disc, "ld_fma discriminant body")
+    fused = _replace_once(
+        fused, "const long double slope = 2.0L * k->a * s + k->b;",
+        "const long double slope = fmal(2.0L * k->a, s, k->b);", "ld_fma slope"
+    )
+    kernels = {
+        "ld_plain": kernel,
+        "ld_fma": fused,
+        "double_fma": transform_to_double(fused),
+        "double_fma_coeff": transform_double_fma_coeff(fused),
+    }
+
+    gen_dir = out_dir / "generated"
+    gen_dir.mkdir(parents=True, exist_ok=True)
+    paths = {}
+    metas = {}
+    for name in VARIANTS:
+        body = kernels[name]
+        banner = (
+            f"/* GENERATED FILE - do not edit.\n"
+            f" * variant: {name}\n"
+            f" * source:  src/asymptotic.c sha256={sha}\n"
+            f" * region lines [{start + 1}, {end}] rewritten; all other code verbatim.\n"
+            f" */\n"
+        )
+        out = banner + prefix + body + suffix
+        path = gen_dir / f"{name}.c"
+        path.write_text(out)
+        paths[name] = path
+        metas[name] = {
+            "path": str(path),
+            "sha256": hashlib.sha256(out.encode()).hexdigest(),
+            "kernel_sha256": hashlib.sha256(body.encode()).hexdigest(),
+        }
+    return {
+        "source": str(src_path),
+        "source_sha256": sha,
+        "kernel_line_range": [start + 1, end],
+        "variants": metas,
+    }
+
+
+def cc_version(cc: str) -> str:
+    try:
+        out = subprocess.run(
+            [cc, "--version"], capture_output=True, text=True, check=False
+        )
+        return (out.stdout or out.stderr).splitlines()[0] if out.stdout or out.stderr else ""
+    except OSError:
+        return ""
+
+
+def environment(repo: Path, cc: str, threads: int) -> dict:
+    src_path = repo / "src" / "asymptotic.c"
+    env = {
+        "repo": str(repo),
+        "cc": cc,
+        "cc_version": cc_version(cc),
+        "base_flags": list(BASE_FLAGS),
+        "threads": threads,
+        "source_sha256": hashlib.sha256(src_path.read_bytes()).hexdigest(),
+        "git_head": _git(repo, "rev-parse", "HEAD"),
+        "git_status_short": _git(repo, "status", "--short"),
+        "uname": os.uname().sysname + " " + os.uname().release + " " + os.uname().machine,
+        "cpu_model": _cpu_model(),
+    }
+    return env
+
+
+def _git(repo: Path, *args: str) -> str:
+    try:
+        out = subprocess.run(
+            ["git", "-C", str(repo), *args], capture_output=True, text=True, check=False
+        )
+        return out.stdout.strip()
+    except OSError:
+        return ""
+
+
+def _cpu_model() -> str:
+    try:
+        for line in Path("/proc/cpuinfo").read_text().splitlines():
+            if line.startswith("model name"):
+                return line.split(":", 1)[1].strip()
+    except OSError:
+        pass
+    return "unknown"
+
+
+def compile_common(repo: Path, out_dir: Path, cc: str, flags: list[str]) -> list[Path]:
+    obj_dir = out_dir / "build" / "common"
+    obj_dir.mkdir(parents=True, exist_ok=True)
+    objs = []
+    for rel in COMMON_SOURCES:
+        src = repo / rel
+        obj = obj_dir / (Path(rel).stem + ".o")
+        cmd = [cc, *flags, f"-I{repo / 'src'}", "-c", str(src), "-o", str(obj)]
+        _run(cmd, out_dir, f"compile common {rel}")
+        objs.append(obj)
+    return objs
+
+
+def compile_probe(
+    repo: Path,
+    out_dir: Path,
+    cc: str,
+    flags: list[str],
+    variant: str,
+    probe_src: Path,
+    cases_dir: Path,
+    common_objs: list[Path],
+) -> Path:
+    """Compile the probe (which #includes the generated variant) and link."""
+    exe_dir = out_dir / "build" / "bin"
+    exe_dir.mkdir(parents=True, exist_ok=True)
+    obj = out_dir / "build" / f"probe_{variant}.o"
+    variant_src = out_dir / "generated" / f"{variant}.c"
+    cmd = [
+        cc,
+        *flags,
+        f"-I{repo / 'src'}",
+        f"-I{cases_dir}",
+        f'-DPROBE_VARIANT_SOURCE="{variant_src}"',
+        f'-DPROBE_VARIANT_NAME="{variant}"',
+        "-c",
+        str(probe_src),
+        "-o",
+        str(obj),
+    ]
+    _run(cmd, out_dir, f"compile probe {variant}")
+    exe = exe_dir / f"probe_{variant}"
+    link = [cc, *flags, str(obj), *[str(o) for o in common_objs], "-lm", "-o", str(exe)]
+    _run(link, out_dir, f"link probe {variant}")
+    return exe
+
+
+def _run(cmd: list[str], out_dir: Path, what: str) -> None:
+    log_dir = out_dir / "logs"
+    log_dir.mkdir(parents=True, exist_ok=True)
+    with open(log_dir / "build_commands.log", "a") as fh:
+        fh.write(what + "\n$ " + " ".join(cmd) + "\n")
+    try:
+        proc = subprocess.run(cmd, capture_output=True, text=True, check=False)
+    except OSError as exc:
+        raise BuildError(f"{what}: {exc}") from exc
+    with open(log_dir / "build_commands.log", "a") as fh:
+        fh.write(f"exit={proc.returncode}\n")
+        if proc.stdout:
+            fh.write(proc.stdout)
+        if proc.stderr:
+            fh.write(proc.stderr)
+    if proc.returncode != 0:
+        raise BuildError(
+            f"{what} failed (exit {proc.returncode}); see logs/build_commands.log\n"
+            + (proc.stderr or "")[-4000:]
+        )
+
+
+def build_all(repo: Path, out_dir: Path, cc: str, flags: list[str], threads: int) -> dict:
+    out_dir.mkdir(parents=True, exist_ok=True)
+    (out_dir / "logs").mkdir(parents=True, exist_ok=True)
+    variant_info = generate_variants(repo, out_dir)
+    common = compile_common(repo, out_dir, cc, flags)
+    probe_src = repo / "benchmarks" / "quadratic_precision" / "fixtures" / "probe.c"
+    cases_dir = out_dir / "cases"
+    exes = {}
+    for name in VARIANTS:
+        exes[name] = compile_probe(
+            repo, out_dir, cc, flags, name, probe_src, cases_dir, common
+        )
+    result = {"variants": variant_info["variants"],
+              "source_sha256": variant_info["source_sha256"],
+               "base_flags": list(flags),
+               "dependencies": dependency_hashes(repo),
+               "compiler": {"command": cc, "version": cc_version(cc)},
+              "probe_sha256": _file_sha(probe_src),
+              "cases_header_sha256": _file_sha(cases_dir / "quadratic_cases.h")}
+    result["executables"] = {k: str(v) for k, v in exes.items()}
+    result["common_objects"] = [str(o) for o in common]
+    (out_dir / "build_manifest.json").write_text(json.dumps(result, indent=2) + "\n")
+    return result
+
+
+def _file_sha(path: Path):
+    if not Path(path).exists():
+        return None
+    return hashlib.sha256(Path(path).read_bytes()).hexdigest()
+
+
+def dependency_hashes(repo: Path):
+    # Include all project headers conservatively, including transitive includes.
+    paths = {repo / rel for rel in COMMON_SOURCES}
+    paths.update((repo / "src").rglob("*.h"))
+    return {str(p.relative_to(repo)): _file_sha(p) for p in sorted(paths)}
+
+
+def verify_manifest(repo: Path, out_dir: Path, flags: list[str], cc: str):
+    """Check that an existing build matches the current source, flags and
+    fixtures.  Returns (manifest, list_of_mismatch_reasons).  Regenerates the
+    variant sources as a side effect so their hashes can be compared."""
+    manifest = json.loads((out_dir / "build_manifest.json").read_text())
+    reasons = []
+    if manifest.get("dependencies") != dependency_hashes(repo):
+        reasons.append("dependencies")
+    if manifest.get("compiler") != {"command": cc, "version": cc_version(cc)}:
+        reasons.append("compiler")
+    current = generate_variants(repo, out_dir)
+    if manifest.get("source_sha256") != current["source_sha256"]:
+        reasons.append("source_sha256")
+    if list(manifest.get("base_flags", [])) != list(flags):
+        reasons.append("base_flags")
+    for name, meta in current["variants"].items():
+        if manifest.get("variants", {}).get(name, {}).get("sha256") != meta["sha256"]:
+            reasons.append(f"variant:{name}")
+    probe_src = repo / "benchmarks" / "quadratic_precision" / "fixtures" / "probe.c"
+    if manifest.get("probe_sha256") != _file_sha(probe_src):
+        reasons.append("probe_sha256")
+    cases_header = out_dir / "cases" / "quadratic_cases.h"
+    if cases_header.exists() and manifest.get("cases_header_sha256") != _file_sha(cases_header):
+        reasons.append("cases_header_sha256")
+    for name, exe in manifest.get("executables", {}).items():
+        if not Path(exe).exists():
+            reasons.append(f"missing_exe:{name}")
+    return manifest, reasons
+
+
+if __name__ == "__main__":
+    import argparse
+
+    ap = argparse.ArgumentParser(description=__doc__)
+    ap.add_argument("--repo", default=str(Path(__file__).resolve().parents[2]))
+    ap.add_argument("--output-dir", default="/tmp/opencode/quadratic-comparison")
+    ap.add_argument("--cc", default=os.environ.get("CC", "cc"))
+    ap.add_argument("--threads", type=int, default=4)
+    args = ap.parse_args()
+    repo = Path(args.repo).resolve()
+    out = Path(args.output_dir).resolve()
+    try:
+        info = build_all(repo, out, args.cc, BASE_FLAGS, args.threads)
+    except BuildError as exc:
+        print(f"build failed: {exc}", file=sys.stderr)
+        sys.exit(1)
+    print(json.dumps(info, indent=2))

benchmarks/quadratic_precision/cases.py

@@ -0,0 +1,441 @@
+#!/usr/bin/env python3
+"""Deterministic fixture generation for the quadratic precision benchmark.
+
+All values are frozen as IEEE-754 hex literals in the generated C header and as
+hex strings in cases.json, so the C probe and the Python oracle see bit-for-bit
+identical inputs.  No external catalog/observer/slab data is required.
+
+Case families (kernel):
+  curated        explicit adversarial / production fixtures
+  fixed          fixed sphere, axis and oblique directions
+  moving         constant sphere velocity (radial/transverse), inward photons
+  growing        rr < 0 (radius grows along the past parameter), both inward
+                 (toward-center) and initially-outward photons, axis + oblique
+  shrinking      rr > 0, kernel only (radius would go negative on the far past)
+  cancellation   translated/far-origin style input cancellation
+
+Route cases use only rr == 0 (positive radius on the whole open past segment)
+plus the two production grazing rows.
+"""
+from __future__ import annotations
+
+import json
+import math
+from pathlib import Path
+
+# Radii spanning subnormal-adjacent to huge scales.
+RADII = [1e-100, 1e-10, 1.0, 1e10, 1e100]
+DIST = [
+    math.nextafter(1.0, math.inf),
+    2.0,
+    10.0,
+    100.0,
+    512.0,
+    1024.0,
+    1e4,
+    1e8,
+    1e10,
+]
+IMPACTS = [
+    0.0,
+    0.5,
+    0.99,
+    1.0 - 1e-6,
+    math.nextafter(1.0, 0.0),
+    1.0,
+    math.nextafter(1.0, math.inf),
+    1.0 + 1e-6,
+    1.1,
+]
+VELS = [0.0, 0.1, 2.0, 10.0]
+RR_GROW = [
+    -1.0,
+    math.nextafter(-1.0, -math.inf),
+    math.nextafter(-1.0, math.inf),
+    -0.999999,
+]
+
+
+def norm3(v):
+    n = math.sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2])
+    return [v[0] / n, v[1] / n, v[2] / n]
+
+
+def cross(a, b):
+    return [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2],
+            a[0] * b[1] - a[1] * b[0]]
+
+
+# Fixed direction frames: (tow, perp).  `tow` points from the sphere centre to
+# the camera; `w = -tow` is the inward (toward-centre) past direction.
+DIRS_AXIS = ([1.0, 0.0, 0.0], [0.0, 1.0, 0.0])
+DIRS_OBLIQUE1 = (
+    norm3([1.0, 2.0, 3.0]),
+    norm3(cross(norm3([1.0, 2.0, 3.0]), [0.0, 0.0, 1.0])),
+)
+DIRS_OBLIQUE2 = (
+    norm3([-2.0, 1.0, 0.7]),
+    norm3(cross(norm3([-2.0, 1.0, 0.7]), [0.0, 1.0, 0.0])),
+)
+DIRS = [DIRS_AXIS, DIRS_OBLIQUE1, DIRS_OBLIQUE2]
+
+
+def _case(cid, category, x, c, w, v, R0, rr):
+    return {
+        "id": cid,
+        "category": category,
+        "x": list(x),
+        "c": list(c),
+        "w": list(w),
+        "v": list(v),
+        "R0": R0,
+        "rr": rr,
+    }
+
+
+def curated_kernel_cases():
+    """Explicit adversarial and production-reproduction fixtures."""
+    cases = []
+    t63 = float.fromhex("0x1.fa8f5c28f5c29p+3")
+    t64 = float.fromhex("0x1.fb17e4b17e4b1p+3")
+    w63 = [-float.fromhex("0x1.d4afba4704cap-2"),
+           float.fromhex("0x1.c7378f8e872d1p-1"),
+           float.fromhex("0x1.15bad4e30e8ddp-8")]
+    w64 = [-float.fromhex("0x1.f98ae1a782104p-2"),
+           float.fromhex("0x1.bd3bb364ac492p-1"),
+           float.fromhex("0x1.102d2a1c6ac74p-7")]
+    # Production Alcubierre grazing rows: x=(0,-24,0), centre=2t, v=2, R=5.
+    cases.append(_case(0, "real63", [0.0, -24.0, 0.0],
+                       [2.0 * t63, 0.0, 0.0], w63, [2.0, 0.0, 0.0], 5.0, 0.0))
+    cases.append(_case(1, "real64", [0.0, -24.0, 0.0],
+                       [2.0 * t64, 0.0, 0.0], w64, [2.0, 0.0, 0.0], 5.0, 0.0))
+
+    D = 10.0
+    R = 5.0
+    cases += [
+        _case(2, "headon_hit", [D, 0.0, 0.0], [0.0] * 3, [-1.0, 0.0, 0.0],
+              [0.0] * 3, R, 0.0),
+        _case(3, "headon_miss", [D, 0.0, 0.0], [0.0] * 3, [1.0, 0.0, 0.0],
+              [0.0] * 3, R, 0.0),
+        _case(4, "clear_miss", [D, 6.0, 0.0], [0.0] * 3, [-1.0, 0.0, 0.0],
+              [0.0] * 3, R, 0.0),
+        _case(5, "grazing_in", [D, math.nextafter(R, 0.0), 0.0], [0.0] * 3,
+              [-1.0, 0.0, 0.0], [0.0] * 3, R, 0.0),
+        _case(6, "exact_tangent", [D, R, 0.0], [0.0] * 3, [-1.0, 0.0, 0.0],
+              [0.0] * 3, R, 0.0),
+        _case(7, "grazing_out", [D, math.nextafter(R, math.inf), 0.0],
+              [0.0] * 3, [-1.0, 0.0, 0.0], [0.0] * 3, R, 0.0),
+        _case(8, "near_boundary_hit",
+              [math.nextafter(R, math.inf), 0.0, 0.0], [0.0] * 3,
+              [-1.0, 0.0, 0.0], [0.0] * 3, R, 0.0),
+        # 1e10 + 0.5, R=1: c loses the transverse term, b*b-4ac rounds to 0.
+        _case(9, "cancel_1e10_05", [1e10, 0.5, 0.0], [0.0] * 3,
+              [-1.0, 0.0, 0.0], [0.0] * 3, 1.0, 0.0),
+        _case(10, "cancel_1e10_1", [1e10, 1.0, 0.0], [0.0] * 3,
+              [-1.0, 0.0, 0.0], [0.0] * 3, 1.0, 0.0),
+    ]
+    # Large-t positive reconstructed minimum (tests/test_asymptotic.c).
+    camera_x = float.fromhex("0x1.6bcc41e901908p+46")
+    t0 = 1e15
+    vx = 0.1
+    cases.append(_case(11, "large_t_01t", [camera_x, 0.0, 0.0],
+                       [vx * t0, 0.0, 0.0], [-1.0, 0.0, 0.0], [vx, 0.0, 0.0],
+                       10.0, 0.0))
+    # Exact linear growing sphere (a == 0), inward.
+    cases.append(_case(12, "linear_grow_in", [100.0, 0.0, 0.0], [0.0] * 3,
+                       [-1.0, 0.0, 0.0], [0.0] * 3, 10.0, -1.0))
+    # Exact linear, initially outward: constant gap, honest miss.
+    cases.append(_case(13, "linear_grow_out", [100.0, 0.0, 0.0], [0.0] * 3,
+                       [1.0, 0.0, 0.0], [0.0] * 3, 10.0, -1.0))
+    # Near-linear oblique, inward vs initially outward, rr just around -1.
+    tow = DIRS_OBLIQUE1[0]
+    perp = DIRS_OBLIQUE1[1]
+    for rr in RR_GROW:
+        for sign, tag in ((-1.0, "in"), (1.0, "out")):
+            w = [sign * tow[i] for i in range(3)]
+            x = [50.0 * tow[i] + 0.25 * perp[i] for i in range(3)]
+            cases.append(_case(len(cases), f"nearlin_{tag}", x, [0.0] * 3, w,
+                               [0.0] * 3, 4.0, rr))
+    # Translated-origin style cancellation: d = x - c with x = c + small.
+    base = 1e10
+    cbase = [base, -base, base * 0.5]
+    xb = [cbase[0] + 10.0, cbase[1] + 0.5, cbase[2] + 0.0]
+    cases.append(_case(len(cases), "translated_origin", xb, cbase,
+                       [-1.0, 0.0, 0.0], [0.0] * 3, 1.0, 0.0))
+    return cases
+
+
+def _broad_fixed(cid):
+    for di, (tow, perp) in enumerate(DIRS):
+        for R in RADII:
+            for dr in DIST:
+                for ir in IMPACTS:
+                    D = dr * R
+                    b = ir * R
+                    x = [D * tow[i] + b * perp[i] for i in range(3)]
+                    w = [-tow[i] for i in range(3)]
+                    yield _case(cid, f"fixed_d{di}", x, [0.0] * 3, w, [0.0] * 3,
+                                R, 0.0)
+                    cid += 1
+    return cid
+
+
+def _broad_moving(cid):
+    for (tow, perp) in DIRS:
+        for R in (1e-10, 1.0, 1e10):
+            for dr in (2.0, 100.0, 1e4):
+                for ir in (0.0, 0.99, 1.0, 1.1):
+                    for vel in VELS:
+                        D = dr * R
+                        b = ir * R
+                        x = [D * tow[i] + b * perp[i] for i in range(3)]
+                        w = [-tow[i] for i in range(3)]
+                        v = [vel * perp[i] for i in range(3)]
+                        yield _case(cid, "moving", x, [0.0] * 3, w, v, R, 0.0)
+                        cid += 1
+    return cid
+
+
+def _broad_growing(cid):
+    for di, (tow, perp) in enumerate(DIRS):
+        for sign, tag in ((-1.0, "in"), (1.0, "out")):
+            for rr in RR_GROW:
+                for dr in (2.0, 10.0, 100.0, 1e4, 1e8):
+                    for ir in (0.0, 0.5, 0.99, 1.0, 1.1):
+                        D = dr * 4.0
+                        b = ir * 4.0
+                        x = [D * tow[i] + b * perp[i] for i in range(3)]
+                        w = [sign * tow[i] for i in range(3)]
+                        yield _case(cid, f"growing_{tag}", x, [0.0] * 3, w,
+                                    [0.0] * 3, 4.0, rr)
+                        cid += 1
+    return cid
+
+
+def _broad_shrinking(cid):
+    for (tow, perp) in (DIRS_AXIS, DIRS_OBLIQUE1):
+        for dr in (2.0, 10.0, 100.0):
+            for ir in (0.0, 0.99, 1.0):
+                D = dr * 10.0
+                b = ir * 10.0
+                x = [D * tow[i] + b * perp[i] for i in range(3)]
+                w = [-tow[i] for i in range(3)]
+                yield _case(cid, "shrinking", x, [0.0] * 3, w, [0.0] * 3, 10.0,
+                            0.1)
+                cid += 1
+    return cid
+
+
+def _materialize(gen, cases):
+    for c in gen:
+        cases.append(c)
+
+
+def _build_all():
+    cases = curated_kernel_cases()
+    cid = 1000
+    for gen in (_broad_fixed, _broad_moving, _broad_growing, _broad_shrinking):
+        gen_cases = []
+        _materialize(gen(cid), gen_cases)
+        if gen_cases:
+            cid = gen_cases[-1]["id"] + 1
+        cases.extend(gen_cases)
+    return cases
+
+
+# ---------------------------------------------------------------------------
+def _route(cid, category, t0, obs, direction, c0, v, R0, rr, model=0,
+           valid_t_min=-1e300):
+    return {
+        "id": cid,
+        "category": category,
+        "t0": t0,
+        "obs": list(obs),
+        "dir": list(direction),
+        "c0": list(c0),
+        "v": list(v),
+        "R0": R0,
+        "rr": rr,
+        "valid_t_min": valid_t_min,
+        "model": model,
+    }
+
+
+def route_cases():
+    cases = []
+    t63 = float.fromhex("0x1.fa8f5c28f5c29p+3")
+    t64 = float.fromhex("0x1.fb17e4b17e4b1p+3")
+    w63 = [-float.fromhex("0x1.d4afba4704cap-2"),
+           float.fromhex("0x1.c7378f8e872d1p-1"),
+           float.fromhex("0x1.15bad4e30e8ddp-8")]
+    w64 = [-float.fromhex("0x1.f98ae1a782104p-2"),
+           float.fromhex("0x1.bd3bb364ac492p-1"),
+           float.fromhex("0x1.102d2a1c6ac74p-7")]
+    # Production rows reproduced through the Alcubierre-style callback (model 1).
+    cases.append(_route(0, "real63", t63, [0.0, -24.0, 0.0], w63,
+                        [2.0 * t63, 0.0, 0.0], [2.0, 0.0, 0.0], 5.0, 0.0,
+                        model=1))
+    cases.append(_route(1, "real64", t64, [0.0, -24.0, 0.0], w64,
+                        [2.0 * t64, 0.0, 0.0], [2.0, 0.0, 0.0], 5.0, 0.0,
+                        model=1))
+    # Same rows through the input-stable callback (model 0).
+    cases.append(_route(2, "real63_stable", t63, [0.0, -24.0, 0.0], w63,
+                        [2.0 * t63, 0.0, 0.0], [2.0, 0.0, 0.0], 5.0, 0.0,
+                        model=0))
+    cases.append(_route(3, "real64_stable", t64, [0.0, -24.0, 0.0], w64,
+                        [2.0 * t64, 0.0, 0.0], [2.0, 0.0, 0.0], 5.0, 0.0,
+                        model=0))
+    R = 5.0
+    cases += [
+        _route(10, "headon_hit", 0.0, [10.0, 0.0, 0.0], [-1.0, 0.0, 0.0],
+               [0.0] * 3, [0.0] * 3, R, 0.0),
+        _route(11, "headon_miss", 0.0, [10.0, 0.0, 0.0], [1.0, 0.0, 0.0],
+               [0.0] * 3, [0.0] * 3, R, 0.0),
+        _route(12, "clear_miss", 0.0, [10.0, 6.0, 0.0], [-1.0, 0.0, 0.0],
+               [0.0] * 3, [0.0] * 3, R, 0.0),
+        _route(13, "grazing_in", 0.0, [10.0, math.nextafter(R, 0.0), 0.0],
+               [-1.0, 0.0, 0.0], [0.0] * 3, [0.0] * 3, R, 0.0),
+        _route(14, "exact_tangent", 0.0, [10.0, R, 0.0], [-1.0, 0.0, 0.0],
+               [0.0] * 3, [0.0] * 3, R, 0.0),
+        _route(15, "grazing_out", 0.0,
+               [10.0, math.nextafter(R, math.inf), 0.0], [-1.0, 0.0, 0.0],
+               [0.0] * 3, [0.0] * 3, R, 0.0),
+        _route(16, "camera_inside", 0.0, [2.0, 0.0, 0.0], [1.0, 0.0, 0.0],
+               [0.0] * 3, [0.0] * 3, R, 0.0),
+        _route(17, "near_boundary_hit", 0.0,
+               [math.nextafter(R, math.inf), 0.0, 0.0], [-1.0, 0.0, 0.0],
+               [0.0] * 3, [0.0] * 3, R, 0.0),
+    ]
+    # Broad fixed-sphere subset through the public route.
+    cid = 100
+    for (tow, perp) in (DIRS_AXIS, DIRS_OBLIQUE1):
+        for R0 in (1e-10, 1.0, 1e10):
+            for dr in (2.0, 10.0, 100.0, 512.0, 1024.0, 1e4):
+                for ir in (0.0, 0.5, 0.99, 1.0, 1.1):
+                    D = dr * R0
+                    b = ir * R0
+                    obs = [D * tow[i] + b * perp[i] for i in range(3)]
+                    direction = [-tow[i] for i in range(3)]
+                    cases.append(_route(cid, f"route_fixed_d{dr:g}", 0.0, obs,
+                                        direction, [0.0] * 3, [0.0] * 3, R0,
+                                        0.0))
+                    cid += 1
+    # Moving spheres (v = 2 tow), model 0.
+    for ir in (0.0, 0.99, 1.1):
+        cases.append(_route(cid, "route_moving", 0.0, [100.0, 0.0, 0.0],
+                            [-1.0, 0.0, 0.0], [0.0] * 3, [2.0, 0.0, 0.0], 5.0,
+                            0.0))
+        cid += 1
+    # Growing worldtubes (rr <= 0, positive radius on the whole past), inward
+    # and initially-outward photons, axis + oblique.  These mirror the
+    # growing_out kernel false-MISS family and exercise the public route where
+    # a kernel MISS bypasses the fallback entirely.  valid_t_min is far below
+    # any sampled time so no artificial history clip is introduced.
+    R0 = 4.0
+    for (tow, perp) in DIRS:
+        for sign, tag in ((-1.0, "in"), (1.0, "out")):
+            for rr in RR_GROW:
+                for dr in (2.0, 10.0, 100.0, 512.0, 1024.0, 1e4):
+                    for ir in (0.0, 1.0, 1.1):
+                        D = dr * R0
+                        b = ir * R0
+                        obs = [D * tow[i] + b * perp[i] for i in range(3)]
+                        direction = [sign * tow[i] for i in range(3)]
+                        cases.append(_route(cid, f"route_grow_{tag}_d{dr:g}",
+                                            0.0, obs, direction, [0.0] * 3,
+                                            [0.0] * 3, R0, rr))
+                        cid += 1
+    return cases
+
+
+# ---------------------------------------------------------------------------
+def _hex(x):
+    return float(x).hex()
+
+
+def _fmt(v):
+    return _hex(v)
+
+
+def write_header(kernel, route, out_path: Path):
+    lines = []
+    lines.append("/* Generated by benchmarks/quadratic_precision/cases.py. */\n")
+    lines.append("#ifndef QUADRATIC_CASES_H\n#define QUADRATIC_CASES_H\n")
+    lines.append("typedef struct {\n  int id;\n  const char *category;\n"
+                 "  double x[3], c[3], w[3], v[3];\n  double R0, rr;\n"
+                 "} QuadKernelCase;\n\n")
+    lines.append("typedef struct {\n  int id;\n  const char *category;\n"
+                 "  double t0;\n  double obs[3];\n  double dir[3];\n"
+                 "  double c0[3];\n  double v[3];\n  double R0, rr;\n"
+                 "  double valid_t_min;\n  int model;\n} QuadRouteCase;\n\n")
+    lines.append("static const QuadKernelCase quad_kernel_cases[] = {\n")
+    for c in kernel:
+        lines.append(
+            "  {.id=%d,.category=\"%s\","
+            ".x={%s,%s,%s},.c={%s,%s,%s},.w={%s,%s,%s},.v={%s,%s,%s},"
+            ".R0=%s,.rr=%s},\n"
+            % (c["id"], c["category"], *[_fmt(z) for z in c["x"]],
+               *[_fmt(z) for z in c["c"]], *[_fmt(z) for z in c["w"]],
+               *[_fmt(z) for z in c["v"]], _fmt(c["R0"]), _fmt(c["rr"])))
+    lines.append("};\n")
+    lines.append("static const int quad_kernel_case_count = %d;\n\n"
+                 % len(kernel))
+    lines.append("static const QuadRouteCase quad_route_cases[] = {\n")
+    for c in route:
+        lines.append(
+            "  {.id=%d,.category=\"%s\",.t0=%s,"
+            ".obs={%s,%s,%s},.dir={%s,%s,%s},.c0={%s,%s,%s},.v={%s,%s,%s},"
+            ".R0=%s,.rr=%s,.valid_t_min=%s,.model=%d},\n"
+            % (c["id"], c["category"], _fmt(c["t0"]),
+               *[_fmt(z) for z in c["obs"]], *[_fmt(z) for z in c["dir"]],
+               *[_fmt(z) for z in c["c0"]], *[_fmt(z) for z in c["v"]],
+               _fmt(c["R0"]), _fmt(c["rr"]), _fmt(c["valid_t_min"]),
+               c["model"]))
+    lines.append("};\n")
+    lines.append("static const int quad_route_case_count = %d;\n"
+                 % len(route))
+    lines.append("#endif\n")
+    out_path.write_text("".join(lines))
+
+
+def write_json(kernel, route, out_path: Path):
+    def enc(c):
+        d = {}
+        for k, v in c.items():
+            if isinstance(v, float):
+                d[k] = v.hex()
+            elif isinstance(v, list):
+                d[k] = [z.hex() for z in v]
+            else:
+                d[k] = v
+        return d
+
+    out_path.write_text(
+        json.dumps(
+            {"deterministic": True, "kernel": [enc(c) for c in kernel],
+             "route": [enc(c) for c in route]},
+            indent=1,
+        )
+        + "\n")
+
+
+def build():
+    kernel = _build_all()
+    route = route_cases()
+    return kernel, route
+
+
+if __name__ == "__main__":
+    k, r = build()
+    print(f"kernel cases: {len(k)}, route cases: {len(r)}")
+    cats = {}
+    for c in k:
+        cats[c["category"]] = cats.get(c["category"], 0) + 1
+    for key in sorted(cats):
+        print(f"  {key}: {cats[key]}")
+    rcats = {}
+    for c in r:
+        rcats[c["category"]] = rcats.get(c["category"], 0) + 1
+    print("route categories:")
+    for key in sorted(rcats):
+        print(f"  {key}: {rcats[key]}")

benchmarks/quadratic_precision/fixtures/probe.c

@@ -0,0 +1,614 @@
+/* Quadratic precision benchmark probe.
+ *
+ * Compiled once per generated variant, with
+ *   -DPROBE_VARIANT_SOURCE=".../ld_fma.c" -DPROBE_VARIANT_NAME="ld_fma"
+ * The probe #includes the generated full module, so the private static
+ * entry_quadratic_coeffs / entry_solve and the public asymptotic_route_camera
+ * are the *actual generated* code.  No production file is modified.
+ *
+ * Modes:
+ *   accuracy   -- run every kernel case through the generated kernel and every
+ *                 route case through the generated public pre-route; emit CSVs.
+ *   microbench -- timed repeated kernel assembly + entry_solve (serial).
+ *   routebench -- timed repeated public pre-route (OpenMP, static schedule).
+ */
+#include PROBE_VARIANT_SOURCE
+
+#include <float.h>
+#include <math.h>
+#include <omp.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <time.h>
+
+#include "quadratic_cases.h"
+
+#ifndef PROBE_VARIANT_NAME
+#define PROBE_VARIANT_NAME "unknown"
+#endif
+
+/* ------------------------------------------------------------------ */
+/* Exact-value coefficient printing: long double uses %La (full mantissa),
+ * double uses %a.  _Generic selects the right printer without knowing the
+ * generated EntryQuadratic field type. */
+static void coeff_ld(FILE *f, long double v) { fprintf(f, "%La", v); }
+static void coeff_d(FILE *f, double v) { fprintf(f, "%a", v); }
+#define PRINT_COEFF(f, v) \
+  _Generic((v), long double : coeff_ld, double : coeff_d)((f), (v))
+
+/* ------------------------------------------------------------------ */
+/* Fixture source: flat identity metric, one Minkowski end whose worldtube is
+ * a constant-velocity (optionally linearly growing/shrinking) sphere.
+ *
+ *   model 0 (input-stable): center = c0 + v*(t - t0), R = R0 + rr*(t - t0)
+ *   model 1 (production-style): center = v*t,        R = R0 + rr*t
+ *
+ * model 0 evaluates exactly at the camera: center(t0) == c0, R(t0) == R0, so
+ * the kernel inputs are the frozen case values.  model 1 mirrors the
+ * Alcubierre-style callback (used only with rr == 0). */
+typedef struct {
+  double c0[3];
+  double v[3];
+  double R0;
+  double rr;
+  double t0;
+  double valid_t_min;
+  int model;
+} FixtureContext;
+
+static double fixture_center(const FixtureContext *ctx, int i, double t) {
+  if (ctx->model == 0)
+    return ctx->c0[i] + ctx->v[i] * (t - ctx->t0);
+  return ctx->v[i] * t;
+}
+
+static double fixture_radius(const FixtureContext *ctx, double t) {
+  if (ctx->model == 0)
+    return ctx->R0 + ctx->rr * (t - ctx->t0);
+  return ctx->R0 + ctx->rr * t;
+}
+
+static SpacetimePointStatus fixture_eval(const SpacetimeSource *source,
+                                         double t, const double x[3],
+                                         MetricData *metric) {
+  (void)source;
+  (void)t;
+  (void)x;
+  *metric = (MetricData){.alpha = 1.0,
+                         .gamma = {{1.0, 0.0, 0.0},
+                                   {0.0, 1.0, 0.0},
+                                   {0.0, 0.0, 1.0}}};
+  return SPACETIME_POINT_OK;
+}
+
+static SpacetimeRayStatus fixture_classify(const SpacetimeSource *source,
+                                           double t, const double x[3]) {
+  (void)source;
+  (void)t;
+  (void)x;
+  return SPACETIME_RAY_ACTIVE;
+}
+
+static size_t fixture_end_count(const SpacetimeSource *source) {
+  (void)source;
+  return 1;
+}
+
+static int fixture_end(const SpacetimeSource *source, size_t index,
+                       SpacetimeAsymptoticEnd *out) {
+  const FixtureContext *ctx = source->context;
+  if (index != 0)
+    return -1;
+  *out = (SpacetimeAsymptoticEnd){
+      .end_id = 0,
+      .exterior_kind = ASYMPTOTIC_EXTERIOR_MINKOWSKI,
+      .mass = 0.0,
+      .frame_origin = {0.0, 0.0, 0.0},
+      .frame_axes = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}};
+  return 0;
+}
+
+static int fixture_worldtube(const SpacetimeSource *source,
+                             SpacetimeEndId end_id, double t,
+                             SpacetimeEscapeWorldtubeSample *out) {
+  const FixtureContext *ctx = source->context;
+  if (end_id != 0)
+    return -1;
+  if (!isfinite(t) || t < ctx->valid_t_min) {
+    *out = (SpacetimeEscapeWorldtubeSample){.valid = 0};
+    return 0;
+  }
+  *out = (SpacetimeEscapeWorldtubeSample){
+      .center = {fixture_center(ctx, 0, t), fixture_center(ctx, 1, t),
+                 fixture_center(ctx, 2, t)},
+      .velocity = {ctx->v[0], ctx->v[1], ctx->v[2]},
+      .radius = fixture_radius(ctx, t),
+      .radius_rate = ctx->rr,
+      .velocity_constant = 1,
+      .valid = 1};
+  return 0;
+}
+
+static double fixture_next_segment(const SpacetimeSource *source,
+                                   SpacetimeEndId end_id, double t) {
+  (void)source;
+  (void)end_id;
+  (void)t;
+  return NAN;
+}
+
+static void fixture_destroy(SpacetimeSource *source) {
+  source->context = NULL;
+  source->ops = NULL;
+}
+
+static const SpacetimeOps fixture_ops = {
+    .eval = fixture_eval,
+    .classify = fixture_classify,
+    .asymptotic_end_count = fixture_end_count,
+    .asymptotic_end = fixture_end,
+    .escape_worldtube_sample = fixture_worldtube,
+    .escape_worldtube_next_segment = fixture_next_segment,
+    .destroy = fixture_destroy,
+};
+
+/* ------------------------------------------------------------------ */
+static void print_environment(void) {
+  fprintf(stdout,
+          "{\"variant\":\"%s\",\"sizeof_long_double\":%zu,\"LDBL_MANT_DIG\":%d,"
+          "\"DBL_MANT_DIG\":%d,\"LDBL_MAX_EXP\":%d,\"hardware_threads\":%d,"
+          "\"omp_max_threads\":%d}\n",
+          PROBE_VARIANT_NAME, sizeof(long double), LDBL_MANT_DIG, DBL_MANT_DIG,
+          LDBL_MAX_EXP, omp_get_num_procs(), omp_get_max_threads());
+}
+
+/* ------------------------------------------------------------------ */
+/* accuracy mode */
+static void mode_accuracy(const char *kernel_out, const char *route_out) {
+  FILE *kf = fopen(kernel_out, "w");
+  if (kf == NULL) {
+    fprintf(stderr, "cannot open %s\n", kernel_out);
+    exit(2);
+  }
+  fprintf(kf,
+          "id,category,status,sigma,x0,x1,x2,c0,c1,c2,w0,w1,w2,v0,v1,v2,R0,rr,"
+          "a,b,c\n");
+  for (int i = 0; i < quad_kernel_case_count; ++i) {
+    const QuadKernelCase *c = &quad_kernel_cases[i];
+    EntryQuadratic k = entry_quadratic_coeffs(c->x, c->c, c->w, c->v, c->R0,
+                                              c->rr);
+    double s = -1.0;
+    EntrySolveResult r = entry_solve(&k, &s);
+    fprintf(kf, "%d,%s,%d,%a", c->id, c->category, (int)r, s);
+    for (int j = 0; j < 3; ++j)
+      fprintf(kf, ",%a", c->x[j]);
+    for (int j = 0; j < 3; ++j)
+      fprintf(kf, ",%a", c->c[j]);
+    for (int j = 0; j < 3; ++j)
+      fprintf(kf, ",%a", c->w[j]);
+    for (int j = 0; j < 3; ++j)
+      fprintf(kf, ",%a", c->v[j]);
+    fprintf(kf, ",%a,%a,", c->R0, c->rr);
+    PRINT_COEFF(kf, k.a);
+    fputc(',', kf);
+    PRINT_COEFF(kf, k.b);
+    fputc(',', kf);
+    PRINT_COEFF(kf, k.c);
+    fputc('\n', kf);
+  }
+  fclose(kf);
+
+  FILE *rf = fopen(route_out, "w");
+  if (rf == NULL) {
+    fprintf(stderr, "cannot open %s\n", route_out);
+    exit(2);
+  }
+  fprintf(rf,
+          "id,category,status,kind,failure_reason,fallback_evals,pi_match,"
+          "lcam_match,F,tol,why,x0,x1,x2,Pi0,Pi1,Pi2,ninf0,ninf1,ninf2,"
+          "canon_t,canonx0,canonx1,canonx2,canonw0,canonw1,canonw2,"
+          "cbx0,cbx1,cbx2,cbr,cbok,ct00,ct01,ct02,rt0,cb0ok,"
+          "kern_ok,kern_status,kern_sigma,kern_a,kern_b,kern_c,"
+          "logcamera,logentry,"
+          "t0,oc0,oc1,oc2,d0,d1,d2,fc0,fc1,fc2,v0,v1,v2,"
+          "R0,rr,model,reason_name\n");
+  for (int i = 0; i < quad_route_case_count; ++i) {
+    const QuadRouteCase *c = &quad_route_cases[i];
+    FixtureContext ctx = {.c0 = {c->c0[0], c->c0[1], c->c0[2]},
+                          .v = {c->v[0], c->v[1], c->v[2]},
+                          .R0 = c->R0,
+                          .rr = c->rr,
+                          .t0 = c->t0,
+                          .valid_t_min = c->valid_t_min,
+                          .model = c->model};
+    SpacetimeSource source = {.ops = &fixture_ops, .context = &ctx};
+    ObserverState obs = {0};
+    obs.coordinate_time = c->t0;
+    for (int j = 0; j < 3; ++j)
+      obs.coordinate_position[j] = c->obs[j];
+    obs.tetrad[0][0] = 1.0;
+    obs.tetrad[1][1] = 1.0;
+    obs.tetrad[2][2] = 1.0;
+    obs.tetrad[3][3] = 1.0;
+
+    int canon_ok = 0;
+    AsymptoticPhotonState canon = {0};
+    GeodesicRayState st = {0};
+    {
+      MetricData metric;
+      if (spacetime_eval(&source, c->t0, obs.coordinate_position, &metric) ==
+              SPACETIME_POINT_OK &&
+          geodesic_initialize_past_ray_metric(&metric, &obs, c->dir, &st) == 0 &&
+          asymptotic_canonical_from_backend(&source, 0, &metric, c->t0, st.x,
+                                            st.Pi, st.log_alpha_p0,
+                                            &canon) == 0)
+        canon_ok = 1;
+    }
+
+    AsymptoticRoute route;
+    AsymptoticStatus status =
+        asymptotic_route_camera(&source, &obs, c->dir, &route);
+
+    int pi_match = 1, lcam_match = 1;
+    if (canon_ok && status == ASYMPTOTIC_OK &&
+        route.kind == ASYMPTOTIC_ROUTE_ENTRY) {
+      for (int j = 0; j < 3; ++j)
+        if (route.Pi[j] != st.Pi[j])
+          pi_match = 0;
+      if (route.log_alpha_p0_camera != st.log_alpha_p0)
+        lcam_match = 0;
+    }
+
+    double F = NAN, tol = NAN;
+    RayReason why = RAY_REASON_NONE;
+    if (status == ASYMPTOTIC_OK && route.kind == ASYMPTOTIC_ROUTE_ENTRY)
+      asymptotic_entry_geometry(&source, route.end_id, route.activate_t,
+                                route.x, &F, &tol, &why);
+
+    SpacetimeEscapeWorldtubeSample cb = {0};
+    int cbok = 0;
+    {
+      const double tt = (status == ASYMPTOTIC_OK &&
+                         route.kind == ASYMPTOTIC_ROUTE_ENTRY)
+                            ? route.activate_t
+                            : c->t0;
+      if (spacetime_escape_worldtube_sample(&source, route.end_id, tt, &cb) ==
+              0 &&
+          cb.valid)
+        cbok = 1;
+    }
+    /* Worldtube sample at the segment start used by the quadratic kernel. */
+    SpacetimeEscapeWorldtubeSample cb0 = {0};
+    int cb0ok = 0;
+    if (spacetime_escape_worldtube_sample(&source, route.end_id, c->t0,
+                                          &cb0) == 0 &&
+        cb0.valid)
+      cb0ok = 1;
+
+    /* Reproduce the first-segment public kernel classification at the actual
+     * canonical inputs, so a reference ENTER / kernel MISS / route ESCAPED is
+     * directly visible and not confused with canonical normalisation. */
+    int kern_ok = 0;
+    int kern_status = -1;
+    double kern_sigma = -1.0;
+    EntryQuadratic kk = {0};
+    SpacetimeAsymptoticEnd end_desc;
+    if (canon_ok && cb0ok &&
+        spacetime_asymptotic_end(&source, 0, &end_desc) == 0) {
+      double c_frame[3], v_frame[3];
+      backend_position_to_frame(&end_desc, cb0.center, c_frame);
+      backend_vector_to_frame(&end_desc, cb0.velocity, v_frame);
+      kk = entry_quadratic_coeffs(canon.x, c_frame, canon.w, v_frame,
+                                  cb0.radius, cb0.radius_rate);
+      kern_status = (int)entry_solve(&kk, &kern_sigma);
+      kern_ok = 1;
+    }
+
+    fprintf(rf, "%d,%s,%d,%d,%d,%u,%d,%d,%a,%a,%d", c->id, c->category,
+            (int)status, (int)route.kind, (int)route.failure_reason,
+            route.entry_fallback_evaluations, pi_match, lcam_match, F, tol,
+            (int)why);
+    for (int j = 0; j < 3; ++j)
+      fprintf(rf, ",%a", route.x[j]);
+    for (int j = 0; j < 3; ++j)
+      fprintf(rf, ",%a", route.Pi[j]);
+    for (int j = 0; j < 3; ++j)
+      fprintf(rf, ",%a", route.n_infinity[j]);
+    fprintf(rf, ",%a", canon.t);
+    for (int j = 0; j < 3; ++j)
+      fprintf(rf, ",%a", canon.x[j]);
+    for (int j = 0; j < 3; ++j)
+      fprintf(rf, ",%a", canon.w[j]);
+    for (int j = 0; j < 3; ++j)
+      fprintf(rf, ",%a", cb.center[j]);
+    fprintf(rf, ",%a,%d", cb.radius, cbok);
+    for (int j = 0; j < 3; ++j)
+      fprintf(rf, ",%a", cb0.center[j]);
+    fprintf(rf, ",%a,%d", cb0.radius, cb0ok);
+    fprintf(rf, ",%d,%d,%a,", kern_ok, kern_status, kern_sigma);
+    PRINT_COEFF(rf, kk.a);
+    fputc(',', rf);
+    PRINT_COEFF(rf, kk.b);
+    fputc(',', rf);
+    PRINT_COEFF(rf, kk.c);
+    fprintf(rf, ",%a,%a", route.log_alpha_p0_camera, route.log_alpha_p0);
+    fprintf(rf, ",%a", c->t0);
+    for (int j = 0; j < 3; ++j)
+      fprintf(rf, ",%a", c->obs[j]);
+    for (int j = 0; j < 3; ++j)
+      fprintf(rf, ",%a", c->dir[j]);
+    for (int j = 0; j < 3; ++j)
+      fprintf(rf, ",%a", c->c0[j]);
+    for (int j = 0; j < 3; ++j)
+      fprintf(rf, ",%a", c->v[j]);
+    fprintf(rf, ",%a,%a,%d,%s\n", c->R0, c->rr, c->model,
+            ray_reason_name(route.failure_reason));
+  }
+  fclose(rf);
+}
+
+/* ------------------------------------------------------------------ */
+/* microbench mode: coefficient assembly + entry_solve, serial. */
+static volatile double g_kernel_sink;
+
+static __attribute__((noinline)) double kernel_batch(const QuadKernelCase *cs,
+                                                     int n, long reps) {
+  double acc = 0.0;
+  for (long r = 0; r < reps; ++r) {
+    for (int i = 0; i < n; ++i) {
+      /* Force the index through an opaque register so the compiler cannot
+       * hoist the pure coefficient solve out of the repetition loop or prove
+       * the loaded fixture invariant.  Identical for every variant. */
+      int idx = i;
+      __asm__ __volatile__("" : "+r"(idx) : : "memory");
+      const QuadKernelCase *c = cs + idx;
+      EntryQuadratic k = entry_quadratic_coeffs(c->x, c->c, c->w, c->v,
+                                                c->R0, c->rr);
+      double s = 0.0;
+      acc += (double)entry_solve(&k, &s) + s * 1e-300;
+    }
+  }
+  return acc;
+}
+
+static void mode_microbench(long target_calls, const char *out) {
+  const int n = quad_kernel_case_count;
+  long reps = target_calls / n;
+  if (reps < 1)
+    reps = 1;
+  /* Warm-up outside the timed region. */
+  g_kernel_sink += kernel_batch(quad_kernel_cases, n, 1);
+  const double t0 = omp_get_wtime();
+  const double acc = kernel_batch(quad_kernel_cases, n, reps);
+  const double t1 = omp_get_wtime();
+  g_kernel_sink += acc;
+  const long calls = (long)n * reps;
+  const double seconds = t1 - t0;
+  FILE *f = fopen(out, "w");
+  if (f == NULL) {
+    fprintf(stderr, "cannot open %s\n", out);
+    exit(2);
+  }
+  fprintf(f,
+          "{\"variant\":\"%s\",\"mode\":\"microbench\",\"cases\":%d,"
+          "\"reps\":%ld,\"calls\":%ld,\"seconds\":%.9f,\"ns_per_call\":%.6f,"
+          "\"sink\":%.17g}\n",
+          PROBE_VARIANT_NAME, n, reps, calls, seconds,
+          seconds * 1e9 / (double)calls, g_kernel_sink);
+  fclose(f);
+  fprintf(stdout, "microbench %s: %ld calls in %.6f s (%.2f ns/call)\n",
+          PROBE_VARIANT_NAME, calls, seconds, seconds * 1e9 / (double)calls);
+}
+
+/* ------------------------------------------------------------------ */
+/* routebench mode: public asymptotic_route_camera, OpenMP static. */
+typedef struct {
+  FixtureContext *ctxs;
+  SpacetimeSource *srcs;
+  ObserverState *obss;
+  const QuadRouteCase *rcs;
+  int n;
+} Preloaded;
+
+typedef struct {
+  long entry, escaped, inside, time_exhausted, invalid, unsupported, other;
+  long fallback_sum;
+} RouteCounts;
+
+static volatile double g_route_sink;
+
+static __attribute__((noinline)) void route_batch(const Preloaded *p, long reps,
+                                                  int threads, double *seconds,
+                                                  RouteCounts *counts) {
+  const int n = p->n;
+  long entry = 0, escaped = 0, inside = 0, texh = 0, inv = 0, unsup = 0,
+       other = 0, fallback = 0;
+  const long total = (long)n * reps;
+  const double t0 = omp_get_wtime();
+#pragma omp parallel num_threads(threads) reduction(+ : entry, escaped, inside, texh, inv, unsup, other, fallback)
+  {
+#pragma omp for schedule(static)
+    for (long k = 0; k < total; ++k) {
+      long kk = k;
+      __asm__ __volatile__("" : "+r"(kk) : : "memory");
+      const int i = (int)(kk % n);
+      AsymptoticRoute route;
+      const AsymptoticStatus st = asymptotic_route_camera(
+          &p->srcs[i], &p->obss[i], p->rcs[i].dir, &route);
+      fallback += (long)route.entry_fallback_evaluations;
+      if (st == ASYMPTOTIC_OK) {
+        switch (route.kind) {
+        case ASYMPTOTIC_ROUTE_ENTRY:
+          ++entry;
+          break;
+        case ASYMPTOTIC_ROUTE_ESCAPED:
+          ++escaped;
+          break;
+        case ASYMPTOTIC_ROUTE_INSIDE:
+          ++inside;
+          break;
+        case ASYMPTOTIC_ROUTE_TIME_RANGE_EXHAUSTED:
+          ++texh;
+          break;
+        default:
+          ++other;
+          break;
+        }
+      } else if (st == ASYMPTOTIC_UNSUPPORTED) {
+        ++unsup;
+      } else {
+        ++inv;
+      }
+    }
+  }
+  *seconds = omp_get_wtime() - t0;
+  counts->entry = entry;
+  counts->escaped = escaped;
+  counts->inside = inside;
+  counts->time_exhausted = texh;
+  counts->invalid = inv;
+  counts->unsupported = unsup;
+  counts->other = other;
+  counts->fallback_sum = fallback;
+}
+
+static void mode_routebench(long target_calls, int threads, double max_seconds,
+                            const char *out) {
+  const int n = quad_route_case_count;
+  Preloaded p;
+  p.n = n;
+  p.rcs = quad_route_cases;
+  p.ctxs = malloc(sizeof(FixtureContext) * (size_t)n);
+  p.srcs = malloc(sizeof(SpacetimeSource) * (size_t)n);
+  p.obss = malloc(sizeof(ObserverState) * (size_t)n);
+  if (!p.ctxs || !p.srcs || !p.obss) {
+    fprintf(stderr, "allocation failure\n");
+    exit(2);
+  }
+  for (int i = 0; i < n; ++i) {
+    const QuadRouteCase *c = &quad_route_cases[i];
+    p.ctxs[i] = (FixtureContext){.c0 = {c->c0[0], c->c0[1], c->c0[2]},
+                                 .v = {c->v[0], c->v[1], c->v[2]},
+                                 .R0 = c->R0,
+                                 .rr = c->rr,
+                                 .t0 = c->t0,
+                                 .valid_t_min = c->valid_t_min,
+                                 .model = c->model};
+    p.srcs[i] = (SpacetimeSource){.ops = &fixture_ops, .context = &p.ctxs[i]};
+    p.obss[i] = (ObserverState){0};
+    p.obss[i].coordinate_time = c->t0;
+    for (int j = 0; j < 3; ++j)
+      p.obss[i].coordinate_position[j] = c->obs[j];
+    p.obss[i].tetrad[0][0] = 1.0;
+    p.obss[i].tetrad[1][1] = 1.0;
+    p.obss[i].tetrad[2][2] = 1.0;
+    p.obss[i].tetrad[3][3] = 1.0;
+  }
+
+  /* Calibrate with one pass, then size reps to the call/time budgets. */
+  double cal_seconds = 0.0;
+  RouteCounts cal_counts;
+  route_batch(&p, 1, threads, &cal_seconds, &cal_counts);
+  const double per_call = cal_seconds / (double)n;
+  long reps = target_calls / n;
+  if (reps < 1)
+    reps = 1;
+  if (per_call > 0.0) {
+    const long by_time = (long)(max_seconds / (per_call * (double)n));
+    if (by_time < 1)
+      reps = 1;
+    else if (reps > by_time)
+      reps = by_time;
+  }
+
+  double seconds = 0.0;
+  RouteCounts counts;
+  route_batch(&p, reps, threads, &seconds, &counts);
+  g_route_sink += (double)counts.entry + (double)counts.escaped;
+
+  FILE *f = fopen(out, "w");
+  if (f == NULL) {
+    fprintf(stderr, "cannot open %s\n", out);
+    exit(2);
+  }
+  fprintf(f,
+          "{\"variant\":\"%s\",\"mode\":\"routebench\",\"cases\":%d,"
+          "\"reps\":%ld,\"calls\":%ld,\"threads\":%d,\"cal_seconds\":%.9f,"
+          "\"seconds\":%.9f,\"ns_per_call\":%.6f,\"entry\":%ld,\"escaped\":%ld,"
+          "\"inside\":%ld,\"time_exhausted\":%ld,\"invalid\":%ld,"
+          "\"unsupported\":%ld,\"other\":%ld,\"fallback_sum\":%ld}\n",
+          PROBE_VARIANT_NAME, n, reps, (long)n * reps, threads, cal_seconds,
+          seconds, seconds * 1e9 / (double)((long)n * reps), counts.entry,
+          counts.escaped, counts.inside, counts.time_exhausted, counts.invalid,
+          counts.unsupported, counts.other, counts.fallback_sum);
+  fclose(f);
+  fprintf(stdout,
+          "routebench %s: %ld calls in %.6f s on %d threads (%.2f ns/call)\n",
+          PROBE_VARIANT_NAME, (long)n * reps, seconds, threads,
+          seconds * 1e9 / (double)((long)n * reps));
+  free(p.ctxs);
+  free(p.srcs);
+  free(p.obss);
+}
+
+/* ------------------------------------------------------------------ */
+static void usage(const char *argv0) {
+  fprintf(stderr,
+          "usage:\n"
+          "  %s accuracy  --kernel-out K.csv --route-out R.csv\n"
+          "  %s microbench --out F.json [--target-calls N]\n"
+          "  %s routebench --out F.json [--target-calls N] [--threads T]"
+          " [--max-seconds S]\n",
+          argv0, argv0, argv0);
+}
+
+static const char *arg_value(int argc, char **argv, const char *flag) {
+  for (int i = 1; i + 1 < argc; ++i)
+    if (strcmp(argv[i], flag) == 0)
+      return argv[i + 1];
+  return NULL;
+}
+
+int main(int argc, char **argv) {
+  print_environment();
+  if (argc < 2) {
+    usage(argv[0]);
+    return 1;
+  }
+  if (strcmp(argv[1], "accuracy") == 0) {
+    const char *k = arg_value(argc, argv, "--kernel-out");
+    const char *r = arg_value(argc, argv, "--route-out");
+    if (!k || !r) {
+      usage(argv[0]);
+      return 1;
+    }
+    mode_accuracy(k, r);
+    return 0;
+  }
+  if (strcmp(argv[1], "microbench") == 0) {
+    const char *out = arg_value(argc, argv, "--out");
+    const char *tc = arg_value(argc, argv, "--target-calls");
+    if (!out) {
+      usage(argv[0]);
+      return 1;
+    }
+    mode_microbench(tc ? atol(tc) : 10000000L, out);
+    return 0;
+  }
+  if (strcmp(argv[1], "routebench") == 0) {
+    const char *out = arg_value(argc, argv, "--out");
+    const char *tc = arg_value(argc, argv, "--target-calls");
+    const char *th = arg_value(argc, argv, "--threads");
+    const char *ms = arg_value(argc, argv, "--max-seconds");
+    if (!out) {
+      usage(argv[0]);
+      return 1;
+    }
+    mode_routebench(tc ? atol(tc) : 1000000L, th ? atoi(th) : 4,
+                    ms ? atof(ms) : 20.0, out);
+    return 0;
+  }
+  usage(argv[0]);
+  return 1;
+}

benchmarks/quadratic_precision/oracle.py

@@ -0,0 +1,571 @@
+#!/usr/bin/env python3
+"""High-precision reference for the quadratic entry benchmark.
+
+Inputs are exact IEEE-754 values, so they are represented exactly as
+``fractions.Fraction``.  The relative-distance quadratic coefficients, the
+discriminant and the polynomial residuals are therefore *exact* rationals; only
+the square root needs ``decimal`` (160 significant digits by default).  This
+avoids ``float.fromhex`` (which would silently drop a long-double coefficient to
+53 bits) and keeps the sign of a near-zero discriminant exact.
+
+Reference classification mirrors the production contract:
+  * c < 0  -> camera already INSIDE (never an entry failure)
+  * c == 0 -> boundary slope b decides; b < 0 enters at once, b == 0 with a < 0
+              enters later, otherwise no crossing
+  * c > 0  -> smallest positive root with inward slope; a missing/tangent root
+              is MISS.  ``a == 0`` is the exact linear branch.
+"""
+from __future__ import annotations
+
+import csv
+import json
+import math
+from decimal import Decimal, localcontext
+from fractions import Fraction
+
+DEFAULT_PREC = 160
+
+
+def parse_hex(s: str) -> Fraction:
+    """Exact rational value of a C ``%a`` / ``%La`` hex float literal."""
+    s = s.strip()
+    low = s.lower()
+    if "nan" in low:
+        raise ValueError("NaN input")
+    if "inf" in low:
+        return Fraction(0)  # only used for flags; never present in fixtures
+    neg = False
+    if s and s[0] in "+-":
+        neg = s[0] == "-"
+        s = s[1:]
+    if s[:2].lower() == "0x":
+        s = s[2:]
+    mant, _, exp_s = s.partition("p")
+    if not exp_s:
+        mant, _, exp_s = s.partition("P")
+    exp = int(exp_s) if exp_s else 0
+    ip, _, fp = mant.partition(".")
+    digits = (ip + fp) or "0"
+    val = Fraction(int(digits, 16), 1)
+    shift = exp - 4 * len(fp)
+    if shift >= 0:
+        val *= Fraction(2) ** shift
+    else:
+        val /= Fraction(2) ** (-shift)
+    return -val if neg else val
+
+
+def frac_to_dec(fr: Fraction, prec: int = DEFAULT_PREC) -> Decimal:
+    with localcontext() as ctx:
+        ctx.prec = prec
+        return Decimal(fr.numerator) / Decimal(fr.denominator)
+
+
+def hex_to_dec(s: str, prec: int = DEFAULT_PREC) -> Decimal:
+    return frac_to_dec(parse_hex(s), prec)
+
+
+def classify(a: Fraction, b: Fraction, c: Fraction, prec: int,
+             R0: Fraction | None = None, rr: Fraction | None = None):
+    """Return (status, root_decimal_or_None, info).
+
+    ``R0``/``rr`` (radius at the segment start and dR/dt) enforce the physical
+    positive-radius domain R(s) = R0 - rr*s > 0 along the past parameter.  A
+    positive root outside that domain is not a physical entry: it is reported
+    as MISS with ``info['domain_clipped']`` set (shrinking worldtubes).
+    """
+    info = {"a_zero": a == 0, "disc_sign": 0, "domain_clipped": False}
+
+    def domain_ok(root: Fraction) -> bool:
+        if R0 is None:
+            return True
+        return (R0 - (rr if rr is not None else Fraction(0)) * root) > 0
+
+    if c < 0:
+        return "INSIDE", None, info
+    if c == 0:
+        if b < 0:
+            return "ENTER", Decimal(0), info
+        if b == 0:
+            if a < 0:
+                return "ENTER", Decimal(0), info
+            return "MISS", None, info
+        if a < 0:
+            root = -b / a
+            if root > 0 and domain_ok(root):
+                return "ENTER", frac_to_dec(root, prec), info
+            if root > 0:
+                info["domain_clipped"] = True
+            return "MISS", None, info
+        return "MISS", None, info
+    if a == 0:
+        if b < 0:
+            root = -c / b
+            if root > 0:
+                if domain_ok(root):
+                    return "ENTER", frac_to_dec(root, prec), info
+                info["domain_clipped"] = True
+        return "MISS", None, info
+    disc = b * b - 4 * a * c
+    info["disc_sign"] = (disc > 0) - (disc < 0)
+    if disc < 0:
+        return "MISS", None, info
+    if disc == 0:
+        return "MISS", None, info  # tangency is not a crossing
+    with localcontext() as ctx:
+        ctx.prec = prec
+        sd = frac_to_dec(disc, prec).sqrt()
+        ad = frac_to_dec(a, prec)
+        bd = frac_to_dec(b, prec)
+        r1 = (-bd - sd) / (2 * ad)
+        r2 = (-bd + sd) / (2 * ad)
+        pos = sorted(r for r in (r1, r2) if r > 0)
+        for r in pos:
+            if 2 * ad * r + bd < 0:  # inward (outside -> inside) slope
+                rfr = _dec_to_frac_snapshot(r)
+                if domain_ok(rfr):
+                    return "ENTER", r, info
+                info["domain_clipped"] = True
+        return "MISS", None, info
+
+
+def _dec_to_frac_snapshot(d: Decimal) -> Fraction:
+    return Fraction(d)
+
+
+def coeffs(x, c, w, v, R0, rr):
+    d = [x[i] - c[i] for i in range(3)]
+    q = [w[i] + v[i] for i in range(3)]
+    qq = sum(q[i] * q[i] for i in range(3))
+    dq = sum(d[i] * q[i] for i in range(3))
+    dd = sum(d[i] * d[i] for i in range(3))
+    a = qq - rr * rr
+    b = 2 * (dq + R0 * rr)
+    cq = dd - R0 * R0
+    return a, b, cq, qq
+
+
+def load_cases(path):
+    data = json.loads(open(path).read())
+    kernel = {}
+    for c in data["kernel"]:
+        kernel[c["id"]] = {
+            "category": c["category"],
+            "x": [parse_hex(z) for z in c["x"]],
+            "c": [parse_hex(z) for z in c["c"]],
+            "w": [parse_hex(z) for z in c["w"]],
+            "v": [parse_hex(z) for z in c["v"]],
+            "R0": parse_hex(c["R0"]),
+            "rr": parse_hex(c["rr"]),
+        }
+    route = {c["id"]: c for c in data["route"]}
+    return kernel, route
+
+
+def kernel_reference(case, prec):
+    a, b, cq, qq = coeffs(case["x"], case["c"], case["w"], case["v"],
+                          case["R0"], case["rr"])
+    status, root, info = classify(a, b, cq, prec, case["R0"], case["rr"])
+    dd = sum(xi * xi for xi in
+             (case["x"][i] - case["c"][i] for i in range(3)))
+    return {
+        "status": status,
+        "root": root,
+        "a": a,
+        "b": b,
+        "c": cq,
+        "qq": qq,
+        "dd": dd,
+        "R0": case["R0"],
+        "rr": case["rr"],
+        "a_zero": info["a_zero"],
+        "disc_sign": info["disc_sign"],
+        "domain_clipped": info["domain_clipped"],
+    }
+
+
+def ulp_dec(cand: float) -> Decimal:
+    return frac_to_dec(Fraction(math.ulp(cand)))
+
+
+def dec_of_float(f: float) -> Decimal:
+    return frac_to_dec(Fraction(f), 300)
+
+
+def poly_resid(a: Fraction, b: Fraction, c: Fraction, s: Fraction) -> Fraction:
+    return a * s * s + b * s + c
+
+
+def poly_scale(a: Fraction, b: Fraction, c: Fraction, s: Fraction) -> Fraction:
+    return abs(a) * s * s + abs(b) * s + abs(c)
+
+
+def analyze_kernel_variant(rows, refs, prec):
+    """Return per-case records for one variant's kernel CSV."""
+    out = []
+    for r in rows:
+        cid = int(r["id"])
+        ref = refs[cid]
+        status = int(r["status"])
+        cand = float.fromhex(r["sigma"]) if r["sigma"] not in ("", "nan") else float("nan")
+        cand_fr = parse_hex(r["sigma"]) if r["sigma"] not in ("", "nan") else None
+        ak = parse_hex(r["a"])
+        bk = parse_hex(r["b"])
+        ck = parse_hex(r["c"])
+        rec = {
+            "id": cid,
+            "category": r["category"],
+            "variant_status": status,
+            "ref_status": ref["status"],
+            "a_zero_ref": ref["a_zero"],
+            "a_zero_kernel": ak == 0,
+            "a_sign_ref": _sign(ref["a"]),
+            "a_sign_kernel": _sign(ak),
+            "disc_sign_ref": ref["disc_sign"],
+            "late": False,
+            "near_linear": False,
+            "near_boundary": False,
+            "ill_conditioned": False,
+            "domain_clipped": ref["domain_clipped"],
+            "root_err_ulps": None,
+            "root_rel_err": None,
+            "resid_kernel_rel": None,
+            "resid_ideal_rel": None,
+            "baseline_resid_rel": None,
+        }
+        # conditioning flags
+        if ref["status"] == "ENTER" and ref["root"] is not None:
+            s = ref["root"]
+            scale = max(abs(ref["R0"]), Fraction(1))
+            rec["late"] = bool(s > frac_to_dec(scale, prec) * (10 ** 6))
+        rec["near_linear"] = abs(ref["a"]) <= Fraction(1, 10**10) * max(
+            ref["qq"], ref["rr"] * ref["rr"], Fraction(1)
+        )
+        r0sq = ref["R0"] * ref["R0"]
+        dd = ref["dd"]
+        rec["near_boundary"] = bool(
+            dd + r0sq != 0
+            and abs(ref["c"]) <= Fraction(1, 10**8) * (dd + r0sq)
+        )
+        rec["ill_conditioned"] = bool(
+            rec["near_linear"] or rec["late"] or rec["near_boundary"]
+            or r["category"].startswith("growing_out")
+            or r["category"] == "shrinking"
+        )
+        # root-level comparison
+        if status == 1 and ref["status"] == "ENTER" and ref["root"] is not None and cand_fr is not None:
+            err = abs(frac_to_dec(cand_fr, prec) - ref["root"])
+            u = ulp_dec(cand)
+            if u > 0:
+                rec["root_err_ulps"] = float(err / u)
+            if ref["root"] != 0:
+                rec["root_rel_err"] = float(err / abs(ref["root"]))
+            rec["resid_kernel_rel"] = _rel(
+                poly_resid(ak, bk, ck, cand_fr), ak, bk, ck, cand_fr
+            )
+            rec["resid_ideal_rel"] = _rel(
+                poly_resid(ref["a"], ref["b"], ref["c"], cand_fr),
+                ref["a"], ref["b"], ref["c"], cand_fr,
+            )
+            try:
+                nearest = float(ref["root"])
+                nfr = Fraction(nearest)
+                rec["baseline_resid_rel"] = _rel(
+                    poly_resid(ref["a"], ref["b"], ref["c"], nfr),
+                    ref["a"], ref["b"], ref["c"], nfr,
+                )
+            except (OverflowError, ValueError):
+                rec["baseline_resid_rel"] = None
+        out.append(rec)
+    return out
+
+
+def _sign(fr: Fraction) -> int:
+    return (fr > 0) - (fr < 0)
+
+
+def _rel(resid: Fraction, a, b, c, s) -> float:
+    scale = poly_scale(a, b, c, s)
+    if scale == 0:
+        return 0.0
+    return float(abs(resid) / scale)
+
+
+def route_reference(row, prec):
+    x = [parse_hex(row[f"canonx{i}"]) for i in range(3)]
+    w = [parse_hex(row[f"canonw{i}"]) for i in range(3)]
+    v = [parse_hex(row[f"v{i}"]) for i in range(3)]
+    center = [parse_hex(row[f"ct0{i}"]) for i in range(3)]
+    rr = parse_hex(row["rr"])
+    R0 = parse_hex(row["rt0"]) if row["rt0"] not in ("", "nan") else parse_hex(row["R0"])
+    a, b, cq, qq = coeffs(x, center, w, v, R0, rr)
+    status, root, info = classify(a, b, cq, prec, R0, rr)
+    return status, root, a, b, cq, qq, info["domain_clipped"]
+
+
+STATUS_NAME = {-1: "INVALID", 0: "MISS", 1: "ENTRY", 2: "UNCERTAIN"}
+KIND_NAME = {
+    0: "INSIDE",
+    1: "ENTRY",
+    2: "ESCAPED",
+    3: "TIME_RANGE_EXHAUSTED",
+    4: "INVALID",
+}
+
+
+def analyze_route_variant(rows, prec, dbl_eps=2.220446049250313e-16):
+    out = []
+    for r in rows:
+        (ref_status, ref_root, a, b, cq, qq,
+         ref_domain_clipped) = route_reference(r, prec)
+        status = int(r["status"])
+        kind = int(r["kind"])
+        F = _maybe_dec(r["F"])
+        tol = _maybe_dec(r["tol"])
+        kern_ok = r.get("kern_ok", "0") == "1"
+        kern_status = int(r["kern_status"]) if kern_ok else None
+        rec = {
+            "id": int(r["id"]),
+            "category": r["category"],
+            "ref_status": ref_status,
+            "ref_domain_clipped": ref_domain_clipped,
+            "status": status,
+            "kind": kind,
+            "kind_name": KIND_NAME.get(kind, "?"),
+            "failure_reason": int(r["failure_reason"]),
+            "fallback": int(r["fallback_evals"]),
+            "pi_match": int(r["pi_match"]),
+            "lcam_match": int(r["lcam_match"]),
+            "kern_status": kern_status if kern_ok else "",
+            "kern_sigma": r.get("kern_sigma", ""),
+            "F_over_tol": None,
+            "inside_mismatch": False,
+            "false_miss": False,
+            "false_candidate": False,
+            "unconfirmed": kind == 4 and r["reason_name"] == "ENTRY_UNCONFIRMED",
+            "kernel_false_miss": False,
+            "kernel_false_miss_escaped": False,
+        }
+        if ref_status == "INSIDE":
+            if kind != 0:
+                rec["inside_mismatch"] = True
+        elif ref_status == "ENTER":
+            if kind == 2:
+                rec["false_miss"] = True
+        elif ref_status == "MISS":
+            if kind == 1:
+                rec["false_candidate"] = True
+        if ref_status == "ENTER" and kern_ok and kern_status == 0:
+            rec["kernel_false_miss"] = True
+            if kind == 2:
+                rec["kernel_false_miss_escaped"] = True
+        if F is not None and tol is not None and tol > 0:
+            rec["F_over_tol"] = float(F / tol)
+        out.append(rec)
+    return out
+
+
+def _maybe_dec(s):
+    if s in ("", "nan"):
+        return None
+    low = s.lower()
+    if "inf" in low:
+        return None
+    return hex_to_dec(s, 300)
+
+
+def summarize_kernel(records_by_variant, refs):
+    """Aggregate counts, ULP distributions, and common-ENTRY intersections."""
+    variants = list(records_by_variant.keys())
+    by_id = {v: {} for v in variants}
+    for v in variants:
+        for rec in records_by_variant[v]:
+            by_id[v][rec["id"]] = rec
+
+    ids = sorted(refs.keys())
+    summary = {"variants": {}}
+    for v in variants:
+        recs = by_id[v]
+        counts = {
+            "enter": 0, "miss": 0, "uncertain": 0, "false_miss": 0,
+            "false_uncertain": 0, "false_candidate": 0,
+            "a_zero_ref": 0, "a_zero_kernel": 0, "a_sign_mismatch": 0,
+            "a_zero_collapse": 0, "a_spurious_nonzero": 0,
+            "domain_clipped_candidate": 0,
+            "ill_conditioned_enter": 0,
+            "near_boundary_enter": 0,
+            "domain_clipped_cases": 0,
+            "late_ref_enter": 0,
+        }
+        ulps = []
+        rels = []
+        for cid in ids:
+            rec = recs[cid]
+            st = rec["variant_status"]
+            if st == 1:
+                counts["enter"] += 1
+            elif st == 0:
+                counts["miss"] += 1
+            elif st == 2:
+                counts["uncertain"] += 1
+            if rec["ref_status"] == "ENTER" and st == 0:
+                counts["false_miss"] += 1
+            if rec["ref_status"] == "ENTER" and st == 2:
+                counts["false_uncertain"] += 1
+            if rec["ref_status"] == "MISS" and st == 1:
+                if rec["domain_clipped"]:
+                    counts["domain_clipped_candidate"] += 1
+                else:
+                    counts["false_candidate"] += 1
+            if rec["a_zero_ref"]:
+                counts["a_zero_ref"] += 1
+            if rec["a_zero_kernel"]:
+                counts["a_zero_kernel"] += 1
+            if not rec["a_zero_ref"] and rec["a_zero_kernel"]:
+                counts["a_zero_collapse"] += 1
+            if rec["a_zero_ref"] and not rec["a_zero_kernel"]:
+                counts["a_spurious_nonzero"] += 1
+            if rec["ill_conditioned"] and rec["ref_status"] == "ENTER":
+                counts["ill_conditioned_enter"] += 1
+            if rec["near_boundary"] and rec["ref_status"] == "ENTER":
+                counts["near_boundary_enter"] += 1
+            if rec["domain_clipped"]:
+                counts["domain_clipped_cases"] += 1
+            if (rec["a_sign_ref"] != 0 and rec["a_sign_kernel"] != 0
+                    and rec["a_sign_ref"] != rec["a_sign_kernel"]):
+                counts["a_sign_mismatch"] += 1
+            if rec["late"]:
+                counts["late_ref_enter"] += 1
+            if rec["root_err_ulps"] is not None:
+                ulps.append(rec["root_err_ulps"])
+            if rec["root_rel_err"] is not None:
+                rels.append(rec["root_rel_err"])
+        summary["variants"][v] = {
+            "counts": counts,
+            "root_ulp_all_ref_enter": _dist(ulps),
+            "root_rel_all_ref_enter": _dist(rels),
+        }
+
+    # Common reference-ENTER subset that every variant solved as ENTRY.
+    common = []
+    for cid in ids:
+        if refs[cid]["status"] != "ENTER":
+            continue
+        if all(by_id[v][cid]["variant_status"] == 1 for v in variants):
+            common.append(cid)
+    summary["common_enter_ids"] = len(common)
+    for v in variants:
+        vals = []
+        vals_nonlate = []
+        rels_nonlate = []
+        for cid in common:
+            rec = by_id[v][cid]
+            e = rec["root_err_ulps"]
+            if e is not None:
+                vals.append(e)
+                if not rec["ill_conditioned"]:
+                    vals_nonlate.append(e)
+                    if rec["root_rel_err"] is not None:
+                        rels_nonlate.append(rec["root_rel_err"])
+        summary["variants"][v]["root_ulp_common_enter"] = _dist(vals)
+        summary["variants"][v]["root_ulp_common_enter_wellcond"] = _dist(vals_nonlate)
+        summary["variants"][v]["root_rel_common_enter_wellcond"] = _dist(rels_nonlate)
+    return summary
+
+
+def _dist(vals):
+    if not vals:
+        return {"n": 0}
+    vs = sorted(vals)
+    n = len(vs)
+    def pct(p):
+        idx = min(n - 1, max(0, int(math.ceil(p * n)) - 1))
+        return vs[idx]
+    return {
+        "n": n,
+        "min": vs[0],
+        "median": pct(0.5),
+        "p95": pct(0.95),
+        "max": vs[-1],
+    }
+
+
+def summarize_route(records_by_variant):
+    summary = {"variants": {}}
+    for v, recs in records_by_variant.items():
+        counts = {
+            "inside_ok": 0, "inside_mismatch": 0, "false_miss": 0,
+            "false_candidate": 0, "unconfirmed": 0, "fallback_used": 0,
+            "accepted_entry": 0, "F_over_tol_gt1": 0, "pi_mismatch": 0,
+            "lcam_mismatch": 0, "kernel_false_miss": 0,
+            "kernel_false_miss_escaped": 0, "ref_domain_clipped": 0,
+        }
+        fmax = 0.0
+        by_cat = {}
+        for rec in recs:
+            cat = by_cat.setdefault(rec["category"], {"cases": 0, "fallback": 0,
+                                                      "entry": 0, "escaped": 0})
+            cat["cases"] += 1
+            if rec["kind"] == 0 and rec["ref_status"] == "INSIDE":
+                counts["inside_ok"] += 1
+            if rec["inside_mismatch"]:
+                counts["inside_mismatch"] += 1
+            if rec["false_miss"]:
+                counts["false_miss"] += 1
+            if rec["false_candidate"]:
+                counts["false_candidate"] += 1
+            if rec["unconfirmed"]:
+                counts["unconfirmed"] += 1
+            if rec["kernel_false_miss"]:
+                counts["kernel_false_miss"] += 1
+            if rec["kernel_false_miss_escaped"]:
+                counts["kernel_false_miss_escaped"] += 1
+            if rec["ref_domain_clipped"]:
+                counts["ref_domain_clipped"] += 1
+            if rec["fallback"] > 0:
+                counts["fallback_used"] += 1
+                cat["fallback"] += 1
+            if rec["kind"] == 1:
+                counts["accepted_entry"] += 1
+                cat["entry"] += 1
+                if rec["F_over_tol"] is not None:
+                    fmax = max(fmax, rec["F_over_tol"])
+                    if rec["F_over_tol"] > 1.0:
+                        counts["F_over_tol_gt1"] += 1
+            if rec["kind"] == 2:
+                cat["escaped"] += 1
+            if not rec["pi_match"]:
+                counts["pi_mismatch"] += 1
+            if not rec["lcam_match"]:
+                counts["lcam_mismatch"] += 1
+        summary["variants"][v] = {"counts": counts, "max_F_over_tol": fmax,
+                                  "by_category": by_cat}
+    return summary
+
+
+def precision_consistency(cases, ids, prec_a, prec_b):
+    """Compare reference classification and nearest-double root at two Decimal
+    precisions.  Coefficients/discriminant are exact rationals, so only the
+    square-root precision can differ.  Returns a list of mismatches."""
+    mismatches = []
+    for cid in sorted(ids):
+        ra = kernel_reference(cases[cid], prec_a)
+        rb = kernel_reference(cases[cid], prec_b)
+        reason = None
+        if ra["status"] != rb["status"]:
+            reason = "status"
+        elif ra["disc_sign"] != rb["disc_sign"]:
+            reason = "disc_sign"
+        elif ra["status"] == "ENTER":
+            try:
+                fa = float(ra["root"])
+                fb = float(rb["root"])
+            except (OverflowError, ValueError):
+                reason = "root_unrepresentable"
+            else:
+                if not (fa == fb or (math.isnan(fa) and math.isnan(fb))):
+                    reason = "nearest_double_root"
+        if reason is not None:
+            mismatches.append({"id": cid, "reason": reason,
+                               "status_a": ra["status"], "status_b": rb["status"]})
+    return mismatches

benchmarks/quadratic_precision/run.py

@@ -0,0 +1,458 @@
+#!/usr/bin/env python3
+"""End-to-end driver for the quadratic precision/performance benchmark.
+
+Self-contained: reads the current working-tree ``src/asymptotic.c``, generates
+and builds four variants, runs the accuracy probe, evaluates the exact-decimal
+reference, and times the kernel and the public pre-route.  All raw artifacts
+are written under ``--output-dir`` (default /tmp/opencode/quadratic-comparison).
+
+No production source, Makefile or git state is modified.
+"""
+from __future__ import annotations
+
+import argparse
+import csv
+import json
+import os
+import re
+import subprocess
+import sys
+
+# Do not leave __pycache__ inside the repository benchmark directory: the
+# imported local modules (build/cases/oracle) must not be cached here.
+sys.dont_write_bytecode = True
+
+from pathlib import Path
+
+HERE = Path(__file__).resolve().parent
+sys.path.insert(0, str(HERE))
+
+import build as buildmod  # noqa: E402
+import cases as casesmod  # noqa: E402
+import oracle  # noqa: E402
+
+VARIANTS = list(buildmod.VARIANTS)
+
+
+def run_cmd(cmd, log_path: Path, env=None):
+    with open(log_path, "a") as fh:
+        fh.write("$ " + " ".join(str(c) for c in cmd) + "\n")
+    proc = subprocess.run(cmd, capture_output=True, text=True, check=False, env=env)
+    with open(log_path, "a") as fh:
+        if proc.stdout:
+            fh.write(proc.stdout)
+        if proc.stderr:
+            fh.write(proc.stderr)
+        fh.write(f"exit={proc.returncode}\n")
+    return proc
+
+
+def read_csv(path: Path):
+    with open(path, newline="") as fh:
+        return list(csv.DictReader(fh))
+
+
+def run_timing(cmd, outp: Path, log_path: Path):
+    # Never ingest output left by a previous invocation, even after a failure.
+    outp.unlink(missing_ok=True)
+    proc = run_cmd(cmd, log_path)
+    if proc.returncode != 0 or not outp.is_file():
+        raise RuntimeError(f"timing failed or produced no fresh output: {log_path}")
+    return json.loads(outp.read_text())
+
+
+def write_metrics_csv(path: Path, records, columns):
+    with open(path, "w", newline="") as fh:
+        w = csv.DictWriter(fh, fieldnames=columns, extrasaction="ignore")
+        w.writeheader()
+        for r in records:
+            w.writerow(r)
+
+
+def write_curated(path: Path, rows, columns):
+    with open(path, "w", newline="") as fh:
+        w = csv.DictWriter(fh, fieldnames=columns, extrasaction="ignore")
+        w.writeheader()
+        for r in rows:
+            w.writerow(r)
+
+
+def build_curated_kernel(kernel_records):
+    ref_ids = [r["id"] for r in kernel_records[VARIANTS[0]] if r["id"] < 1000]
+    by = {v: {r["id"]: r for r in kernel_records[v]} for v in VARIANTS}
+    rows = []
+    for cid in sorted(ref_ids):
+        row = {"id": cid, "category": by[VARIANTS[0]][cid]["category"],
+               "ref": by[VARIANTS[0]][cid]["ref_status"]}
+        for v in VARIANTS:
+            row[f"{v}_status"] = by[v][cid]["variant_status"]
+            row[f"{v}_ulp"] = by[v][cid]["root_err_ulps"]
+            row[f"{v}_rel"] = by[v][cid]["root_rel_err"]
+        rows.append(row)
+    return rows
+
+
+def build_curated_route(route_records):
+    ref_ids = [r["id"] for r in route_records[VARIANTS[0]] if r["id"] < 100]
+    by = {v: {r["id"]: r for r in route_records[v]} for v in VARIANTS}
+    rows = []
+    for cid in sorted(ref_ids):
+        row = {"id": cid, "category": by[VARIANTS[0]][cid]["category"],
+               "ref": by[VARIANTS[0]][cid]["ref_status"]}
+        for v in VARIANTS:
+            row[f"{v}_kind"] = by[v][cid]["kind"]
+            row[f"{v}_fallback"] = by[v][cid]["fallback"]
+            row[f"{v}_F_over_tol"] = by[v][cid]["F_over_tol"]
+        rows.append(row)
+    return rows
+
+
+def main():
+    ap = argparse.ArgumentParser(description=__doc__)
+    ap.add_argument("--repo", default=str(HERE.parents[1]))
+    ap.add_argument("--output-dir", default="/tmp/opencode/quadratic-comparison")
+    ap.add_argument("--cc", default=os.environ.get("CC", "cc"))
+    ap.add_argument("--threads", type=int, default=4)
+    ap.add_argument("--precision", type=int, default=160)
+    ap.add_argument("--rounds", type=int, default=3)
+    ap.add_argument("--kernel-target-calls", type=int, default=10_000_000)
+    ap.add_argument("--route-target-calls", type=int, default=1_000_000)
+    ap.add_argument("--skip-build", action="store_true")
+    args = ap.parse_args()
+
+    repo = Path(args.repo).resolve()
+    out = Path(args.output_dir).resolve()
+    (out / "raw").mkdir(parents=True, exist_ok=True)
+    (out / "cases").mkdir(parents=True, exist_ok=True)
+    (out / "logs").mkdir(parents=True, exist_ok=True)
+    import time as _time
+    t_start = _time.time()
+    phase_times = {}
+
+    # ---------------------------------------------------------------- cases
+    t0 = _time.time()
+    kernel_cases, route_cases = casesmod.build()
+    casesmod.write_header(kernel_cases, route_cases, out / "cases" / "quadratic_cases.h")
+    casesmod.write_json(kernel_cases, route_cases, out / "cases" / "cases.json")
+    phase_times["case_gen"] = _time.time() - t0
+    print(f"cases: kernel={len(kernel_cases)} route={len(route_cases)}")
+
+    # ------------------------------------------------------------ environment
+    env_info = buildmod.environment(repo, args.cc, args.threads)
+    env_info["python_version"] = sys.version.split()[0]
+    env_info["invocation"] = " ".join([sys.executable, *sys.argv])
+    (out / "environment.json").write_text(json.dumps(env_info, indent=2) + "\n")
+    (out / "raw" / "invocation.txt").write_text(
+        " ".join([sys.executable, *sys.argv]) + "\n"
+    )
+    print(f"source sha256: {env_info['source_sha256'][:16]}  cc: {env_info['cc_version']}")
+
+    # ---------------------------------------------------------------- build
+    t0 = _time.time()
+    if not args.skip_build:
+        info = buildmod.build_all(repo, out, args.cc, buildmod.BASE_FLAGS, args.threads)
+        exes = {k: Path(v) for k, v in info["executables"].items()}
+    else:
+        manifest, reasons = buildmod.verify_manifest(repo, out, buildmod.BASE_FLAGS, args.cc)
+        if reasons:
+            print("refusing --skip-build: existing build does not match the "
+                  f"current source/flags/fixtures ({', '.join(reasons)}); "
+                  "rerun without --skip-build", file=sys.stderr)
+            sys.exit(1)
+        exes = {k: Path(v) for k, v in manifest["executables"].items()}
+    phase_times["build"] = _time.time() - t0
+    for name, exe in exes.items():
+        if not exe.exists():
+            print(f"missing executable for {name}: {exe}", file=sys.stderr)
+            sys.exit(1)
+
+    # -------------------------------------------------------------- accuracy
+    t0 = _time.time()
+    env_log = out / "logs" / "probe_environment.log"
+    env_log.write_text("")
+    for name in VARIANTS:
+        kout = out / "raw" / f"kernel_{name}.csv"
+        rout = out / "raw" / f"route_{name}.csv"
+        proc = run_cmd(
+            [str(exes[name]), "accuracy", "--kernel-out", str(kout),
+             "--route-out", str(rout)],
+            out / "logs" / f"accuracy_{name}.log",
+        )
+        with open(env_log, "a") as fh:
+            fh.write(proc.stdout)
+        if proc.returncode != 0:
+            print(f"accuracy probe failed for {name}", file=sys.stderr)
+            sys.exit(1)
+    phase_times["accuracy"] = _time.time() - t0
+
+    # ---------------------------------------------------------------- oracle
+    t0 = _time.time()
+    kernel_ref_cases, route_cases_json = oracle.load_cases(out / "cases" / "cases.json")
+    refs = {
+        cid: oracle.kernel_reference(case, args.precision)
+        for cid, case in kernel_ref_cases.items()
+    }
+    kernel_records = {}
+    route_records = {}
+    for name in VARIANTS:
+        krows = read_csv(out / "raw" / f"kernel_{name}.csv")
+        kernel_records[name] = oracle.analyze_kernel_variant(krows, refs, args.precision)
+        rrows = read_csv(out / "raw" / f"route_{name}.csv")
+        route_records[name] = oracle.analyze_route_variant(rrows, args.precision)
+
+    kcols = ["id", "category", "ref_status", "variant_status", "root_err_ulps",
+             "root_rel_err",
+             "resid_kernel_rel", "resid_ideal_rel", "baseline_resid_rel",
+             "a_zero_ref", "a_zero_kernel", "a_sign_ref", "a_sign_kernel",
+             "disc_sign_ref", "near_linear", "near_boundary",
+             "ill_conditioned", "domain_clipped", "late"]
+    for name in VARIANTS:
+        write_metrics_csv(out / "raw" / f"kernel_metrics_{name}.csv",
+                          kernel_records[name], kcols)
+    rcols = ["id", "category", "ref_status", "ref_domain_clipped", "status",
+             "kind", "kind_name", "failure_reason", "fallback", "pi_match",
+             "lcam_match", "kern_status", "kern_sigma", "F_over_tol",
+             "inside_mismatch", "false_miss", "false_candidate", "unconfirmed",
+             "kernel_false_miss", "kernel_false_miss_escaped"]
+    for name in VARIANTS:
+        write_metrics_csv(out / "raw" / f"route_metrics_{name}.csv",
+                          route_records[name], rcols)
+
+    ksummary = oracle.summarize_kernel(kernel_records, refs)
+    rsummary = oracle.summarize_route(route_records)
+    consistency_ids = {
+        rec["id"] for rec in kernel_records[VARIANTS[0]]
+        if rec["id"] < 1000 or rec["near_linear"] or rec["late"]
+        or rec["near_boundary"] or rec["domain_clipped"]
+    }
+    precision_mismatches = oracle.precision_consistency(
+        kernel_ref_cases, consistency_ids, args.precision, 240
+    )
+    curated_kernel = build_curated_kernel(kernel_records)
+    curated_route = build_curated_route(route_records)
+    write_curated(out / "raw" / "curated_kernel.csv", curated_kernel,
+                  ["id", "category", "ref"] + [f"{v}_{f}" for v in VARIANTS
+                                               for f in ("status", "ulp", "rel")])
+    write_curated(out / "raw" / "curated_route.csv", curated_route,
+                  ["id", "category", "ref"] + [f"{v}_{f}" for v in VARIANTS
+                                               for f in ("kind", "fallback",
+                                                         "F_over_tol")])
+
+    # ------------------------------------------------------------- hardware FMA
+    phase_times["oracle"] = _time.time() - t0
+    fma_info = verify_fma(out, args.threads)
+
+    # ---------------------------------------------------------------- timing
+    t0 = _time.time()
+    timings = timed_runs(exes, out, args, kernel_target_calls=args.kernel_target_calls,
+                         route_target_calls=args.route_target_calls)
+    lin_target = max(100_000, args.kernel_target_calls // 20)
+    timings["linearity"] = linearity_check(exes, out, lin_target)
+    phase_times["timing"] = _time.time() - t0
+
+    # ---------------------------------------------------------------- summary
+    summary = {
+        "source_sha256": env_info["source_sha256"],
+        "flags": buildmod.BASE_FLAGS,
+        "precision": args.precision,
+        "precision_consistency": {
+            "ids_checked": len(consistency_ids),
+            "mismatches": precision_mismatches,
+        },
+        "threads_timing": args.threads,
+        "cases": {"kernel": len(kernel_cases), "route": len(route_cases)},
+        "kernel_accuracy": ksummary,
+        "route_accuracy": rsummary,
+        "curated_kernel": curated_kernel,
+        "fma": fma_info,
+        "timing": timings,
+        "phase_seconds": phase_times,
+        "total_seconds": _time.time() - t_start,
+    }
+    (out / "summary.json").write_text(json.dumps(summary, indent=2) + "\n")
+    print_summary(summary)
+    print(f"\nartifacts under {out}")
+
+
+def verify_fma(out: Path, threads: int):
+    info = {"objects": {}}
+    for name in VARIANTS:
+        obj = out / "build" / f"probe_{name}.o"
+        if not obj.exists():
+            continue
+        proc = subprocess.run(["objdump", "-d", str(obj)], capture_output=True,
+                              text=True, check=False)
+        n = sum(1 for line in proc.stdout.splitlines()
+                if "vfmadd" in line or "vfmsub" in line or "fmadd" in line)
+        nm = subprocess.run(["nm", "-u", str(obj)], capture_output=True,
+                            text=True, check=False)
+        undef = sorted({tok for line in nm.stdout.splitlines()
+                        for tok in line.split()
+                        if tok.startswith("fma")})
+        kd = subprocess.run(["objdump", "-dr", str(obj)], capture_output=True,
+                            text=True, check=False)
+        kb_lines = []
+        inside = False
+        for line in kd.stdout.splitlines():
+            if re.match(r"^[0-9a-f]+ <entry_solve", line):
+                inside = True
+                kb_lines.append(line)
+                continue
+            if inside:
+                if re.match(r"^[0-9a-f]+ <", line):
+                    break
+                kb_lines.append(line)
+        kb = "\n".join(kb_lines)
+        x87 = sum(1 for line in kb.splitlines()
+                  if re.search(r"\b(fld|fstp|fmul|fadd|fsub|fdiv|fcom)\b", line))
+        vfma = sum(1 for line in kb.splitlines()
+                   if "vfmadd" in line or "vfmsub" in line)
+        calls_fmal = sum(1 for line in kb.splitlines() if "fmal" in line)
+        info["objects"][name] = {
+            "fma_instructions": n,
+            "lowered_hardware": n > 0,
+            "undefined_fma_symbols": undef,
+            "entry_solve_x87": x87,
+            "entry_solve_vfmadd": vfma,
+            "entry_solve_fmal_calls": calls_fmal,
+        }
+        print(f"arith verify {name}: fused={n} undef={undef} "
+              f"entry_solve[x87={x87} vfmadd={vfma} fmal_calls={calls_fmal}]")
+    (out / "raw" / "fma_verify.json").write_text(json.dumps(info, indent=2) + "\n")
+    return info
+
+
+def linearity_check(exes, out: Path, target_calls: int):
+    """Confirm the microbench time scales with the call count (no hoisting)."""
+    res = {}
+    for name in ("ld_fma", "double_fma"):
+        secs = []
+        for mult in (1, 2):
+            outp = out / "raw" / f"linearity_{name}_{mult}.json"
+            d = run_timing([str(exes[name]), "microbench", "--out", str(outp),
+                     "--target-calls", str(target_calls * mult)],
+                    outp, out / "logs" / f"linearity_{name}_{mult}.log")
+            secs.append(d["seconds"])
+        res[name] = {"t1": secs[0], "t2": secs[1],
+                     "ratio": secs[1] / secs[0] if secs[0] > 0 else None}
+    return res
+
+
+def timed_runs(exes, out: Path, args, kernel_target_calls, route_target_calls):
+    schedules = []
+    fwd = list(VARIANTS)
+    schedules.append(fwd)
+    schedules.append(list(reversed(fwd)))
+    for i in range(2, args.rounds):
+        schedules.append(fwd if i % 2 == 0 else list(reversed(fwd)))
+    schedules = schedules[: args.rounds]
+
+    results = {"microbench": [], "routebench": []}
+
+    def order(seq):
+        return [exes[n] for n in seq]
+
+    for rnd, seq in enumerate(schedules):
+        for name in seq:
+            outp = out / "raw" / f"microbench_{name}_r{rnd}.json"
+            d = run_timing([str(exes[name]), "microbench", "--out", str(outp),
+                     "--target-calls", str(kernel_target_calls)],
+                    outp, out / "logs" / f"microbench_{name}_r{rnd}.log")
+            d["round"] = rnd
+            results["microbench"].append(d)
+    for rnd, seq in enumerate(schedules):
+        for name in seq:
+            outp = out / "raw" / f"routebench_{name}_r{rnd}.json"
+            d = run_timing([str(exes[name]), "routebench", "--out", str(outp),
+                     "--target-calls", str(route_target_calls),
+                     "--threads", str(args.threads), "--max-seconds", "20"],
+                    outp, out / "logs" / f"routebench_{name}_r{rnd}.log")
+            d["round"] = rnd
+            results["routebench"].append(d)
+    return results
+
+
+def print_summary(summary):
+    print("\n=== kernel accuracy (per variant) ===")
+    print(f"{'variant':<18}{'ENTER':>7}{'MISS':>7}{'UNC':>6}{'falseMiss':>10}"
+          f"{'falseCand':>10}{'domClipCand':>12}{'a0ref':>7}{'a0k':>6}{'asign':>6}"
+          f"{'illcond':>8}{'nearBnd':>8}{'medULP':>10}{'p95ULP':>11}{'maxULP':>11}")
+    for v, s in summary["kernel_accuracy"]["variants"].items():
+        c = s["counts"]
+        d = s["root_ulp_all_ref_enter"]
+        print(f"{v:<18}{c['enter']:>7}{c['miss']:>7}{c['uncertain']:>6}"
+              f"{c['false_miss']:>10}{c['false_candidate']:>10}"
+              f"{c['domain_clipped_candidate']:>12}"
+              f"{c['a_zero_ref']:>7}{c['a_zero_kernel']:>6}"
+              f"{c['a_sign_mismatch']:>6}{c['ill_conditioned_enter']:>8}"
+              f"{c['near_boundary_enter']:>8}"
+              f"{d.get('median', float('nan')):>10.4g}{d.get('p95', float('nan')):>11.4g}"
+              f"{d.get('max', float('nan')):>11.4g}")
+    print("well-conditioned common-ENTER ULP (excludes late / near-linear /"
+          " growing-out / shrinking):")
+    for v, s in summary["kernel_accuracy"]["variants"].items():
+        d = s["root_ulp_common_enter_wellcond"]
+        r = s["root_rel_common_enter_wellcond"]
+        print(f"  {v:<18}n={d.get('n',0):>5} median={d.get('median', float('nan')):>9.3g}"
+              f" p95={d.get('p95', float('nan')):>9.3g}"
+              f" max={d.get('max', float('nan')):>9.3g}"
+              f"   |rel err median={r.get('median', float('nan')):>9.3g}"
+              f" p95={r.get('p95', float('nan')):>9.3g}")
+    print(f"common reference-ENTER ids solved ENTRY by all variants: "
+          f"{summary['kernel_accuracy']['common_enter_ids']}")
+    print("\n=== curated kernel cases (status/ULP) ===")
+    print(f"{'id':>4} {'category':<20} {'ref':<6} "
+          + " ".join(f"{v:>16}" for v in VARIANTS))
+    for row in summary.get("curated_kernel", []):
+        cells = []
+        for v in VARIANTS:
+            st = row.get(f"{v}_status")
+            u = row.get(f"{v}_ulp")
+            cells.append(f"{st}/{float(u):.3g}" if u not in (None, "") else f"{st}/-")
+        print(f"{row['id']:>4} {row['category']:<20} {row.get('ref',''):<6} "
+              + " ".join(f"{c:>16}" for c in cells))
+    print("\n=== route accuracy (per variant) ===")
+    print(f"{'variant':<18}{'insideOK':>9}{'insideBad':>10}{'falseMiss':>10}"
+          f"{'falseCand':>10}{'unconf':>8}{'kernFM':>8}{'kernFMesc':>10}"
+          f"{'fallback':>9}{'F>tol':>7}{'maxF/tol':>10}")
+    for v, s in summary["route_accuracy"]["variants"].items():
+        c = s["counts"]
+        print(f"{v:<18}{c['inside_ok']:>9}{c['inside_mismatch']:>10}"
+              f"{c['false_miss']:>10}{c['false_candidate']:>10}"
+              f"{c['unconfirmed']:>8}{c['kernel_false_miss']:>8}"
+              f"{c['kernel_false_miss_escaped']:>10}"
+              f"{c['fallback_used']:>9}"
+              f"{c['F_over_tol_gt1']:>7}{s['max_F_over_tol']:>10.3g}")
+    pc = summary.get("precision_consistency", {})
+    print(f"reference precision {summary['precision']} vs 240: checked "
+          f"{pc.get('ids_checked')} curated/near-linear/late ids, "
+          f"mismatches={len(pc.get('mismatches', []))}")
+    print("\nroute fallback / entry / escaped by category "
+          "(ld_plain | ld_fma | double_fma):")
+    ra = summary["route_accuracy"]["variants"]
+    cats = sorted({c for s in ra.values() for c in s.get("by_category", {})})
+    for cat in cats:
+        cells = []
+        for v in ("ld_plain", "ld_fma", "double_fma"):
+            d = ra[v].get("by_category", {}).get(cat, {})
+            cells.append(f"fb={d.get('fallback',0)} e={d.get('entry',0)} "
+                         f"x={d.get('escaped',0)}")
+        print(f"  {cat:<20} " + " | ".join(cells))
+    print("\n=== timing (median ns/call over rounds) ===")
+    for mode in ("microbench", "routebench"):
+        per = {}
+        for rec in summary["timing"][mode]:
+            per.setdefault(rec["variant"], []).append(rec["ns_per_call"])
+        for v, vals in per.items():
+            vals.sort()
+            med = vals[len(vals) // 2]
+            print(f"{mode:<12}{v:<18}median={med:>10.2f} ns/call  "
+                  f"rounds={len(vals)}")
+    lin = summary["timing"].get("linearity", {})
+    for v, d in lin.items():
+        print(f"linearity  {v:<18}t1={d['t1']:.4f}s t2={d['t2']:.4f}s "
+              f"ratio={d['ratio']:.3f} (expect ~2)")
+
+
+if __name__ == "__main__":
+    main()

nr_spacetime_movie_renderer_design.md

@@ -1007,6 +1007,35 @@ residual、Chebyshev 表或解析主项；运行期不得建表。
 - Minkowski entry quadratic 用稳定根公式（$q=-\tfrac12(b+\mathrm{copysign}
   (\sqrt\Delta,b))$，取最小正根）；$c=0$（相机在边界）时按 $b=\mathrm dF/
   \mathrm ds$ 分类。worldtube sample 必须有限且 $R>0$，否则 `INVALID`。
+  Alcubierre 等平直外区共用此解法；相对位置/速度、系数、判别式和根采用
+  `long double` 中间量，最后才转换为 backend 的 double 状态。某些平台的
+  `long double` 不提供额外有效位，因此精度提升不能替代入口验证。
+  送入根公式前，以系数最大绝对值的二进制指数共同缩放 $a,b,c$（最大值进入
+  $[1/2,1)$），避免判别式乘积溢出；不改变传播参数及根，也不覆盖系数构造前
+  已发生的误差。非有限系数或缩放后非零系数落入 subnormal 范围，转入不确定
+  路径而非当作退化线性式。判别式 $b^2-4ac$ 和入口斜率采用普通 `long double`
+  运算，不使用显式 FMA；精度/成本对照见 `benchmarks/quadratic_precision/`。
+  本机软件 `fmal` 成本明显，测试未显示足以抵偿成本的精度收益；不确定带、
+  几何验证及后备定位仍保留，不能把此样本结论当作全参数可靠性保证。
+- **入口快路径与后备定位分离**：解析相交只提供候选入口；采用前须在实际 backend
+  坐标、实际入口时间重新检查 worldtube 残差。通过原有几何舍入容差的候选沿快路径
+  激活；未通过的候选沿同一外区 geodesic，从原相机事件重新求值并用确定在外/
+  严格在内的区间定位首次入口。公共 localizer 不依赖具体 metric backend，
+  各已支持的外区提供轨迹求值与首次入口 bracket；不能把离散采样无交点当作 miss。
+  miss 判定必须排除真实入口；擦边或求根病态造成的数值不确定性不能作为
+  `ESCAPED` 的依据，允许保守地生成候选入口并进行后备验证。
+  浮点判别式等于零不构成精确擦边证明；重建最小点的单次正残差也不构成
+  miss 证明（大时间/坐标的舍入会移动该最小点）。精确擦边可由独立几何证据
+  排除入口，例如固定半径段内某个不变坐标的精确距离已不小于球半径；
+  无法获得这种证据且无法构造可信 bracket 时仍返回 `ENTRY_UNCONFIRMED`。
+  后备路径返回 bracket 收敛后的可表示内侧轨迹状态，不做径向位置投影，也不放宽
+  原 worldtube 检查容差。未能确认首次入口须报告 `INCOMPLETE/ENTRY_UNCONFIRMED`，
+  callback、历史及几何错误仍传播各自具体 reason，不能改写为 capture 或 escape。
+  定位在 pre-route 中完成，不在 slab sweep 中倒回相机时间；不重置相机 `L0`，
+  内区 accepted-step/lookback 预算仍从最终激活事件起算。成功的后备入口由
+  `entry_fallback_evaluations` 记录 localizer 的轨迹求值次数（不含构造 bracket
+  的探测），不混入内区 RHS 计数；临时状态由
+  调用者独占，不引入共享可变缓存。此标记不改变 lens-map 二进制布局。
 - **构造期验证优先**：`spacetime_create_*()` 成功即承诺该 source 已可安全光追。
   每个 constructor 在安装 ops/context 后调用公共 `spacetime_source_finalize()`：
   检查 ops/context 完整、`end_id` 唯一且非 `NONE`、exterior kind 受支持、

src/asymptotic.c

@@ -1,5 +1,6 @@
 #include "asymptotic.h"
 
+#include "asymptotic_entry.h"
 #include "asymptotic_schwarzschild.h"
 
 #include <float.h>
@@ -188,73 +189,190 @@ int asymptotic_backend_from_canonical(const SpacetimeSource *source,
   return 0;
 }
 
-/* Solve |d + q s|^2 = (R0 - rr s)^2 for the smallest s >= 0 with outside ->
+/* Long-double coefficients of the relative-distance quadratic
- * inside crossing.  Returns 1 on entry (sets s), 0 on miss, -1 on error. */
+ * F(sigma) = |d + q sigma|^2 - (R0 - rr sigma)^2 = a sigma^2 + b sigma + c,
-static int solve_entry_quadratic(const double d[3], const double q[3],
+ * where d = x_cur - c_frame and q = w_frame + v_frame.  Keeping them in long
-                                 double R0, double rr, double *s_out) {
+ * double preserves the cancellation-prone grazing entries; the caller hands
-  const double qq = dot3(q, q);
+ * the coefficients to both the root solve and the conservative fallback. */
-  const double a = qq - rr * rr;
+typedef struct {
-  const double b = 2.0 * (dot3(d, q) + R0 * rr);
+  long double a, b, c;
-  const double c = dot3(d, d) - R0 * R0;
+} EntryQuadratic;
-  if (c < 0.0) {
+
+static EntryQuadratic entry_quadratic_coeffs(
+    const double x_cur[3], const double c_frame[3], const double w_frame[3],
+    const double v_frame[3], double R0, double rr) {
+  long double d[3], q[3];
+  for (int i = 0; i < 3; ++i) {
+    d[i] = (long double)x_cur[i] - (long double)c_frame[i];
+    q[i] = (long double)w_frame[i] + (long double)v_frame[i];
+  }
+  long double qq = 0.0L, dot_dq = 0.0L, dot_dd = 0.0L;
+  for (int i = 0; i < 3; ++i) {
+    qq += q[i] * q[i];
+    dot_dq += d[i] * q[i];
+    dot_dd += d[i] * d[i];
+  }
+  const long double R0_ld = (long double)R0;
+  const long double rr_ld = (long double)rr;
+  EntryQuadratic k;
+  k.a = qq - rr_ld * rr_ld;
+  k.b = 2.0L * (dot_dq + R0_ld * rr_ld);
+  k.c = dot_dd - R0_ld * R0_ld;
+  return k;
+}
+
+typedef enum {
+  ENTRY_SOLVE_MISS = 0,     /* proven no outside->inside crossing */
+  ENTRY_SOLVE_ENTRY = 1,    /* a first-crossing candidate parameter was found */
+  ENTRY_SOLVE_UNCERTAIN = 2 /* discriminant/degeneracy at the resolution floor */
+} EntrySolveResult;
+
+/* For normalized coefficients, compute b*b - 4*a*c without overflowing.
+ * Near-cancellation is handled conservatively by the uncertainty band. */
+static long double entry_discriminant(const EntryQuadratic *k,
+                                      long double *scale) {
+  const long double four_a = 4.0L * k->a;
+  const long double ac = four_a * k->c;
+  *scale = k->b * k->b + fabsl(ac);
+  return k->b * k->b - ac;
+}
+
+/* Solve the quadratic for the smallest sigma >= 0 with an outside -> inside
+ * crossing.  The stable roots are evaluated in long double.  A discriminant
+ * that is negative but within its own rounding bound, an exact double root, or
+ * a root that does not move inward is not a proof: it is reported as
+ * ENTRY_SOLVE_UNCERTAIN so the caller attempts a strict-inside bracket instead
+ * of fabricating a miss.  Exact algebra (a == 0) stays the only linear branch;
+ * a tiny-but-nonzero `a` always goes through the discriminant, so a future
+ * entry at a huge parameter is never discarded by an approximate threshold. */
+static EntrySolveResult entry_solve(const EntryQuadratic *k, double *s_out) {
+  if (!isfinite(k->a) || !isfinite(k->b) || !isfinite(k->c))
+    return ENTRY_SOLVE_UNCERTAIN;
+  /* A common power-of-two scale preserves roots and coefficient signs without
+   * adding division rounding.  Normalize before squaring/products, retaining
+   * the original coefficients in the caller for geometric fallback.  Do not
+   * silently turn an underflowed coefficient into a linear/boundary case. */
+  EntryQuadratic scaled = *k;
+  const long double magnitude = fmaxl(fabsl(k->a),
+                                      fmaxl(fabsl(k->b), fabsl(k->c)));
+  if (magnitude > 0.0L) {
+    int exponent;
+    (void)frexpl(magnitude, &exponent);
+    scaled.a = scalbnl(k->a, -exponent);
+    scaled.b = scalbnl(k->b, -exponent);
+    scaled.c = scalbnl(k->c, -exponent);
+    if ((k->a != 0.0L && fabsl(scaled.a) < LDBL_MIN) ||
+        (k->b != 0.0L && fabsl(scaled.b) < LDBL_MIN) ||
+        (k->c != 0.0L && fabsl(scaled.c) < LDBL_MIN))
+      return ENTRY_SOLVE_UNCERTAIN;
+  }
+  k = &scaled;
+  if (k->c < 0.0L) {
     /* Strictly inside; the lifecycle normally handles this as INSIDE. */
     *s_out = 0.0;
-    return 1;
+    return ENTRY_SOLVE_ENTRY;
   }
-  if (c == 0.0) {
+  if (k->c == 0.0L) {
     /* On the boundary: classify by dF/ds = b.  Past-inward enters at once;
      * outward/tangent rays may still re-enter later when the sphere shrinks
      * (a < 0), so do not declare a permanent miss on the zero root. */
-    if (b < 0.0) {
+    if (k->b < 0.0L) {
       *s_out = 0.0;
-      return 1;
+      return ENTRY_SOLVE_ENTRY;
     }
-    if (b == 0.0) {
+    if (k->b == 0.0L) {
-      if (a < 0.0) {
+      if (k->a < 0.0L) {
         *s_out = 0.0;
-        return 1;
+        return ENTRY_SOLVE_ENTRY;
       }
-      return 0;
+      return ENTRY_SOLVE_MISS;
     }
-    if (a < 0.0) {
+    if (k->a < 0.0L) {
-      *s_out = -b / a;
+      *s_out = (double)(-k->b / k->a);
-      return 1;
+      return ENTRY_SOLVE_ENTRY;
     }
-    return 0;
+    return ENTRY_SOLVE_MISS;
   }
-  /* Compare the quadratic coefficient against the velocity-squared scale it
+  /* c > 0: the ray starts outside. */
-   * is built from; mixing in R0^2 would let a large radius misclassify a
+  if (k->a == 0.0L) {
-   * genuinely quadratic entry as linear. */
+    /* Exact linear branch only.  b >= 0 never crosses for sigma > 0. */
-  const double scale = qq + rr * rr;
+    if (!(k->b < 0.0L))
-  if (fabs(a) <= 32.0 * DBL_EPSILON * scale) {
+      return ENTRY_SOLVE_MISS;
-    if (!isfinite(b) || b >= 0.0)
+    const long double s = -k->c / k->b;
-      return 0;
+    if (!(s > 0.0L))
-    const double s = -c / b;
+      return ENTRY_SOLVE_MISS;
-    if (s <= 0.0)
+    *s_out = (double)s;
-      return 0;
+    return ENTRY_SOLVE_ENTRY;
-    *s_out = s;
-    return 1;
   }
-  const double disc = b * b - 4.0 * a * c;
+  /* Ordinary long-double products; near-zero differences need fallback. */
-  if (!isfinite(disc) || disc <= 0.0)
+  long double disc_scale;
-    return 0;
+  const long double disc = entry_discriminant(k, &disc_scale);
-  /* Numerically stable quadratic roots: q avoids cancellation in the root
+  /* Conservative discriminant resolution: the long-double evaluation error
+   * plus the rounding the double inputs already carry through the frame
+   * rotation/translation into d and q.  The double term dominates and keeps a
+   * near-tangent discriminant from being read as a proven miss. */
+  const long double disc_err =
+      64.0L * ((long double)DBL_EPSILON + (long double)LDBL_EPSILON) *
+      disc_scale;
+  if (!isfinite(disc))
+    return ENTRY_SOLVE_UNCERTAIN;
+  if (disc < -disc_err)
+    return ENTRY_SOLVE_MISS;
+  if (fabsl(disc) <= disc_err)
+    return ENTRY_SOLVE_UNCERTAIN;
+  /* Numerically stable quadratic roots: qq2 avoids cancellation in the root
    * with the same sign as b, which is exactly the small entry root when the
    * camera sits just outside a large sphere. */
-  const double root = sqrt(disc);
+  const long double root = sqrtl(disc);
-  const double qq2 = -0.5 * (b + copysign(root, b));
+  const long double qq2 = -0.5L * (k->b + copysignl(root, k->b));
-  const double r1 = qq2 / a;
+  if (qq2 == 0.0L)
-  const double r2 = c / qq2;
+    return ENTRY_SOLVE_UNCERTAIN;
+  const long double r1 = qq2 / k->a;
+  const long double r2 = k->c / qq2;
   /* The first outside->inside crossing is the smallest positive root. */
-  double s = INFINITY;
+  long double s = INFINITY;
-  if (r1 > 0.0)
+  if (r1 > 0.0L)
     s = r1;
-  if (r2 > 0.0 && r2 < s)
+  if (r2 > 0.0L && r2 < s)
     s = r2;
   if (!(s < INFINITY))
-    return 0;
+    return (r1 > 0.0L || r2 > 0.0L) ? ENTRY_SOLVE_UNCERTAIN : ENTRY_SOLVE_MISS;
-  *s_out = s;
+  /* First crossing must move inward (dF/dsigma < 0).  A nonnegative slope
-  return 1;
+   * means the stable-root selection picked the exit root or the roots merged;
+   * that is uncertain, not a proof of a miss. */
+  const long double slope = 2.0L * k->a * s + k->b;
+  if (!(slope < 0.0L))
+    return ENTRY_SOLVE_UNCERTAIN;
+  *s_out = (double)s;
+  return ENTRY_SOLVE_ENTRY;
+}
+
+static void minkowski_route_entry(const SpacetimeAsymptoticEnd *end, double t0,
+                                  const double x_frame[3],
+                                  const double w_frame[3], double s_entry,
+                                  SpacetimeEndId end_id,
+                                  AsymptoticRoute *route);
+
+/* Opaque evaluator context: repropagate constant-velocity motion from the
+ * ORIGINAL camera state to a total past parameter, never from a nearby root. */
+typedef struct {
+  const SpacetimeAsymptoticEnd *end;
+  double t0;
+  const double *x_frame;
+  const double *w_frame;
+  double log_alpha_p0;
+  SpacetimeEndId end_id;
+} MinkowskiEntryContext;
+
+static AsymptoticStatus minkowski_entry_evaluate(void *opaque, double parameter,
+                                                 AsymptoticRoute *state) {
+  const MinkowskiEntryContext *ctx = opaque;
+  *state = (AsymptoticRoute){0};
+  minkowski_route_entry(ctx->end, ctx->t0, ctx->x_frame, ctx->w_frame,
+                        parameter, ctx->end_id, state);
+  state->log_alpha_p0 = ctx->log_alpha_p0;
+  state->log_alpha_p0_camera = ctx->log_alpha_p0;
+  return ASYMPTOTIC_OK;
 }
 
 static void minkowski_route_escaped(const SpacetimeAsymptoticEnd *end,
@@ -287,13 +405,248 @@ static void minkowski_route_entry(const SpacetimeAsymptoticEnd *end, double t0,
     route->Pi[i] = -w_backend[i];
 }
 
+/* Recover a first-entry bracket inside the CURRENT constant-motion segment
+ * when the closed-form candidate fails geometric validation or the
+ * discriminant is uncertain.  The outside endpoint is the segment start (the
+ * previous segments produced no entry), and the inside endpoint is either the
+ * convex minimum or a modest, geometrically grown step past the candidate
+ * root.  Both are confirmed with the actual worldtube callback.  On success
+ * `route` carries the localized ENTRY and a nonzero fallback evaluation count.
+ * Any failure leaves `route->failure_reason` set and never reports an escape;
+ * one reconstructed probe cannot certify a miss. */
+static AsymptoticStatus minkowski_fallback_entry(
+    const SpacetimeSource *source, const SpacetimeAsymptoticEnd *end, double t0,
+    const double x_frame[3], const double w_frame[3], double log_alpha_p0,
+    SpacetimeEndId end_id, double s_base, double s_segment,
+    const EntryQuadratic *k, EntrySolveResult solve, double candidate_sigma,
+    AsymptoticRoute *route) {
+  MinkowskiEntryContext ctx = {.end = end,
+                               .t0 = t0,
+                               .x_frame = x_frame,
+                               .w_frame = w_frame,
+                               .log_alpha_p0 = log_alpha_p0,
+                               .end_id = end_id};
+
+  /* Outside endpoint: the segment start must still be outside. */
+  AsymptoticRoute outside_state;
+  AsymptoticStatus status =
+      minkowski_entry_evaluate(&ctx, s_base, &outside_state);
+  if (status != ASYMPTOTIC_OK) {
+    route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED;
+    route->end_id = end_id;
+    return status;
+  }
+  RayReason why = RAY_REASON_NONE;
+  double f_out = 0.0, tol_out = 0.0;
+  status = asymptotic_entry_geometry(source, end_id, outside_state.activate_t,
+                                     outside_state.x, &f_out, &tol_out, &why);
+  if (status != ASYMPTOTIC_OK) {
+    route->failure_reason = why;
+    route->end_id = end_id;
+    return status;
+  }
+  if (!(f_out >= 0.0)) {
+    /* Already inside at the segment start: an earlier segment missed the
+     * crossing.  Do not fabricate a bracket from it. */
+    route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED;
+    route->end_id = end_id;
+    return ASYMPTOTIC_INVALID;
+  }
+
+  int have_inside = 0;
+  double inside_sigma = 0.0;
+
+  if (k->a > 0.0L) {
+    /* Convex: only the quadratic minimum can be strictly inside, and F is
+     * monotonically decreasing from the segment start to that minimum, so the
+     * bracket still straddles the first crossing. */
+    const long double sigma_min_ld = -k->b / (2.0L * k->a);
+    if (!(sigma_min_ld > 0.0L)) {
+      route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED;
+      route->end_id = end_id;
+      return ASYMPTOTIC_INVALID;
+    }
+    double probe = (double)sigma_min_ld;
+    if (probe > s_segment)
+      probe = s_segment;
+    if (!(probe > 0.0)) {
+      route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED;
+      route->end_id = end_id;
+      return ASYMPTOTIC_INVALID;
+    }
+    AsymptoticRoute probe_state;
+    status = minkowski_entry_evaluate(&ctx, s_base + probe, &probe_state);
+    if (status != ASYMPTOTIC_OK) {
+      route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED;
+      route->end_id = end_id;
+      return status;
+    }
+    double f_probe = 0.0, tol_probe = 0.0;
+    status = asymptotic_entry_geometry(source, end_id, probe_state.activate_t,
+                                       probe_state.x, &f_probe, &tol_probe,
+                                       &why);
+    if (status != ASYMPTOTIC_OK) {
+      route->failure_reason = why;
+      route->end_id = end_id;
+      return status;
+    }
+    if (f_probe < 0.0) {
+      have_inside = 1;
+      inside_sigma = probe;
+    } else {
+      /* One reconstructed probe is not a miss proof: coordinate-time rounding
+       * can shift the minimum and hide an inside point at a nearby parameter. */
+      route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED;
+      route->end_id = end_id;
+      return ASYMPTOTIC_INVALID;
+    }
+  } else {
+    /* Monotone (a == 0 linear) or concave after the first crossing: step
+     * modestly past the candidate root and grow geometrically, staying inside
+     * the declared segment and the positive-radius domain. */
+    if (solve != ENTRY_SOLVE_ENTRY || !(candidate_sigma >= 0.0)) {
+      route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED;
+      route->end_id = end_id;
+      return ASYMPTOTIC_INVALID;
+    }
+    double sigma = candidate_sigma;
+    if (sigma > s_segment)
+      sigma = s_segment;
+    AsymptoticRoute probe_state;
+    status = minkowski_entry_evaluate(&ctx, s_base + sigma, &probe_state);
+    if (status != ASYMPTOTIC_OK) {
+      route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED;
+      route->end_id = end_id;
+      return status;
+    }
+    double f_sigma = 0.0, tol_sigma = 0.0;
+    status = asymptotic_entry_geometry(source, end_id, probe_state.activate_t,
+                                       probe_state.x, &f_sigma, &tol_sigma,
+                                       &why);
+    if (status != ASYMPTOTIC_OK) {
+      route->failure_reason = why;
+      route->end_id = end_id;
+      return status;
+    }
+    if (f_sigma < 0.0) {
+      have_inside = 1;
+      inside_sigma = sigma;
+    } else {
+      const long double slope =
+          2.0L * k->a * (long double)sigma + k->b; /* < 0 for an entry */
+      double delta = 0.0;
+      if (slope < 0.0L)
+        delta = 2.0 * fabs(f_sigma) / fabs((double)slope);
+      const double ulp_term =
+          16.0 * DBL_EPSILON * fmax(1.0, fabs(s_base + sigma));
+      if (!(delta > ulp_term))
+        delta = ulp_term;
+      for (int attempt = 0; attempt < 64 && !have_inside; ++attempt) {
+        const double probe = sigma + delta;
+        if (!(probe > sigma) || probe > s_segment)
+          break;
+        AsymptoticRoute grown_state;
+        status = minkowski_entry_evaluate(&ctx, s_base + probe, &grown_state);
+        if (status != ASYMPTOTIC_OK) {
+          route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED;
+          route->end_id = end_id;
+          return status;
+        }
+        double f_grown = 0.0, tol_grown = 0.0;
+        status = asymptotic_entry_geometry(source, end_id,
+                                           grown_state.activate_t,
+                                           grown_state.x, &f_grown, &tol_grown,
+                                           &why);
+        if (status != ASYMPTOTIC_OK) {
+          route->failure_reason = why;
+          route->end_id = end_id;
+          return status;
+        }
+        if (f_grown < 0.0) {
+          have_inside = 1;
+          inside_sigma = probe;
+        } else {
+          delta *= 2.0;
+        }
+      }
+      if (!have_inside) {
+        route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED;
+        route->end_id = end_id;
+        return ASYMPTOTIC_INVALID;
+      }
+    }
+  }
+
+  /* Bracket confirmed: hand it to the common, exterior-independent driver. */
+  AsymptoticRoute out;
+  unsigned int evaluations = 0;
+  RayReason localize_reason = RAY_REASON_NONE;
+  status = asymptotic_entry_localize(
+      source, end_id, minkowski_entry_evaluate, &ctx, s_base,
+      s_base + inside_sigma, &out, &evaluations, &localize_reason);
+  if (status == ASYMPTOTIC_OK) {
+    *route = out;
+    route->failure_reason = RAY_REASON_NONE;
+    route->entry_fallback_evaluations = evaluations;
+    return ASYMPTOTIC_OK;
+  }
+  if (status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED ||
+      status == ASYMPTOTIC_UNSUPPORTED) {
+    route->failure_reason = localize_reason;
+    route->end_id = end_id;
+    return status;
+  }
+  route->failure_reason =
+      (localize_reason == RAY_REASON_ESCAPE_LOCALIZATION_FAILED ||
+       localize_reason == RAY_REASON_NONE ||
+       localize_reason == RAY_REASON_PROTOCOL_ERROR)
+          ? RAY_REASON_ENTRY_UNCONFIRMED
+          : localize_reason;
+  route->end_id = end_id;
+  return ASYMPTOTIC_INVALID;
+}
+
+/* A fixed coordinate outside the sphere is an independent algebraic miss
+ * certificate, including exact tangency.  Restrict this cheap certificate to
+ * identity axes and zero origin so frame reconstruction cannot change the
+ * original camera component.  Sterbenz's lemma certifies the subtraction when
+ * both nonzero operands have the same sign and are within a factor of two. */
+static int minkowski_coordinate_miss(
+    const SpacetimeAsymptoticEnd *end, const double x_cur[3],
+    const double w_frame[3], const SpacetimeEscapeWorldtubeSample *sample) {
+  if (sample->radius_rate != 0.0)
+    return 0;
+  for (int i = 0; i < 3; ++i)
+    if (end->frame_origin[i] != 0.0)
+      return 0;
+  for (int i = 0; i < 3; ++i)
+    for (int j = 0; j < 3; ++j)
+      if (end->frame_axes[i][j] != (i == j ? 1.0 : 0.0))
+        return 0;
+  for (int i = 0; i < 3; ++i) {
+    if (w_frame[i] != 0.0 || sample->velocity[i] != 0.0)
+      continue;
+    const double x = x_cur[i] + end->frame_origin[i];
+    const double c = sample->center[i];
+    const int exact = x == 0.0 || c == 0.0 ||
+        (signbit(x) == signbit(c) && fabs(x) * 0.5 <= fabs(c) &&
+         fabs(c) * 0.5 <= fabs(x));
+    const double d = x - c;
+    if (isfinite(x) && isfinite(d) && exact && fabs(d) >= sample->radius)
+      return 1;
+  }
+  return 0;
+}
+
 static AsymptoticStatus minkowski_preroute(
     const SpacetimeSource *source, const SpacetimeAsymptoticEnd *end, double t0,
-    const double x_frame[3], const double w_frame[3], SpacetimeEndId end_id,
+    const double x_frame[3], const double w_frame[3], double log_alpha_p0,
-    AsymptoticRoute *route) {
+    SpacetimeEndId end_id, AsymptoticRoute *route) {
   /* Walk constant-velocity motion segments.  A quadratic root is only valid
    * inside the current segment and while the radius stays positive; otherwise
    * advance to the next segment boundary and re-sample. */
+  route->end_id = end_id;
+  route->failure_reason = RAY_REASON_NONE;
   double s = 0.0;
   for (int segment = 0; segment < 1000000; ++segment) {
     const double t = t0 - s;
@@ -313,32 +666,68 @@ static AsymptoticStatus minkowski_preroute(
        * path below. */
       return ASYMPTOTIC_UNSUPPORTED;
     }
-    double c_frame[3], v_frame[3], d[3], q[3];
+    double c_frame[3], v_frame[3];
     backend_position_to_frame(end, sample.center, c_frame);
     backend_vector_to_frame(end, sample.velocity, v_frame);
-    for (int i = 0; i < 3; ++i) {
-      d[i] = x_cur[i] - c_frame[i];
-      q[i] = w_frame[i] + v_frame[i];
-    }
     const double boundary = spacetime_escape_worldtube_next_segment(
         source, end->end_id, t);
     const double s_segment = isfinite(boundary) ? (t - boundary) : INFINITY;
     if (!(s_segment >= 0.0))
       return ASYMPTOTIC_INVALID;
-    double sigma;
+    const EntryQuadratic k = entry_quadratic_coeffs(
-    const int hit = solve_entry_quadratic(d, q, sample.radius,
+        x_cur, c_frame, w_frame, v_frame, sample.radius, sample.radius_rate);
-                                          sample.radius_rate, &sigma);
+    double sigma = 0.0;
+    EntrySolveResult solve = entry_solve(&k, &sigma);
+    if (solve == ENTRY_SOLVE_UNCERTAIN &&
+        minkowski_coordinate_miss(end, x_cur, w_frame, &sample))
+      solve = ENTRY_SOLVE_MISS;
     /* The backend constructor guarantees R > 0 throughout every segment, so
      * a root inside the segment is a real entry.  A root past the segment
      * boundary is not adopted here; the next segment is sampled instead.
      * The cheap R > 0 test at the root guards against a backend that
      * bypasses its constructor. */
-    if (hit && sigma >= 0.0 && sigma <= s_segment) {
+    if (solve == ENTRY_SOLVE_ENTRY && sigma >= 0.0 && sigma <= s_segment) {
       if (sample.radius - sample.radius_rate * sigma <= 0.0)
         return ASYMPTOTIC_INVALID;
+      AsymptoticRoute candidate = {.kind = ASYMPTOTIC_ROUTE_INVALID};
       minkowski_route_entry(end, t0, x_frame, w_frame, s + sigma, end_id,
-                            route);
+                            &candidate);
-      return ASYMPTOTIC_OK;
+      candidate.log_alpha_p0 = log_alpha_p0;
+      candidate.log_alpha_p0_camera = log_alpha_p0;
+      candidate.failure_reason = RAY_REASON_NONE;
+      int valid = 0;
+      RayReason why = RAY_REASON_NONE;
+      status = asymptotic_entry_validate(source, end_id, &candidate, &valid,
+                                         &why);
+      if (status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) {
+        route->failure_reason = why;
+        return status;
+      }
+      if (status != ASYMPTOTIC_OK) {
+        route->failure_reason = why;
+        route->end_id = end_id;
+        return status;
+      }
+      if (valid) {
+        /* Fast path: the reconstructed entry is already on the boundary
+         * within the geometric ULP band. */
+        *route = candidate;
+        return ASYMPTOTIC_OK;
+      }
+      /* The closed-form candidate lands off the reconstructed boundary:
+       * fall back to the common numerical localizer inside this segment. */
+      status = minkowski_fallback_entry(source, end, t0, x_frame, w_frame,
+                                        log_alpha_p0, end_id, s, s_segment, &k,
+                                        solve, sigma, route);
+      return status;
+    } else if (solve == ENTRY_SOLVE_UNCERTAIN) {
+      /* Near-tangent/degenerate discriminant: not a proof of a miss.  Attempt
+       * a strict-inside bracket; a positive reconstructed probe cannot prove
+       * that the continuous trajectory misses. */
+      status = minkowski_fallback_entry(source, end, t0, x_frame, w_frame,
+                                        log_alpha_p0, end_id, s, s_segment, &k,
+                                        solve, 0.0, route);
+      return status;
     }
     if (!isfinite(s_segment)) {
       /* Open final segment with no entry: a genuine miss. */
@@ -355,6 +744,127 @@ static AsymptoticStatus minkowski_preroute(
   return ASYMPTOTIC_INVALID;
 }
 
+/* Evaluate the exact Schwarzschild inward orbit from the original camera at
+ * the radius parameter p = -rho.  Recomputed per call (rotate + integral), not
+ * projected from an earlier state, and using no backend metric. */
+typedef struct {
+  const SpacetimeAsymptoticEnd *end;
+  const SchwarzschildCanonical *camera;
+  double log_alpha_p0_camera;
+  SpacetimeEndId end_id;
+} SchwarzschildEntryContext;
+
+static AsymptoticStatus schwarzschild_entry_evaluate(void *opaque,
+                                                     double parameter,
+                                                     AsymptoticRoute *state) {
+  const SchwarzschildEntryContext *ctx = opaque;
+  const double rho = -parameter;
+  double x[3], Pi[3], log_alpha_p0 = 0.0, activate_t = 0.0;
+  if (asymptotic_schwarzschild_inward_state_at_radius(
+          ctx->end, ctx->camera, rho, x, Pi, &log_alpha_p0, &activate_t))
+    return ASYMPTOTIC_INVALID;
+  *state = (AsymptoticRoute){0};
+  state->kind = ASYMPTOTIC_ROUTE_ENTRY;
+  state->end_id = ctx->end_id;
+  state->activate_t = activate_t;
+  for (int i = 0; i < 3; ++i) {
+    state->x[i] = x[i];
+    state->Pi[i] = Pi[i];
+  }
+  state->log_alpha_p0 = log_alpha_p0;
+  state->log_alpha_p0_camera = ctx->log_alpha_p0_camera;
+  state->failure_reason = RAY_REASON_NONE;
+  state->entry_fallback_evaluations = 0;
+  return ASYMPTOTIC_OK;
+}
+
+/* Bracket the first inward radius crossing when the closed-form entry state
+ * lands off the reconstructed worldtube boundary.  The parameter p = -rho
+ * increases inward from the original camera radius; the inside end is nudged
+ * just below the worldtube radius, and expanded inward only while staying
+ * above the turning radius and the rho > 2 state domain.  Failure is an
+ * explicit unconfirmed entry, never an escape. */
+static AsymptoticStatus schwarzschild_fallback_entry(
+    const SpacetimeSource *source, const SpacetimeAsymptoticEnd *end,
+    const SchwarzschildCanonical *camera, double worldtube_radius,
+    double log_alpha_p0_camera, AsymptoticRoute *route) {
+  SchwarzschildEntryContext ctx = {.end = end,
+                                   .camera = camera,
+                                   .log_alpha_p0_camera = log_alpha_p0_camera,
+                                   .end_id = end->end_id};
+  double floor = 2.0 + 1e-12 * fmax(1.0, worldtube_radius);
+  const double rho_turn = asymptotic_schwarzschild_turning_rho(camera->beta);
+  if (isfinite(rho_turn) && rho_turn > floor)
+    floor = rho_turn;
+  if (!(floor < worldtube_radius)) {
+    route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED;
+    route->end_id = end->end_id;
+    return ASYMPTOTIC_INVALID;
+  }
+  double rho_inside = nextafter(worldtube_radius, -INFINITY);
+  if (!(rho_inside > floor))
+    rho_inside = 0.5 * (worldtube_radius + floor);
+  double decrement = 0.0;
+  int found = 0;
+  for (int attempt = 0; attempt < 64; ++attempt) {
+    if (!(rho_inside > floor))
+      break;
+    double x[3], Pi[3], log_alpha_p0 = 0.0, activate_t = 0.0;
+    if (asymptotic_schwarzschild_inward_state_at_radius(
+            end, camera, rho_inside, x, Pi, &log_alpha_p0, &activate_t))
+      break;
+    double F = 0.0, tol = 0.0;
+    RayReason why = RAY_REASON_NONE;
+    const AsymptoticStatus st = asymptotic_entry_geometry(
+        source, end->end_id, activate_t, x, &F, &tol, &why);
+    if (st != ASYMPTOTIC_OK) {
+      route->failure_reason = why;
+      route->end_id = end->end_id;
+      return st;
+    }
+    if (F < 0.0) {
+      found = 1;
+      break;
+    }
+    if (decrement == 0.0)
+      decrement = (worldtube_radius - rho_inside) * 2.0;
+    else
+      decrement *= 2.0;
+    rho_inside = worldtube_radius - decrement;
+  }
+  if (!found) {
+    route->failure_reason = RAY_REASON_ENTRY_UNCONFIRMED;
+    route->end_id = end->end_id;
+    return ASYMPTOTIC_INVALID;
+  }
+  AsymptoticRoute out;
+  unsigned int evaluations = 0;
+  RayReason localize_reason = RAY_REASON_NONE;
+  const AsymptoticStatus status = asymptotic_entry_localize(
+      source, end->end_id, schwarzschild_entry_evaluate, &ctx, -camera->rho,
+      -rho_inside, &out, &evaluations, &localize_reason);
+  if (status == ASYMPTOTIC_OK) {
+    *route = out;
+    route->failure_reason = RAY_REASON_NONE;
+    route->entry_fallback_evaluations = evaluations;
+    return ASYMPTOTIC_OK;
+  }
+  if (status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED ||
+      status == ASYMPTOTIC_UNSUPPORTED) {
+    route->failure_reason = localize_reason;
+    route->end_id = end->end_id;
+    return status;
+  }
+  route->failure_reason =
+      (localize_reason == RAY_REASON_ESCAPE_LOCALIZATION_FAILED ||
+       localize_reason == RAY_REASON_NONE ||
+       localize_reason == RAY_REASON_PROTOCOL_ERROR)
+          ? RAY_REASON_ENTRY_UNCONFIRMED
+          : localize_reason;
+  route->end_id = end->end_id;
+  return ASYMPTOTIC_INVALID;
+}
+
 static AsymptoticStatus schwarzschild_route(
     const SpacetimeSource *source, const SpacetimeAsymptoticEnd *end,
     const MetricData *metric, const GeodesicRayState *state,
@@ -397,14 +907,38 @@ static AsymptoticStatus schwarzschild_route(
   }
   route->end_id = end->end_id;
   if (kind == SCH_ROUTE_ENTRY) {
-    route->kind = ASYMPTOTIC_ROUTE_ENTRY;
+    AsymptoticRoute candidate = {.kind = ASYMPTOTIC_ROUTE_ENTRY,
-    route->activate_t = activate_t;
+                                 .end_id = end->end_id,
+                                 .activate_t = activate_t,
+                                 .log_alpha_p0 = log_alpha_p0,
+                                 .log_alpha_p0_camera = state->log_alpha_p0,
+                                 .failure_reason = RAY_REASON_NONE};
     for (int i = 0; i < 3; ++i) {
-      route->x[i] = x[i];
+      candidate.x[i] = x[i];
-      route->Pi[i] = Pi[i];
+      candidate.Pi[i] = Pi[i];
     }
-    route->log_alpha_p0 = log_alpha_p0;
+    int valid = 0;
-    return ASYMPTOTIC_OK;
+    RayReason why = RAY_REASON_NONE;
+    const AsymptoticStatus validate_status = asymptotic_entry_validate(
+        source, end->end_id, &candidate, &valid, &why);
+    if (validate_status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) {
+      route->failure_reason = why;
+      route->kind = ASYMPTOTIC_ROUTE_TIME_RANGE_EXHAUSTED;
+      return validate_status;
+    }
+    if (validate_status != ASYMPTOTIC_OK) {
+      route->failure_reason = why;
+      return validate_status;
+    }
+    if (valid) {
+      *route = candidate;
+      return ASYMPTOTIC_OK;
+    }
+    /* The closed-form entry is off the reconstructed boundary: localize it in
+     * the radius parameter against the same exact inward transfer. */
+    return schwarzschild_fallback_entry(source, end, &camera,
+                                        sample.radius / end->mass,
+                                        state->log_alpha_p0, route);
   }
   route->kind = ASYMPTOTIC_ROUTE_ESCAPED;
   for (int i = 0; i < 3; ++i)
@@ -521,14 +1055,20 @@ AsymptoticStatus asymptotic_route_camera(const SpacetimeSource *source,
     AsymptoticRoute candidate = {.kind = ASYMPTOTIC_ROUTE_INVALID};
     const AsymptoticStatus status = minkowski_preroute(
         source, &end, state.coordinate_time, canonical.x, canonical.w,
-        end.end_id, &candidate);
+        state.log_alpha_p0, end.end_id, &candidate);
     if (status == ASYMPTOTIC_TIME_RANGE_EXHAUSTED) {
       route->kind = ASYMPTOTIC_ROUTE_TIME_RANGE_EXHAUSTED;
       route->end_id = end.end_id;
+      route->failure_reason = candidate.failure_reason;
       return ASYMPTOTIC_TIME_RANGE_EXHAUSTED;
     }
-    if (status != ASYMPTOTIC_OK)
+    if (status != ASYMPTOTIC_OK) {
+      /* Propagate the specific validation/fallback failure the segment walk
+       * refused with, not a generic preroute error. */
+      route->failure_reason = candidate.failure_reason;
+      route->end_id = candidate.end_id;
       return status;
+    }
     if (candidate.kind == ASYMPTOTIC_ROUTE_ENTRY) {
       const double s = state.coordinate_time - candidate.activate_t;
       if (!have_entry || s < best_s) {

src/asymptotic.h

@@ -48,6 +48,16 @@ typedef struct {
   /* Terminal infinity endpoint for ESCAPED. */
   double n_infinity[3];
   double frequency_ratio;
+  /* Diagnostic reason for an ASYMPTOTIC_INVALID return: set specifically by the
+   * validation/fallback failure that refused the route, so the lifecycle can
+   * report the concrete cause instead of a generic preroute failure.  NONE on
+   * success. */
+  RayReason failure_reason;
+  /* Nonzero when this route was produced by the generic bracketed first-entry
+   * localizer rather than a closed-form/fast entry solve.  It records the
+   * number of evaluator calls the localizer spent; zero on the fast path.  This
+   * is private in-memory provenance only and is not serialized. */
+  unsigned int entry_fallback_evaluations;
 } AsymptoticRoute;
 
 /* Pre-route one camera ray against every declared end's worldtube. */

src/asymptotic_entry.c

@@ -0,0 +1,326 @@
+#include "asymptotic_entry.h"
+
+#include <float.h>
+#include <math.h>
+#include <stddef.h>
+
+/* See asymptotic_entry.h for the contract.  This module deliberately keeps no
+ * global mutable state: every cache/scratch value lives on the stack of the
+ * calling trace, so it stays thread-safe under the coarse-grained OpenMP ray
+ * parallelism of the renderer. */
+
+static void set_failure(RayReason *failure, RayReason reason) {
+  if (failure != NULL)
+    *failure = reason;
+}
+
+/* Left-to-right double accumulation, matching the geodesic event layer's
+ * worldtube F. */
+static double dot3(const double a[3], const double b[3]) {
+  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
+}
+
+static int entry_state_finite(const AsymptoticRoute *state) {
+  if (!isfinite(state->activate_t) || !isfinite(state->log_alpha_p0) ||
+      !isfinite(state->log_alpha_p0_camera))
+    return 0;
+  for (int i = 0; i < 3; ++i)
+    if (!isfinite(state->x[i]) || !isfinite(state->Pi[i]))
+      return 0;
+  return 1;
+}
+
+AsymptoticStatus asymptotic_entry_geometry(const SpacetimeSource *source,
+                                           SpacetimeEndId end_id, double t,
+                                           const double x[3], double *F,
+                                           double *tol, RayReason *failure) {
+  set_failure(failure, RAY_REASON_PROTOCOL_ERROR);
+  if (source == NULL || x == NULL || F == NULL || tol == NULL) {
+    set_failure(failure, RAY_REASON_INVALID_ARGUMENT);
+    return ASYMPTOTIC_INVALID;
+  }
+  if (!isfinite(t)) {
+    set_failure(failure, RAY_REASON_WORLDTUBE_GEOMETRY_INVALID);
+    return ASYMPTOTIC_INVALID;
+  }
+  SpacetimeEscapeWorldtubeSample sample;
+  if (spacetime_escape_worldtube_sample(source, end_id, t, &sample)) {
+    set_failure(failure, RAY_REASON_WORLDTUBE_SAMPLE_FAILED);
+    return ASYMPTOTIC_INVALID;
+  }
+  if (!sample.valid) {
+    /* The backend cannot describe the worldtube at this time; this is history
+     * exhaustion, never a miss. */
+    set_failure(failure, RAY_REASON_TIME_RANGE_EXHAUSTED);
+    return ASYMPTOTIC_TIME_RANGE_EXHAUSTED;
+  }
+  if (!(sample.radius > 0.0) || !isfinite(sample.radius) ||
+      !isfinite(sample.radius_rate)) {
+    set_failure(failure, RAY_REASON_WORLDTUBE_GEOMETRY_INVALID);
+    return ASYMPTOTIC_INVALID;
+  }
+  for (int i = 0; i < 3; ++i) {
+    if (!isfinite(sample.center[i]) || !isfinite(sample.velocity[i])) {
+      set_failure(failure, RAY_REASON_WORLDTUBE_GEOMETRY_INVALID);
+      return ASYMPTOTIC_INVALID;
+    }
+  }
+  double d[3];
+  for (int i = 0; i < 3; ++i) {
+    d[i] = x[i] - sample.center[i];
+    if (!isfinite(d[i])) {
+      set_failure(failure, RAY_REASON_WORLDTUBE_GEOMETRY_INVALID);
+      return ASYMPTOTIC_INVALID;
+    }
+  }
+  /* Exact same grouping as the geodesic event layer: F uses
+   * dot3(d,d) - radius^2, while the geometric ULP band accumulates d2 with an
+   * explicit left-to-right loop.  Keeping both identical means a candidate that
+   * passes this validator at the tolerance threshold is grouped exactly like
+   * the geodesic's own f_before. */
+  const double value = dot3(d, d) - sample.radius * sample.radius;
+  const double r2 = sample.radius * sample.radius;
+  double d2 = 0.0;
+  for (int i = 0; i < 3; ++i)
+    d2 += d[i] * d[i];
+  const double geometry_tol = 128.0 * DBL_EPSILON * fmax(r2, d2);
+  /* Overflow to inf and NaN propagation both land here. */
+  if (!isfinite(value) || !isfinite(geometry_tol)) {
+    set_failure(failure, RAY_REASON_WORLDTUBE_GEOMETRY_INVALID);
+    return ASYMPTOTIC_INVALID;
+  }
+  *F = value;
+  *tol = geometry_tol;
+  set_failure(failure, RAY_REASON_NONE);
+  return ASYMPTOTIC_OK;
+}
+
+AsymptoticStatus asymptotic_entry_validate(const SpacetimeSource *source,
+                                           SpacetimeEndId end_id,
+                                           const AsymptoticRoute *candidate,
+                                           int *valid, RayReason *failure) {
+  set_failure(failure, RAY_REASON_PROTOCOL_ERROR);
+  if (source == NULL || candidate == NULL || valid == NULL) {
+    set_failure(failure, RAY_REASON_INVALID_ARGUMENT);
+    return ASYMPTOTIC_INVALID;
+  }
+  *valid = 0;
+  if (candidate->kind != ASYMPTOTIC_ROUTE_ENTRY) {
+    set_failure(failure, RAY_REASON_PROTOCOL_ERROR);
+    return ASYMPTOTIC_INVALID;
+  }
+  if (!entry_state_finite(candidate)) {
+    set_failure(failure, RAY_REASON_PROTOCOL_ERROR);
+    return ASYMPTOTIC_INVALID;
+  }
+  double value = 0.0, tol = 0.0;
+  const AsymptoticStatus status = asymptotic_entry_geometry(
+      source, end_id, candidate->activate_t, candidate->x, &value, &tol,
+      failure);
+  if (status != ASYMPTOTIC_OK)
+    return status;
+  /* Boundary-near means within the geometric ULP band on either side.  Outside
+   * and arbitrary-deep-inside are both non-candidates, not protocol errors. */
+  *valid = fabs(value) <= tol;
+  set_failure(failure, RAY_REASON_NONE);
+  return ASYMPTOTIC_OK;
+}
+
+/* Map an evaluator's own failure to a diagnostic reason.  The worldtube
+ * callback/history/geometry reasons come from asymptotic_entry_geometry; this
+ * only covers the case where the evaluator itself refuses to produce a state. */
+static AsymptoticStatus entry_evaluator_failure(AsymptoticStatus status,
+                                                RayReason *failure) {
+  switch (status) {
+  case ASYMPTOTIC_TIME_RANGE_EXHAUSTED:
+    set_failure(failure, RAY_REASON_TIME_RANGE_EXHAUSTED);
+    break;
+  case ASYMPTOTIC_UNSUPPORTED:
+    set_failure(failure, RAY_REASON_UNSUPPORTED);
+    break;
+  default:
+    /* The evaluator refused to produce a state, so the first entry is not
+     * confirmed. */
+    set_failure(failure, RAY_REASON_ENTRY_UNCONFIRMED);
+    break;
+  }
+  return status;
+}
+
+AsymptoticStatus asymptotic_entry_localize(
+    const SpacetimeSource *source, SpacetimeEndId end_id,
+    AsymptoticEntryEvaluator evaluate, void *context,
+    double outside_parameter, double inside_parameter, AsymptoticRoute *out,
+    unsigned int *evaluations, RayReason *failure) {
+  set_failure(failure, RAY_REASON_ENTRY_UNCONFIRMED);
+  if (evaluations != NULL)
+    *evaluations = 0;
+  if (source == NULL || evaluate == NULL || out == NULL) {
+    set_failure(failure, RAY_REASON_INVALID_ARGUMENT);
+    return ASYMPTOTIC_INVALID;
+  }
+  if (!isfinite(outside_parameter) || !isfinite(inside_parameter) ||
+      !(outside_parameter < inside_parameter)) {
+    set_failure(failure, RAY_REASON_INVALID_ARGUMENT);
+    return ASYMPTOTIC_INVALID;
+  }
+
+  unsigned int count = 0;
+  AsymptoticRoute lo_state, hi_state;
+  double f_lo = 0.0, f_hi = 0.0, tol_lo = 0.0, tol_hi = 0.0;
+
+  AsymptoticStatus status = evaluate(context, outside_parameter, &lo_state);
+  ++count;
+  if (status != ASYMPTOTIC_OK) {
+    if (evaluations != NULL)
+      *evaluations = count;
+    return entry_evaluator_failure(status, failure);
+  }
+  if (!entry_state_finite(&lo_state)) {
+    if (evaluations != NULL)
+      *evaluations = count;
+    set_failure(failure, RAY_REASON_ENTRY_UNCONFIRMED);
+    return ASYMPTOTIC_INVALID;
+  }
+  status = asymptotic_entry_geometry(source, end_id, lo_state.activate_t,
+                                     lo_state.x, &f_lo, &tol_lo, failure);
+  if (status != ASYMPTOTIC_OK) {
+    if (evaluations != NULL)
+      *evaluations = count;
+    return status;
+  }
+
+  status = evaluate(context, inside_parameter, &hi_state);
+  ++count;
+  if (status != ASYMPTOTIC_OK) {
+    if (evaluations != NULL)
+      *evaluations = count;
+    return entry_evaluator_failure(status, failure);
+  }
+  if (!entry_state_finite(&hi_state)) {
+    if (evaluations != NULL)
+      *evaluations = count;
+    set_failure(failure, RAY_REASON_ENTRY_UNCONFIRMED);
+    return ASYMPTOTIC_INVALID;
+  }
+  status = asymptotic_entry_geometry(source, end_id, hi_state.activate_t,
+                                     hi_state.x, &f_hi, &tol_hi, failure);
+  if (status != ASYMPTOTIC_OK) {
+    if (evaluations != NULL)
+      *evaluations = count;
+    return status;
+  }
+
+  if (!(f_lo >= 0.0) || !(f_hi < 0.0)) {
+    /* The caller owns the first-entry bracket; a bracket that does not straddle
+     * the boundary is an unconfirmed entry, never an escape. */
+    set_failure(failure, RAY_REASON_ENTRY_UNCONFIRMED);
+    if (evaluations != NULL)
+      *evaluations = count;
+    return ASYMPTOTIC_INVALID;
+  }
+
+  /* Parameter increases backward, so the inside end must not be later in
+   * coordinate time than the outside end. */
+  if (!(hi_state.activate_t <= lo_state.activate_t)) {
+    set_failure(failure, RAY_REASON_ENTRY_UNCONFIRMED);
+    if (evaluations != NULL)
+      *evaluations = count;
+    return ASYMPTOTIC_INVALID;
+  }
+
+  if (f_lo == 0.0) {
+    /* Exact boundary at the outside end with a strictly inside other end: this
+     * is the legitimate inward first-entry state already on the worldtube, so
+     * no bisection is needed. */
+    *out = lo_state;
+    out->kind = ASYMPTOTIC_ROUTE_ENTRY;
+    out->end_id = end_id;
+    if (evaluations != NULL)
+      *evaluations = count;
+    set_failure(failure, RAY_REASON_NONE);
+    return ASYMPTOTIC_OK;
+  }
+
+  double lo = outside_parameter;
+  double hi = inside_parameter;
+  unsigned int iteration = 0;
+  for (; iteration < ASYMPTOTIC_ENTRY_BISECTION_LIMIT; ++iteration) {
+    const double span = hi - lo;
+    const double mid = isfinite(span) ? lo + span * 0.5
+                                      : lo * 0.5 + hi * 0.5;
+    if (!(mid > lo && mid < hi))
+      break; /* Parameter midpoint cannot be represented; bracket is adjacent. */
+    AsymptoticRoute mid_state;
+    status = evaluate(context, mid, &mid_state);
+    ++count;
+    if (status != ASYMPTOTIC_OK) {
+      if (evaluations != NULL)
+        *evaluations = count;
+      return entry_evaluator_failure(status, failure);
+    }
+    if (!entry_state_finite(&mid_state) ||
+        !(hi_state.activate_t <= mid_state.activate_t &&
+          mid_state.activate_t <= lo_state.activate_t)) {
+      if (evaluations != NULL)
+        *evaluations = count;
+      set_failure(failure, RAY_REASON_ENTRY_UNCONFIRMED);
+      return ASYMPTOTIC_INVALID;
+    }
+    double f_mid = 0.0, tol_mid = 0.0;
+    status = asymptotic_entry_geometry(source, end_id, mid_state.activate_t,
+                                       mid_state.x, &f_mid, &tol_mid, failure);
+    if (status != ASYMPTOTIC_OK) {
+      if (evaluations != NULL)
+        *evaluations = count;
+      return status;
+    }
+    if (f_mid >= 0.0) {
+      lo = mid;
+      lo_state = mid_state;
+    } else {
+      hi = mid;
+      hi_state = mid_state;
+      f_hi = f_mid;
+      tol_hi = tol_mid;
+    }
+  }
+
+  if (iteration >= ASYMPTOTIC_ENTRY_BISECTION_LIMIT) {
+    /* The parameter interval never contracted to adjacent doubles within the
+     * implementation guard; do not fabricate an entry. */
+    set_failure(failure, RAY_REASON_ENTRY_UNCONFIRMED);
+    if (evaluations != NULL)
+      *evaluations = count;
+    return ASYMPTOTIC_INVALID;
+  }
+
+  /* Final consistency: keep the strict-inside adjacent endpoint.  Its negative
+   * residual need not fit the fast-path band at coarse coordinate resolution;
+   * the bracket, rather than a radial displacement, establishes the entry. */
+  if (!(f_hi < 0.0) ||
+      !(hi_state.activate_t <= lo_state.activate_t) ||
+      !isfinite(hi_state.activate_t) || !isfinite(hi_state.log_alpha_p0) ||
+      !isfinite(hi_state.log_alpha_p0_camera)) {
+    set_failure(failure, RAY_REASON_ENTRY_UNCONFIRMED);
+    if (evaluations != NULL)
+      *evaluations = count;
+    return ASYMPTOTIC_INVALID;
+  }
+  for (int i = 0; i < 3; ++i) {
+    if (!isfinite(hi_state.x[i]) || !isfinite(hi_state.Pi[i])) {
+      set_failure(failure, RAY_REASON_ENTRY_UNCONFIRMED);
+      if (evaluations != NULL)
+        *evaluations = count;
+      return ASYMPTOTIC_INVALID;
+    }
+  }
+
+  *out = hi_state;
+  out->kind = ASYMPTOTIC_ROUTE_ENTRY;
+  out->end_id = end_id;
+  if (evaluations != NULL)
+    *evaluations = count;
+  set_failure(failure, RAY_REASON_NONE);
+  return ASYMPTOTIC_OK;
+}

src/asymptotic_entry.h

@@ -0,0 +1,136 @@
+#ifndef ASYMPTOTIC_ENTRY_H
+#define ASYMPTOTIC_ENTRY_H
+
+#include "asymptotic.h"
+
+/* Backend-independent numerical entry localizer for the common asymptotic
+ * exterior.
+ *
+ * The camera pre-route must decide, without touching any backend metric outside
+ * a worldtube, whether a past-directed camera ray crosses an escape worldtube
+ * from outside to inside and where the FIRST such entry lies.  A supported
+ * exterior model may provide a closed-form (quadratic/analytic) entry; when it
+ * cannot, the caller supplies a path-parameter bracket that is already known to
+ * straddle the first entry and this module refines it numerically.
+ *
+ * The module never evaluates a metric, never loads a slab, never sweeps the
+ * movie in reverse and never snaps a position onto a radial shell.  It only
+ * consumes the worldtube `escape_worldtube_sample` callback through
+ * `asymptotic_entry_geometry`, and an opaque evaluator callback that
+ * repropagates the exact supported exterior geodesic from its camera state to a
+ * path parameter.  The driver is therefore independent of the exterior model
+ * and does not solve the entry equation itself.
+ *
+ * Parameter and time convention (18A.5/18A.6): the evaluator parameter
+ * increases along the renderer's backward propagation; the returned state's
+ * `activate_t` is the coordinate time at that parameter, so it is
+ * nonincreasing as the parameter increases.  `outside_parameter` is the
+ * smaller-parameter end (worldtube outside, or exactly on the boundary) and
+ * `inside_parameter` is the larger-parameter end that is strictly inside.
+ *
+ * The caller is responsible for providing a correct first-entry bracket.  This
+ * driver does not search arbitrary samples for a crossing; a bracket that does
+ * not straddle the boundary is reported as an unconfirmed entry
+ * (RAY_REASON_ENTRY_UNCONFIRMED), never as an escape. */
+typedef AsymptoticStatus (*AsymptoticEntryEvaluator)(void *context,
+                                                     double parameter,
+                                                     AsymptoticRoute *state);
+
+/* Convenience bundle for callers that want to keep the callback and its opaque
+ * context together.  Not required by any entry point. */
+typedef struct {
+  AsymptoticEntryEvaluator evaluate;
+  void *context;
+} AsymptoticEntryPropagator;
+
+/* Hard implementation guard on the number of bisection refinements.  This
+ * bounds the double-parameter bisection; it is not a physical parameter.  Any
+ * bracket near a finite nonzero entry contracts in roughly 60 halvings;
+ * very wide exponent ranges may instead exhaust the guard explicitly. */
+#define ASYMPTOTIC_ENTRY_BISECTION_LIMIT 256u
+
+/* Evaluate the worldtube function
+ *
+ *   F(t, x) = |x - center(t)|^2 - radius(t)^2
+ *
+ * at one state through the worldtube sample callback alone (no metric
+ * evaluation).  `*F` is accumulated with the same left-to-right double
+ * arithmetic as the geodesic event layer, and `*tol` is the matching geometric
+ * ULP band 128 * DBL_EPSILON * max(radius^2, |x-center|^2).
+ *
+ * Status / failure reason:
+ *   ASYMPTOTIC_OK                       -> *failure = RAY_REASON_NONE
+ *   ASYMPTOTIC_TIME_RANGE_EXHAUSTED     -> RAY_REASON_TIME_RANGE_EXHAUSTED
+ *   ASYMPTOTIC_INVALID (callback failed)-> RAY_REASON_WORLDTUBE_SAMPLE_FAILED
+ *   ASYMPTOTIC_INVALID (bad geometry)   -> RAY_REASON_WORLDTUBE_GEOMETRY_INVALID
+ *
+ * A non-positive/non-finite radius, non-finite radius_rate, centers or
+ * velocities, and any non-finite (overflow/NaN) F or tolerance are geometry
+ * failures.  `*failure` may be NULL. */
+AsymptoticStatus asymptotic_entry_geometry(const SpacetimeSource *source,
+                                           SpacetimeEndId end_id, double t,
+                                           const double x[3], double *F,
+                                           double *tol, RayReason *failure);
+
+/* Validate one caller-produced entry candidate.  The candidate must be finite
+ * (activate_t, x, Pi, log_alpha_p0, log_alpha_p0_camera), carry kind
+ * ASYMPTOTIC_ROUTE_ENTRY, and sit
+ * on the worldtube boundary within the geometric ULP band:
+ *
+ *   |F| <= tol  ->  *valid = 1
+ *
+ * A candidate outside the worldtube (F > tol) or arbitrarily deep inside
+ * (F < -tol) is NOT an entry candidate: it yields *valid = 0 but still returns
+ * ASYMPTOTIC_OK, because the caller owns the exterior solve and a non-candidate
+ * is not a protocol error.  A wrong kind or non-finite state is a protocol
+ * error.  Worldtube callback/history/geometry failures propagate with the same
+ * reasons as asymptotic_entry_geometry.
+ *
+ * This function does not fabricate a velocity or re-derive the entry: the
+ * caller already guarantees that its exterior solve produced an inward entry. */
+AsymptoticStatus asymptotic_entry_validate(const SpacetimeSource *source,
+                                           SpacetimeEndId end_id,
+                                           const AsymptoticRoute *candidate,
+                                           int *valid, RayReason *failure);
+
+/* Numerically locate the first outside -> inside entry inside a known bracket.
+ *
+ * `evaluate(context, parameter, state)` repropagates the exact exterior geodesic
+ * from its camera state to `parameter`; it must fill `state->activate_t` (the
+ * coordinate time at that parameter), `state->x`, `state->Pi`,
+ * `state->log_alpha_p0`, `state->log_alpha_p0_camera` and `state->end_id`.
+ * Every call is counted into `*evaluations` (may be NULL).
+ *
+ * The two bracket endpoints are evaluated first.  The driver requires
+ * F(outside_parameter) >= 0 and F(inside_parameter) < 0; a bracket with an
+ * exact boundary at the outside end (F == 0) and a strictly inside other end is
+ * returned directly as the legitimate inward entry.  Otherwise it bisects the
+ * parameter.  Each midpoint reconstructs a fresh state through the evaluator
+ * (no propagation from a prior midpoint, so no accumulated rounding or
+ * projection error), keeps the low end outside (F >= 0) and the high end
+ * strictly inside (F < 0), and never turns an F >= 0 midpoint into an escape.
+ *
+ * Bisection stops only when the parameter midpoint can no longer be represented
+ * strictly between the two ends (adjacent doubles), not on an F tolerance, so a
+ * curved trajectory is not stopped early by a coarser coordinate-time
+ * resolution.  On success the strictly-inside endpoint adjacent to the entry in
+ * path parameter is returned: its Pi/L/L_camera are copied from the evaluator
+ * state unchanged, with kind forced to ENTRY and end_id set to `end_id`.  The
+ * result is a representable bracketing of the entry (the true entry lies
+ * between the final outside and inside endpoints), not an absolute positional
+ * error claim; the returned inside state may legitimately have F < -tol.
+ *
+ * `ASYMPTOTIC_ENTRY_BISECTION_LIMIT` is an implementation guard on parameter
+ * representability, not a physical parameter.  If the bracket never contracts,
+ * if the endpoints do not straddle the boundary, or if the endpoint coordinate
+ * times are not ordered (inside time <= outside time), the result is
+ * ASYMPTOTIC_INVALID with RAY_REASON_ENTRY_UNCONFIRMED.  Callback, history and
+ * geometry failures from the evaluator/geometry propagate their exact status
+ * and reason; nothing is silently reported as a successful entry or escape. */
+AsymptoticStatus asymptotic_entry_localize(
+    const SpacetimeSource *source, SpacetimeEndId end_id,
+    AsymptoticEntryEvaluator evaluate, void *context,
+    double outside_parameter, double inside_parameter, AsymptoticRoute *out,
+    unsigned int *evaluations, RayReason *failure);
+
+#endif

src/asymptotic_schwarzschild.c

@@ -408,6 +408,47 @@ int asymptotic_schwarzschild_state_from_canonical(
   return 0;
 }
 
+int asymptotic_schwarzschild_inward_state_at_radius(
+    const SpacetimeAsymptoticEnd *end, const SchwarzschildCanonical *camera,
+    double rho, double x[3], double Pi[3], double *log_alpha_p0,
+    double *activate_t) {
+  if (end == NULL || camera == NULL || x == NULL || Pi == NULL)
+    return -1;
+  if (!(rho > 2.0) || !(rho <= camera->rho))
+    return -1;
+  /* The inward branch only exists while the orbit has not turned before the
+   * requested radius.  Allow a few ULP at the grazing limit so rounding in the
+   * turning root does not reject a legitimate boundary radius; a genuine
+   * inside-the-turning radius still fails through the Q >= 0 check below. */
+  const double rho_turn = asymptotic_schwarzschild_turning_rho(camera->beta);
+  if (isfinite(rho_turn) &&
+      rho < rho_turn - 16.0 * DBL_EPSILON * fmax(1.0, rho_turn))
+    return -1;
+  const double dphi = asymptotic_schwarzschild_phi(rho, camera->beta) -
+                      asymptotic_schwarzschild_phi(camera->rho, camera->beta);
+  if (!isfinite(dphi))
+    return -1;
+  double rhat_rho[3];
+  if (camera->beta > 0.0)
+    rotate_axis(camera->rhat, camera->Lhat, -dphi, rhat_rho);
+  else
+    for (int i = 0; i < 3; ++i)
+      rhat_rho[i] = camera->rhat[i];
+  SchwarzschildCanonical state = *camera;
+  state.rho = rho;
+  for (int i = 0; i < 3; ++i)
+    state.rhat[i] = rhat_rho[i];
+  state.radial_sign = -1;
+  if (asymptotic_schwarzschild_state_from_canonical(end, &state, x, Pi,
+                                                    log_alpha_p0))
+    return -1;
+  if (activate_t != NULL) {
+    const double T = sch_time_transfer(camera->rho, rho, camera->beta);
+    *activate_t = camera->t - end->mass * T;
+  }
+  return 0;
+}
+
 int asymptotic_schwarzschild_finish(const SpacetimeAsymptoticEnd *end,
                                     const SchwarzschildCanonical *canonical,
                                     double n_infinity[3],
@@ -476,25 +517,10 @@ int asymptotic_schwarzschild_preroute(
   }
 
   if (camera->beta < beta_R) {
-    const double dphi = asymptotic_schwarzschild_phi(R, camera->beta) -
+    /* Shared exact inward transfer, also used by the common entry fallback. */
-                        asymptotic_schwarzschild_phi(camera->rho, camera->beta);
+    if (asymptotic_schwarzschild_inward_state_at_radius(
-    double rhat_entry[3];
+            end, camera, R, x, Pi, log_alpha_p0, activate_t))
-    if (camera->beta > 0.0)
-      rotate_axis(camera->rhat, camera->Lhat, -dphi, rhat_entry);
-    else
-      for (int i = 0; i < 3; ++i)
-        rhat_entry[i] = camera->rhat[i];
-    SchwarzschildCanonical entry = *camera;
-    entry.rho = R;
-    for (int i = 0; i < 3; ++i)
-      entry.rhat[i] = rhat_entry[i];
-    entry.radial_sign = -1;
-    if (asymptotic_schwarzschild_state_from_canonical(end, &entry, x, Pi,
-                                                      log_alpha_p0))
       return -1;
-    const double T =
-        sch_time_transfer(camera->rho, R, camera->beta);
-    *activate_t = camera->t - end->mass * T;
     *kind = SCH_ROUTE_ENTRY;
     return 0;
   }

src/asymptotic_schwarzschild.h

@@ -42,6 +42,19 @@ int asymptotic_schwarzschild_state_from_canonical(
     const SpacetimeAsymptoticEnd *end, const SchwarzschildCanonical *canonical,
     double x[3], double Pi[3], double *log_alpha_p0);
 
+/* Repropagate the exact inward orbit from the ORIGINAL camera canonical state
+ * to `rho >= turning_rho` (and > 2) by recomputing the swept angle and the
+ * coordinate-time integral -- never by projecting a nearby state.  Fills the
+ * backend entry state, its local L = ln(alpha p^0), and the activation
+ * coordinate time.  Returns -1 when `rho` lies outside the reachable inward
+ * domain.  This is the radius-parameter evaluator used by the common first
+ * entry localizer; `R` may dip slightly below the worldtube radius because the
+ * common driver only needs a strictly-inside bracket. */
+int asymptotic_schwarzschild_inward_state_at_radius(
+    const SpacetimeAsymptoticEnd *end, const SchwarzschildCanonical *camera,
+    double rho, double x[3], double Pi[3], double *log_alpha_p0,
+    double *activate_t);
+
 /* Infinity endpoint for an outward crossing at the canonical radius. */
 int asymptotic_schwarzschild_finish(const SpacetimeAsymptoticEnd *end,
                                     const SchwarzschildCanonical *canonical,

src/geodesic.c

@@ -83,6 +83,8 @@ const char *ray_reason_name(RayReason reason) {
     return "THRESHOLD_EVENT_UNCONFIRMED";
   case RAY_REASON_SLAB_LOAD_FAILED:
     return "SLAB_LOAD_FAILED";
+  case RAY_REASON_ENTRY_UNCONFIRMED:
+    return "ENTRY_UNCONFIRMED";
   case RAY_REASON_COUNT:
     break;
   }
@@ -136,6 +138,7 @@ RayReason ray_reason_category(RayReason reason) {
   case RAY_REASON_SUBINTEGRATION_TARGET_MISSED:
   case RAY_REASON_ESCAPE_EVENT_UNCONFIRMED:
   case RAY_REASON_THRESHOLD_EVENT_UNCONFIRMED:
+  case RAY_REASON_ENTRY_UNCONFIRMED:
     return RAY_REASON_INTEGRATION_ERROR;
   case RAY_REASON_SLAB_LOAD_FAILED:
     return RAY_REASON_IO_ERROR;
@@ -2169,7 +2172,10 @@ RayEndpoint geodesic_trace_past(const SpacetimeSource *source,
   }
   if (route_status != ASYMPTOTIC_OK) {
     out.outcome = RAY_OUTCOME_INCOMPLETE;
-    out.reason = RAY_REASON_CAMERA_PREROUTE_FAILED;
+    out.reason = ray_reason_valid(route.failure_reason) &&
+                         route.failure_reason != RAY_REASON_NONE
+                     ? route.failure_reason : RAY_REASON_CAMERA_PREROUTE_FAILED;
+    out.end_id = route.end_id;
     return out;
   }
   if (route.kind == ASYMPTOTIC_ROUTE_ESCAPED) {

src/geodesic.h

@@ -64,6 +64,7 @@ typedef enum {
   RAY_REASON_THRESHOLD_EVENT_UNCONFIRMED,  /* threshold event retry exhausted */
   /* I/O-derived. */
   RAY_REASON_SLAB_LOAD_FAILED,             /* spacetime_load_slab failed */
+  RAY_REASON_ENTRY_UNCONFIRMED,            /* no trustworthy first-entry bracket */
   RAY_REASON_COUNT                        /* sentinel: valid ids are < COUNT */
 } RayReason;
 

src/lens_map.c

@@ -297,6 +297,7 @@ int lens_map_read(const char *path, LensMapProvenance *provenance,
           read_u32(file, &end_id, &crc) || read_u32(file, &outcome, &crc) ||
           read_u32(file, &reason, &crc) || outcome > RAY_OUTCOME_INCOMPLETE ||
           !ray_reason_valid((RayReason)reason);
+      if (failed) break;
       v->end_id = (SpacetimeEndId)end_id;
       v->outcome = (RayOutcome)outcome;
       v->reason = (RayReason)reason;

src/ray.c

@@ -149,8 +149,11 @@ void ray_pool_preroute(RayPool *p, const SpacetimeSource *source) {
     if (status == ASYMPTOTIC_UNSUPPORTED || status == ASYMPTOTIC_INVALID) {
       p->endpoint[i].outcome = RAY_OUTCOME_INCOMPLETE;
       p->endpoint[i].reason = status == ASYMPTOTIC_UNSUPPORTED
-                                  ? RAY_REASON_UNSUPPORTED
+                                   ? RAY_REASON_UNSUPPORTED
-                                  : RAY_REASON_CAMERA_PREROUTE_FAILED;
+                                   : (ray_reason_valid(route.failure_reason) &&
+                                              route.failure_reason != RAY_REASON_NONE
+                                          ? route.failure_reason
+                                          : RAY_REASON_CAMERA_PREROUTE_FAILED);
       p->endpoint[i].end_id = route.end_id;
       p->status[i] = RAY_POOL_FAILED;
       continue;

tests/test_asymptotic.c

@@ -3,6 +3,7 @@
 #include "ray.h"
 #include "spacetime.h"
 
+#include <float.h>
 #include <math.h>
 #include <stdio.h>
 
@@ -50,6 +51,7 @@ typedef struct {
   int sample_nonpositive_radius;
   int fail_on_sample_call; /* 1-based callback invocation to fail. */
   int sample_call_count;
+  double frame_origin[3];
 } SyntheticContext;
 
 static SpacetimePointStatus synthetic_eval(const SpacetimeSource *source,
@@ -95,7 +97,8 @@ static int synthetic_end(const SpacetimeSource *source, size_t index,
       .end_id = 0,
       .exterior_kind = kind,
       .mass = mass,
-      .frame_origin = {0.0, 0.0, 0.0},
+      .frame_origin = {context->frame_origin[0], context->frame_origin[1],
+                        context->frame_origin[2]},
       .frame_axes = {{1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}};
   return 0;
 }
@@ -134,12 +137,17 @@ static int synthetic_worldtube(const SpacetimeSource *source,
     *out = (SpacetimeEscapeWorldtubeSample){.valid = 0};
     return 0;
   }
+  /* The sample contract is the worldtube value at this coordinate time, so the
+   * reported radius must carry its own time dependence: R(t) = R0 + rr t, with
+   * dR/dt = radius_rate.  A constant `radius` with a nonzero rate would make
+   * the closed-form segment model and the callback geometry disagree. */
+  const double radius_t = context->radius + context->radius_rate * t;
   if (context->has_segment && t < context->segment_t) {
     /* Second segment: center moves toward +x as t decreases. */
     *out = (SpacetimeEscapeWorldtubeSample){
         .center = {context->segment_t - t, 0.0, 0.0},
         .velocity = {-1.0, 0.0, 0.0},
-        .radius = context->radius,
+        .radius = radius_t,
         .radius_rate = context->radius_rate,
         .velocity_constant = context->constant,
         .valid = 1};
@@ -148,7 +156,7 @@ static int synthetic_worldtube(const SpacetimeSource *source,
   *out = (SpacetimeEscapeWorldtubeSample){
       .center = {context->vx * t + 0.5 * context->accel * t * t, 0.0, 0.0},
       .velocity = {context->vx + context->accel * t, 0.0, 0.0},
-      .radius = context->radius,
+      .radius = radius_t,
       .radius_rate = context->radius_rate,
       .velocity_constant = context->constant,
       .valid = 1};
@@ -710,6 +718,166 @@ static void test_moving_sphere(void) {
         "co-moving ray misses");
 }
 
+/* Fixed-observer tetrad used by the production Alcubierre observer-track rows
+ * 63/64, with spatial axes (e1, e2, e3) = (y-hat, z-hat, x-hat).  The literal
+ * values are embedded here so this regression does not depend on the
+ * untracked observer CSV. */
+static ObserverState track_observer(double coordinate_time) {
+  ObserverState o = {0};
+  o.coordinate_time = coordinate_time;
+  o.coordinate_position[0] = 0.0;
+  o.coordinate_position[1] = -24.0;
+  o.coordinate_position[2] = 0.0;
+  o.tetrad[0][0] = 1.0;
+  o.tetrad[1][2] = 1.0;
+  o.tetrad[2][3] = 1.0;
+  o.tetrad[3][1] = 1.0;
+  return o;
+}
+
+/* Two exact production RayPool pre-route samples (frame 63 sample 12315 and
+ * frame 64 sample 3994).  They are grazing (disc/b^2 ~ 1e-4), so the plain
+ * double root solve left the reconstructed entry state at F ~ 1.0-1.2 x
+ * geom_tol, which the event layer rejected as OUTSIDE_WORLDTUBE.  The moving
+ * sphere fixture reproduces the Alcubierre worldtube (center = 2 t, radius 5);
+ * the observer time/position/tetrad and the camera direction are the exact raw
+ * production values.  The route must land inside the geometric tolerance with
+ * the entry direction (Pi) unchanged. */
+static void test_grazing_production_entries(void) {
+  SyntheticContext context = {.vx = 2.0,
+                              .accel = 0.0,
+                              .radius = 5.0,
+                              .radius_rate = 0.0,
+                              .valid_t_min = -1.0e30,
+                              .constant = 1};
+  SpacetimeSource source = {.ops = &synthetic_ops, .context = &context};
+  struct {
+    double time;
+    double direction[3];
+    double pi[3];
+  } cases[2] = {
+      {0x1.fa8f5c28f5c29p+3,
+       {0x1.c7378f8e872d1p-1, 0x1.15bad4e30e8ddp-8, -0x1.d4afba4704cap-2},
+       {0x1.d4afba4704cap-2, -0x1.c7378f8e872d1p-1, -0x1.15bad4e30e8ddp-8}},
+      {0x1.fb17e4b17e4b1p+3,
+       {0x1.bd3bb364ac492p-1, 0x1.102d2a1c6ac74p-7, -0x1.f98ae1a782104p-2},
+       {0x1.f98ae1a782104p-2, -0x1.bd3bb364ac492p-1, -0x1.102d2a1c6ac74p-7}},
+  };
+  for (int c = 0; c < 2; ++c) {
+    ObserverState observer = track_observer(cases[c].time);
+    AsymptoticRoute route;
+    CHECK(asymptotic_route_camera(&source, &observer, cases[c].direction,
+                                  &route) == ASYMPTOTIC_OK &&
+              route.kind == ASYMPTOTIC_ROUTE_ENTRY,
+          "grazing production entry found");
+    CHECK(route.entry_fallback_evaluations == 0,
+          "grazing production entry stays on the fast path");
+    for (int i = 0; i < 3; ++i)
+      CHECK(route.Pi[i] == cases[c].pi[i],
+            "grazing entry direction is unchanged");
+    SpacetimeEscapeWorldtubeSample sample;
+    CHECK(spacetime_escape_worldtube_sample(&source, route.end_id,
+                                            route.activate_t, &sample) == 0 &&
+              sample.valid && sample.radius > 0.0,
+          "grazing entry sample valid");
+    double d2 = 0.0;
+    for (int k = 0; k < 3; ++k) {
+      const double dk = route.x[k] - sample.center[k];
+      d2 += dk * dk;
+    }
+    const double r2 = sample.radius * sample.radius;
+    double value;
+    CHECK(asymptotic_worldtube_value(&source, route.end_id, route.activate_t,
+                                     route.x, &value) == 0,
+          "grazing entry worldtube value");
+    const double geom_tol = 128.0 * DBL_EPSILON * fmax(r2, d2);
+    /* The event layer rejects the entry (OUTSIDE_WORLDTUBE) exactly when
+     * F > geom_tol; require the residual to sit inside the tolerance band
+     * rather than accepting an arbitrary sign. */
+    CHECK(value <= geom_tol && value >= -geom_tol,
+          "grazing entry F within geometric tolerance");
+  }
+}
+
+/* Linear (a == 0) entry: a growing sphere whose radius rate cancels the
+ * relative closing speed, so qq == rr^2 and the quadratic degenerates.  The
+ * stable solver must still take the smallest positive root.  The fixture's
+ * sample() reports the consistent radius R(t) = 10 - t, i.e. R(s) = 10 + s
+ * along the past parameter s = -t, so the contact point is on the true
+ * ruled-surface boundary. */
+static void test_linear_a_zero_entry(void) {
+  SyntheticContext context = {.vx = 0.0,
+                              .accel = 0.0,
+                              .radius = 10.0,
+                              .radius_rate = -1.0, /* R(t) = 10 - t */
+                              .valid_t_min = -1.0e30,
+                              .constant = 1};
+  SpacetimeSource source = {.ops = &synthetic_ops, .context = &context};
+  const ObserverState camera = flat_observer(100.0, 0.0, 0.0);
+  AsymptoticRoute route;
+  CHECK(asymptotic_route_camera(&source, &camera, (double[]){-1.0, 0.0, 0.0},
+                                &route) == ASYMPTOTIC_OK &&
+            route.kind == ASYMPTOTIC_ROUTE_ENTRY,
+        "linear a~0 entry");
+  CHECK(fabs(route.activate_t + 45.0) < 1e-12, "linear entry time");
+  CHECK(fabs(route.x[0] - 55.0) < 1e-12 && fabs(route.x[1]) < 1e-12 &&
+            fabs(route.x[2]) < 1e-12,
+        "linear entry position");
+}
+
+/* Large-coordinate-time cancellation: the long-double closed-form root is
+ * accurate in the frame, but the entry state reconstructed in double at a huge
+ * t0 loses the sub-ULP part of the event and lands far outside the geometric
+ * ULP band (F ~ 0.3 >> tol).  The common fallback must repropagate from the
+ * original camera, localize the first entry numerically, and return a
+ * strict-inside endpoint (F < 0) while preserving the camera direction Pi and
+ * reference L exactly.  The fixture is a legitimate constant-velocity
+ * worldtube, not a nonlinearity injection. */
+static void test_fallback_reconstruction_cancellation(void) {
+  SyntheticContext context = {.vx = 0x1.999999999999ap-4, /* 0.1 */
+                              .accel = 0.0,
+                              .radius = 10.0,
+                              .radius_rate = 0.0,
+                              .valid_t_min = -1.0e30,
+                              .constant = 1};
+  SpacetimeSource source = {.ops = &synthetic_ops, .context = &context};
+  const double t0 = 1.0e15;
+  /* Exactly 0.1 * 1e15 + 100.123456789, pinned as a hex literal.  The offset
+   * is not aligned to the t0 ULP, so activate_t = t0 - s rounds and the
+   * reconstructed boundary residual exceeds the tolerance band. */
+  const double camera_x = 0x1.6bcc41e901908p+46;
+  ObserverState observer = flat_observer(camera_x, 0.0, 0.0);
+  observer.coordinate_time = t0;
+  const double direction[3] = {-1.0, 0.0, 0.0};
+  AsymptoticRoute route;
+  CHECK(asymptotic_route_camera(&source, &observer, direction, &route) ==
+                ASYMPTOTIC_OK &&
+            route.kind == ASYMPTOTIC_ROUTE_ENTRY,
+        "cancellation fallback still produces an entry");
+  CHECK(route.entry_fallback_evaluations > 0,
+        "cancellation entry used the common fallback");
+  CHECK(route.failure_reason == RAY_REASON_NONE,
+        "successful fallback has no failure reason");
+  double value;
+  CHECK(asymptotic_worldtube_value(&source, route.end_id, route.activate_t,
+                                   route.x, &value) == 0 &&
+            value <= 0.0,
+        "fallback entry state is inside the worldtube");
+
+  MetricData metric;
+  GeodesicRayState camera_state;
+  CHECK(spacetime_eval(&source, t0, observer.coordinate_position, &metric) ==
+                0 &&
+            geodesic_initialize_past_ray_metric(&metric, &observer, direction,
+                                                &camera_state) == 0,
+        "cancellation camera state");
+  for (int i = 0; i < 3; ++i)
+    CHECK(route.Pi[i] == camera_state.Pi[i],
+          "fallback preserves the entry direction exactly");
+  CHECK(route.log_alpha_p0_camera == camera_state.log_alpha_p0,
+        "fallback preserves the camera reference L exactly");
+}
+
 static void test_accelerated_worldtube_unsupported(void) {
   /* A genuinely accelerating (non-constant velocity) worldtube has no strict
    * relative-motion interval bound, so the route is explicitly unsupported.
@@ -789,11 +957,98 @@ static void test_ray_pool_lifecycle(void) {
   spacetime_destroy(&source);
 }
 
+static void test_zero_discriminant_is_not_miss(void) {
+  SpacetimeSource source;
+  CHECK(spacetime_create_minkowski(&source, 1.0) == 0,
+        "zero-discriminant source");
+  const ObserverState observer = flat_observer(1e10, 0.5, 0.0);
+  AsymptoticRoute route;
+  /* Forming c = 1e20 + 0.25 - 1 loses the transverse contribution even in
+   * 80-bit arithmetic; b*b - 4*a*c then rounds to zero despite a real entry. */
+  CHECK(asymptotic_route_camera(&source, &observer,
+                                (double[]){-1.0, 0.0, 0.0}, &route) ==
+            ASYMPTOTIC_OK && route.kind == ASYMPTOTIC_ROUTE_ENTRY &&
+            route.entry_fallback_evaluations > 0,
+        "rounded zero discriminant uses entry fallback, not escape");
+  double value = 0.0;
+  CHECK(asymptotic_worldtube_value(&source, 0, route.activate_t, route.x,
+                                   &value) == ASYMPTOTIC_OK && value < 0.0,
+        "zero-discriminant fallback produces actual inside state");
+  const ObserverState tangent = flat_observer(1e10, 1.0, 0.0);
+  CHECK(asymptotic_route_camera(&source, &tangent,
+                                (double[]){-1.0, 0.0, 0.0}, &route) ==
+            ASYMPTOTIC_OK && route.kind == ASYMPTOTIC_ROUTE_ESCAPED,
+        "fixed transverse coordinate independently certifies exact tangency");
+  spacetime_destroy(&source);
+}
+
+static void test_positive_reconstructed_minimum_is_not_miss(void) {
+  SyntheticContext ctx = {.vx = 0.1, .radius = 10.0, .constant = 1};
+  SpacetimeSource source = {.ops = &synthetic_ops, .context = &ctx};
+  ObserverState observer = flat_observer(0x1.6bcc41e901904p+46,
+                                         0x1.3ffffde7210bfp+3, 0.0);
+  observer.coordinate_time = 1e15;
+  const double parameter = 0x1.bca8814065f1ep+6;
+  const double witness[3] = {observer.coordinate_position[0] - parameter,
+                             observer.coordinate_position[1], 0.0};
+  double value = 0.0;
+  CHECK(asymptotic_worldtube_value(&source, 0,
+                                   observer.coordinate_time - parameter,
+                                   witness, &value) == ASYMPTOTIC_OK && value < 0,
+        "strict-inside witness exists despite positive reconstructed minimum");
+  AsymptoticRoute route;
+  CHECK(asymptotic_route_camera(&source, &observer,
+                                (double[]){-1.0, 0.0, 0.0}, &route) ==
+            ASYMPTOTIC_INVALID &&
+            route.failure_reason == RAY_REASON_ENTRY_UNCONFIRMED,
+        "positive probe without a miss certificate fails explicitly");
+  RayPool pool;
+  CHECK(ray_pool_init(&pool, 1) == 0, "ambiguous entry pool");
+  CHECK(ray_pool_append(&pool, &observer, (double[]){-1.0, 0.0, 0.0},
+                        0, 0) == 0, "ambiguous entry ray");
+  ray_pool_preroute(&pool, &source);
+  CHECK(pool.status[0] == RAY_POOL_FAILED &&
+            pool.endpoint[0].outcome == RAY_OUTCOME_INCOMPLETE &&
+            pool.endpoint[0].reason == RAY_REASON_ENTRY_UNCONFIRMED,
+        "pool propagates unconfirmed entry instead of fabricating escape");
+  ray_pool_destroy(&pool);
+}
+
+static void test_translated_frame_cannot_certify_miss(void) {
+  SyntheticContext ctx = {.radius = 1.0, .constant = 1,
+                           .valid_t_min = -1e100,
+                           .frame_origin = {0.0, 1e10, 0.0}};
+  SpacetimeSource source = {.ops = &synthetic_ops, .context = &ctx};
+  const ObserverState observer = flat_observer(1e10, 1.0 - 0x1p-22, 0.0);
+  double value = 0.0;
+  const double witness[3] = {0.0, observer.coordinate_position[1], 0.0};
+  CHECK(asymptotic_worldtube_value(&source, 0, -1e10, witness, &value) ==
+            ASYMPTOTIC_OK && value < 0.0,
+        "original untranslated trajectory has an inside witness");
+  AsymptoticRoute route;
+  CHECK(asymptotic_route_camera(&source, &observer,
+                                (double[]){-1.0, 0.0, 0.0}, &route) ==
+            ASYMPTOTIC_INVALID &&
+            route.failure_reason == RAY_REASON_ENTRY_UNCONFIRMED,
+        "lossy translated frame must not certify a miss");
+  ctx.frame_origin[1] = 0.0;
+  CHECK(asymptotic_route_camera(&source, &observer,
+                                (double[]){-1.0, 0.0, 0.0}, &route) ==
+            ASYMPTOTIC_OK && route.kind == ASYMPTOTIC_ROUTE_ENTRY,
+        "same unshifted trajectory confirms an entry");
+}
+
 int main(void) {
   test_fixed_sphere();
   test_large_radius_quadratic();
   test_round_trip();
   test_moving_sphere();
+  test_grazing_production_entries();
+  test_linear_a_zero_entry();
+  test_fallback_reconstruction_cancellation();
+  test_zero_discriminant_is_not_miss();
+  test_positive_reconstructed_minimum_is_not_miss();
+  test_translated_frame_cannot_certify_miss();
   test_accelerated_worldtube_unsupported();
   test_piecewise_segment_entry();
   test_boundary_semantics_minkowski();

tests/test_asymptotic_entry.c

@@ -0,0 +1,781 @@
+#include "asymptotic_entry.h"
+#include "spacetime.h"
+
+#include <float.h>
+#include <math.h>
+#include <stdio.h>
+
+/* Backend-independent core regression for the numerical entry localizer.  It
+ * deliberately links no analytic backend and no geodesic integrator: the fake
+ * SpacetimeSource exposes only `escape_worldtube_sample` (plus a deliberately
+ * trapped `eval`), and the evaluator is an analytic path-parameter callback.
+ *
+ * Nothing here depends on an untracked production track, CSV or binary. */
+
+static int failures = 0;
+
+#define CHECK(condition, message)                                        \
+  do {                                                                   \
+    if (!(condition)) {                                                  \
+      fprintf(stderr, "FAIL %s:%d: %s\n", __FILE__, __LINE__, message);  \
+      ++failures;                                                        \
+    }                                                                    \
+  } while (0)
+
+#ifndef TEST_PI
+#define TEST_PI 3.14159265358979323846
+#endif
+
+/* ------------------------------------------------------------------ */
+/* Fake worldtube source                                              */
+/* ------------------------------------------------------------------ */
+
+typedef struct {
+  double center0[3];
+  double center_vel[3]; /* dc/dt */
+  double radius0;
+  double radius_rate; /* dR/dt */
+  double valid_t_min; /* sample is valid for t >= valid_t_min */
+  int callback_fails; /* always return -1 */
+  int fail_at_call;   /* 1-based sample-call index to fail, 0 disabled */
+  int nan_radius;
+  int zero_radius;
+  double hole_center;    /* isolated invalid time window */
+  double hole_halfwidth; /* 0 disables the window */
+  int call_count;
+  int eval_calls; /* trap: how often the metric eval callback ran */
+} EntryWorldtube;
+
+static SpacetimePointStatus entry_eval_trap(const SpacetimeSource *source,
+                                            double t, const double x[3],
+                                            MetricData *metric) {
+  EntryWorldtube *wt = source->context;
+  ++wt->eval_calls;
+  (void)t;
+  (void)x;
+  (void)metric;
+  /* This source is deliberately outside the metric domain.  The localizer must
+   * never reach here because it does no metric evaluation. */
+  return SPACETIME_POINT_OUT_OF_DOMAIN;
+}
+
+static int entry_worldtube_cb(const SpacetimeSource *source,
+                              SpacetimeEndId end_id, double t,
+                              SpacetimeEscapeWorldtubeSample *out) {
+  EntryWorldtube *wt = source->context;
+  if (end_id != 0)
+    return -1;
+  ++wt->call_count;
+  if (wt->fail_at_call > 0 && wt->call_count == wt->fail_at_call)
+    return -1;
+  if (wt->callback_fails)
+    return -1;
+  if (!isfinite(t) || t < wt->valid_t_min) {
+    *out = (SpacetimeEscapeWorldtubeSample){.valid = 0};
+    return 0;
+  }
+  if (wt->hole_halfwidth > 0.0 &&
+      fabs(t - wt->hole_center) <= wt->hole_halfwidth) {
+    *out = (SpacetimeEscapeWorldtubeSample){.valid = 0};
+    return 0;
+  }
+  if (wt->nan_radius) {
+    *out = (SpacetimeEscapeWorldtubeSample){.radius = NAN, .valid = 1};
+    return 0;
+  }
+  if (wt->zero_radius) {
+    *out = (SpacetimeEscapeWorldtubeSample){.radius = 0.0, .valid = 1};
+    return 0;
+  }
+  *out = (SpacetimeEscapeWorldtubeSample){
+      .center = {wt->center0[0] + wt->center_vel[0] * t,
+                 wt->center0[1] + wt->center_vel[1] * t,
+                 wt->center0[2] + wt->center_vel[2] * t},
+      .velocity = {wt->center_vel[0], wt->center_vel[1], wt->center_vel[2]},
+      .radius = wt->radius0 + wt->radius_rate * t,
+      .radius_rate = wt->radius_rate,
+      .velocity_constant = 1,
+      .valid = 1};
+  return 0;
+}
+
+static const SpacetimeOps entry_ops = {
+    .eval = entry_eval_trap,
+    .escape_worldtube_sample = entry_worldtube_cb,
+};
+
+static SpacetimeSource entry_source(EntryWorldtube *wt) {
+  return (SpacetimeSource){.ops = &entry_ops, .context = wt};
+}
+
+/* ------------------------------------------------------------------ */
+/* Analytic path-parameter evaluator                                  */
+/* ------------------------------------------------------------------ */
+
+typedef struct {
+  double camera_t;
+  double camera_x[3];
+  double w[3]; /* unit past direction (straight mode) */
+  int arc_mode;
+  double arc_center[3];
+  double arc_radius;
+  double arc_theta0;
+  double L0;
+  double L0camera;
+  int evaluator_fails_at;
+  AsymptoticStatus fail_status;
+  int evaluator_call_count;
+  int nonfinite_at;
+  int reversed_time_at;
+} EntryEvaluator;
+
+static void entry_trajectory(const EntryEvaluator *c, double parameter,
+                             double x[3], double w[3], double *t) {
+  if (c->arc_mode) {
+    /* Circular analytic arc: parameter is arc length.  Not a physical
+     * geodesic, but a generic curved callback that exercises the driver beyond
+     * straight lines. */
+    const double theta = c->arc_theta0 + parameter / c->arc_radius;
+    x[0] = c->arc_center[0] + c->arc_radius * cos(theta);
+    x[1] = c->arc_center[1] + c->arc_radius * sin(theta);
+    x[2] = c->arc_center[2];
+    w[0] = -sin(theta);
+    w[1] = cos(theta);
+    w[2] = 0.0;
+  } else {
+    for (int i = 0; i < 3; ++i) {
+      x[i] = c->camera_x[i] + parameter * c->w[i];
+      w[i] = c->w[i];
+    }
+  }
+  *t = c->camera_t - parameter;
+}
+
+static AsymptoticStatus entry_evaluator_cb(void *context, double parameter,
+                                           AsymptoticRoute *state) {
+  EntryEvaluator *c = context;
+  ++c->evaluator_call_count;
+  if (c->evaluator_fails_at > 0 &&
+      c->evaluator_call_count == c->evaluator_fails_at)
+    return c->fail_status;
+  double x[3], w[3], t;
+  entry_trajectory(c, parameter, x, w, &t);
+  *state = (AsymptoticRoute){0};
+  state->kind = ASYMPTOTIC_ROUTE_ENTRY;
+  state->end_id = 0;
+  state->activate_t = t;
+  for (int i = 0; i < 3; ++i) {
+    state->x[i] = x[i];
+    state->Pi[i] = -w[i];
+  }
+  state->log_alpha_p0 = c->L0;
+  state->log_alpha_p0_camera = c->L0camera;
+  if (c->evaluator_call_count == c->nonfinite_at)
+    state->log_alpha_p0_camera = NAN;
+  if (c->evaluator_call_count == c->reversed_time_at)
+    state->activate_t = c->camera_t + 1.0;
+  return ASYMPTOTIC_OK;
+}
+
+/* Independent test-side oracle: the same worldtube F the driver sees, but
+ * computed directly from the analytic trajectory.  Used only to find the true
+ * first entry for comparison. */
+typedef struct {
+  const EntryEvaluator *ev;
+  const EntryWorldtube *wt;
+} EntryOracle;
+
+static double entry_oracle_F(void *context, double parameter) {
+  const EntryOracle *o = context;
+  double x[3], w[3], t;
+  entry_trajectory(o->ev, parameter, x, w, &t);
+  double d[3];
+  for (int i = 0; i < 3; ++i)
+    d[i] = x[i] - (o->wt->center0[i] + o->wt->center_vel[i] * t);
+  const double d2 = d[0] * d[0] + d[1] * d[1] + d[2] * d[2];
+  const double radius = o->wt->radius0 + o->wt->radius_rate * t;
+  return d2 - radius * radius;
+}
+
+static double entry_oracle_root(const EntryEvaluator *ev,
+                                const EntryWorldtube *wt, double lo,
+                                double hi) {
+  EntryOracle o = {.ev = ev, .wt = wt};
+  if (!(entry_oracle_F(&o, lo) >= 0.0 && entry_oracle_F(&o, hi) < 0.0))
+    return NAN;
+  for (int i = 0; i < 200; ++i) {
+    const double mid = 0.5 * (lo + hi);
+    if (!(mid > lo && mid < hi))
+      break;
+    if (entry_oracle_F(&o, mid) >= 0.0)
+      lo = mid;
+    else
+      hi = mid;
+  }
+  return 0.5 * (lo + hi);
+}
+
+static double path_parameter(const EntryEvaluator *ev,
+                             const AsymptoticRoute *state) {
+  return ev->camera_t - state->activate_t;
+}
+
+/* ------------------------------------------------------------------ */
+/* Tests                                                              */
+/* ------------------------------------------------------------------ */
+
+static void test_geometry_contract(void) {
+  EntryWorldtube wt = {.radius0 = 10.0, .valid_t_min = -1.0e300};
+  SpacetimeSource source = entry_source(&wt);
+  double F = NAN, tol = NAN;
+  RayReason reason = RAY_REASON_COUNT;
+
+  CHECK(asymptotic_entry_geometry(&source, 0, 0.0, (double[]){10.0, 0.0, 0.0},
+                                  &F, &tol, &reason) == ASYMPTOTIC_OK,
+        "boundary geometry status");
+  CHECK(F == 0.0, "boundary F is exactly zero");
+  CHECK(tol > 0.0 && isfinite(tol), "boundary tolerance finite positive");
+  CHECK(reason == RAY_REASON_NONE, "boundary reason none");
+
+  CHECK(asymptotic_entry_geometry(&source, 0, 0.0, (double[]){20.0, 0.0, 0.0},
+                                  &F, &tol, &reason) == ASYMPTOTIC_OK &&
+            F == 300.0,
+        "outside F is positive 300");
+  CHECK(asymptotic_entry_geometry(&source, 0, 0.0, (double[]){5.0, 0.0, 0.0},
+                                  &F, &tol, &reason) == ASYMPTOTIC_OK &&
+            F == -75.0,
+        "inside F is negative 75");
+
+  /* History exhaustion beats a miss. */
+  EntryWorldtube hole = {.radius0 = 10.0, .valid_t_min = 0.0};
+  source = entry_source(&hole);
+  CHECK(asymptotic_entry_geometry(&source, 0, -1.0, (double[]){10.0, 0.0, 0.0},
+                                  &F, &tol, &reason) ==
+                ASYMPTOTIC_TIME_RANGE_EXHAUSTED &&
+            reason == RAY_REASON_TIME_RANGE_EXHAUSTED,
+        "valid=0 is history exhaustion");
+
+  /* Callback failure is distinct from an invalid geometry. */
+  EntryWorldtube fail = {.radius0 = 10.0, .callback_fails = 1};
+  source = entry_source(&fail);
+  CHECK(asymptotic_entry_geometry(&source, 0, 0.0, (double[]){10.0, 0.0, 0.0},
+                                  &F, &tol, &reason) == ASYMPTOTIC_INVALID &&
+            reason == RAY_REASON_WORLDTUBE_SAMPLE_FAILED,
+        "callback failure reason");
+
+  EntryWorldtube nanr = {.radius0 = 10.0, .nan_radius = 1};
+  source = entry_source(&nanr);
+  CHECK(asymptotic_entry_geometry(&source, 0, 0.0, (double[]){10.0, 0.0, 0.0},
+                                  &F, &tol, &reason) == ASYMPTOTIC_INVALID &&
+            reason == RAY_REASON_WORLDTUBE_GEOMETRY_INVALID,
+        "NaN radius is invalid geometry");
+
+  EntryWorldtube zeror = {.radius0 = 10.0, .zero_radius = 1};
+  source = entry_source(&zeror);
+  CHECK(asymptotic_entry_geometry(&source, 0, 0.0, (double[]){10.0, 0.0, 0.0},
+                                  &F, &tol, &reason) == ASYMPTOTIC_INVALID &&
+            reason == RAY_REASON_WORLDTUBE_GEOMETRY_INVALID,
+        "non-positive radius is invalid geometry");
+
+  source = entry_source(&wt);
+  CHECK(asymptotic_entry_geometry(&source, 0, NAN, (double[]){10.0, 0.0, 0.0},
+                                  &F, &tol, &reason) == ASYMPTOTIC_INVALID &&
+            reason == RAY_REASON_WORLDTUBE_GEOMETRY_INVALID,
+        "NaN time is invalid geometry");
+  CHECK(asymptotic_entry_geometry(NULL, 0, 0.0, (double[]){10.0, 0.0, 0.0},
+                                  &F, &tol, &reason) == ASYMPTOTIC_INVALID &&
+            reason == RAY_REASON_INVALID_ARGUMENT,
+        "NULL source rejected");
+  CHECK(asymptotic_entry_geometry(&source, 0, 0.0, NULL, &F, &tol,
+                                  &reason) == ASYMPTOTIC_INVALID,
+        "NULL position rejected");
+}
+
+static void test_validate_contract(void) {
+  EntryWorldtube wt = {.radius0 = 10.0, .valid_t_min = -1.0e300};
+  SpacetimeSource source = entry_source(&wt);
+  int valid = -1;
+  RayReason reason = RAY_REASON_COUNT;
+  AsymptoticRoute candidate = {.kind = ASYMPTOTIC_ROUTE_ENTRY,
+                               .end_id = 0,
+                               .activate_t = 0.0,
+                               .x = {10.0, 0.0, 0.0},
+                               .Pi = {-1.0, 0.0, 0.0},
+                               .log_alpha_p0 = 0.5,
+                               .log_alpha_p0_camera = 0.25};
+
+  CHECK(asymptotic_entry_validate(&source, 0, &candidate, &valid, &reason) ==
+                ASYMPTOTIC_OK &&
+            valid == 1,
+        "boundary candidate is valid");
+
+  candidate.x[0] = 11.0; /* F = 21 > tol */
+  CHECK(asymptotic_entry_validate(&source, 0, &candidate, &valid, &reason) ==
+                ASYMPTOTIC_OK &&
+            valid == 0 && reason == RAY_REASON_NONE,
+        "outside candidate is valid=0 with OK status");
+
+  candidate.x[0] = 5.0; /* F = -75, far inside */
+  CHECK(asymptotic_entry_validate(&source, 0, &candidate, &valid, &reason) ==
+                ASYMPTOTIC_OK &&
+            valid == 0,
+        "deep-inside candidate is valid=0 with OK status");
+
+  candidate.x[0] = 10.0;
+  candidate.kind = ASYMPTOTIC_ROUTE_ESCAPED;
+  CHECK(asymptotic_entry_validate(&source, 0, &candidate, &valid, &reason) ==
+                ASYMPTOTIC_INVALID &&
+            valid == 0 && reason == RAY_REASON_PROTOCOL_ERROR,
+        "wrong candidate kind is a protocol error");
+
+  candidate.kind = ASYMPTOTIC_ROUTE_ENTRY;
+  candidate.x[0] = NAN;
+  CHECK(asymptotic_entry_validate(&source, 0, &candidate, &valid, &reason) ==
+                ASYMPTOTIC_INVALID &&
+            reason == RAY_REASON_PROTOCOL_ERROR,
+        "non-finite candidate is a protocol error");
+
+  /* History exhaustion propagates through validation. */
+  candidate.x[0] = 10.0;
+  candidate.activate_t = -1.0;
+  EntryWorldtube hole = {.radius0 = 10.0, .valid_t_min = 0.0};
+  source = entry_source(&hole);
+  CHECK(asymptotic_entry_validate(&source, 0, &candidate, &valid, &reason) ==
+                ASYMPTOTIC_TIME_RANGE_EXHAUSTED &&
+            reason == RAY_REASON_TIME_RANGE_EXHAUSTED,
+        "validation propagates history exhaustion");
+}
+
+static void test_fixed_sphere_localize(void) {
+  EntryWorldtube wt = {.radius0 = 10.0, .valid_t_min = -1.0e300};
+  SpacetimeSource source = entry_source(&wt);
+  EntryEvaluator ev = {.camera_t = 0.0,
+                       .camera_x = {50.0, 0.0, 0.0},
+                       .w = {-1.0, 0.0, 0.0},
+                       .L0 = 0.75,
+                       .L0camera = 0.5};
+  AsymptoticRoute out;
+  unsigned int evaluations = 0;
+  RayReason reason = RAY_REASON_COUNT;
+
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0,
+                                  45.0, &out, &evaluations,
+                                  &reason) == ASYMPTOTIC_OK,
+        "fixed-sphere localize succeeds");
+  CHECK(out.kind == ASYMPTOTIC_ROUTE_ENTRY && out.end_id == 0,
+        "localized kind and end");
+  const double s = path_parameter(&ev, &out);
+  const double s_true = entry_oracle_root(&ev, &wt, 30.0, 45.0);
+  CHECK(isfinite(s_true), "oracle found the same bracket");
+  CHECK(s >= s_true && s - s_true <= 1e-9,
+        "localized just past first entry");
+  CHECK(fabs(s - 40.0) <= 1e-9, "fixed-sphere entry at s=40");
+  CHECK(out.Pi[0] == 1.0 && out.Pi[1] == 0.0 && out.Pi[2] == 0.0,
+        "direction preserved exactly");
+  CHECK(out.log_alpha_p0 == 0.75 && out.log_alpha_p0_camera == 0.5,
+        "L and camera L preserved exactly");
+  CHECK(evaluations >= 2 && evaluations <= 260, "evaluation count bounded");
+  CHECK(ev.evaluator_call_count == (int)evaluations,
+        "evaluator calls counted once each");
+  CHECK(wt.eval_calls == 0, "no metric evaluation outside the worldtube");
+}
+
+static void test_too_early_hint(void) {
+  /* Outside endpoint is the camera (a deliberately too-early, corrupted
+   * bracket); the localizer still returns the true first entry, not the
+   * inside hint and not the camera. */
+  EntryWorldtube wt = {.radius0 = 10.0, .valid_t_min = -1.0e300};
+  SpacetimeSource source = entry_source(&wt);
+  EntryEvaluator ev = {.camera_t = 0.0,
+                       .camera_x = {50.0, 0.0, 0.0},
+                       .w = {-1.0, 0.0, 0.0},
+                       .L0 = 0.1,
+                       .L0camera = 0.2};
+  AsymptoticRoute out;
+  unsigned int evaluations = 0;
+  RayReason reason = RAY_REASON_COUNT;
+
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 0.0,
+                                  45.0, &out, &evaluations,
+                                  &reason) == ASYMPTOTIC_OK,
+        "too-early hint still localizes");
+  const double s = path_parameter(&ev, &out);
+  CHECK(fabs(s - 40.0) <= 1e-9, "returns actual first entry, not the hint");
+  CHECK(s > 1.0 && s < 45.0, "not the camera and not the inside hint");
+  CHECK(evaluations <= 260, "hint evaluation budget");
+}
+
+static void test_moving_sphere_localize(void) {
+  EntryWorldtube wt = {.radius0 = 10.0,
+                       .center_vel = {0.5, 0.0, 0.0},
+                       .valid_t_min = -1.0e300};
+  SpacetimeSource source = entry_source(&wt);
+  EntryEvaluator ev = {.camera_t = 0.0,
+                       .camera_x = {100.0, 0.0, 0.0},
+                       .w = {-1.0, 0.0, 0.0},
+                       .L0 = 0.3,
+                       .L0camera = 0.4};
+  AsymptoticRoute out;
+  unsigned int evaluations = 0;
+  RayReason reason = RAY_REASON_COUNT;
+
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 150.0,
+                                  200.0, &out, &evaluations,
+                                  &reason) == ASYMPTOTIC_OK,
+        "translated+ moving sphere localize");
+  const double s = path_parameter(&ev, &out);
+  const double s_true = entry_oracle_root(&ev, &wt, 150.0, 200.0);
+  CHECK(fabs(s - s_true) <= 1e-8 && fabs(s - 180.0) <= 1e-8,
+        "moving-sphere entry at s=180");
+  CHECK(evaluations <= 260, "moving-sphere evaluation budget");
+}
+
+static void test_radius_rate_localize(void) {
+  /* radius(t) = radius0 + radius_rate * t with radius_rate = -1 and t = -s, so
+   * R grows as 10 + s; the entry is at s = 45. */
+  EntryWorldtube wt = {.radius0 = 10.0,
+                       .radius_rate = -1.0,
+                       .valid_t_min = -1.0e300};
+  SpacetimeSource source = entry_source(&wt);
+  EntryEvaluator ev = {.camera_t = 0.0,
+                       .camera_x = {100.0, 0.0, 0.0},
+                       .w = {-1.0, 0.0, 0.0},
+                       .L0 = 0.6,
+                       .L0camera = 0.6};
+  AsymptoticRoute out;
+  unsigned int evaluations = 0;
+  RayReason reason = RAY_REASON_COUNT;
+
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0,
+                                  60.0, &out, &evaluations,
+                                  &reason) == ASYMPTOTIC_OK,
+        "linear radius-rate localize");
+  const double s = path_parameter(&ev, &out);
+  const double s_true = entry_oracle_root(&ev, &wt, 30.0, 60.0);
+  CHECK(fabs(s - s_true) <= 1e-8 && fabs(s - 45.0) <= 1e-8,
+        "linear radius-rate entry at s=45");
+  CHECK(evaluations <= 260, "radius-rate evaluation budget");
+}
+
+static void test_rotated_frame_localize(void) {
+  /* Camera on a rotated axis: (40,30,0), past direction toward the origin. */
+  EntryWorldtube wt = {.radius0 = 10.0, .valid_t_min = -1.0e300};
+  SpacetimeSource source = entry_source(&wt);
+  EntryEvaluator ev = {.camera_t = 0.0,
+                       .camera_x = {40.0, 30.0, 0.0},
+                       .w = {-0.8, -0.6, 0.0},
+                       .L0 = 0.2,
+                       .L0camera = 0.1};
+  AsymptoticRoute out;
+  unsigned int evaluations = 0;
+  RayReason reason = RAY_REASON_COUNT;
+
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0,
+                                  45.0, &out, &evaluations,
+                                  &reason) == ASYMPTOTIC_OK,
+        "rotated flat frame localize");
+  const double s = path_parameter(&ev, &out);
+  const double s_true = entry_oracle_root(&ev, &wt, 30.0, 45.0);
+  CHECK(fabs(s - s_true) <= 1e-9 && fabs(s - 40.0) <= 1e-9,
+        "rotated-frame entry at s=40");
+  CHECK(fabs(out.x[1] - 6.0) <= 1e-6, "rotated entry position on sphere");
+}
+
+static void test_grazing_first_entry(void) {
+  /* Grazing pass: the camera is offset by 9.9 from the sphere axis.  The first
+   * entry at s ~ 48.589 is inside the bracket; the exit at s ~ 51.410 is not.
+   * Bisection must return the first entry, not the later exit. */
+  EntryWorldtube wt = {.radius0 = 10.0, .valid_t_min = -1.0e300};
+  SpacetimeSource source = entry_source(&wt);
+  EntryEvaluator ev = {.camera_t = 0.0,
+                       .camera_x = {50.0, 9.9, 0.0},
+                       .w = {-1.0, 0.0, 0.0},
+                       .L0 = 0.0,
+                       .L0camera = 0.0};
+  AsymptoticRoute out;
+  unsigned int evaluations = 0;
+  RayReason reason = RAY_REASON_COUNT;
+
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 48.0,
+                                  50.0, &out, &evaluations,
+                                  &reason) == ASYMPTOTIC_OK,
+        "grazing first entry localizes");
+  const double s = path_parameter(&ev, &out);
+  const double first = 50.0 - sqrt(100.0 - 9.9 * 9.9);
+  CHECK(fabs(s - first) <= 1e-8, "grazing entry is the first crossing");
+  CHECK(s < 51.4, "not the later exit crossing");
+  CHECK(fabs(out.x[1] - 9.9) <= 1e-9, "grazing impact parameter preserved");
+  CHECK(evaluations <= 260, "grazing evaluation budget");
+}
+
+static void test_curved_arc_localize(void) {
+  /* Circular analytic arc of radius 30 and worldtube centered at (25,0,0)
+   * radius 8; entry at arc length ~ 87.40. */
+  EntryWorldtube wt = {.radius0 = 8.0,
+                       .center0 = {25.0, 0.0, 0.0},
+                       .valid_t_min = -1.0e300};
+  SpacetimeSource source = entry_source(&wt);
+  EntryEvaluator ev = {.camera_t = 0.0,
+                       .arc_mode = 1,
+                       .arc_center = {0.0, 0.0, 0.0},
+                       .arc_radius = 30.0,
+                       .arc_theta0 = TEST_PI,
+                       .L0 = 0.9,
+                       .L0camera = 0.8};
+  AsymptoticRoute out;
+  unsigned int evaluations = 0;
+  RayReason reason = RAY_REASON_COUNT;
+
+  const double lo = 30.0 * (5.9 - TEST_PI);
+  const double hi = 30.0 * (6.2 - TEST_PI);
+  const double s_true = entry_oracle_root(&ev, &wt, lo, hi);
+  CHECK(isfinite(s_true), "curved oracle bracket");
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, lo, hi,
+                                  &out, &evaluations,
+                                  &reason) == ASYMPTOTIC_OK,
+        "curved arc localize");
+  const double s = path_parameter(&ev, &out);
+  CHECK(s >= s_true - 1e-9 && s - s_true <= 1e-8,
+        "curved arc entry matches the oracle");
+  /* d^2(theta) = 1525 - 1500 cos(theta) = 8^2 on the arc. */
+  const double expected =
+      30.0 * (2.0 * TEST_PI - acos((1525.0 - 64.0) / 1500.0) - TEST_PI);
+  CHECK(fabs(s - expected) <= 1e-8, "curved arc entry matches analytic root");
+  CHECK(evaluations <= 260, "curved arc evaluation budget");
+}
+
+static void test_boundary_entry_exact(void) {
+  EntryWorldtube wt = {.radius0 = 10.0, .valid_t_min = -1.0e300};
+  SpacetimeSource source = entry_source(&wt);
+  EntryEvaluator ev = {.camera_t = 0.0,
+                       .camera_x = {50.0, 0.0, 0.0},
+                       .w = {-1.0, 0.0, 0.0},
+                       .L0 = 0.4,
+                       .L0camera = 0.4};
+  AsymptoticRoute out;
+  unsigned int evaluations = 0;
+  RayReason reason = RAY_REASON_COUNT;
+
+  /* F == 0 exactly at the outside endpoint and strictly inside at 45. */
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 40.0,
+                                  45.0, &out, &evaluations,
+                                  &reason) == ASYMPTOTIC_OK,
+        "exact boundary entry localize");
+  CHECK(out.kind == ASYMPTOTIC_ROUTE_ENTRY && out.activate_t == -40.0 &&
+            out.x[0] == 10.0,
+        "boundary endpoint returned directly");
+  CHECK(evaluations == 2, "boundary path needs no bisection");
+}
+
+static void test_unconfirmed_bracket(void) {
+  EntryWorldtube wt = {.radius0 = 10.0, .valid_t_min = -1.0e300};
+  SpacetimeSource source = entry_source(&wt);
+  EntryEvaluator ev = {.camera_t = 0.0,
+                       .camera_x = {50.0, 0.0, 0.0},
+                       .w = {-1.0, 0.0, 0.0}};
+  AsymptoticRoute out;
+  unsigned int evaluations = 0;
+  RayReason reason = RAY_REASON_COUNT;
+
+  /* Both endpoints outside: no strict-inside bracket. */
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 10.0,
+                                  20.0, &out, &evaluations,
+                                  &reason) == ASYMPTOTIC_INVALID &&
+            reason == RAY_REASON_ENTRY_UNCONFIRMED,
+        "false candidate outside bracket is unconfirmed, not escaped");
+
+  /* Both endpoints strictly inside: also no entry bracket. */
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 45.0,
+                                  50.0, &out, &evaluations,
+                                  &reason) == ASYMPTOTIC_INVALID &&
+            reason == RAY_REASON_ENTRY_UNCONFIRMED,
+        "both-inside bracket is unconfirmed");
+
+  /* Reversed bracket ordering. */
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 45.0,
+                                  30.0, &out, &evaluations,
+                                  &reason) == ASYMPTOTIC_INVALID &&
+            reason == RAY_REASON_INVALID_ARGUMENT,
+        "reversed bracket rejected");
+}
+
+static void test_callback_failure_propagation(void) {
+  EntryEvaluator ev = {.camera_t = 0.0,
+                       .camera_x = {50.0, 0.0, 0.0},
+                       .w = {-1.0, 0.0, 0.0}};
+  AsymptoticRoute out;
+  unsigned int evaluations = 0;
+  RayReason reason = RAY_REASON_COUNT;
+
+  EntryWorldtube fail = {.radius0 = 10.0,
+                         .valid_t_min = -1.0e300,
+                         .callback_fails = 1};
+  SpacetimeSource source = entry_source(&fail);
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0,
+                                  45.0, &out, &evaluations,
+                                  &reason) == ASYMPTOTIC_INVALID &&
+            reason == RAY_REASON_WORLDTUBE_SAMPLE_FAILED,
+        "endpoint callback failure propagates");
+
+  fail = (EntryWorldtube){.radius0 = 10.0,
+                          .valid_t_min = -1.0e300,
+                          .fail_at_call = 2};
+  source = entry_source(&fail);
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0,
+                                  45.0, &out, &evaluations,
+                                  &reason) == ASYMPTOTIC_INVALID &&
+            reason == RAY_REASON_WORLDTUBE_SAMPLE_FAILED,
+        "inside endpoint callback failure propagates");
+
+  fail = (EntryWorldtube){.radius0 = 10.0,
+                          .valid_t_min = -1.0e300,
+                          .fail_at_call = 3};
+  source = entry_source(&fail);
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0,
+                                  45.0, &out, &evaluations,
+                                  &reason) == ASYMPTOTIC_INVALID &&
+            reason == RAY_REASON_WORLDTUBE_SAMPLE_FAILED,
+        "midpoint callback failure propagates");
+}
+
+static void test_history_hole_propagation(void) {
+  /* The inside endpoint falls past the valid history: the driver must report
+   * TIME_RANGE_EXHAUSTED, never a miss or a fabricated entry. */
+  EntryWorldtube wt = {.radius0 = 10.0, .valid_t_min = -40.0};
+  SpacetimeSource source = entry_source(&wt);
+  EntryEvaluator ev = {.camera_t = 0.0,
+                       .camera_x = {50.0, 0.0, 0.0},
+                       .w = {-1.0, 0.0, 0.0}};
+  AsymptoticRoute out;
+  unsigned int evaluations = 0;
+  RayReason reason = RAY_REASON_COUNT;
+
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0,
+                                  45.0, &out, &evaluations,
+                                  &reason) == ASYMPTOTIC_TIME_RANGE_EXHAUSTED &&
+            reason == RAY_REASON_TIME_RANGE_EXHAUSTED,
+        "endpoint history hole propagates");
+
+  /* A midpoint-only history hole: both endpoints are valid, but the first
+   * bisection midpoint (t = -37.5) falls in an isolated invalid window. */
+  wt = (EntryWorldtube){.radius0 = 10.0,
+                        .valid_t_min = -1.0e300,
+                        .hole_center = -37.5,
+                        .hole_halfwidth = 0.5};
+  source = entry_source(&wt);
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0,
+                                  45.0, &out, &evaluations,
+                                  &reason) == ASYMPTOTIC_TIME_RANGE_EXHAUSTED &&
+            reason == RAY_REASON_TIME_RANGE_EXHAUSTED,
+        "midpoint history hole propagates, never a miss");
+}
+
+static void test_evaluator_failure_propagation(void) {
+  EntryWorldtube wt = {.radius0 = 10.0, .valid_t_min = -1.0e300};
+  SpacetimeSource source = entry_source(&wt);
+  AsymptoticRoute out;
+  unsigned int evaluations = 0;
+  RayReason reason = RAY_REASON_COUNT;
+
+  EntryEvaluator ev = {.camera_t = 0.0,
+                       .camera_x = {50.0, 0.0, 0.0},
+                       .w = {-1.0, 0.0, 0.0},
+                       .evaluator_fails_at = 1,
+                       .fail_status = ASYMPTOTIC_TIME_RANGE_EXHAUSTED};
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0,
+                                  45.0, &out, &evaluations,
+                                  &reason) == ASYMPTOTIC_TIME_RANGE_EXHAUSTED &&
+            reason == RAY_REASON_TIME_RANGE_EXHAUSTED,
+        "evaluator history failure propagates");
+
+  ev = (EntryEvaluator){.camera_t = 0.0,
+                        .camera_x = {50.0, 0.0, 0.0},
+                        .w = {-1.0, 0.0, 0.0},
+                        .evaluator_fails_at = 1,
+                        .fail_status = ASYMPTOTIC_INVALID};
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0,
+                                  45.0, &out, &evaluations,
+                                  &reason) == ASYMPTOTIC_INVALID &&
+            reason == RAY_REASON_ENTRY_UNCONFIRMED,
+        "evaluator invalid failure maps to unconfirmed");
+
+  ev = (EntryEvaluator){.camera_t = 0.0,
+                        .camera_x = {50.0, 0.0, 0.0},
+                        .w = {-1.0, 0.0, 0.0},
+                        .evaluator_fails_at = 2,
+                        .fail_status = ASYMPTOTIC_INVALID};
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 30.0,
+                                  45.0, &out, &evaluations,
+                                  &reason) == ASYMPTOTIC_INVALID,
+        "inside endpoint evaluator failure propagates");
+}
+
+static AsymptoticStatus wide_parameter_path(void *context, double parameter,
+                                            AsymptoticRoute *state) {
+  (void)context;
+  *state = (AsymptoticRoute){.kind = ASYMPTOTIC_ROUTE_ENTRY,
+                            .end_id = 0,
+                            .activate_t = -parameter,
+                            .x = {2.0 - parameter * 1e-308, 0.0, 0.0},
+                            .Pi = {1.0, 0.0, 0.0}};
+  return ASYMPTOTIC_OK;
+}
+
+static void test_representability_and_state_checks(void) {
+  EntryWorldtube wt = {.radius0 = 1.0, .valid_t_min = -DBL_MAX};
+  SpacetimeSource source = entry_source(&wt);
+  AsymptoticRoute out;
+  unsigned int evaluations = 0;
+  RayReason reason = RAY_REASON_COUNT;
+  /* Both endpoints are finite, but subtracting them overflows.  This must not
+   * be mistaken for an adjacent bracket and return the far-inside endpoint. */
+  CHECK(asymptotic_entry_localize(&source, 0, wide_parameter_path, NULL,
+                                  -1.6e308, 1.6e308, &out, &evaluations,
+                                  &reason) == ASYMPTOTIC_OK,
+         "overflow-safe parameter midpoint");
+  CHECK(fabs(out.x[0] - 1.0) < 1e-14 && evaluations > 2,
+         "wide bracket contracts to entry, not initial inside endpoint");
+  wt.radius0 = 10.0;
+  EntryEvaluator ev = {.camera_x = {50.0, 0.0, 0.0},
+                        .w = {-1.0, 0.0, 0.0}, .nonfinite_at = 1};
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 40.0,
+                                  45.0, &out, &evaluations, &reason) ==
+            ASYMPTOTIC_INVALID && reason == RAY_REASON_ENTRY_UNCONFIRMED,
+         "boundary shortcut rejects nonfinite camera energy reference");
+  ev.evaluator_call_count = 0;
+  ev.nonfinite_at = 3;
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 35.0,
+                                  45.0, &out, &evaluations, &reason) ==
+            ASYMPTOTIC_INVALID && reason == RAY_REASON_ENTRY_UNCONFIRMED,
+         "midpoint rejects nonfinite camera energy reference");
+  ev.evaluator_call_count = 0;
+  ev.nonfinite_at = 0;
+  ev.reversed_time_at = 3;
+  CHECK(asymptotic_entry_localize(&source, 0, entry_evaluator_cb, &ev, 35.0,
+                                  45.0, &out, &evaluations, &reason) ==
+            ASYMPTOTIC_INVALID && reason == RAY_REASON_ENTRY_UNCONFIRMED,
+         "midpoint cannot reverse coordinate time");
+}
+
+int main(void) {
+  test_geometry_contract();
+  test_validate_contract();
+  test_fixed_sphere_localize();
+  test_too_early_hint();
+  test_moving_sphere_localize();
+  test_radius_rate_localize();
+  test_rotated_frame_localize();
+  test_grazing_first_entry();
+  test_curved_arc_localize();
+  test_boundary_entry_exact();
+  test_unconfirmed_bracket();
+  test_callback_failure_propagation();
+  test_history_hole_propagation();
+  test_evaluator_failure_propagation();
+  test_representability_and_state_checks();
+  if (failures == 0)
+    puts("asymptotic entry regression passed");
+  else
+    fprintf(stderr, "%d asymptotic entry regression failures\n", failures);
+  return failures == 0 ? 0 : 1;
+}

tests/test_asymptotic_quadratic.c

@@ -0,0 +1,65 @@
+/* Exercise the private numerical kernel directly, including coefficient
+ * ranges that cannot be represented by a public double worldtube fixture.
+ * The build rule omits the separately compiled asymptotic.c. */
+#include "../src/asymptotic.c"
+
+#include <stdio.h>
+
+static int failures;
+
+#define CHECK(condition, message) do { \
+  if (!(condition)) { \
+    fprintf(stderr, "FAIL %s:%d: %s\n", __FILE__, __LINE__, message); \
+    ++failures; \
+  } \
+} while (0)
+
+static void check_scaled(int exponent) {
+  const long double scale = scalbnl(1.0L, exponent);
+  double root = -1.0;
+  EntryQuadratic k = {scale, -3.0L * scale, 2.0L * scale};
+  CHECK(entry_solve(&k, &root) == ENTRY_SOLVE_ENTRY && root == 1.0,
+        "common scale preserves smallest inward root");
+  k = (EntryQuadratic){scale, -scale, scale};
+  CHECK(entry_solve(&k, &root) == ENTRY_SOLVE_MISS,
+        "common scale preserves a clear miss");
+  k = (EntryQuadratic){0.0L, -scale, 2.0L * scale};
+  CHECK(entry_solve(&k, &root) == ENTRY_SOLVE_ENTRY && root == 2.0,
+        "common scale preserves linear entry");
+  k = (EntryQuadratic){scale, -scale, 0.0L};
+  CHECK(entry_solve(&k, &root) == ENTRY_SOLVE_ENTRY && root == 0.0,
+        "common scale preserves boundary entry");
+  k = (EntryQuadratic){-scale, scale, 2.0L * scale};
+  CHECK(entry_solve(&k, &root) == ENTRY_SOLVE_ENTRY && root == 2.0,
+        "common scale preserves concave entry");
+}
+
+static void test_product_cancellation(void) {
+  const long double u = scalbnl(1.0L, 1 - LDBL_MANT_DIG);
+  const EntryQuadratic k = {1.0L + u, -2.0L, 1.0L - 0.5L * u};
+  long double scale;
+  (void)entry_discriminant(&k, &scale);
+  /* Exact dyadic oracle: 4 - 4(1+u)(1-u/2) = -2u + 2u^2.
+   * A separately rounded 4*a*c is 4 and loses this nonzero discriminant. */
+  double root;
+  CHECK(entry_solve(&k, &root) == ENTRY_SOLVE_UNCERTAIN,
+         "product cancellation must remain uncertain, not a proven miss");
+}
+
+int main(void) {
+  check_scaled(0);
+  check_scaled(LDBL_MAX_EXP - 4);
+  check_scaled(LDBL_MIN_EXP + 4);
+  check_scaled(LDBL_MIN_EXP - LDBL_MANT_DIG + 2);
+  test_product_cancellation();
+  double root;
+  EntryQuadratic k = {LDBL_MIN, LDBL_MAX / 8.0L, 1.0L};
+  CHECK(entry_solve(&k, &root) == ENTRY_SOLVE_UNCERTAIN,
+        "scaling cannot silently erase a nonzero coefficient");
+  k = (EntryQuadratic){1.0L, 2.0L, -INFINITY};
+  CHECK(entry_solve(&k, &root) == ENTRY_SOLVE_UNCERTAIN,
+        "nonfinite coefficient is not a normal entry");
+  if (!failures)
+    puts("asymptotic quadratic regression passed");
+  return failures ? 1 : 0;
+}

tests/test_asymptotic_schwarzschild.c

@@ -6,6 +6,7 @@
 
 #include <math.h>
 #include <stdio.h>
+#include <stdlib.h>
 
 static int failures = 0;
 #define CHECK(condition, message)                                        \
@@ -189,6 +190,8 @@ static void test_preroute_entry(void) {
                 ASYMPTOTIC_OK &&
             route.kind == ASYMPTOTIC_ROUTE_ENTRY,
         "outside camera enters");
+  CHECK(route.entry_fallback_evaluations == 0,
+        "analytic entry stays on the fast path");
   double value;
   CHECK(asymptotic_worldtube_value(&source, route.end_id, route.activate_t,
                                    route.x, &value) == 0 &&
@@ -521,6 +524,122 @@ static void test_preroute_branches(void) {
   spacetime_destroy(&source);
 }
 
+/* Translated-origin wrapper around the analytic Schwarzschild KS source: the
+ * inner metric is evaluated at x - origin and the worldtube/end are shifted by
+ * the same origin.  This models a black hole at a large coordinate offset; the
+ * double reconstruction of the boundary state loses the sub-ULP offset and
+ * trips the common entry fallback, while the exact inward orbit transfer stays
+ * valid.  It exercises the production fallback path, not a synthetic
+ * nonlinearity. */
+typedef struct {
+  SpacetimeSource inner;
+  double origin[3];
+} ShiftedOriginContext;
+
+static SpacetimePointStatus shifted_origin_eval(const SpacetimeSource *source,
+                                                double t, const double x[3],
+                                                MetricData *metric) {
+  const ShiftedOriginContext *ctx = source->context;
+  const double local[3] = {x[0] - ctx->origin[0], x[1] - ctx->origin[1],
+                           x[2] - ctx->origin[2]};
+  return spacetime_eval(&ctx->inner, t, local, metric);
+}
+
+static SpacetimeRayStatus shifted_origin_classify(const SpacetimeSource *source,
+                                                  double t,
+                                                  const double x[3]) {
+  const ShiftedOriginContext *ctx = source->context;
+  const double local[3] = {x[0] - ctx->origin[0], x[1] - ctx->origin[1],
+                           x[2] - ctx->origin[2]};
+  return spacetime_classify(&ctx->inner, t, local);
+}
+
+static size_t shifted_origin_end_count(const SpacetimeSource *source) {
+  const ShiftedOriginContext *ctx = source->context;
+  return spacetime_asymptotic_end_count(&ctx->inner);
+}
+
+static int shifted_origin_end(const SpacetimeSource *source, size_t index,
+                              SpacetimeAsymptoticEnd *out) {
+  const ShiftedOriginContext *ctx = source->context;
+  if (spacetime_asymptotic_end(&ctx->inner, index, out))
+    return -1;
+  for (int i = 0; i < 3; ++i)
+    out->frame_origin[i] = ctx->origin[i];
+  return 0;
+}
+
+static int shifted_origin_worldtube(const SpacetimeSource *source,
+                                    SpacetimeEndId end_id, double t,
+                                    SpacetimeEscapeWorldtubeSample *out) {
+  const ShiftedOriginContext *ctx = source->context;
+  if (spacetime_escape_worldtube_sample(&ctx->inner, end_id, t, out))
+    return -1;
+  for (int i = 0; i < 3; ++i)
+    out->center[i] += ctx->origin[i];
+  return 0;
+}
+
+static void shifted_origin_destroy(SpacetimeSource *source) {
+  ShiftedOriginContext *ctx = source->context;
+  if (ctx != NULL) {
+    spacetime_destroy(&ctx->inner);
+    free(ctx);
+  }
+  source->context = NULL;
+  source->ops = NULL;
+}
+
+static const SpacetimeOps shifted_origin_ops = {
+    .eval = shifted_origin_eval,
+    .classify = shifted_origin_classify,
+    .asymptotic_end_count = shifted_origin_end_count,
+    .asymptotic_end = shifted_origin_end,
+    .escape_worldtube_sample = shifted_origin_worldtube,
+    .destroy = shifted_origin_destroy};
+
+static void test_translated_origin_fallback(void) {
+  ShiftedOriginContext *ctx = malloc(sizeof *ctx);
+  CHECK(ctx != NULL, "shifted-origin context");
+  if (ctx == NULL)
+    return;
+  ctx->origin[0] = 1.0e6;
+  ctx->origin[1] = 2.0e6;
+  ctx->origin[2] = -3.0e6;
+  CHECK(spacetime_create_schwarzschild_ks(&ctx->inner, 1.0, 256.0) == 0,
+        "shifted-origin inner source");
+  SpacetimeSource source = {.ops = &shifted_origin_ops, .context = ctx};
+
+  ObserverCamera cam = {.look_ra_deg = 180.0, .look_dec_deg = 0.0};
+  for (int i = 0; i < 3; ++i)
+    cam.position[i] = ctx->origin[i];
+  cam.position[0] += 500.0;
+  MetricData metric;
+  CHECK(spacetime_eval(&source, 0.0, cam.position, &metric) == 0,
+        "shifted-origin camera metric");
+  ObserverState observer;
+  CHECK(observer_from_coordinate_camera(&metric, &cam, &observer, NULL) ==
+            OBSERVER_BUILD_OK,
+        "shifted-origin camera observer");
+  const double direction[3] = {cos(0.3), sin(0.3), 0.0};
+  AsymptoticRoute route;
+  CHECK(asymptotic_route_camera(&source, &observer, direction, &route) ==
+                ASYMPTOTIC_OK &&
+            route.kind == ASYMPTOTIC_ROUTE_ENTRY,
+        "shifted-origin entry found");
+  CHECK(route.entry_fallback_evaluations > 0,
+        "shifted-origin entry used the common fallback");
+  CHECK(route.failure_reason == RAY_REASON_NONE,
+        "shifted-origin fallback has no failure reason");
+  double value;
+  CHECK(asymptotic_worldtube_value(&source, route.end_id, route.activate_t,
+                                   route.x, &value) == 0 &&
+            value <= 0.0,
+        "shifted-origin fallback state is inside the worldtube");
+  CHECK(route.activate_t < 0.0, "shifted-origin entry is in the past");
+  spacetime_destroy(&source);
+}
+
 int main(void) {
   test_round_trip();
   test_finish_matches_integration();
@@ -532,6 +651,7 @@ int main(void) {
   test_time_reference();
   test_grazing_reference();
   test_preroute_branches();
+  test_translated_origin_fallback();
   if (failures == 0)
     puts("asymptotic schwarzschild regression passed");
   else

tests/test_frame.c

@@ -1469,6 +1469,31 @@ int main(void) {
       lens_map_destroy(&dloaded); unlink(dp_path); goto done;
     }
     lens_map_destroy(&dloaded);
+    /* Truncate within the first v3 vertex, including each terminal field.
+     * Failed reads must reject the map and release its partially read mesh. */
+    {
+      unsigned char prefix[316];
+      FILE *fixture = fopen(dp_path, "rb");
+      int fixture_failed = fixture == NULL ||
+          fread(prefix, 1, sizeof prefix, fixture) != sizeof prefix;
+      if (fixture != NULL && fclose(fixture)) fixture_failed = 1;
+      if (fixture_failed) {
+        fputs("lens-map truncation fixture read failed\n", stderr);
+        unlink(dp_path); goto done;
+      }
+      const size_t cuts[] = {232, 303, 304, 307, 308, 311, 312, 315};
+      for (size_t c = 0; c < sizeof cuts / sizeof cuts[0]; ++c) {
+        FILE *short_file = fopen(dp_path, "wb");
+        int short_failed = short_file == NULL ||
+            fwrite(prefix, 1, cuts[c], short_file) != cuts[c];
+        if (short_file != NULL && fclose(short_file)) short_failed = 1;
+        if (short_failed || !lens_map_read(dp_path, NULL, &dloaded) ||
+            dloaded.frames != NULL || dloaded.frame_count != 0) {
+          fputs("lens-map truncated vertex rejection regression failed\n", stderr);
+          lens_map_destroy(&dloaded); unlink(dp_path); goto done;
+        }
+      }
+    }
     /* Unknown wire code and non-finite/out-of-bounds DP fields must be rejected
      * by the shared schema validator, not accepted as a usable map. */
     {

tests/test_geodesic.c

@@ -53,6 +53,8 @@ static int check_reason_names(void) {
           RAY_REASON_INTEGRATION_ERROR ||
       ray_reason_category(RAY_REASON_INVALID_ESCAPE_DIRECTION) !=
           RAY_REASON_INTEGRATION_ERROR ||
+      ray_reason_category(RAY_REASON_ENTRY_UNCONFIRMED) !=
+          RAY_REASON_INTEGRATION_ERROR ||
       ray_reason_category(RAY_REASON_SLAB_LOAD_FAILED) != RAY_REASON_IO_ERROR) {
     fputs("detail reasons map to the wrong coarse category\n", stderr);
     failed = 1;